[Axiom-developer] 20080119.03.tpd.patch

[daly]

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