PEP 8 fixes:
E101 - mix of tabs and spaces W191 - use of tabs W291 - trailing whitespace W293 - trailing whitespace E111 - 4 spaces per indentation level E112 - 4 spaces per indentation level E113 - 4 spaces per indentation level W292 - no newline at end of file W391 - blank line at end of file see http://pep8.readthedocs.org/
This commit is contained in:
+1
-1
@@ -33,4 +33,4 @@ def main():
|
||||
print "Plotting..."
|
||||
plot(functions, xlim=xlim)
|
||||
|
||||
main()
|
||||
main()
|
||||
|
||||
@@ -2,4 +2,4 @@
|
||||
f = airyai
|
||||
f_diff = lambda z: airyai(z, derivative=1)
|
||||
f_int = lambda z: airyai(z, derivative=-1)
|
||||
plot([f, f_diff, f_int], [-10,5])
|
||||
plot([f, f_diff, f_int], [-10,5])
|
||||
|
||||
@@ -1,2 +1,2 @@
|
||||
# Airy function Ai(z) in the complex plane
|
||||
cplot(airyai, [-8,8], [-8,8], points=50000)
|
||||
cplot(airyai, [-8,8], [-8,8], points=50000)
|
||||
|
||||
@@ -3,4 +3,4 @@ f0 = lambda x: ber(0,x)
|
||||
f1 = lambda x: bei(0,x)
|
||||
f2 = lambda x: ber(2,x)
|
||||
f3 = lambda x: bei(2,x)
|
||||
plot([f0,f1,f2,f3],[0,10],[-10,10])
|
||||
plot([f0,f1,f2,f3],[0,10],[-10,10])
|
||||
|
||||
@@ -3,4 +3,4 @@ i0 = lambda x: besseli(0,x)
|
||||
i1 = lambda x: besseli(1,x)
|
||||
i2 = lambda x: besseli(2,x)
|
||||
i3 = lambda x: besseli(3,x)
|
||||
plot([i0,i1,i2,i3],[0,5],[0,5])
|
||||
plot([i0,i1,i2,i3],[0,5],[0,5])
|
||||
|
||||
@@ -1,2 +1,2 @@
|
||||
# Modified Bessel function I_n(z) in the complex plane
|
||||
cplot(lambda z: besseli(1,z), [-8,8], [-8,8], points=50000)
|
||||
cplot(lambda z: besseli(1,z), [-8,8], [-8,8], points=50000)
|
||||
|
||||
@@ -3,4 +3,4 @@ j0 = lambda x: besselj(0,x)
|
||||
j1 = lambda x: besselj(1,x)
|
||||
j2 = lambda x: besselj(2,x)
|
||||
j3 = lambda x: besselj(3,x)
|
||||
plot([j0,j1,j2,j3],[0,14])
|
||||
plot([j0,j1,j2,j3],[0,14])
|
||||
|
||||
@@ -1,2 +1,2 @@
|
||||
# Bessel function J_n(z) in the complex plane
|
||||
cplot(lambda z: besselj(1,z), [-8,8], [-8,8], points=50000)
|
||||
cplot(lambda z: besselj(1,z), [-8,8], [-8,8], points=50000)
|
||||
|
||||
@@ -3,4 +3,4 @@ k0 = lambda x: besselk(0,x)
|
||||
k1 = lambda x: besselk(1,x)
|
||||
k2 = lambda x: besselk(2,x)
|
||||
k3 = lambda x: besselk(3,x)
|
||||
plot([k0,k1,k2,k3],[0,8],[0,5])
|
||||
plot([k0,k1,k2,k3],[0,8],[0,5])
|
||||
|
||||
@@ -1,2 +1,2 @@
|
||||
# Modified Bessel function of 2nd kind K_n(z) in the complex plane
|
||||
cplot(lambda z: besselk(1,z), [-8,8], [-8,8], points=50000)
|
||||
cplot(lambda z: besselk(1,z), [-8,8], [-8,8], points=50000)
|
||||
|
||||
@@ -3,4 +3,4 @@ y0 = lambda x: bessely(0,x)
|
||||
y1 = lambda x: bessely(1,x)
|
||||
y2 = lambda x: bessely(2,x)
|
||||
y3 = lambda x: bessely(3,x)
|
||||
plot([y0,y1,y2,y3],[0,10],[-4,1])
|
||||
plot([y0,y1,y2,y3],[0,10],[-4,1])
|
||||
|
||||
@@ -1,2 +1,2 @@
|
||||
# Bessel function of 2nd kind Y_n(z) in the complex plane
|
||||
cplot(lambda z: bessely(1,z), [-8,8], [-8,8], points=50000)
|
||||
cplot(lambda z: bessely(1,z), [-8,8], [-8,8], points=50000)
|
||||
|
||||
@@ -2,4 +2,4 @@
|
||||
f = airybi
|
||||
f_diff = lambda z: airybi(z, derivative=1)
|
||||
f_int = lambda z: airybi(z, derivative=-1)
|
||||
plot([f, f_diff, f_int], [-10,2], [-1,2])
|
||||
plot([f, f_diff, f_int], [-10,2], [-1,2])
|
||||
|
||||
@@ -1,2 +1,2 @@
|
||||
# Airy function Bi(z) in the complex plane
|
||||
cplot(airybi, [-8,8], [-8,8], points=50000)
|
||||
cplot(airybi, [-8,8], [-8,8], points=50000)
|
||||
|
||||
@@ -4,4 +4,4 @@ f1 = lambda x: chebyt(1,x)
|
||||
f2 = lambda x: chebyt(2,x)
|
||||
f3 = lambda x: chebyt(3,x)
|
||||
f4 = lambda x: chebyt(4,x)
|
||||
plot([f0,f1,f2,f3,f4],[-1,1])
|
||||
plot([f0,f1,f2,f3,f4],[-1,1])
|
||||
|
||||
@@ -4,4 +4,4 @@ f1 = lambda x: chebyu(1,x)
|
||||
f2 = lambda x: chebyu(2,x)
|
||||
f3 = lambda x: chebyu(3,x)
|
||||
f4 = lambda x: chebyu(4,x)
|
||||
plot([f0,f1,f2,f3,f4],[-1,1])
|
||||
plot([f0,f1,f2,f3,f4],[-1,1])
|
||||
|
||||
@@ -4,4 +4,4 @@ F2 = lambda x: coulombf(0,1,x)
|
||||
F3 = lambda x: coulombf(0,5,x)
|
||||
F4 = lambda x: coulombf(0,10,x)
|
||||
F5 = lambda x: coulombf(0,x/2,x)
|
||||
plot([F1,F2,F3,F4,F5], [0,25], [-1.2,1.6])
|
||||
plot([F1,F2,F3,F4,F5], [0,25], [-1.2,1.6])
|
||||
|
||||
@@ -1,2 +1,2 @@
|
||||
# Regular Coulomb wave function in the complex plane
|
||||
cplot(lambda z: coulombf(1,1,z), points=50000)
|
||||
cplot(lambda z: coulombf(1,1,z), points=50000)
|
||||
|
||||
@@ -4,4 +4,4 @@ F2 = lambda x: coulombg(0,1,x)
|
||||
F3 = lambda x: coulombg(0,5,x)
|
||||
F4 = lambda x: coulombg(0,10,x)
|
||||
F5 = lambda x: coulombg(0,x/2,x)
|
||||
plot([F1,F2,F3,F4,F5], [0,30], [-2,2])
|
||||
plot([F1,F2,F3,F4,F5], [0,30], [-2,2])
|
||||
|
||||
@@ -1,2 +1,2 @@
|
||||
# Irregular Coulomb wave function in the complex plane
|
||||
cplot(lambda z: coulombg(1,1,z), points=50000)
|
||||
cplot(lambda z: coulombg(1,1,z), points=50000)
|
||||
|
||||
@@ -1,2 +1,2 @@
|
||||
# Complete elliptic integrals K(m) and E(m)
|
||||
plot([ellipk, ellipe], [-2,1], [0,3], points=600)
|
||||
plot([ellipk, ellipe], [-2,1], [0,3], points=600)
|
||||
|
||||
@@ -1,2 +1,2 @@
|
||||
# Scorer function Gi(x) and Gi'(x) on the real line
|
||||
plot([scorergi, diffun(scorergi)], [-10,10])
|
||||
plot([scorergi, diffun(scorergi)], [-10,10])
|
||||
|
||||
@@ -1,2 +1,2 @@
|
||||
# Scorer function Gi(z) in the complex plane
|
||||
cplot(scorergi, [-8,8], [-8,8], points=50000)
|
||||
cplot(scorergi, [-8,8], [-8,8], points=50000)
|
||||
|
||||
@@ -3,4 +3,4 @@ h0 = lambda x: hankel1(0,x)
|
||||
h1 = lambda x: hankel1(1,x)
|
||||
h2 = lambda x: hankel1(2,x)
|
||||
h3 = lambda x: hankel1(3,x)
|
||||
plot([h0,h1,h2,h3],[0,6],[-2,1])
|
||||
plot([h0,h1,h2,h3],[0,6],[-2,1])
|
||||
|
||||
@@ -1,2 +1,2 @@
|
||||
# Hankel function H1_n(z) in the complex plane
|
||||
cplot(lambda z: hankel1(1,z), [-8,8], [-8,8], points=50000)
|
||||
cplot(lambda z: hankel1(1,z), [-8,8], [-8,8], points=50000)
|
||||
|
||||
@@ -3,4 +3,4 @@ h0 = lambda x: hankel2(0,x)
|
||||
h1 = lambda x: hankel2(1,x)
|
||||
h2 = lambda x: hankel2(2,x)
|
||||
h3 = lambda x: hankel2(3,x)
|
||||
plot([h0,h1,h2,h3],[0,6],[-1,2])
|
||||
plot([h0,h1,h2,h3],[0,6],[-1,2])
|
||||
|
||||
@@ -1,2 +1,2 @@
|
||||
# Hankel function H2_n(z) in the complex plane
|
||||
cplot(lambda z: hankel2(1,z), [-8,8], [-8,8], points=50000)
|
||||
cplot(lambda z: hankel2(1,z), [-8,8], [-8,8], points=50000)
|
||||
|
||||
@@ -4,4 +4,4 @@ f1 = lambda x: hermite(1,x)
|
||||
f2 = lambda x: hermite(2,x)
|
||||
f3 = lambda x: hermite(3,x)
|
||||
f4 = lambda x: hermite(4,x)
|
||||
plot([f0,f1,f2,f3,f4],[-2,2],[-25,25])
|
||||
plot([f0,f1,f2,f3,f4],[-2,2],[-25,25])
|
||||
|
||||
@@ -1,2 +1,2 @@
|
||||
# Scorer function Hi(x) and Hi'(x) on the real line
|
||||
plot([scorerhi, diffun(scorerhi)], [-10,2], [0,2])
|
||||
plot([scorerhi, diffun(scorerhi)], [-10,2], [0,2])
|
||||
|
||||
@@ -1,2 +1,2 @@
|
||||
# Scorer function Hi(z) in the complex plane
|
||||
cplot(scorerhi, [-8,8], [-8,8], points=50000)
|
||||
cplot(scorerhi, [-8,8], [-8,8], points=50000)
|
||||
|
||||
@@ -3,4 +3,4 @@ f0 = lambda x: ker(0,x)
|
||||
f1 = lambda x: kei(0,x)
|
||||
f2 = lambda x: ker(2,x)
|
||||
f3 = lambda x: kei(2,x)
|
||||
plot([f0,f1,f2,f3],[0,5],[-1,4])
|
||||
plot([f0,f1,f2,f3],[0,5],[-1,4])
|
||||
|
||||
@@ -1,2 +1,2 @@
|
||||
# Klein J-function as function of the number-theoretic nome
|
||||
fp.cplot(lambda q: fp.kleinj(qbar=q), [-1,1], [-1,1], points=50000)
|
||||
fp.cplot(lambda q: fp.kleinj(qbar=q), [-1,1], [-1,1], points=50000)
|
||||
|
||||
@@ -1,2 +1,2 @@
|
||||
# Klein J-function as function of the half-period ratio
|
||||
fp.cplot(lambda t: fp.kleinj(tau=t), [-1,2], [0,1.5], points=50000)
|
||||
fp.cplot(lambda t: fp.kleinj(tau=t), [-1,2], [0,1.5], points=50000)
|
||||
|
||||
@@ -4,4 +4,4 @@ f1 = lambda x: laguerre(1,0,x)
|
||||
f2 = lambda x: laguerre(2,0,x)
|
||||
f3 = lambda x: laguerre(3,0,x)
|
||||
f4 = lambda x: laguerre(4,0,x)
|
||||
plot([f0,f1,f2,f3,f4],[0,10],[-10,10])
|
||||
plot([f0,f1,f2,f3,f4],[0,10],[-10,10])
|
||||
|
||||
@@ -1,2 +1,2 @@
|
||||
# Branches 0 and -1 of the Lambert W function
|
||||
plot([lambertw, lambda x: lambertw(x,-1)], [-2,2], [-5,2], points=2000)
|
||||
plot([lambertw, lambda x: lambertw(x,-1)], [-2,2], [-5,2], points=2000)
|
||||
|
||||
@@ -1,2 +1,2 @@
|
||||
# Principal branch of the Lambert W function W(z)
|
||||
cplot(lambertw, [-1,1], [-1,1], points=50000)
|
||||
cplot(lambertw, [-1,1], [-1,1], points=50000)
|
||||
|
||||
@@ -4,4 +4,4 @@ f1 = lambda x: legendre(1,x)
|
||||
f2 = lambda x: legendre(2,x)
|
||||
f3 = lambda x: legendre(3,x)
|
||||
f4 = lambda x: legendre(4,x)
|
||||
plot([f0,f1,f2,f3,f4],[-1,1])
|
||||
plot([f0,f1,f2,f3,f4],[-1,1])
|
||||
|
||||
@@ -3,4 +3,4 @@ f1 = lambda x: lommels1(-1,2.5,x)
|
||||
f2 = lambda x: lommels1(0,0.5,x)
|
||||
f3 = lambda x: lommels1(0,6,x)
|
||||
f4 = lambda x: lommels1(0.5,3,x)
|
||||
plot([f1,f2,f3,f4], [0,20])
|
||||
plot([f1,f2,f3,f4], [0,20])
|
||||
|
||||
@@ -3,4 +3,4 @@ f1 = lambda x: lommels2(-1,2.5,x)
|
||||
f2 = lambda x: lommels2(1.5,2,x)
|
||||
f3 = lambda x: lommels2(2.5,1,x)
|
||||
f4 = lambda x: lommels2(3.5,-0.5,x)
|
||||
plot([f1,f2,f3,f4], [0,8], [-8,8])
|
||||
plot([f1,f2,f3,f4], [0,8], [-8,8])
|
||||
|
||||
@@ -12,4 +12,4 @@ fp.splot(Y(4,0), [0,fp.pi], [0,2*fp.pi], points=300)
|
||||
# fp.splot(Y(4,0), [0,fp.pi], [0,2*fp.pi], points=300)
|
||||
# fp.splot(Y(4,1), [0,fp.pi], [0,2*fp.pi], points=300)
|
||||
# fp.splot(Y(4,2), [0,fp.pi], [0,2*fp.pi], points=300)
|
||||
# fp.splot(Y(4,3), [0,fp.pi], [0,2*fp.pi], points=300)
|
||||
# fp.splot(Y(4,3), [0,fp.pi], [0,2*fp.pi], points=300)
|
||||
|
||||
@@ -8,4 +8,4 @@ def Y(l,m):
|
||||
return [x,y,z]
|
||||
return g
|
||||
|
||||
fp.splot(Y(4,1), [0,fp.pi], [0,2*fp.pi], points=300)
|
||||
fp.splot(Y(4,1), [0,fp.pi], [0,2*fp.pi], points=300)
|
||||
|
||||
@@ -8,4 +8,4 @@ def Y(l,m):
|
||||
return [x,y,z]
|
||||
return g
|
||||
|
||||
fp.splot(Y(4,2), [0,fp.pi], [0,2*fp.pi], points=300)
|
||||
fp.splot(Y(4,2), [0,fp.pi], [0,2*fp.pi], points=300)
|
||||
|
||||
@@ -8,4 +8,4 @@ def Y(l,m):
|
||||
return [x,y,z]
|
||||
return g
|
||||
|
||||
fp.splot(Y(4,3), [0,fp.pi], [0,2*fp.pi], points=300)
|
||||
fp.splot(Y(4,3), [0,fp.pi], [0,2*fp.pi], points=300)
|
||||
|
||||
@@ -458,4 +458,3 @@ def doctests(filter=[]):
|
||||
|
||||
if __name__ == '__main__':
|
||||
doctests()
|
||||
|
||||
|
||||
@@ -717,9 +717,9 @@ class levin_class:
|
||||
return value,err
|
||||
|
||||
def levin(ctx, method = "levin", variant = "u"):
|
||||
L = levin_class(method = method, variant = variant)
|
||||
L.ctx = ctx
|
||||
return L
|
||||
L = levin_class(method = method, variant = variant)
|
||||
L.ctx = ctx
|
||||
return L
|
||||
|
||||
levin.__doc__ = levin_class.__doc__
|
||||
defun(levin)
|
||||
@@ -1217,7 +1217,7 @@ def adaptive_extrapolation(ctx, update, emfun, kwargs):
|
||||
if verbose:
|
||||
print("%s error: %s" % (L.name, ctx.nstr(lerror)))
|
||||
if lerror <= tol:
|
||||
return est
|
||||
return est
|
||||
if lerror < error:
|
||||
error = lerror
|
||||
best = est
|
||||
|
||||
@@ -6776,7 +6776,7 @@ Evaluation is supported for real and complex arguments::
|
||||
-0.7363054628673177346778998
|
||||
>>> siegelz(3+4j)
|
||||
(-0.1852895764366314976003936 - 0.2773099198055652246992479j)
|
||||
|
||||
|
||||
The first four derivatives are supported, using the
|
||||
optional *derivative* keyword argument::
|
||||
|
||||
@@ -6784,7 +6784,7 @@ optional *derivative* keyword argument::
|
||||
56.89689348495089294249178
|
||||
>>> diff(siegelz, 1234567, n=3)
|
||||
56.89689348495089294249178
|
||||
|
||||
|
||||
|
||||
The Z-function has a Maclaurin expansion::
|
||||
|
||||
@@ -9933,7 +9933,7 @@ of each other::
|
||||
... for k in range(5):
|
||||
... A[n,k] = stirling1(n,k)
|
||||
... B[n,k] = stirling2(n,k)
|
||||
...
|
||||
...
|
||||
>>> A * B
|
||||
[1.0 0.0 0.0 0.0 0.0]
|
||||
[0.0 1.0 0.0 0.0 0.0]
|
||||
|
||||
@@ -47,7 +47,7 @@ def besselj(ctx, n, z, derivative=0, **kwargs):
|
||||
return T
|
||||
v = ctx.hypercomb(h, [n,d], **kwargs)
|
||||
else:
|
||||
# Fast case: J_n(x), n int, appropriate magnitude for fixed-point calculation
|
||||
# Fast case: J_n(x), n int, appropriate magnitude for fixed-point calculation
|
||||
if (not derivative) and n_isint and abs(M) < 10 and abs(n) < 20:
|
||||
try:
|
||||
return ctx._besselj(n, z)
|
||||
|
||||
@@ -687,7 +687,7 @@ def elliprc(ctx, x, y, pv=True):
|
||||
|
||||
.. math ::
|
||||
|
||||
R_C(x,y) =
|
||||
R_C(x,y) =
|
||||
\begin{cases}
|
||||
\dfrac{1}{\sqrt{y-x}}
|
||||
\cos^{-1}\left(\sqrt{\dfrac{x}{y}}\right), & x < y \\
|
||||
|
||||
@@ -826,9 +826,9 @@ def _hyp1f2(ctx, a_s, b_s, z, **kwargs):
|
||||
# Asymptotic series is in terms of 3F0
|
||||
can_use_asymptotic = (not kwargs.get('force_series')) and \
|
||||
(ctx.mag(absz) > 19) and \
|
||||
(ctx.sqrt(absz) > 1.5*orig) #and \
|
||||
#ctx._hyp_check_convergence([a1, a1-b1+1, a1-b2+1], [],
|
||||
# 1/absz, orig+40+asymp_extraprec)
|
||||
(ctx.sqrt(absz) > 1.5*orig) # and \
|
||||
# ctx._hyp_check_convergence([a1, a1-b1+1, a1-b2+1], [],
|
||||
# 1/absz, orig+40+asymp_extraprec)
|
||||
|
||||
# TODO: much of the following could be shared with 2F3 instead of
|
||||
# copypasted
|
||||
|
||||
@@ -208,8 +208,8 @@ def qhyper(ctx, a_s, b_s, q, z, **kwargs):
|
||||
|
||||
.. math ::
|
||||
|
||||
\,_r\phi_s \left[\begin{matrix}
|
||||
a_1 & a_2 & \ldots & a_r \\
|
||||
\,_r\phi_s \left[\begin{matrix}
|
||||
a_1 & a_2 & \ldots & a_r \\
|
||||
b_1 & b_2 & \ldots & b_s
|
||||
\end{matrix} ; q,z \right] =
|
||||
\sum_{n=0}^\infty
|
||||
|
||||
+14
-17
@@ -51,7 +51,7 @@ formula can not compute to the wanted precision.
|
||||
import math
|
||||
|
||||
class RSCache(object):
|
||||
def __init__(ctx):
|
||||
def __init__(ctx):
|
||||
ctx._rs_cache = [0, 10, {}, {}]
|
||||
|
||||
from .functions import defun
|
||||
@@ -105,8 +105,7 @@ def _coef(ctx, J, eps):
|
||||
"""
|
||||
|
||||
newJ = J+2 # compute more coefficients that are needed
|
||||
neweps6 = eps/2. # compute with a slight more precision
|
||||
# that are needed
|
||||
neweps6 = eps/2. # compute with a slight more precision that are needed
|
||||
|
||||
# PREPARATION FOR THE COMPUTATION OF V(N) AND W(N)
|
||||
# See II Section 3.16
|
||||
@@ -412,9 +411,8 @@ def Rzeta_simul(ctx, s, der=0):
|
||||
cc = {}
|
||||
cont = {}
|
||||
cont, pipowers = coef(ctx, J, eps6)
|
||||
cc=cont.copy() # we need a copy since we have
|
||||
# to change his values.
|
||||
Fp={} # this is the adequate locus of this
|
||||
cc=cont.copy() # we need a copy since we have to change his values.
|
||||
Fp={} # this is the adequate locus of this
|
||||
for n in range(M, 3*L-2):
|
||||
Fp[n] = 0
|
||||
Fp={}
|
||||
@@ -560,11 +558,11 @@ def Rzeta_simul(ctx, s, der=0):
|
||||
xfortcoef[mu,k,ell]=xfortcoef[mu,k,ell]/((2*ctx.j)**ell)
|
||||
|
||||
def trunc_a(t):
|
||||
wp = ctx.prec
|
||||
ctx.prec = wp + 2
|
||||
aa = ctx.sqrt(t/(2*ctx.pi))
|
||||
ctx.prec = wp
|
||||
return aa
|
||||
wp = ctx.prec
|
||||
ctx.prec = wp + 2
|
||||
aa = ctx.sqrt(t/(2*ctx.pi))
|
||||
ctx.prec = wp
|
||||
return aa
|
||||
|
||||
# computing the tcoef[k,ell]
|
||||
xtcoef={}
|
||||
@@ -1012,11 +1010,11 @@ def Rzeta_set(ctx, s, derivatives=[0]):
|
||||
fortcoef[mu,k,ell]=fortcoef[mu,k,ell]/((2*ctx.j)**ell)
|
||||
|
||||
def trunc_a(t):
|
||||
wp = ctx.prec
|
||||
ctx.prec = wp + 2
|
||||
aa = ctx.sqrt(t/(2*ctx.pi))
|
||||
ctx.prec = wp
|
||||
return aa
|
||||
wp = ctx.prec
|
||||
ctx.prec = wp + 2
|
||||
aa = ctx.sqrt(t/(2*ctx.pi))
|
||||
ctx.prec = wp
|
||||
return aa
|
||||
|
||||
# computing the tcoef[chi,k,ell]
|
||||
tcoef={}
|
||||
@@ -1403,4 +1401,3 @@ def rs_z(ctx, w, derivative=0):
|
||||
return z_half(ctx, w, derivative)
|
||||
else:
|
||||
return z_offline(ctx, w, derivative)
|
||||
|
||||
|
||||
@@ -157,7 +157,6 @@ def siegelz(ctx, t, **kwargs):
|
||||
if d > 4:
|
||||
h = lambda x: ctx.siegelz(x, derivative=4)
|
||||
return ctx.diff(h, t, n=d-4)
|
||||
|
||||
|
||||
|
||||
_zeta_zeros = [
|
||||
@@ -883,7 +882,7 @@ def secondzeta_exp_term(ctx, s, a):
|
||||
mg = abs(term)
|
||||
v = a**(0.5*s)*totsum/ctx.gamma(0.5*s)
|
||||
return v
|
||||
|
||||
|
||||
def secondzeta_singular_term(ctx, s, a, **kwargs):
|
||||
factor = a**(0.5*(s-1))/(4*ctx.sqrt(ctx.pi)*ctx.gamma(0.5*s))
|
||||
extraprec = ctx.mag(factor)
|
||||
@@ -903,7 +902,7 @@ def secondzeta_singular_term(ctx, s, a, **kwargs):
|
||||
term = f(n)
|
||||
totsum += term
|
||||
n +=1
|
||||
term = f(n)
|
||||
term = f(n)
|
||||
mg1 = mg2
|
||||
mg2 = abs(term)
|
||||
totsum += term
|
||||
|
||||
@@ -141,17 +141,17 @@ def separate_my_zero(ctx, my_zero_number, zero_number_block, T, V, prec):
|
||||
for k in range(1,len(V)):
|
||||
v1 = V[k]
|
||||
if v0*v1 < 0:
|
||||
variations +=1
|
||||
if variations == my_zero_number:
|
||||
k0 = k
|
||||
leftv = v0
|
||||
rightv = v1
|
||||
variations +=1
|
||||
if variations == my_zero_number:
|
||||
k0 = k
|
||||
leftv = v0
|
||||
rightv = v1
|
||||
v0 = v1
|
||||
t1 = T[k0]
|
||||
t0 = T[k0-1]
|
||||
ctx.prec = prec
|
||||
wpz = wpzeros(my_zero_number*ctx.log(my_zero_number))
|
||||
|
||||
|
||||
guard = 4*ctx.mag(my_zero_number)
|
||||
precs = [ctx.prec+4]
|
||||
index=0
|
||||
@@ -160,7 +160,7 @@ def separate_my_zero(ctx, my_zero_number, zero_number_block, T, V, prec):
|
||||
precs = [precs[0] // 2 +3+2*index] + precs
|
||||
ctx.prec = precs[0] + guard
|
||||
r = ctx.findroot(lambda x:ctx.siegelz(x), (t0,t1), solver ='illinois', verbose=False)
|
||||
#print "first step at", ctx.dps, "digits"
|
||||
#print "first step at", ctx.dps, "digits"
|
||||
z=ctx.mpc(0.5,r)
|
||||
for prec in precs[1:]:
|
||||
ctx.prec = prec + guard
|
||||
@@ -190,7 +190,7 @@ def sure_number_block(ctx, n):
|
||||
def compute_triple_tvb(ctx, n):
|
||||
t = ctx.grampoint(n)
|
||||
v = ctx._fp.siegelz(t)
|
||||
if ctx.mag(abs(v))<ctx.mag(t)-45:
|
||||
if ctx.mag(abs(v))<ctx.mag(t)-45:
|
||||
v = ctx.siegelz(t)
|
||||
b = v*(-1)**n
|
||||
return t,v,b
|
||||
@@ -384,7 +384,7 @@ def zetazero(ctx, n, info=False, round=True):
|
||||
(0.5 + 727690.9069482075392389420041147142092708393819935j)
|
||||
>>> chop(zeta(_)/_)
|
||||
0.0
|
||||
|
||||
|
||||
with *info=True*, :func:`~mpmath.zetazero` gives additional information::
|
||||
|
||||
>>> mp.dps = 15
|
||||
@@ -395,7 +395,7 @@ def zetazero(ctx, n, info=False, round=True):
|
||||
it is the 6-th zero between them. Finally (01311110) is the pattern
|
||||
of zeros in this interval. The numbers indicate the number of zeros
|
||||
in each Gram interval (Rosser blocks between parenthesis). In this case
|
||||
there is only one Rosser block of length nine.
|
||||
there is only one Rosser block of length nine.
|
||||
"""
|
||||
n = int(n)
|
||||
if n < 0:
|
||||
|
||||
@@ -2379,4 +2379,3 @@ def mpf_gamma_int(n, prec, rnd=round_fast):
|
||||
if n < SMALL_FACTORIAL_CACHE_SIZE:
|
||||
return mpf_pos(small_factorial_cache[n-1], prec, rnd)
|
||||
return mpf_gamma(from_int(n), prec, rnd)
|
||||
|
||||
|
||||
@@ -910,15 +910,15 @@ def mpf_atan2(y, x, prec, rnd=round_fast):
|
||||
return mpf_pos(tquo, prec, rnd)
|
||||
|
||||
def mpf_asin(x, prec, rnd=round_fast):
|
||||
sign, man, exp, bc = x
|
||||
if bc+exp > 0 and x not in (fone, fnone):
|
||||
raise ComplexResult("asin(x) is real only for -1 <= x <= 1")
|
||||
# asin(x) = 2*atan(x/(1+sqrt(1-x**2)))
|
||||
wp = prec + 15
|
||||
a = mpf_mul(x, x)
|
||||
b = mpf_add(fone, mpf_sqrt(mpf_sub(fone, a, wp), wp), wp)
|
||||
c = mpf_div(x, b, wp)
|
||||
return mpf_shift(mpf_atan(c, prec, rnd), 1)
|
||||
sign, man, exp, bc = x
|
||||
if bc+exp > 0 and x not in (fone, fnone):
|
||||
raise ComplexResult("asin(x) is real only for -1 <= x <= 1")
|
||||
# asin(x) = 2*atan(x/(1+sqrt(1-x**2)))
|
||||
wp = prec + 15
|
||||
a = mpf_mul(x, x)
|
||||
b = mpf_add(fone, mpf_sqrt(mpf_sub(fone, a, wp), wp), wp)
|
||||
c = mpf_div(x, b, wp)
|
||||
return mpf_shift(mpf_atan(c, prec, rnd), 1)
|
||||
|
||||
def mpf_acos(x, prec, rnd=round_fast):
|
||||
# acos(x) = 2*atan(sqrt(1-x**2)/(1+x))
|
||||
|
||||
+1
-1
@@ -194,7 +194,7 @@ _lanczos_g = 7
|
||||
_lanczos_p = (0.99999999999980993, 676.5203681218851, -1259.1392167224028,
|
||||
771.32342877765313, -176.61502916214059, 12.507343278686905,
|
||||
-0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7)
|
||||
|
||||
|
||||
def _gamma_real(x):
|
||||
_intx = int(x)
|
||||
if _intx == x:
|
||||
|
||||
@@ -55,7 +55,7 @@ class MatrixCalculusMethods(object):
|
||||
Computes the matrix exponential of a square matrix `A`, which is defined
|
||||
by the power series
|
||||
|
||||
.. math ::
|
||||
.. math ::
|
||||
|
||||
\exp(A) = I + A + \frac{A^2}{2!} + \frac{A^3}{3!} + \ldots
|
||||
|
||||
|
||||
+15
-15
@@ -49,7 +49,7 @@ def hessenberg_reduce_0(ctx, A, T):
|
||||
Q' A Q = H and Q' Q = Q Q' = 1
|
||||
|
||||
where H is an upper Hessenberg matrix, meaning that it only contains zeros
|
||||
below the first subdiagonal. Here ' denotes the hermitian transpose (i.e.
|
||||
below the first subdiagonal. Here ' denotes the hermitian transpose (i.e.
|
||||
transposition and conjugation).
|
||||
|
||||
parameters:
|
||||
@@ -192,7 +192,7 @@ def hessenberg(ctx, A, overwrite_a = False):
|
||||
|
||||
Q' A Q = H and Q' Q = Q Q' = 1
|
||||
|
||||
where H is an upper right Hessenberg matrix. Here ' denotes the hermitian
|
||||
where H is an upper right Hessenberg matrix. Here ' denotes the hermitian
|
||||
transpose (i.e. transposition and conjugation).
|
||||
|
||||
input:
|
||||
@@ -380,7 +380,7 @@ def hessenberg_qr(ctx, A, Q):
|
||||
|
||||
Q' A Q = R and Q' Q = Q Q' = 1
|
||||
|
||||
where R is an upper right triangular matrix. Here ' denotes the hermitian
|
||||
where R is an upper right triangular matrix. Here ' denotes the hermitian
|
||||
transpose (i.e. transposition and conjugation).
|
||||
|
||||
parameters:
|
||||
@@ -593,12 +593,12 @@ def eig_tr_r(ctx, A):
|
||||
rmax = max(rmax, abs(r))
|
||||
if rmax > simin:
|
||||
for k in xrange(j, i+1):
|
||||
ER[k,i] /= rmax
|
||||
ER[k,i] /= rmax
|
||||
rmax = 1
|
||||
|
||||
if rmax != 1:
|
||||
for k in xrange(0, i + 1):
|
||||
ER[k,i] /= rmax
|
||||
ER[k,i] /= rmax
|
||||
|
||||
return ER
|
||||
|
||||
@@ -651,12 +651,12 @@ def eig_tr_l(ctx, A):
|
||||
rmax = max(rmax, abs(r))
|
||||
if rmax > simin:
|
||||
for k in xrange(i, j + 1):
|
||||
EL[i,k] /= rmax
|
||||
EL[i,k] /= rmax
|
||||
rmax = 1
|
||||
|
||||
if rmax != 1:
|
||||
for k in xrange(i, n):
|
||||
EL[i,k] /= rmax
|
||||
EL[i,k] /= rmax
|
||||
|
||||
return EL
|
||||
|
||||
@@ -816,14 +816,14 @@ def eig_sort(ctx, E, EL = False, ER = False, f = "real"):
|
||||
"""
|
||||
|
||||
if isinstance(f, str):
|
||||
if f == "real":
|
||||
f = ctx.re
|
||||
elif f == "imag":
|
||||
f = ctx.im
|
||||
elif cmp == "abs":
|
||||
f = abs
|
||||
else:
|
||||
raise RuntimeError("unknown function %s" % f)
|
||||
if f == "real":
|
||||
f = ctx.re
|
||||
elif f == "imag":
|
||||
f = ctx.im
|
||||
elif cmp == "abs":
|
||||
f = abs
|
||||
else:
|
||||
raise RuntimeError("unknown function %s" % f)
|
||||
|
||||
n = len(E)
|
||||
|
||||
|
||||
@@ -79,8 +79,7 @@ def r_sy_tridiag(ctx, A, D, E, calc_ev = True):
|
||||
|
||||
n = A.rows
|
||||
for i in xrange(n - 1, 0, -1):
|
||||
|
||||
# scale the vector
|
||||
# scale the vector
|
||||
|
||||
scale = 0
|
||||
for k in xrange(0, i):
|
||||
@@ -90,14 +89,14 @@ def r_sy_tridiag(ctx, A, D, E, calc_ev = True):
|
||||
if scale != 0:
|
||||
scale_inv = 1/scale
|
||||
|
||||
# sadly there are floating point numbers not equal to zero whose reciprocal is infinity
|
||||
# sadly there are floating point numbers not equal to zero whose reciprocal is infinity
|
||||
|
||||
if i == 1 or scale == 0 or ctx.isinf(scale_inv):
|
||||
E[i] = A[i-1,i] # nothing to do
|
||||
D[i] = 0
|
||||
continue
|
||||
|
||||
# calculate parameters for housholder transformation
|
||||
# calculate parameters for housholder transformation
|
||||
|
||||
H = 0
|
||||
for k in xrange(0, i):
|
||||
@@ -113,7 +112,7 @@ def r_sy_tridiag(ctx, A, D, E, calc_ev = True):
|
||||
A[i-1,i] = F - G
|
||||
F = 0
|
||||
|
||||
# apply housholder transformation
|
||||
# apply housholder transformation
|
||||
|
||||
for j in xrange(0, i):
|
||||
if calc_ev:
|
||||
@@ -206,7 +205,7 @@ def c_he_tridiag_0(ctx, A, D, E, T):
|
||||
T[n-1] = 1
|
||||
for i in xrange(n - 1, 0, -1):
|
||||
|
||||
# scale the vector
|
||||
# scale the vector
|
||||
|
||||
scale = 0
|
||||
for k in xrange(0, i):
|
||||
@@ -216,7 +215,7 @@ def c_he_tridiag_0(ctx, A, D, E, T):
|
||||
if scale != 0:
|
||||
scale_inv = 1 / scale
|
||||
|
||||
# sadly there are floating point numbers not equal to zero whose reciprocal is infinity
|
||||
# sadly there are floating point numbers not equal to zero whose reciprocal is infinity
|
||||
|
||||
if scale == 0 or ctx.isinf(scale_inv):
|
||||
E[i] = 0
|
||||
@@ -235,7 +234,7 @@ def c_he_tridiag_0(ctx, A, D, E, T):
|
||||
T[i-1] = T[i]
|
||||
continue
|
||||
|
||||
# calculate parameters for housholder transformation
|
||||
# calculate parameters for housholder transformation
|
||||
|
||||
H = 0
|
||||
for k in xrange(0, i):
|
||||
@@ -258,7 +257,7 @@ def c_he_tridiag_0(ctx, A, D, E, T):
|
||||
A[i-1,i] += G
|
||||
F = 0
|
||||
|
||||
# apply housholder transformation
|
||||
# apply housholder transformation
|
||||
|
||||
for j in xrange(0, i):
|
||||
A[i,j] = A[j,i] / H
|
||||
@@ -281,8 +280,8 @@ def c_he_tridiag_0(ctx, A, D, E, T):
|
||||
|
||||
for k in xrange(0, j + 1):
|
||||
A[k,j] -= ctx.conj(F) * T[k] + ctx.conj(G) * A[k,i]
|
||||
# as we use the lower left part for storage
|
||||
# we have to use the transpose of the normal formula
|
||||
# as we use the lower left part for storage
|
||||
# we have to use the transpose of the normal formula
|
||||
|
||||
T[i-1] = TZ
|
||||
D[i] = H
|
||||
@@ -887,7 +886,7 @@ def gauss_quadrature(ctx, n, qtype = "legendre", alpha = 0, beta = 0):
|
||||
tridiag_eigen(ctx, d, e, z)
|
||||
|
||||
for i in xrange(len(z)):
|
||||
z[i] *= z[i]
|
||||
z[i] *= z[i]
|
||||
|
||||
z = z.transpose()
|
||||
return (d, w * z)
|
||||
|
||||
+19
-19
@@ -639,7 +639,7 @@ class LinearAlgebraMethods(object):
|
||||
[0.0 1.0]
|
||||
|
||||
"""
|
||||
|
||||
|
||||
# check values before continuing
|
||||
assert isinstance(A, ctx.matrix)
|
||||
m = A.rows
|
||||
@@ -654,22 +654,22 @@ class LinearAlgebraMethods(object):
|
||||
# temporarily increase the precision and initialize
|
||||
with ctx.extradps(edps):
|
||||
tau = ctx.matrix(n,1)
|
||||
A = A.copy()
|
||||
A = A.copy()
|
||||
|
||||
# ---------------
|
||||
# FACTOR MATRIX A
|
||||
# FACTOR MATRIX A
|
||||
# ---------------
|
||||
if cmplx:
|
||||
one = ctx.mpc('1.0', '0.0')
|
||||
zero = ctx.mpc('0.0', '0.0')
|
||||
rzero = ctx.mpf('0.0')
|
||||
|
||||
# main loop to factor A (complex)
|
||||
|
||||
# main loop to factor A (complex)
|
||||
for j in xrange(0, n):
|
||||
alpha = A[j,j]
|
||||
alphr = ctx.re(alpha)
|
||||
alphi = ctx.im(alpha)
|
||||
|
||||
|
||||
if (m-j) >= 2:
|
||||
xnorm = ctx.fsum( A[i,j]*ctx.conj(A[i,j]) for i in xrange(j+1, m) )
|
||||
xnorm = ctx.re( ctx.sqrt(xnorm) )
|
||||
@@ -688,7 +688,7 @@ class LinearAlgebraMethods(object):
|
||||
tau[j] = ctx.mpc( (beta - alphr) / beta, -alphi / beta )
|
||||
t = -ctx.conj(tau[j])
|
||||
za = one / (alpha - beta)
|
||||
|
||||
|
||||
for i in xrange(j+1, m):
|
||||
A[i,j] *= za
|
||||
|
||||
@@ -724,29 +724,29 @@ class LinearAlgebraMethods(object):
|
||||
beta = ctx.sqrt(alpha**2 + xnorm**2)
|
||||
else:
|
||||
beta = -ctx.sqrt(alpha**2 + xnorm**2)
|
||||
|
||||
|
||||
tau[j] = (beta - alpha) / beta
|
||||
t = -tau[j]
|
||||
da = one / (alpha - beta)
|
||||
|
||||
|
||||
for i in xrange(j+1, m):
|
||||
A[i,j] *= da
|
||||
|
||||
|
||||
A[j,j] = one
|
||||
for k in xrange(j+1, n):
|
||||
y = ctx.fsum( A[i,j] * A[i,k] for i in xrange(j, m) )
|
||||
temp = t * y
|
||||
for i in xrange(j,m):
|
||||
A[i,k] += A[i,j] * temp
|
||||
|
||||
A[i,k] += A[i,j] * temp
|
||||
|
||||
A[j,j] = beta
|
||||
|
||||
# return factorization in same internal format as LAPACK
|
||||
if (mode == 'raw') or (mode == 'RAW'):
|
||||
return A, tau
|
||||
|
||||
|
||||
# ----------------------------------
|
||||
# FORM Q USING BACKWARD ACCUMULATION
|
||||
# FORM Q USING BACKWARD ACCUMULATION
|
||||
# ----------------------------------
|
||||
|
||||
# form R before the values are overwritten
|
||||
@@ -766,7 +766,7 @@ class LinearAlgebraMethods(object):
|
||||
A[j,j] = one
|
||||
for i in xrange(0, j):
|
||||
A[i,j] = zero
|
||||
|
||||
|
||||
# main loop to form Q
|
||||
for j in xrange(n-1, -1, -1):
|
||||
t = -tau[j]
|
||||
@@ -774,18 +774,18 @@ class LinearAlgebraMethods(object):
|
||||
|
||||
for k in xrange(j+1, p):
|
||||
if cmplx:
|
||||
y = ctx.fsum(A[i,j] * ctx.conj(A[i,k]) for i in xrange(j+1, m))
|
||||
y = ctx.fsum(A[i,j] * ctx.conj(A[i,k]) for i in xrange(j+1, m))
|
||||
temp = t * ctx.conj(y)
|
||||
else:
|
||||
y = ctx.fsum(A[i,j] * A[i,k] for i in xrange(j+1, m))
|
||||
temp = t * y
|
||||
temp = t * y
|
||||
A[j,k] = temp
|
||||
for i in xrange(j+1, m):
|
||||
A[i,k] += A[i,j] * temp
|
||||
|
||||
|
||||
for i in xrange(j+1, m):
|
||||
A[i, j] *= t
|
||||
|
||||
|
||||
return A, R[0:p,0:n]
|
||||
|
||||
# ------------------
|
||||
|
||||
+11
-12
@@ -406,7 +406,7 @@ class _matrix(object):
|
||||
|
||||
def __get_element(self, key):
|
||||
'''
|
||||
Fast extraction of the i,j element from the matrix
|
||||
Fast extraction of the i,j element from the matrix
|
||||
This function is for private use only because is unsafe:
|
||||
1. Does not check on the value of key it expects key to be a integer tuple (i,j)
|
||||
2. Does not check bounds
|
||||
@@ -431,11 +431,11 @@ class _matrix(object):
|
||||
|
||||
|
||||
def __getitem__(self, key):
|
||||
'''
|
||||
'''
|
||||
Getitem function for mp matrix class with slice index enabled
|
||||
it allows the following assingments
|
||||
scalar to a slice of the matrix
|
||||
B = A[:,2:6]
|
||||
B = A[:,2:6]
|
||||
'''
|
||||
# Convert vector to matrix indexing
|
||||
if isinstance(key, int) or isinstance(key,slice):
|
||||
@@ -446,7 +446,7 @@ class _matrix(object):
|
||||
key = (key, 0)
|
||||
else:
|
||||
raise IndexError('insufficient indices for matrix')
|
||||
|
||||
|
||||
if isinstance(key[0],slice) or isinstance(key[1],slice):
|
||||
|
||||
#Rows
|
||||
@@ -461,7 +461,7 @@ class _matrix(object):
|
||||
else:
|
||||
# Single row
|
||||
rows = [key[0]]
|
||||
|
||||
|
||||
# Columns
|
||||
if isinstance(key[1],slice):
|
||||
# Check bounds
|
||||
@@ -471,21 +471,21 @@ class _matrix(object):
|
||||
columns = xrange(*key[1].indices(self.__cols))
|
||||
else:
|
||||
raise IndexError('Column index out of bounds')
|
||||
|
||||
|
||||
else:
|
||||
# Single column
|
||||
columns = [key[1]]
|
||||
|
||||
# Create matrix slice
|
||||
m = self.ctx.matrix(len(rows),len(columns))
|
||||
m = self.ctx.matrix(len(rows),len(columns))
|
||||
|
||||
# Assign elements to the output matrix
|
||||
for i,x in enumerate(rows):
|
||||
for j,y in enumerate(columns):
|
||||
m.__set_element((i,j),self.__get_element((x,y)))
|
||||
|
||||
|
||||
return m
|
||||
|
||||
|
||||
else:
|
||||
# single element extraction
|
||||
if key[0] >= self.__rows or key[1] >= self.__cols:
|
||||
@@ -511,7 +511,7 @@ class _matrix(object):
|
||||
key = (key, 0)
|
||||
else:
|
||||
raise IndexError('insufficient indices for matrix')
|
||||
# Slice indexing
|
||||
# Slice indexing
|
||||
if isinstance(key[0],slice) or isinstance(key[1],slice):
|
||||
# Rows
|
||||
if isinstance(key[0],slice):
|
||||
@@ -563,11 +563,10 @@ class _matrix(object):
|
||||
self.__data[key] = value
|
||||
elif key in self.__data:
|
||||
del self.__data[key]
|
||||
|
||||
|
||||
if self._LU:
|
||||
self._LU = None
|
||||
return
|
||||
|
||||
|
||||
def __iter__(self):
|
||||
for i in xrange(self.__rows):
|
||||
|
||||
@@ -156,4 +156,3 @@ if __name__ == '__main__':
|
||||
r.write_results(show_missing=True, summary=True, coverdir="/tmp")
|
||||
else:
|
||||
testit(importdir, testdir)
|
||||
|
||||
|
||||
@@ -125,10 +125,10 @@ def test_complex_misc():
|
||||
|
||||
def test_complex_zeros():
|
||||
for a in [0,2]:
|
||||
for b in [0,3]:
|
||||
for c in [0,4]:
|
||||
for d in [0,5]:
|
||||
assert mpc(a,b)*mpc(c,d) == complex(a,b)*complex(c,d)
|
||||
for b in [0,3]:
|
||||
for c in [0,4]:
|
||||
for d in [0,5]:
|
||||
assert mpc(a,b)*mpc(c,d) == complex(a,b)*complex(c,d)
|
||||
|
||||
def test_hash():
|
||||
for i in range(-256, 256):
|
||||
@@ -197,23 +197,23 @@ def test_exact_integer_arithmetic():
|
||||
# XXX: re-fix this so that all operations are tested with all rounding modes
|
||||
random.seed(0)
|
||||
for prec in [6, 10, 25, 40, 100, 250, 725]:
|
||||
for rounding in ['d', 'u', 'f', 'c', 'n']:
|
||||
mp.dps = prec
|
||||
M = 10**(prec-2)
|
||||
M2 = 10**(prec//2-2)
|
||||
for i in range(10):
|
||||
a = random.randint(-M, M)
|
||||
b = random.randint(-M, M)
|
||||
assert mpf(a, rounding=rounding) == a
|
||||
assert int(mpf(a, rounding=rounding)) == a
|
||||
assert int(mpf(str(a), rounding=rounding)) == a
|
||||
assert mpf(a) + mpf(b) == a + b
|
||||
assert mpf(a) - mpf(b) == a - b
|
||||
assert -mpf(a) == -a
|
||||
a = random.randint(-M2, M2)
|
||||
b = random.randint(-M2, M2)
|
||||
assert mpf(a) * mpf(b) == a*b
|
||||
assert mpf_mul(from_int(a), from_int(b), mp.prec, rounding) == from_int(a*b)
|
||||
for rounding in ['d', 'u', 'f', 'c', 'n']:
|
||||
mp.dps = prec
|
||||
M = 10**(prec-2)
|
||||
M2 = 10**(prec//2-2)
|
||||
for i in range(10):
|
||||
a = random.randint(-M, M)
|
||||
b = random.randint(-M, M)
|
||||
assert mpf(a, rounding=rounding) == a
|
||||
assert int(mpf(a, rounding=rounding)) == a
|
||||
assert int(mpf(str(a), rounding=rounding)) == a
|
||||
assert mpf(a) + mpf(b) == a + b
|
||||
assert mpf(a) - mpf(b) == a - b
|
||||
assert -mpf(a) == -a
|
||||
a = random.randint(-M2, M2)
|
||||
b = random.randint(-M2, M2)
|
||||
assert mpf(a) * mpf(b) == a*b
|
||||
assert mpf_mul(from_int(a), from_int(b), mp.prec, rounding) == from_int(a*b)
|
||||
mp.dps = 15
|
||||
|
||||
def test_odd_int_bug():
|
||||
|
||||
@@ -189,4 +189,3 @@ def test_mpmathify():
|
||||
assert mpmathify('(1.0+1.0j)') == mpc(1, 1)
|
||||
assert mpmathify('(1.2e-10 - 3.4e5j)') == mpc('1.2e-10', '-3.4e5')
|
||||
assert mpmathify('1j') == mpc(1j)
|
||||
|
||||
|
||||
+26
-26
@@ -146,34 +146,34 @@ def test_eig_dyn():
|
||||
run_eig(A, verbose = v)
|
||||
|
||||
def test_eig():
|
||||
v = 0
|
||||
AS = []
|
||||
v = 0
|
||||
AS = []
|
||||
|
||||
A = mp.matrix([[2, 1, 0], # jordan block of size 3
|
||||
[0, 2, 1],
|
||||
[0, 0, 2]])
|
||||
AS.append(A)
|
||||
AS.append(A.transpose())
|
||||
A = mp.matrix([[2, 1, 0], # jordan block of size 3
|
||||
[0, 2, 1],
|
||||
[0, 0, 2]])
|
||||
AS.append(A)
|
||||
AS.append(A.transpose())
|
||||
|
||||
A = mp.matrix([[2, 0, 0], # jordan block of size 2
|
||||
[0, 2, 1],
|
||||
[0, 0, 2]])
|
||||
AS.append(A)
|
||||
AS.append(A.transpose())
|
||||
A = mp.matrix([[2, 0, 0], # jordan block of size 2
|
||||
[0, 2, 1],
|
||||
[0, 0, 2]])
|
||||
AS.append(A)
|
||||
AS.append(A.transpose())
|
||||
|
||||
A = mp.matrix([[2, 0, 1], # jordan block of size 2
|
||||
[0, 2, 0],
|
||||
[0, 0, 2]])
|
||||
AS.append(A)
|
||||
AS.append(A.transpose())
|
||||
A = mp.matrix([[2, 0, 1], # jordan block of size 2
|
||||
[0, 2, 0],
|
||||
[0, 0, 2]])
|
||||
AS.append(A)
|
||||
AS.append(A.transpose())
|
||||
|
||||
A= mp.matrix([[0, 0, 1], # cyclic
|
||||
[1, 0, 0],
|
||||
[0, 1, 0]])
|
||||
AS.append(A)
|
||||
AS.append(A.transpose())
|
||||
A= mp.matrix([[0, 0, 1], # cyclic
|
||||
[1, 0, 0],
|
||||
[0, 1, 0]])
|
||||
AS.append(A)
|
||||
AS.append(A.transpose())
|
||||
|
||||
for A in AS:
|
||||
run_hessenberg(A, verbose = v)
|
||||
run_schur(A, verbose = v)
|
||||
run_eig(A, verbose = v)
|
||||
for A in AS:
|
||||
run_hessenberg(A, verbose = v)
|
||||
run_schur(A, verbose = v)
|
||||
run_eig(A, verbose = v)
|
||||
|
||||
@@ -321,37 +321,37 @@ def test_gauss_quadrature_static():
|
||||
run_gauss("laguerre", a , b)
|
||||
|
||||
def test_gauss_quadrature_dynamic(verbose = False):
|
||||
n = 5
|
||||
n = 5
|
||||
|
||||
A = mp.randmatrix(2 * n, 1)
|
||||
A = mp.randmatrix(2 * n, 1)
|
||||
|
||||
def F(x):
|
||||
r = 0
|
||||
for i in xrange(len(A) - 1, -1, -1):
|
||||
r = r * x + A[i]
|
||||
return r
|
||||
def F(x):
|
||||
r = 0
|
||||
for i in xrange(len(A) - 1, -1, -1):
|
||||
r = r * x + A[i]
|
||||
return r
|
||||
|
||||
def run(qtype, FW, R, alpha = 0, beta = 0):
|
||||
X, W = mp.gauss_quadrature(n, qtype, alpha = alpha, beta = beta)
|
||||
def run(qtype, FW, R, alpha = 0, beta = 0):
|
||||
X, W = mp.gauss_quadrature(n, qtype, alpha = alpha, beta = beta)
|
||||
|
||||
a = 0
|
||||
for i in xrange(len(X)):
|
||||
a += W[i] * F(X[i])
|
||||
a = 0
|
||||
for i in xrange(len(X)):
|
||||
a += W[i] * F(X[i])
|
||||
|
||||
b = mp.quad(lambda x: FW(x) * F(x), R)
|
||||
b = mp.quad(lambda x: FW(x) * F(x), R)
|
||||
|
||||
c = mp.fabs(a - b)
|
||||
c = mp.fabs(a - b)
|
||||
|
||||
if verbose:
|
||||
print(qtype, c, a, b)
|
||||
if verbose:
|
||||
print(qtype, c, a, b)
|
||||
|
||||
assert c < 1e-5
|
||||
assert c < 1e-5
|
||||
|
||||
run("legendre", lambda x: 1, [-1, 1])
|
||||
run("legendre01", lambda x: 1, [0, 1])
|
||||
run("hermite", lambda x: mp.exp(-x*x), [-mp.inf, mp.inf])
|
||||
run("laguerre", lambda x: mp.exp(-x), [0, mp.inf])
|
||||
run("glaguerre", lambda x: mp.sqrt(x)*mp.exp(-x), [0, mp.inf], alpha = 1 / mp.mpf(2))
|
||||
run("chebyshev1", lambda x: 1/mp.sqrt(1-x*x), [-1, 1])
|
||||
run("chebyshev2", lambda x: mp.sqrt(1-x*x), [-1, 1])
|
||||
run("jacobi", lambda x: (1-x)**(1/mp.mpf(3)) * (1+x)**(1/mp.mpf(5)), [-1, 1], alpha = 1 / mp.mpf(3), beta = 1 / mp.mpf(5) )
|
||||
run("legendre", lambda x: 1, [-1, 1])
|
||||
run("legendre01", lambda x: 1, [0, 1])
|
||||
run("hermite", lambda x: mp.exp(-x*x), [-mp.inf, mp.inf])
|
||||
run("laguerre", lambda x: mp.exp(-x), [0, mp.inf])
|
||||
run("glaguerre", lambda x: mp.sqrt(x)*mp.exp(-x), [0, mp.inf], alpha = 1 / mp.mpf(2))
|
||||
run("chebyshev1", lambda x: 1/mp.sqrt(1-x*x), [-1, 1])
|
||||
run("chebyshev2", lambda x: mp.sqrt(1-x*x), [-1, 1])
|
||||
run("jacobi", lambda x: (1-x)**(1/mp.mpf(3)) * (1+x)**(1/mp.mpf(5)), [-1, 1], alpha = 1 / mp.mpf(3), beta = 1 / mp.mpf(5) )
|
||||
|
||||
@@ -914,7 +914,7 @@ def test_gammainc_expint_n():
|
||||
# These tests are intended to check all cases of the low-level code
|
||||
# for upper gamma and expint with small integer index.
|
||||
# Need to cover positive/negative arguments; small/large/huge arguments
|
||||
# for both positive and negative indices, as well as indices 0 and 1
|
||||
# for both positive and negative indices, as well as indices 0 and 1
|
||||
# which may be special-cased
|
||||
mp.dps = 15
|
||||
assert expint(-3,3.5).ae(0.021456366563296693987)
|
||||
|
||||
@@ -289,4 +289,3 @@ def test_exp_hp():
|
||||
# IntegerPart[N[Exp[1/10] * 10^4000, 4000]]
|
||||
# ...92167105162069688129
|
||||
assert int(r * 10**mp.dps) % 10**20 == 92167105162069688129
|
||||
|
||||
|
||||
+13
-13
@@ -101,23 +101,23 @@ def test_levin_3():
|
||||
assert err < eps
|
||||
|
||||
def test_levin_nsum():
|
||||
mp.dps = 17
|
||||
mp.dps = 17
|
||||
|
||||
with mp.extraprec(mp.prec):
|
||||
z = mp.mpf(10) ** (-10)
|
||||
a = mp.nsum(lambda n: n**(-(1+z)), [1, mp.inf], method = "l") - 1 / z
|
||||
assert abs(a - mp.euler) < 1e-10
|
||||
with mp.extraprec(mp.prec):
|
||||
z = mp.mpf(10) ** (-10)
|
||||
a = mp.nsum(lambda n: n**(-(1+z)), [1, mp.inf], method = "l") - 1 / z
|
||||
assert abs(a - mp.euler) < 1e-10
|
||||
|
||||
eps = mp.exp(0.8 * mp.log(mp.eps))
|
||||
eps = mp.exp(0.8 * mp.log(mp.eps))
|
||||
|
||||
a = mp.nsum(lambda n: (-1)**(n-1) / n, [1, mp.inf], method = "sidi")
|
||||
assert abs(a - mp.log(2)) < eps
|
||||
a = mp.nsum(lambda n: (-1)**(n-1) / n, [1, mp.inf], method = "sidi")
|
||||
assert abs(a - mp.log(2)) < eps
|
||||
|
||||
z = 2 + 1j
|
||||
f = lambda n: mp.rf(2 / mp.mpf(3), n) * mp.rf(4 / mp.mpf(3), n) * z**n / (mp.rf(1 / mp.mpf(3), n) * mp.fac(n))
|
||||
v = mp.nsum(f, [0, mp.inf], method = "levin", steps = [10 for x in xrange(1000)])
|
||||
exact = mp.hyp2f1(2 / mp.mpf(3), 4 / mp.mpf(3), 1 / mp.mpf(3), z)
|
||||
assert abs(exact - v) < eps
|
||||
z = 2 + 1j
|
||||
f = lambda n: mp.rf(2 / mp.mpf(3), n) * mp.rf(4 / mp.mpf(3), n) * z**n / (mp.rf(1 / mp.mpf(3), n) * mp.fac(n))
|
||||
v = mp.nsum(f, [0, mp.inf], method = "levin", steps = [10 for x in xrange(1000)])
|
||||
exact = mp.hyp2f1(2 / mp.mpf(3), 4 / mp.mpf(3), 1 / mp.mpf(3), z)
|
||||
assert abs(exact - v) < eps
|
||||
|
||||
def test_cohen_alt_0():
|
||||
mp.dps = 17
|
||||
|
||||
@@ -302,7 +302,7 @@ def test_qr():
|
||||
if n1 > maxnorm:
|
||||
raise ValueError('Excessive norm value')
|
||||
|
||||
if dtype == 'real':
|
||||
if dtype == 'real':
|
||||
n1 = norm(eye(m) - Q.T * Q)
|
||||
#print ' Norm of I - Q.T * Q = ', n1
|
||||
if n1 > maxnorm:
|
||||
@@ -313,7 +313,7 @@ def test_qr():
|
||||
if n1 > maxnorm:
|
||||
raise ValueError('Excessive norm value')
|
||||
|
||||
if dtype == 'complex':
|
||||
if dtype == 'complex':
|
||||
n1 = norm(eye(m) - Q.T * Q.conjugate())
|
||||
#print ' Norm of I - Q.T * Q.conjugate() = ', n1
|
||||
if n1 > maxnorm:
|
||||
|
||||
@@ -57,17 +57,17 @@ def test_matrix_basic():
|
||||
A9[0,0] = -100
|
||||
assert A9 != A10
|
||||
assert nstr(A9)
|
||||
|
||||
|
||||
def test_matrix_slices():
|
||||
A = matrix([ [1, 2, 3],
|
||||
[4, 5 ,6],
|
||||
A = matrix([ [1, 2, 3],
|
||||
[4, 5 ,6],
|
||||
[7, 8 ,9]])
|
||||
V = matrix([1,2,3,4,5])
|
||||
|
||||
# Get slice
|
||||
assert A[:,:] == A
|
||||
assert A[:,:] == A
|
||||
assert A[:,1] == matrix([[2],[5],[8]])
|
||||
assert A[2,:] == matrix([[7, 8 ,9]])
|
||||
assert A[2,:] == matrix([[7, 8 ,9]])
|
||||
assert A[1:3,1:3] == matrix([[5,6],[8,9]])
|
||||
assert V[2:4] == matrix([3,4])
|
||||
try:
|
||||
@@ -75,24 +75,24 @@ def test_matrix_slices():
|
||||
assert False
|
||||
except IndexError:
|
||||
pass
|
||||
|
||||
|
||||
# Assign slice with matrix
|
||||
A1 = matrix(3)
|
||||
A1[:,:] = A
|
||||
assert A1[:,:] == matrix([[1, 2, 3],
|
||||
[4, 5 ,6],
|
||||
assert A1[:,:] == matrix([[1, 2, 3],
|
||||
[4, 5 ,6],
|
||||
[7, 8 ,9]])
|
||||
A1[0,:] = matrix([[10, 11, 12]])
|
||||
assert A1 == matrix([ [10, 11, 12],
|
||||
[4, 5 ,6],
|
||||
[7, 8 ,9]])
|
||||
assert A1 == matrix([ [10, 11, 12],
|
||||
[4, 5 ,6],
|
||||
[7, 8 ,9]])
|
||||
A1[:,2] = matrix([[13], [14], [15]])
|
||||
assert A1 == matrix([ [10, 11, 13],
|
||||
[4, 5 ,14],
|
||||
[7, 8 ,15]])
|
||||
assert A1 == matrix([ [10, 11, 13],
|
||||
[4, 5 ,14],
|
||||
[7, 8 ,15]])
|
||||
A1[:2,:2] = matrix([[16, 17], [18 , 19]])
|
||||
assert A1 == matrix([ [16, 17, 13],
|
||||
[18, 19 ,14],
|
||||
assert A1 == matrix([ [16, 17, 13],
|
||||
[18, 19 ,14],
|
||||
[7, 8 ,15]])
|
||||
V[1:3] = 10
|
||||
assert V == matrix([1,10,10,4,5])
|
||||
@@ -107,12 +107,12 @@ def test_matrix_slices():
|
||||
assert False
|
||||
except IndexError:
|
||||
pass
|
||||
|
||||
|
||||
# Assign slice with scalar
|
||||
A1[:,2] = 10
|
||||
assert A1 == matrix([ [16, 17, 10],
|
||||
[18, 19 ,10],
|
||||
[7, 8 ,10]])
|
||||
assert A1 == matrix([ [16, 17, 10],
|
||||
[18, 19 ,10],
|
||||
[7, 8 ,10]])
|
||||
A1[:,:] = 40
|
||||
for x in A1:
|
||||
assert x == 40
|
||||
@@ -198,4 +198,3 @@ def test_matrix_numpy():
|
||||
l = [[1, 2], [3, 4], [5, 6]]
|
||||
a = numpy.matrix(l)
|
||||
assert matrix(l) == matrix(a)
|
||||
|
||||
|
||||
@@ -75,5 +75,3 @@ def test_trivial():
|
||||
assert findroot(lambda x: 0, 1) == 1
|
||||
assert findroot(lambda x: x, 0) == 0
|
||||
#assert findroot(lambda x, y: x + y, (1, -1)) == (1, -1)
|
||||
|
||||
|
||||
|
||||
@@ -12,4 +12,3 @@ def test_nstr():
|
||||
'''[ 0.75 0.1909 -0.02992]
|
||||
[ 0.1909 0.6563 0.2057]
|
||||
[-0.02992 0.2057 6.445e-21]'''
|
||||
|
||||
|
||||
@@ -227,4 +227,3 @@ if __name__ == '__main__':
|
||||
mapf(testit, tasks)
|
||||
t2 = clock()
|
||||
print("Cumulative wall time:", t2-t1)
|
||||
|
||||
|
||||
Reference in New Issue
Block a user