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)

(1) -> 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

Let's try that again with more precision:

digits(100)

(1)

Gam(0,1)

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

Let us 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
  -- near the smallest representable number, can not make it smaller because
  -- otherwise 1/FPMIN overflows
  FPMIN ==> 1.0e-308::DoubleFloat

  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 + 1))$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/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/52052831343972396-25px004.spad
      using old system compiler.
   SFX abbreviates package SpecialFunctionsExtended 
------------------------------------------------------------------------
   initializing NRLIB SFX for SpecialFunctionsExtended 
   compiling into NRLIB SFX 
   compiling exported Gamma : (DoubleFloat,DoubleFloat) -> DoubleFloat
Time: 0 SEC.

(time taken in buildFunctor:  0)

;;;     ***       |SpecialFunctionsExtended| REDEFINED

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

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

   Cumulative Statistics for Constructor SpecialFunctionsExtended
      Time: 0.01 seconds

   finalizing NRLIB SFX 
   Processing SpecialFunctionsExtended for Browser database:
--->-->SpecialFunctionsExtended(constructor): Not documented!!!!
--->-->SpecialFunctionsExtended((Gamma ((DoubleFloat) (DoubleFloat) (DoubleFloat)))): Not documented!!!!
--->-->SpecialFunctionsExtended(): Missing Description
; compiling file "/var/aw/var/LatexWiki/SFX.NRLIB/SFX.lsp" (written 22 APR 2025 09:04:22 AM):

; wrote /var/aw/var/LatexWiki/SFX.NRLIB/SFX.fasl
; compilation finished in 0:00:00.024
------------------------------------------------------------------------
   SpecialFunctionsExtended is now explicitly exposed in frame initial 
   SpecialFunctionsExtended will be automatically loaded when needed 
      from /var/aw/var/LatexWiki/SFX.NRLIB/SFX

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)

(2)

Gamma(0,1)$SFX

(3)

Gamma(0,2)$SFX

(4)

Gamma(1,1)$SFX

(5)

Gamma(1,1.1)$SFX

(6)

Gamma(5,10)$SFX

(7)

Gamma(5,11)$SFX

(8)

Gamma(7,0)$SFX

(9)

Note: we need package call because otherwise we get builtin function, either from Expression(Integer) or (unimplemented) from DoubleFloat?.

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-308::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/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/sfx2.as using
      Aldor compiler and options 
-O -Fasy -Fao -Flsp -lfricas -Mno-ALDOR_W_WillObsolete -DFriCAS -Y $FRICAS/algebra -I $FRICAS/algebra
      Use the system command )set compiler args to change these 
      options.
"/var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/sfx2.as", line 1: 
#include "axiom.as";
^
[L1 C1] #1 (Error) Could not open file `axiom.as'.

   The )library system command was not called after compilation.

Gamma(a,x)

(10)

Gamma(0,1)

(11)

Gamma(0,2)

(12)

Gamma(1,1)

(13)

Gamma(1,1.1)

(14)

Gamma(5,10)

(15)

Gamma(5,11)

(16)

Gamma(1,0)

(17)

Gamma(7,0)

(18)

  Gam(7,0)
              0.0
                                             Type: Float

Gam(7,0)

(19)

Gam(7,0.1)

(20)

Gam(7,0.2)

(21)

)set functions compile on

j:=120

(22)

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 31 exposed and 40 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.

(23)

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

(24)

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

(25)

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

(26)

Fraction Fraction Integer is also nonesense.