Ticket #14247: trac_14247-bell_number_mpmath-ts.patch

File trac_14247-bell_number_mpmath-ts.patch, 3.0 KB (added by Travis Scrimshaw, 10 years ago)
  • sage/combinat/combinat.py

    # HG changeset patch
    # User Travis Scrimshaw <tscrim@ucdavis.edu>
    # Date 1363028736 25200
    # Node ID fb896d8263fd66316d4aebf991bfc89558e4799a
    # Parent  7c4771346653db738f1aa04d2c6d0da3d1d1368e
    #14247: Give bell_number the option to temporarily set mpmath's precision.
    
    diff --git a/sage/combinat/combinat.py b/sage/combinat/combinat.py
    a b from sage.misc.superseded import depreca 
    203203from combinat_cython import _stirling_number2
    204204######### combinatorial sequences
    205205
    206 def bell_number(n, algorithm='dobinski'):
     206def bell_number(n, algorithm='dobinski', **options):
    207207    r"""
    208208    Return the `n`-th Bell number (the number of ways to partition a set
    209209    of n elements into pairwise disjoint nonempty subsets). If `n \leq 0`,
    def bell_number(n, algorithm='dobinski') 
    221221      - ``'gap'`` -- Wrap GAP's ``Bell``
    222222      - ``'mpmath'`` -- Wrap mpmath's ``bell``
    223223
     224    .. WARNING::
     225
     226        When using the mpmath algorithm to compute bell numbers larger than 20,
     227        you will need to (temporarily) raise the precision using the option
     228        ``prec`` since it returns incorrect results due to low precision.
     229
    224230    Let `B_n` denote the `n`-th Bell number. Dobinski's formula is:
    225231
    226232    .. MATH::
    def bell_number(n, algorithm='dobinski') 
    320326        ...
    321327        TypeError: no conversion of this rational to integer
    322328
     329    Using mpmath on 30 produces incorrect results due to low precision::
     330
     331        sage: k = bell_number(30, 'mpmath'); k
     332        846749014511809388871680
     333        sage: k == bell_number(30)
     334        False
     335
     336    However we must (temporarily) set mpmath's precision higher by using the
     337    ``prec`` option. To guarentee correct results, one needs at least
     338    `\log_2(B_n)` bits::
     339
     340        sage: k = bell_number(30, 'mpmath', prec=30); k
     341        846749014511809332450147
     342        sage: k == bell_number(30)
     343        True
     344
    323345    TESTS::
    324346
    325347        sage: all([bell_number(n) == bell_number(n,'gap') for n in range(200, 220)])
    326348        True
    327         sage: all([bell_number(n,'gap') == bell_number(n,'mpmath') for n in range(1, 20)])
     349        sage: all([bell_number(n) == bell_number(n,'mpmath', prec=500) for n in range(200, 220)])
    328350        True
    329351
    330     .. WARNING::
    331 
    332         Using ``'mpmath'`` currently returns incorrect results for
    333         `n \geq 30`::
    334 
    335             sage: bell_number(30)   
    336             846749014511809332450147
    337             sage: bell_number(30, 'mpmath')
    338             846749014511809388871680
    339 
    340352    AUTHORS:
    341353
    342354    - Robert Gerbicz
    def bell_number(n, algorithm='dobinski') 
    349361    """
    350362    if algorithm == 'mpmath':
    351363        from sage.libs.mpmath.all import bell
     364        if 'prec' in options:
     365            from sage.libs.mpmath.all import mp
     366            old_prec = mp.dps
     367            mp.dps = options['prec']
     368            ret = ZZ(int(bell(n)))
     369            mp.dps = old_prec
     370            return ret
    352371        return ZZ(int(bell(n)))
    353372    if n < 200 or algorithm == 'gap':
    354373        return ZZ(gap.eval("Bell(%s)"%ZZ(n)))