"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "mpmath/libmp/gammazeta.py" between
mpmath-1.0.0.tar.gz and mpmath-1.1.0.tar.gz

About: mpmath is a Python library for arbitrary-precision floating-point arithmetic.

gammazeta.py  (mpmath-1.0.0):gammazeta.py  (mpmath-1.1.0)
skipping to change at line 353 skipping to change at line 353
mpf_twinprime = def_mpf_constant(twinprime_fixed) mpf_twinprime = def_mpf_constant(twinprime_fixed)
#-----------------------------------------------------------------------# #-----------------------------------------------------------------------#
# # # #
# Bernoulli numbers # # Bernoulli numbers #
# # # #
#-----------------------------------------------------------------------# #-----------------------------------------------------------------------#
MAX_BERNOULLI_CACHE = 3000 MAX_BERNOULLI_CACHE = 3000
""" r"""
Small Bernoulli numbers and factorials are used in numerous summations, Small Bernoulli numbers and factorials are used in numerous summations,
so it is critical for speed that sequential computation is fast and that so it is critical for speed that sequential computation is fast and that
values are cached up to a fairly high threshold. values are cached up to a fairly high threshold.
On the other hand, we also want to support fast computation of isolated On the other hand, we also want to support fast computation of isolated
large numbers. Currently, no such acceleration is provided for integer large numbers. Currently, no such acceleration is provided for integer
factorials (though it is for large floating-point factorials, which are factorials (though it is for large floating-point factorials, which are
computed via gamma if the precision is low enough). computed via gamma if the precision is low enough).
For sequential computation of Bernoulli numbers, we use Ramanujan's formula For sequential computation of Bernoulli numbers, we use Ramanujan's formula
skipping to change at line 794 skipping to change at line 794
t = mpc_sub(mpc_mul(logxpa, (reph, im), wp), (repa, im), wp) t = mpc_sub(mpc_mul(logxpa, (reph, im), wp), (repa, im), wp)
t = mpc_mul(mpc_exp(t, wp), s, prec, rounding) t = mpc_mul(mpc_exp(t, wp), s, prec, rounding)
return t return t
#-----------------------------------------------------------------------# #-----------------------------------------------------------------------#
# # # #
# Polygamma functions # # Polygamma functions #
# # # #
#-----------------------------------------------------------------------# #-----------------------------------------------------------------------#
""" r"""
For all polygamma (psi) functions, we use the Euler-Maclaurin summation For all polygamma (psi) functions, we use the Euler-Maclaurin summation
formula. It looks slightly different in the m = 0 and m > 0 cases. formula. It looks slightly different in the m = 0 and m > 0 cases.
For m = 0, we have For m = 0, we have
oo oo
___ B ___ B
(0) 1 \ 2 k -2 k (0) 1 \ 2 k -2 k
psi (z) ~ log z + --- - ) ------ z psi (z) ~ log z + --- - ) ------ z
2 z /___ (2 k)! 2 z /___ (2 k)!
k = 1 k = 1
skipping to change at line 882 skipping to change at line 882
Computation of the digamma function (psi function of order 0) Computation of the digamma function (psi function of order 0)
of a real argument. of a real argument.
""" """
sign, man, exp, bc = x sign, man, exp, bc = x
wp = prec + 10 wp = prec + 10
if not man: if not man:
if x == finf: return x if x == finf: return x
if x == fninf or x == fnan: return fnan if x == fninf or x == fnan: return fnan
if x == fzero or (exp >= 0 and sign): if x == fzero or (exp >= 0 and sign):
raise ValueError("polygamma pole") raise ValueError("polygamma pole")
# Near 0 -- fixed-point arithmetic becomes bad
if exp+bc < -5:
v = mpf_psi0(mpf_add(x, fone, prec, rnd), prec, rnd)
return mpf_sub(v, mpf_div(fone, x, wp, rnd), prec, rnd)
# Reflection formula # Reflection formula
if sign and exp+bc > 3: if sign and exp+bc > 3:
c, s = mpf_cos_sin_pi(x, wp) c, s = mpf_cos_sin_pi(x, wp)
q = mpf_mul(mpf_div(c, s, wp), mpf_pi(wp), wp) q = mpf_mul(mpf_div(c, s, wp), mpf_pi(wp), wp)
p = mpf_psi0(mpf_sub(fone, x, wp), wp) p = mpf_psi0(mpf_sub(fone, x, wp), wp)
return mpf_sub(p, q, prec, rnd) return mpf_sub(p, q, prec, rnd)
# The logarithmic term is accurate enough # The logarithmic term is accurate enough
if (not sign) and bc + exp > wp: if (not sign) and bc + exp > wp:
return mpf_log(mpf_sub(x, fone, wp), prec, rnd) return mpf_log(mpf_sub(x, fone, wp), prec, rnd)
# Initial recurrence to obtain a large enough x # Initial recurrence to obtain a large enough x
skipping to change at line 1054 skipping to change at line 1058
if not (m & 1): if not (m & 1):
v = mpf_neg(v[0]), mpf_neg(v[1]) v = mpf_neg(v[0]), mpf_neg(v[1])
return v return v
#-----------------------------------------------------------------------# #-----------------------------------------------------------------------#
# # # #
# Riemann zeta function # # Riemann zeta function #
# # # #
#-----------------------------------------------------------------------# #-----------------------------------------------------------------------#
""" r"""
We use zeta(s) = eta(s) / (1 - 2**(1-s)) and Borwein's approximation We use zeta(s) = eta(s) / (1 - 2**(1-s)) and Borwein's approximation
n-1 n-1
___ k ___ k
-1 \ (-1) (d_k - d_n) -1 \ (-1) (d_k - d_n)
eta(s) ~= ---- ) ------------------ eta(s) ~= ---- ) ------------------
d_n /___ s d_n /___ s
k = 0 (k + 1) k = 0 (k + 1)
where where
k k
 End of changes. 4 change blocks. 
3 lines changed or deleted 7 lines changed or added

Home  |  About  |  Features  |  All  |  Newest  |  Dox  |  Diffs  |  RSS Feeds  |  Screenshots  |  Comments  |  Imprint  |  Privacy  |  HTTP(S)