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
-~----------~----~----~----~------~----~------~--~---
