[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Axiom-developer] 20080221.02.wxh.patch (236)
From: |
daly |
Subject: |
[Axiom-developer] 20080221.02.wxh.patch (236) |
Date: |
Thu, 21 Feb 2008 00:33:53 -0600 |
Increase the total digits in Pi representation in special functions (Waldek)
=======================================================================
diff --git a/changelog b/changelog
index c7165ea..eea576d 100644
--- a/changelog
+++ b/changelog
@@ -1,3 +1,4 @@
+20080221 wxh src/interp/sfsfun.boot increase precision of PI (236)
20080221 tpd src/input/gamma.input investigate complex gamma issues
20080220 tpd src/hyper/bookvol11 add additional hyperdoc page translations
20080219 tpd src/hyper/Makefile handle hyperdoc bitmaps properly
diff --git a/src/interp/sfsfun.boot.pamphlet b/src/interp/sfsfun.boot.pamphlet
index 50bd7b5..6536e2e 100644
--- a/src/interp/sfsfun.boot.pamphlet
+++ b/src/interp/sfsfun.boot.pamphlet
@@ -105,6 +105,11 @@ negintp(x) ==
else
false
+-- Lisp PI is a long float and causes type errors, here we give
+-- enough digits to have double accuracy even after conversion
+-- to binary
+DEFCONSTANT(dfPi,3.14159265358979323846264338328)
+
--- Small float implementation of Gamma function. Valid for
--- real arguments. Maximum error (relative) seems to be
--- 2-4 ulps for real x 2<x<9, and up to ten ulps for larger x
@@ -132,7 +137,7 @@ gammaStirling(x) ==
EXP(lnrgamma(x))
lnrgammaRatapprox(x) ==
- (x-.5)*LOG(x) - x + LOG(SQRT(2.0*PI)) + phiRatapprox(x)
+ (x-.5)*LOG(x) - x + LOG(SQRT(2.0*dfPi)) + phiRatapprox(x)
phiRatapprox(x) ==
arg := 1/(x**2)
@@ -167,7 +172,7 @@ gammaRatapprox (x) ==
prod := */[x+i for i in 0..n-1]
result := gammaRatapprox(reducedarg)/prod
else
- Pi := PI
+ Pi := dfPi
lx := MULTIPLE_-VALUE_-LIST(FLOOR(x))
intpartx := CAR(lx)+1
restx := CADR(lx)
@@ -247,12 +252,12 @@ PiMinusLogSinPi(z1,z2,z) ==
cgammaG(z1,z2) - logH(z1,z2,z)
cgammaG(z1,z2) ==
- LOG(2*PI) + PI*z2 - COMPLEX(0.0,1.0)*PI*(z1-.5)
+ LOG(2*dfPi) + dfPi*z2 - COMPLEX(0.0,1.0)*dfPi*(z1-.5)
logH(z1,z2,z) ==
z1bar := CADR(MULTIPLE_-VALUE_-LIST(FLOOR(z1))) ---frac part of z1
- piz1bar := PI*z1bar
- piz2 := PI*z2
+ piz1bar := dfPi*z1bar
+ piz2 := dfPi*z2
twopiz2 := 2.0*piz2
i := COMPLEX(0.0,1.0)
part2 := EXP(twopiz2)*(2.0*SIN(piz1bar)**2 + SIN(2.0*piz1bar)*i)
@@ -284,7 +289,7 @@ logS(z1,z2,z,n,zpn) ==
--- adjusted by 2 Pi if it is negative.
cgammaAdjust(z) ==
if IMAGPART(z)<0.0
- then result := z + COMPLEX(0.0, 2.0*PI)
+ then result := z + COMPLEX(0.0, 2.0*dfPi)
else result := z
result
@@ -292,7 +297,7 @@ clngammacase3(z) ==
(z- .5)*LOG(z) - z + cgammaBernsum(z)
cgammaBernsum (z) ==
- sum := LOG(2.0*PI)/2.0
+ sum := LOG(2.0*dfPi)/2.0
zterm := z
zsquaredinv := 1.0/(z*z)
l:= [.083333333333333333333, -.0027777777777777777778,_
@@ -365,7 +370,7 @@ rPsi(n,x) ==
skipit := 1
else
skipit := 0
- sign*((PI**(n+1))*cotdiffeval(n,PI*(-x),skipit) +
rPsi(n,1.0-x))
+ sign*((dfPi**(n+1))*cotdiffeval(n,dfPi*(-x),skipit) +
rPsi(n,1.0-x))
else if n=0
then
- rPsiW(n,x)
@@ -508,7 +513,7 @@ PsiXotic(n,result) ==
cpsireflect(n,z) ==
m := MOD(n,2)
sign := (-1)**m
- sign*PI**(n+1)*cotdiffeval(n,PI*z,0) + cPsi(n,1.0-z)
+ sign*dfPi**(n+1)*cotdiffeval(n,dfPi*z,0) + cPsi(n,1.0-z)
--- c parameter to 0F1, possibly complex
--- z argument to 0F1
@@ -600,7 +605,7 @@ f01(c,z)==
--- t := SQRT(-z)
--- c1 := c-1.0
--- p := PHASE(c)
---- if ABS(c)>10.0*ABS(t) and p>=0.0 and PHASE(c)<.90*PI
+--- if ABS(c)>10.0*ABS(t) and p>=0.0 and PHASE(c)<.90*dfPi
--- then BesselJAsymptOrder(c1,2*t)*cgamma(c/(t**(c1)))
--- else if ABS(t)>10.0*ABS(c) and ABS(t)>50.0
--- then BesselJAsympt(c1,2*t)*cgamma(c/(t**(c1)))
@@ -904,7 +909,7 @@ BesselasymptB(mu,z,zsqr,zfth) ==
--- Asymptotic series only works when |v| < |z|.
BesselJAsympt (v,z) ==
- pi := PI
+ pi := dfPi
mu := 4.0*v*v
zsqr := z*z
zfth := zsqr*zsqr
@@ -931,9 +936,9 @@ BesselIAsympt(v,z,n) ==
(float(r)*eight*z)
sum1 := sum1 + term1
sum2 := sum2 + ABS(term1)
- sqrttwopiz := SQRT(two*PI*z)
+ sqrttwopiz := SQRT(two*dfPi*z)
EXP(z)/sqrttwopiz*(1.0 + sum1 ) +_
- EXP(-(float(n)+.5)*PI*i)*EXP(-z)/sqrttwopiz*(1.0+ sum2)
+ EXP(-(float(n)+.5)*dfPi*i)*EXP(-z)/sqrttwopiz*(1.0+ sum2)
---Asymptotic formula for BesselJ when order is large comes from
@@ -951,7 +956,7 @@ BesselJAsymptOrder(v,z) ==
-- cothalpha := 1.0/tanhalpha
ca := 1.0/tanhalpha
- Pi := PI
+ Pi := dfPi
ca2:=ca*ca
ca4:=ca2*ca2
ca8:=ca4*ca4
@@ -974,7 +979,7 @@ BesselJAsymptOrder(v,z) ==
--- See Olver, p. 376-382.
BesselIAsymptOrder(v,vz) ==
z := vz/v
- Pi := PI
+ Pi := dfPi
--- Use reflection formula (Atlas, p. 492) if v not in right half plane;
Is this always accurate?
if REALPART(v)<0.0
then return BesselIAsymptOrder(-v,vz) +
2.0/Pi*SIN(-v*Pi)*BesselKAsymptOrder(-v,vz)
@@ -1015,7 +1020,7 @@ BesselKAsymptOrder (v,vz) ==
u4p :=
(3675.0/32768.0+(-96833.0/40960.0+(144001.0/16384.0+(-7436429.0/663552.0+37182145.0/7962624.0*p2)*p2)*p2)*p2)*p4
u5p :=
((59535.0/262144.0+(-67608983.0/9175040.0+(250881631.0/5898240.0+(-108313205.0/1179648.0+(5391411025.0/63700992.0-5391411025.0/191102976.0*p2)*p2)*p2)*p2)*p2)*p4*p)*(-1.0)
hornerresult := horner([u5p,u4p,u3p,u2p,u1p,u0p],vinv)
- SQRT(PI/(2.0*v))*EXP(-v*eta)/(SQRT(opzsqroh))*hornerresult
+ SQRT(dfPi/(2.0*v))*EXP(-v*eta)/(SQRT(opzsqroh))*hornerresult
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Axiom-developer] 20080221.02.wxh.patch (236),
daly <=