axiom-developer
[Top][All Lists]

## [Axiom-developer] 20080119.03.tpd.patch

 From: daly Subject: [Axiom-developer] 20080119.03.tpd.patch Date: Sat, 19 Jan 2008 20:48:29 -0600

```E1(0.0) is infinity. This fixes the special case to give correct results.
e1.input has a new regression test for infinity at 0.0

Because of the infinity the range of E1 is OnePointCompletion(DoubleFloat)
and we need to add some explicit coercions.

Tim

===========================================================================
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
index 241d23e..0744b97 100644
@@ -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

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

<<package DFSFUN DoubleFloatSpecialFunctions>>=
-        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)

```