$SPAD/src/input numericgamma.input

Bill Page December 3, 2016

Abstract

Contents

2

— * —

)set break resume
)sys rm -f numericgamma.output
)spool numericgamma.output
)set message test on
)set message auto off
)clear all
)sys cp $AXIOM/../../src/input/numericgamma.input.pamphlet .
)lisp (tangle "numericgamma.input.pamphlet" "sfx.spad" "sfx.spad")
)co sfx.spad
--S 1 of 36
Gam(a:Float,x:Float):Float ==
  if x < 0.0 or a < 0.0 then error "Invalid arguments"
  if x = 0.0 then return Gamma(a)
  ITMAX ==> 100        -- Maximum allowed number of iterations
  FPMIN ==> 1.0e-1000  -- near the smallest representable number
                       -- (there is no smallest representable float)
  EPS := (10.0^(-digits()$Float+1))$Float  -- Relative accuracy
  an: Float
  del: Float
  b:Float:=x+1.0-a     -- Set up for evaluating continued fractions
  c:Float:=1.0/FPMIN   -- by modified Lentz's method
  d:Float:=1.0/b       -- with b_0 = 0
  h:Float:=d
  i:=1
  repeat               -- iterate to convergence
    an:=-i*(i-a)
    b:=b+2.0
    d:=an*d+b
    if abs(d) < FPMIN then d:=FPMIN
    c:=b+an/c;
    if abs(c) < FPMIN then c:=FPMIN
    d:=1.0/d
    del:=d*c
    h:=h*del
    if i > ITMAX or abs(del-1.0) < EPS then break
    i:=i+1
  if i > ITMAX then error("a too large, ITMAX too small")
  exp(-x)*x^a*h        -- put factors in front
--R
--R   Function declaration Gam : (Float,Float) -> Float has been added to
--R      workspace.
--R                                                                   Type: Void

--E 1
--S 2 of 36
Gam(0,1)
--R
--R   Compiling function Gam with type (Float,Float) -> Float
--R
--RDaly Bug
--R   Error signalled from user code in function Gam:
--R      a too large, ITMAX too small
--E 2
--S 3 of 36
Gam(1.1.1)
--R
--R   There are 1 exposed and 1 unexposed library operations named elt
--R      having 1 argument(s) but none was determined to be applicable.
--R      Use HyperDoc Browse, or issue
--R                               )display op elt
--R      to learn more about the available operations. Perhaps
--R      package-calling the operation or using coercions on the arguments
--R      will allow you to apply the operation.
--R
--RDaly Bug
--R   Cannot find application of object of type Float to argument(s) of
--R      type(s)
--R                                    Float
--R
--E 3
--S 4 of 36
Gam(5,10)
--R
--R
--R   (2)  0.7020645138 4706574415
--R                                                                  Type: Float
--E 4
--S 5 of 36
Gam(5,11)
--R
--R
--R   (3)  0.3625104156 5228203538
--R                                                                  Type: Float
--E 5
--S 6 of 36
Gam(7,0)
--R
--R

--R   (4)  720.0000000000 0011369
--R                                                                  Type: Float
--E 6
--S 7 of 36
digits 100
--R
--R
--R   (5)  20
--R                                                        Type: PositiveInteger
--E 7
--S 8 of 36
Gam(0,1)
--R
--R
--RDaly Bug
--R   Error signalled from user code in function Gam:
--R      a too large, ITMAX too small
--E 8
--S 9 of 36
Gam(1,1.1)
--R
--R
--R   (6)
--R  0.3328710836 9807955328 8846906431 3155216124 7952156921 2491793331 386750747
--R  0 8541284431 1612617072 7005478519
--R                                                                  Type: Float
--E 9
--S 10 of 36
Gam(1,1)
--R
--R
--R   (7)
--R  0.3678794411 7144232159 5523770161 4608674458 1113103176 7834507836 801697461
--R  4 9574489980 3357147274 3459196437
--R                                                                  Type: Float
--E 10
--S 11 of 36
Gam(1,1.1)
--R
--R
--R   (8)
--R  0.3328710836 9807955328 8846906431 3155216124 7952156921 2491793331 386750747
--R  0 8541284431 1612617072 7005478519
--R                                                                  Type: Float
--E 11

--S 12 of 36
Gam(5,10)
--R
--R
--R   (9)
--R  0.7020645138 4706574414 6387196628 3546367191 6532623256 0684622278 670587055
--R  0 5584357048 3474646670 2985365058
--R                                                                  Type: Float
--E 12
--S 13 of 36
Gam(5,11)
--R
--R
--R   (10)
--R  0.3625104156 5228203538 0753904311 4079803866 4530925132 7036797697 419049037
--R  4 2658968752 0305953551 1648548436
--R                                                                  Type: Float
--E 13

———-

Note that in the case of Gam(7,0) the function Gamma(a) is called in the above code but it is actually implemented in Axiom as Gamma(a)$DoubleFloat with a fixed precision so the precision of the result is not affected by the setting of digits().

— * —

--S 14 of 36
Gam(7,0)
--R
--R
--R   (11)  720.0000000000 0011368683 7721616029 7393798828 125
--R                                                                  Type: Float
--E 14
--S 15 of 36
Gam(7,0.1)
--R
--R
--R   (12)
--R  719.9999999869 1035963050 9717349089 5137595484 2683243460 6577519316 5312727
--R  417 6619922456 9102294155 8764196
--R                                                                  Type: Float
--E 15
--S 16 of 36
Gam(7,0.2)
--R

--R
--R   (13)
--R  719.9999984646 1597708521 8246915701 3222705579 4693807497 2229652513 6047137
--R  980 7138425860 0596921944 0451807
--R                                                                  Type: Float
--E 16

———-

— sfx.spad —

)abbrev package SFX SpecialFunctionExtended
SpecialFunctionExtended: Exports == Implementation where
ITMAX ==> 100.0::DoubleFloat    -- Maximum allowed number of iterations
FPMIN ==> 1.0e-323::DoubleFloat -- near the smallest representable number
Exports ==> with
NGamma:(DoubleFloat,DoubleFloat)->DoubleFloat
Implementation ==> add
NGamma(a:DoubleFloat,x:DoubleFloat):DoubleFloat ==
if x < 0 or a < 0 then error "Invalid arguments"
if x = 0 then return Gamma(a)
EPS := (10.0::DoubleFloat**(-digits()$DoubleFloat))$DoubleFloat
b:DoubleFloat:=x+1-a     -- Set up for evaluating continued fractions
c:DoubleFloat:=1/FPMIN   -- by modified Lentz's method
d:DoubleFloat:=1/b       -- with b_0 = 0
h:DoubleFloat:=d
i:DoubleFloat:=1
repeat               -- iterate to convergence
  an:DoubleFloat:=-i*(i-a)
  b:=b+2.0::DoubleFloat
  d:=an*d+b
  if abs(d) < FPMIN then d:=FPMIN
  c:=b+an/c
  if abs(c) < FPMIN then c:=FPMIN
  d:=1/d
  del:DoubleFloat:=d*c
  h:=h*del
  if i > ITMAX or abs(del-1) < EPS then break
  i:=i+1
if i > ITMAX then error("a too large, ITMAX too small")
exp(-x+log(x)*a)*h        -- put factors in front

———-

— * —

--S 17 of 36
NGamma(a,x)
--R
--R   There are 1 exposed and 0 unexposed library operations named NGamma
--R      having 2 argument(s) but none was determined to be applicable.
--R      Use HyperDoc Browse, or issue
--R                             )display op NGamma
--R      to learn more about the available operations. Perhaps
--R      package-calling the operation or using coercions on the arguments
--R      will allow you to apply the operation.
--R
--RDaly Bug
--R   Cannot find a definition or applicable library operation named
--R      NGamma with argument type(s)
--R                                 Variable(a)
--R                                 Variable(x)
--R
--R      Perhaps you should use "@" to indicate the required return type,
--R      or "$" to specify which version of the function you need.
--E 17
--S 18 of 36
machineFraction(NGamma(0,1))
--R
--R          7904139241557877
--R   (14)  -----------------
--R         36028797018963968
--R                                                      Type: Fraction(Integer)
--E 18
--S 19 of 36
machineFraction(NGamma(0,2))
--R
--R          7047306297697619
--R   (15)  ------------------
--R         144115188075855872
--R                                                      Type: Fraction(Integer)
--E 19
--S 20 of 36
machineFraction(NGamma(1,1))
--R
--R          828390857088487
--R   (16)  ----------------
--R         2251799813685248
--R                                                      Type: Fraction(Integer)
--E 20
--S 21 of 36
machineFraction(NGamma(1,1.1))

--R
--R
--R          5996472354020337
--R   (17)  -----------------
--R         18014398509481984
--R                                                      Type: Fraction(Integer)
--E 21
--S 22 of 36
machineFraction(NGamma(5,10))
--R
--R         197613592684481
--R   (18)  ---------------
--R         281474976710656
--R                                                      Type: Fraction(Integer)
--E 22
--S 23 of 36
machineFraction(NGamma(5,11))
--R
--R         3265203545699083
--R   (19)  ----------------
--R         9007199254740992
--R                                                      Type: Fraction(Integer)
--E 23
--S 24 of 36
machineFraction(NGamma(7,0))
--R
--R         6333186975989761
--R   (20)  ----------------
--R           8796093022208
--R                                                      Type: Fraction(Integer)
--E 24
--S 25 of 36
machineFraction(NGamma(7,0.1))
--R
--R
--R         6333190201121973
--R   (21)  ----------------
--R           8796093022208
--R                                                      Type: Fraction(Integer)
--E 25
--S 26 of 36
machineFraction(NGamma(7,0.2))
--R
--R
--R         6333187094743041

--R   (22)  ----------------
--R           8796093022208
--R                                                      Type: Fraction(Integer)
--E 26
)set functions compile on
--S 27 of 36
j:=120
--R
--R
--R   (23)  120
--R                                                        Type: PositiveInteger
--E 27
--S 28 of 36
nume(a) == cons(1::Float,[((a-i)*i)::Float for i in 1..])
--R
--R                                                                   Type: Void
--E 28
--S 29 of 36
dene(a,x) == [(x+2*i+1-a)::Float for i in 0..]
--R
--R                                                                   Type: Void
--E 29
--S 30 of 36
cfe(a,x) == continuedFraction(0,nume(a),dene(a,x))
--R
--R                                                                   Type: Void
--E 30
--S 31 of 36
ccfe(a,x) == convergents cfe(a,x)
--R
--R                                                                   Type: Void
--E 31
--S 32 of 36
gamcfe(a,x) == exp(-x)*x^a*(ccfe(a,x).j)::Float
--R
--R                                                                   Type: Void
--E 32
--S 33 of 36
gamcfe(2,3)
--R
--R   Compiling function nume with type PositiveInteger -> Stream(Float)
--R   Compiling function dene with type (PositiveInteger,PositiveInteger)

--R       -> Stream(Float)
--R   Compiling function cfe with type (PositiveInteger,PositiveInteger)
--R       -> ContinuedFraction(Float)
--R   Compiling function ccfe with type (PositiveInteger,PositiveInteger)
--R       -> Stream(Fraction(Float))
--R   Compiling function gamcfe with type (PositiveInteger,PositiveInteger
--R      ) -> Expression(Float)
--R
--R   (29)
--R  0.1991482734 7145577191 7369662600 2471065267 9836875369 2862270510 910424242
--R  6 7092079820 0616216976 9465333782
--R                                                      Type: Expression(Float)
--E 33
--S 34 of 36
E1fun(x) == gamcfe(0,x)
--R
--R                                                                   Type: Void
--E 34
--S 35 of 36
E1fun(2.0)
--R
--R   Compiling function nume with type NonNegativeInteger -> Stream(Float
--R      )
--R   Compiling function dene with type (NonNegativeInteger,Float) ->
--R      Stream(Float)
--R   Compiling function cfe with type (NonNegativeInteger,Float) ->
--R      ContinuedFraction(Float)
--R   Compiling function ccfe with type (NonNegativeInteger,Float) ->
--R      Stream(Fraction(Float))
--R   Compiling function gamcfe with type (NonNegativeInteger,Float) ->
--R      Float
--R   Compiling function E1fun with type Float -> Float
--R
--R   (31)
--R  0.0489005107 0806111956 7239826914 3472898212 1544510421 3277251841 716377988
--R  0 9149832755 9949235928 1965882172 4
--R                                                                  Type: Float
--E 35
--S 36 of 36
E1fun(2.0)-E1(2.0)
--R
--R
--R   (32)  1.1102230246251565E-16
--R                                        Type: OnePointCompletion(DoubleFloat)
--E 36
)spool

)lisp (bye) ———-

12

References

[1] http://axiom-wiki.newsynthesis.org/SandBoxGamma