diff --git a/changelog b/changelog index e2767b4..2a361bf 100644 --- a/changelog +++ b/changelog @@ -1,3 +1,5 @@ +20080119 tpd src/input/e1.input regression test E1(0.0) +20080119 tpd src/algebra/special.spad handle E1(0.0) properly 20080119 tpd changelog credit FreeAbelianGroup fix to Franz Lehner 20080119 tpd src/input/Makefile add e1.input regression test file 20080119 tpd src/input/e1.input create A&S reference regression for E1 diff --git a/src/algebra/special.spad.pamphlet b/src/algebra/special.spad.pamphlet index 241d23e..0744b97 100644 --- a/src/algebra/special.spad.pamphlet +++ b/src/algebra/special.spad.pamphlet @@ -30,6 +30,7 @@ DoubleFloatSpecialFunctions(): Exports == Impl where NNI ==> NonNegativeInteger R ==> DoubleFloat C ==> Complex DoubleFloat + OPR ==> OnePointCompletion R Exports ==> with Gamma: R -> R @@ -39,7 +40,7 @@ DoubleFloatSpecialFunctions(): Exports == Impl where ++ Gamma(x) is the Euler gamma function, \spad{Gamma(x)}, defined by ++ \spad{Gamma(x) = integrate(t^(x-1)*exp(-t), t=0..%infinity)}. - E1: R -> R + E1: R -> OPR ++ E1(x) is the Exponential Integral function ++ The current implementation is a piecewise approximation ++ involving one poly from -4..4 and a second poly for x > 4 @@ -377,8 +378,8 @@ rewritten the polynomial using Horner's method so the large powers of $x$ are only computed once. <>= - E1(x:R):R == - x = 0.0::R => error "E1 undefined at zero" + E1(x:R):OPR == + x = 0.0::R => infinity() x > 4.0::R => t1:R:=0.14999948967737774608E-15::R t2:R:=0.9999999999993112::R @@ -426,7 +427,11 @@ powers of $x$ are only computed once. t23:R:=15474250491.067253436::R tv:R:=(tu*x-t23) tw:R:=(-1.0::R*x) - (tv * exp(tw)::R)/x**22 + tx:R:=exp(tw) + ty:R:=tv*tx + tz:R:=x**22 + taz:R:=ty/tz + taz::OPR x > -4.0::R => a1:R:=0.476837158203125E-22::R a2:R:=0.10967254638671875E-20::R @@ -473,7 +478,8 @@ powers of $x$ are only computed once. au:R:=(at*x+a22) a23:R:=0.5772156649015328::R av:R:=au*x-a23 - - 1.0::R*log(abs(x)) + av + aw:R:=- 1.0::R*log(abs(x)) + av + aw::OPR error "E1: no approximation available" polygamma(k,z) == CPSI(k, z)$Lisp diff --git a/src/input/e1.input.pamphlet b/src/input/e1.input.pamphlet index 0c24a64..329cb01 100644 --- a/src/input/e1.input.pamphlet +++ b/src/input/e1.input.pamphlet @@ -19,7 +19,7 @@ This is a regression test for E1(x) @ Here we enter the value of Gamma <<*>>= ---S 1 of 6 +--S 1 of 7 G:DFLOAT:=0.577215664901532860606512::DFLOAT --R --R (1) 0.57721566490153287 @@ -34,8 +34,8 @@ listed in Abramowitz and Stegun, ``Handbook of Mathematical Functions'', Dover Publications, Inc. New York 1965. pp238 <<*>>= ---S 2 of 6 -f(x)==x^-1 * (E1(x) + log(x) + G) +--S 2 of 7 +f(x)==x^-1 * (E1(x)::DFLOAT + log(x) + G) --R Type: Void --E 2 @ @@ -47,7 +47,7 @@ Abramowitz and Stegun, ``Handbook of Mathematical Functions'', Dover Publications, Inc. New York 1965. pp238 <<*>>= ---S 3 of 6 +--S 3 of 7 [[0.01,0.9975055452, f(0.01), f(0.01)-0.9975055452],_ [0.02,0.9950221392, f(0.02), f(0.02)-0.9950221392],_ [0.03,0.9925497201, f(0.03), f(0.03)-0.9925497201],_ @@ -97,7 +97,7 @@ Dover Publications, Inc. New York 1965. pp238 [0.47,0.8937670423, f(0.47), f(0.47)-0.8937670423],_ [0.48,0.8917309048, f(0.48), f(0.48)-0.8917309048],_ [0.49,0.8897032920, f(0.49), f(0.49)-0.8897032920],_ -[0.50,0.8876841584, f(0.50), f(0.50)-0.8876841584]] +[0.50,0.8876841584, f(0.50), f(0.50)-0.8876841584]]::LIST(LIST(DFLOAT)) --R --R Compiling function f with type Float -> DoubleFloat --R @@ -270,7 +270,7 @@ is the reference value of $E1(x)$ from the book Abramowitz and Stegun, ``Handbook of Mathematical Functions'', Dover Publications, Inc. New York 1965. pp239-241 <<*>>= ---S 4 of 6 +--S 4 of 7 [[0.50, 0.559773595, E1(0.50), E1(0.50)-0.559773595],_ [0.51, 0.547822352, E1(0.51), E1(0.51)-0.547822352],_ [0.52, 0.536219798, E1(0.52), E1(0.52)-0.536219798],_ @@ -421,7 +421,7 @@ Dover Publications, Inc. New York 1965. pp239-241 [1.97, 0.050976988, E1(1.97), E1(1.97)-0.050976988],_ [1.98, 0.050274392, E1(1.98), E1(1.98)-0.050274392],_ [1.99, 0.049582291, E1(1.99), E1(1.99)-0.049582291],_ -[2.00, 0.048900511, E1(2.00), E1(2.00)-0.048900511]] +[2.00, 0.048900511, E1(2.00), E1(2.00)-0.048900511]]::LIST(LIST(DFLOAT)) --R --R --R (4) @@ -835,8 +835,8 @@ Dover Publications, Inc. New York 1965. pp239-241 Now that we move into larger numbers we need a new scaling function. <<*>>= ---S 5 of 6 -g(x)==x * %e^x * E1(x) +--S 5 of 7 +g(x)==x * %e^x * E1(x)::DFLOAT --R Type: Void --E 5 @ @@ -845,7 +845,7 @@ And we compute the scaled value of E1(x) in the range 2.0 to 10.0 from Abramowitz and Stegun,``Handbook of Mathematical Functions'', Dover Publications, Inc. New York 1965. pp242-243 <<*>>= ---S 6 of 6 +--S 6 of 7 [[2.0,0.722657234,g(2.0),g(2.0)-0.722657234],_ [2.1,0.730791502,g(2.1),g(2.1)-0.730791502],_ [2.2,0.738431132,g(2.2),g(2.2)-0.738431132],_ @@ -1204,6 +1204,12 @@ Dover Publications, Inc. New York 1965. pp242-243 --R Type: List List Expression DoubleFloat --E 6 +--S 7 of 7 +E1(0.0) +--R +--R (7) infinity +--R Type: OnePointCompletion DoubleFloat +--E 7 )spool )lisp (bye)