details: http://hg.sympy.org/sympy/rev/3a9219d09389 changeset: 1848:3a9219d09389 user: Fredrik Johansson <[EMAIL PROTECTED]> date: Mon Oct 20 17:19:52 2008 +0200 description: Improved calculation of Bernoulli numbers
Use mpmath to compute large Bernoulli numbers (bernoulli(10000) is now fast) Support bernoulli(10**10, evaluate=False).evalf() via mpmath Merge a recent mpmath bugfix without which bernoulli(10**10, evaluate=False).evalf() would give an inaccurate result diffs (168 lines): diff -r 8d7fc8ac251a -r 3a9219d09389 sympy/core/evalf.py --- a/sympy/core/evalf.py Wed Oct 22 22:39:07 2008 +0200 +++ b/sympy/core/evalf.py Mon Oct 20 17:19:52 2008 +0200 @@ -17,6 +17,7 @@ from sympy.mpmath.libmpf import MP_BASE, from_man_exp from sympy.mpmath.calculus import shanks_extrapolation, richardson_extrapolation +from sympy.mpmath.gammazeta import mpf_bernoulli import math @@ -616,6 +617,15 @@ # We still have undefined symbols raise NotImplementedError +def evalf_bernoulli(expr, prec, options): + arg = expr.args[0] + if not arg.is_Integer: + raise ValueError("Bernoulli number index must be an integer") + n = int(arg) + b = mpf_bernoulli(n, prec, round_nearest) + if b == fzero: + return None, None, None, None + return b, None, prec, None #----------------------------------------------------------------------------# # # @@ -943,6 +953,8 @@ C.Integral : evalf_integral, C.Sum : evalf_sum, C.Piecewise : evalf_piecewise, + + C.bernoulli : evalf_bernoulli, } def evalf(x, prec, options): diff -r 8d7fc8ac251a -r 3a9219d09389 sympy/functions/combinatorial/numbers.py --- a/sympy/functions/combinatorial/numbers.py Wed Oct 22 22:39:07 2008 +0200 +++ b/sympy/functions/combinatorial/numbers.py Mon Oct 20 17:19:52 2008 +0200 @@ -10,6 +10,8 @@ from sympy.core import * from sympy.core.function import Function from sympy.core.basic import S + +from sympy.mpmath import bernfrac def _product(a, b): p = 1 @@ -236,6 +238,10 @@ if n.is_odd: return S.Zero n = int(n) + # Use mpmath for enormous Bernoulli numbers + if n > 500: + p, q = bernfrac(n) + return Rational(int(p), q) case = n % 6 highest_cached = cls._highest[case] if n <= highest_cached: diff -r 8d7fc8ac251a -r 3a9219d09389 sympy/functions/combinatorial/tests/test_comb_numbers.py --- a/sympy/functions/combinatorial/tests/test_comb_numbers.py Wed Oct 22 22:39:07 2008 +0200 +++ b/sympy/functions/combinatorial/tests/test_comb_numbers.py Mon Oct 20 17:19:52 2008 +0200 @@ -20,6 +20,15 @@ assert bernoulli(1, x) == x-Rational(1,2) assert bernoulli(2, x) == x**2-x+Rational(1,6) assert bernoulli(3, x) == x**3 - (3*x**2)/2 + x/2 + + # Should be fast; computed with mpmath + b = bernoulli(1000) + assert b.p % 10**10 == 7950421099 + assert b.q == 342999030 + + b = bernoulli(10**6, evaluate=False).evalf() + assert str(b) == '-2.23799235765713e+4767529' + def test_fibonacci(): assert [fibonacci(n) for n in range(-3, 5)] == [2, -1, 1, 0, 1, 1, 2, 3] diff -r 8d7fc8ac251a -r 3a9219d09389 sympy/mpmath/gammazeta.py --- a/sympy/mpmath/gammazeta.py Wed Oct 22 22:39:07 2008 +0200 +++ b/sympy/mpmath/gammazeta.py Mon Oct 20 17:19:52 2008 +0200 @@ -24,7 +24,7 @@ from libmpf import (\ lshift, sqrt_fixed, fzero, fone, fnone, fhalf, ftwo, finf, fninf, fnan, - from_int, to_int, to_fixed, from_man_exp, + from_int, to_int, to_fixed, from_man_exp, from_rational, mpf_pos, mpf_neg, mpf_abs, mpf_add, mpf_sub, mpf_mul, mpf_mul_int, mpf_div, mpf_sqrt, mpf_pow_int, mpf_rdiv_int, @@ -423,6 +423,8 @@ lgn = math.log(n,2) return int(2.326 + 0.5*lgn + n*(lgn - 4.094)) +BERNOULLI_PREC_CUTOFF = bernoulli_size(MAX_BERNOULLI_CACHE) + def mpf_bernoulli(n, prec, rnd=None): """Computation of Bernoulli numbers (numerically)""" if n < 2: @@ -432,11 +434,19 @@ return fone if n == 1: return mpf_neg(fhalf) + # For odd n > 1, the Bernoulli numbers are zero if n & 1: return fzero - wp = prec + 30 + # If precision is extremely high, we can save time by computing + # the Bernoulli number at a lower precision that is sufficient to + # obtain the exact fraction, round to the exact fraction, and + # convert the fraction back to an mpf value at the original precision + if prec > BERNOULLI_PREC_CUTOFF and prec > bernoulli_size(n)*1.1 + 1000: + p, q = bernfrac(n) + return from_rational(p, q, prec, rnd or round_floor) if n > MAX_BERNOULLI_CACHE: return mpf_bernoulli_huge(n, prec, rnd) + wp = prec + 30 # Reuse nearby precisions wp += 32 - (prec & 31) cached = bernoulli_cache.get(wp) @@ -493,7 +503,7 @@ def mpf_bernoulli_huge(n, prec, rnd=None): wp = prec + 10 piprec = wp + int(math.log(n,2)) - v = mpf_gamma_int(n-1, wp) + v = mpf_gamma_int(n+1, wp) v = mpf_mul(v, mpf_zeta_int(n, wp), wp) v = mpf_mul(v, mpf_pow_int(mpf_pi(piprec), -n, wp)) v = mpf_shift(v, 1-n) @@ -648,11 +658,11 @@ # back here def mpf_gamma_int(n, prec, rounding=round_fast): if n < 1000: - return from_int(int_fac(n+1), prec, rounding) + return from_int(int_fac(n-1), prec, rounding) # XXX: choose the cutoff less arbitrarily size = int(n*math.log(n,2)) if prec > size/20.0: - return from_int(int_fac(n+1), prec, rounding) + return from_int(int_fac(n-1), prec, rounding) return mpf_gamma(from_int(n), prec, rounding) def mpf_factorial(x, prec, rounding=round_fast): diff -r 8d7fc8ac251a -r 3a9219d09389 sympy/mpmath/tests/test_gammazeta.py --- a/sympy/mpmath/tests/test_gammazeta.py Wed Oct 22 22:39:07 2008 +0200 +++ b/sympy/mpmath/tests/test_gammazeta.py Mon Oct 20 17:19:52 2008 +0200 @@ -38,15 +38,15 @@ assert bernoulli(6).ae(1./42) assert str(bernoulli(10)) == '0.0757575757575758' assert str(bernoulli(234)) == '7.62772793964344e+267' - assert str(bernoulli(10**5)) == '-5.82235253813873e+376745' - assert str(bernoulli(10**8+2)) == '1.19570351452842e+676752568' - assert str(bernoulli(10**100)) == '-2.58183325604736e+987675256497386331227838638980680030172857347883537824464410652557820800494271520411283004120790908423' + assert str(bernoulli(10**5)) == '-5.82229431461335e+376755' + assert str(bernoulli(10**8+2)) == '1.19570355039953e+676752584' + assert str(bernoulli(10**100)) == '-2.58183325604736e+987675256497386331227838638980680030172857347883537824464410652557820800494271520411283004120790908623' mp.dps = 50 assert str(bernoulli(10)) == '0.075757575757575757575757575757575757575757575757576' assert str(bernoulli(234)) == '7.6277279396434392486994969020496121553385863373331e+267' - assert str(bernoulli(10**5)) == '-5.8223525381387322109718142542038307938400075005058e+376745' - assert str(bernoulli(10**8+2)) == '1.1957035145284240522328803157892147712077025377042e+676752568' - assert str(bernoulli(10**100)) == '-2.5818332560473632073252488656039475548106223822913e+987675256497386331227838638980680030172857347883537824464410652557820800494271520411283004120790908423' + assert str(bernoulli(10**5)) == '-5.8222943146133508236497045360612887555320691004308e+376755' + assert str(bernoulli(10**8+2)) == '1.1957035503995297272263047884604346914602088317782e+676752584' + assert str(bernoulli(10**100)) == '-2.5818332560473632073252488656039475548106223822913e+987675256497386331227838638980680030172857347883537824464410652557820800494271520411283004120790908623' mp.dps = 1000 assert bernoulli(10).ae(mpf(5)/66) mp.dps = 15 --~--~---------~--~----~------------~-------~--~----~ You received this message because you are subscribed to the Google Groups "sympy-commits" group. To post to this group, send email to sympy-commits@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sympy-commits?hl=en -~----------~----~----~----~------~----~------~--~---