From Numerical Recipes in C, 2nd edition

Section 6.2 Incomplete Gamma Function

http://www.library.cornell.edu/nr/bookcpdf/c6-2.pdf

This Axiom interpreter function implements equation (6.2.7)

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 not smallest number for 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 fraction
  c:Float := 1.0/FPMIN      --   by modified Lentz' method (section 5.2)
  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

   Function declaration Gam : (Float,Float) -> Float has been added to
      workspace.

Gam(0,1)

Compiling function Gam with type (Float,Float) -> Float 

   Error signalled from user code in function Gam: 
      a too large, ITMAX too small
Gam(0,1.1)

   Error signalled from user code in function Gam: 
      a too large, ITMAX too small
Gam(1,1)

(1)

Gam(1,1.1)

(2)

Gam(5,10)

(3)

Gam(5,11)

(4)

Gam(7,0)

(5)

Let's try that again with more precision:

digits(100)

(6)

Gam(0,1)

   Error signalled from user code in function Gam:
      a too large, ITMAX too small
Gam(0,1.1)

   Error signalled from user code in function Gam:
      a too large, ITMAX too small
Gam(1,1)

(7)

Gam(1,1.1)

(8)

Gam(5,10)

(9)

Gam(5,11)

(10)

Gam(7,0)

(11)

Notice 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 and so the precison of result in this case is not affected by the setting of digits().

So lets do this using machine floating point (DoubleFloat?) in Axiom's library programming language SPAD

)abbrev package SFX SpecialFunctionsExtended
SpecialFunctionsExtended: Exports == Implementation where
  ITMAX  ==> 100.0::DoubleFloat   -- Maximum allowed number of iterations
  FPMIN ==> 1.0e-323::DoubleFloat -- near the smallest representable number

  Exports ==> with
    Gamma:(DoubleFloat,DoubleFloat)->DoubleFloat

  Implementation ==> add

    Gamma(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  -- Relative accuracy

      b:DoubleFloat := x+1-a        -- Set up for evaluating continued fraction
      c:DoubleFloat := 1/FPMIN      --   by modified Lentz' method (section 5.2)
      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")
      return exp(-x+log(x)*a) * h       -- Put factors in front

   Compiling FriCAS source code from file 
      /var/zope2/var/LatexWiki/6599990542159375921-25px004.spad using 
      old system compiler.
   SFX abbreviates package SpecialFunctionsExtended 
   processing macro definition ITMAX ==> ::((elt (Float) float)(100,Zero,10),DoubleFloat) 
   processing macro definition FPMIN ==> ::((elt (Float) float)(One,-323,10),DoubleFloat) 
   processing macro definition Exports ==> -- the constructor category 
   processing macro definition Implementation ==> -- the constructor capsule 
------------------------------------------------------------------------
   initializing NRLIB SFX for SpecialFunctionsExtended 
   compiling into NRLIB SFX 
   compiling exported Gamma : (DoubleFloat,DoubleFloat) -> DoubleFloat
Time: 0.10 SEC.

(time taken in buildFunctor:  0)

;;;     ***       |SpecialFunctionsExtended| REDEFINED

;;;     ***       |SpecialFunctionsExtended| REDEFINED
Time: 0.05 SEC.

   Warnings: 
      [1] Gamma:  d has no value
      [2] Gamma:  c has no value

   Cumulative Statistics for Constructor SpecialFunctionsExtended
      Time: 0.15 seconds

   finalizing NRLIB SFX 
   Processing SpecialFunctionsExtended for Browser database:
--->-->SpecialFunctionsExtended((Gamma ((DoubleFloat) (DoubleFloat) (DoubleFloat)))): Not documented!!!!
--->-->SpecialFunctionsExtended(constructor): Not documented!!!!
--->-->SpecialFunctionsExtended(): Missing Description
------------------------------------------------------------------------
   SpecialFunctionsExtended is now explicitly exposed in frame initial 
   SpecialFunctionsExtended will be automatically loaded when needed 
      from /var/zope2/var/LatexWiki/SFX.NRLIB/code

Apparently in SPAD we must write 0 and 1 instead of 0.0 and 1.0 because 0 and 1 are functions exported by DoubleFloat but 0.0 and 1.0 are by default, constructs of Float. We must also be careful to declare the type of other constants.

Gamma(a,x)

(12)

Gamma(0,1)

(13)

Gamma(0,2)

(14)

Gamma(1,1)

(15)

Gamma(1,1.1)

(16)

Gamma(5,10)

(17)

Gamma(5,11)

(18)

Gamma(7,0)

(19)

Exactly the same code can also be compiled using Aldor (Axiom library compiler version 2):

#include "axiom.as";
#pile
SpecialFunctionsExtended2: Exports == Implementation where
  ITMAX  ==> 100.0::DoubleFloat   -- Maximum allowed number of iterations
  FPMIN ==> 1.0e-323::DoubleFloat -- near the smallest representable number

  Exports ==> with
    Gamma:(DoubleFloat,DoubleFloat)->DoubleFloat

  Implementation ==> add

    Gamma(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  -- Relative accuracy

      b:DoubleFloat := x+1-a        -- Set up for evaluating continued fraction
      c:DoubleFloat := 1/FPMIN      --   by modified Lentz' method (section 5.2)
      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")
      return exp(-x+log(x)*a) * h       -- Put factors in front

   Compiling FriCAS source code from file 
      /var/zope2/var/LatexWiki/sfx2.as using AXIOM-XL compiler and 
      options 
-O -Fasy -Fao -Flsp -laxiom -Mno-AXL_W_WillObsolete -DAxiom -Y $AXIOM/algebra
      Use the system command )set compiler args to change these 
      options.
#1 (Warning) Deprecated message prefix: use `ALDOR_' instead of `_AXL'
"/var/zope2/var/LatexWiki/sfx2.as", line 16: 
      EPS := (10.0::DoubleFloat**(-digits()$DoubleFloat))$DoubleFloat  -- Relative accuracy
......^............................^
[L16 C7] #3 (Error) (After Macro Expansion) No meaning for identifier `EPS'.
Expanded expression was: EPS
[L16 C36] #4 (Error) (After Macro Expansion) Argument 1 of `-' did not match any possible parameter type.
    The rejected type is PositiveInteger.
    Expected one of:
      -- DoubleFloat
      -- DoubleFloat, DoubleFloat
Expanded expression was: digits$DoubleFloat()

"/var/zope2/var/LatexWiki/sfx2.as", line 33: 
        if i>ITMAX or abs(del-1) < EPS then break
...................................^
[L33 C36] #2 (Error) (After Macro Expansion) No meaning for identifier `EPS'.
Expanded expression was: EPS

   The )library system command was not called after compilation.

Gamma(a,x)

(20)

Gamma(0,1)

(21)

Gamma(0,2)

(22)

Gamma(1,1)

(23)

Gamma(1,1.1)

(24)

Gamma(5,10)

(25)

Gamma(5,11)

(26)

Gamma(1,0)

(27)

Gamma(7,0)

(28)

  Gam(7,0)
              0.0
                                             Type: Float

Gam(7,0)

(29)

Gam(7,0.1)

(30)

Gam(7,0.2)

(31)

)set functions compile on

j:=120

(32)

nume(a) == cons(1 :: Float,[((a-i)*i) :: Float for i in 1..]);

dene(a,x) == [(x+2*i+1-a) :: Float for i in 0..];

cfe(a,x) == continuedFraction(0,nume(a),dene(a,x));

ccfe(a,x) == convergents cfe(a,x);

gamcfe(a,x) == exp(-x)*x^a*ccfe(a,x).j;

gamcfe(2,3)

Compiling function nume with type PositiveInteger -> Stream Float

Compiling function dene with type (PositiveInteger,PositiveInteger)
       -> Stream Float

Compiling function cfe with type (PositiveInteger,PositiveInteger)
       -> ContinuedFraction Float

Compiling function ccfe with type (PositiveInteger,PositiveInteger)
       -> Stream Fraction Float 
   There are 34 exposed and 23 unexposed library operations named * 
      having 2 argument(s) but none was determined to be applicable. 
      Use HyperDoc Browse, or issue
                                )display op *
      to learn more about the available operations. Perhaps 
      package-calling the operation or using coercions on the arguments
      will allow you to apply the operation.
   Cannot find a definition or applicable library operation named * 
      with argument type(s) 
                             Expression Integer
                               Fraction Float

      Perhaps you should use "@" to indicate the required return type, 
      or "$" to specify which version of the function you need.
   FriCAS will attempt to step through and interpret the code.

(33)

E1(x)== gamcfe(0,x)

E1(2.0)

   Cannot compile map: gamcfe
   We will attempt to interpret the code.

Compiling function nume with type NonNegativeInteger -> Stream Float

Compiling function dene with type (NonNegativeInteger,Float) -> 
      Stream Float

Compiling function cfe with type (NonNegativeInteger,Float) -> 
      ContinuedFraction Float

Compiling function ccfe with type (NonNegativeInteger,Float) -> 
      Stream Fraction Float

(34)

ff:Fraction Float

   Fraction Float is not a valid type.

Something similar happens if the argument is Fraction Integer

nume(a) == cons(1,[((a-i)*i) for i in 1..]);

   Compiled code for nume has been cleared.
   Compiled code for cfe has been cleared.
   Compiled code for ccfe has been cleared.
   Compiled code for gamcfe has been cleared.
   Compiled code for E1 has been cleared.
   1 old definition(s) deleted for function or rule nume

dene(a,x) == [(x+2*i+1-a) for i in 0..];

   Compiled code for dene has been cleared.
   1 old definition(s) deleted for function or rule dene

cfe(a,x) == continuedFraction(0,nume(a),dene(a,x));

   1 old definition(s) deleted for function or rule cfe

ccfe(a,x) == convergents cfe(a,x);

   1 old definition(s) deleted for function or rule ccfe

ccfe(0,2::Float)

Compiling function nume with type NonNegativeInteger -> Stream 
      Integer

Compiling function dene with type (NonNegativeInteger,Float) -> 
      Stream Float

Compiling function cfe with type (NonNegativeInteger,Float) -> 
      ContinuedFraction Float

Compiling function ccfe with type (NonNegativeInteger,Float) -> 
      Stream Fraction Float

(35)

ccfe(0,2::Fraction Integer)

Compiling function dene with type (NonNegativeInteger,Fraction 
      Integer) -> Stream Fraction Integer

Compiling function cfe with type (NonNegativeInteger,Fraction 
      Integer) -> ContinuedFraction Fraction Integer

Compiling function ccfe with type (NonNegativeInteger,Fraction 
      Integer) -> Stream Fraction Fraction Integer

(36)

Fraction Fraction Integer is also nonesense.