/usr/share/axiom-20120501/input/numericgamma.input

)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


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


--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)
axiom-test 20120501-8 / usr / share / axiom-20120501 / input / numericgamma.input
