# 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 from sage.misc.superseded import depreca from combinat_cython import _stirling_number2 ######### combinatorial sequences def bell_number(n, algorithm='dobinski'): def bell_number(n, algorithm='dobinski', **options): r""" Return the n-th Bell number (the number of ways to partition a set of n elements into pairwise disjoint nonempty subsets). If n \leq 0, def bell_number(n, algorithm='dobinski') - 'gap' -- Wrap GAP's Bell - 'mpmath' -- Wrap mpmath's bell .. WARNING:: When using the mpmath algorithm to compute bell numbers larger than 20, you will need to (temporarily) raise the precision using the option prec since it returns incorrect results due to low precision. Let B_n denote the n-th Bell number. Dobinski's formula is: .. MATH:: def bell_number(n, algorithm='dobinski') ... TypeError: no conversion of this rational to integer Using mpmath on 30 produces incorrect results due to low precision:: sage: k = bell_number(30, 'mpmath'); k 846749014511809388871680 sage: k == bell_number(30) False However we must (temporarily) set mpmath's precision higher by using the prec option. To guarentee correct results, one needs at least \log_2(B_n) bits:: sage: k = bell_number(30, 'mpmath', prec=30); k 846749014511809332450147 sage: k == bell_number(30) True TESTS:: sage: all([bell_number(n) == bell_number(n,'gap') for n in range(200, 220)]) True sage: all([bell_number(n,'gap') == bell_number(n,'mpmath') for n in range(1, 20)]) sage: all([bell_number(n) == bell_number(n,'mpmath', prec=500) for n in range(200, 220)]) True .. WARNING:: Using 'mpmath' currently returns incorrect results for n \geq 30:: sage: bell_number(30) 846749014511809332450147 sage: bell_number(30, 'mpmath') 846749014511809388871680 AUTHORS: - Robert Gerbicz def bell_number(n, algorithm='dobinski') """ if algorithm == 'mpmath': from sage.libs.mpmath.all import bell if 'prec' in options: from sage.libs.mpmath.all import mp old_prec = mp.dps mp.dps = options['prec'] ret = ZZ(int(bell(n))) mp.dps = old_prec return ret return ZZ(int(bell(n))) if n < 200 or algorithm == 'gap': return ZZ(gap.eval("Bell(%s)"%ZZ(n)))