login  home  contents  what's new  discussion  bug reports     help  links  subscribe  changes  refresh  edit

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)

fricas
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.
Type: Void

fricas
Gam(0,1)
fricas
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:

fricas
digits(100)

\label{eq1}20(1)
Type: PositiveInteger?
fricas
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

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
spad
   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.02 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.02 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 04 APR 2022 05:20:29 PM):
; /var/aw/var/LatexWiki/SFX.NRLIB/SFX.fasl written ; compilation finished in 0:00:00.026 ------------------------------------------------------------------------ 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.

fricas
Gamma(a,x)

\label{eq2}\Gamma \left({a , \: x}\right)(2)
Type: Expression(Integer)
fricas
Gamma(0,1)$SFX

\label{eq3}0.21938393439551213(3)
Type: DoubleFloat?
fricas
Gamma(0,2)$SFX

\label{eq4}0.04890051070806015(4)
Type: DoubleFloat?
fricas
Gamma(1,1)$SFX

\label{eq5}0.36787944117144233(5)
Type: DoubleFloat?
fricas
Gamma(1,1.1)$SFX

\label{eq6}0.33287108369807955(6)
Type: DoubleFloat?
fricas
Gamma(5,10)$SFX

\label{eq7}0.7020645138470669(7)
Type: DoubleFloat?
fricas
Gamma(5,11)$SFX

\label{eq8}0.36251041565228215(8)
Type: DoubleFloat?
fricas
Gamma(7,0)$SFX

\label{eq9}720.0000000000001(9)
Type: DoubleFloat?

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

aldor
#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
aldor
   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.
   The )library system command was not called after compilation.

fricas
Gamma(a,x)

\label{eq10}\Gamma \left({a , \: x}\right)(10)
Type: Expression(Integer)
fricas
Gamma(0,1)

\label{eq11}\Gamma \left({0, \: 1}\right)(11)
Type: Expression(Integer)
fricas
Gamma(0,2)

\label{eq12}\Gamma \left({0, \: 2}\right)(12)
Type: Expression(Integer)
fricas
Gamma(1,1)

\label{eq13}\Gamma \left({1, \: 1}\right)(13)
Type: Expression(Integer)
fricas
Gamma(1,1.1)

\label{eq14}0.33287108369807955(14)
Type: DoubleFloat?
fricas
Gamma(5,10)

\label{eq15}\Gamma \left({5, \:{10}}\right)(15)
Type: Expression(Integer)
fricas
Gamma(5,11)

\label{eq16}\Gamma \left({5, \:{11}}\right)(16)
Type: Expression(Integer)
fricas
Gamma(1,0)

\label{eq17}\Gamma \left({1, \: 0}\right)(17)
Type: Expression(Integer)
fricas
Gamma(7,0)

\label{eq18}\Gamma \left({7, \: 0}\right)(18)
Type: Expression(Integer)

This is false (result is 720):
  Gam(7,0)
              0.0
                                             Type: Float

Sorry I forgot the limiting case. Actually the algoritm in Recipes is a little more sophisticated. It uses the series expansion of \gamma eq. (6.25) for x<a+1 to compute \Gamma since it converges more rapidly.

fricas
Gam(7,0)

\label{eq19}720.0000000000 \<u> 0000000000 \</u> 0008(19)
Type: Float
fricas
Gam(7,0.1)

\label{eq20}719.9999999869 \<u> 1035963050 \</u> 9717349012 \<u> 8710490499 \</u> 2
980729057 \<u> 2468299528 \</u> 7297035908 \<u> 3400303838 \</u> 758284
2912 \<u> 4885996(20)
Type: Float
fricas
Gam(7,0.2)

\label{eq21}719.9999984646 \<u> 1597708521 \</u> 8246906824 \<u> 6256588795 \</u> 3
666381154 \<u> 4681930285 \</u> 0299928361 \<u> 6981388256 \</u> 179859
9100 \<u> 3380769(21)
Type: Float

The command )set functions compile on is necessary in order to avoid a seg fault with certain function definitions in the Axiom interpreter. See: #196
fricas
)set functions compile on

fricas
j:=120

\label{eq22}120(22)
Type: PositiveInteger?
fricas
nume(a) == cons(1 :: Float,[((a-i)*i) :: Float for i in 1..]);
Type: Void
fricas
dene(a,x) == [(x+2*i+1-a) :: Float for i in 0..];
Type: Void
fricas
cfe(a,x) == continuedFraction(0,nume(a),dene(a,x));
Type: Void
fricas
ccfe(a,x) == convergents cfe(a,x);
Type: Void
fricas
gamcfe(a,x) == exp(-x)*x^a*ccfe(a,x).j;
Type: Void
fricas
gamcfe(2,3)
fricas
Compiling function nume with type PositiveInteger -> Stream(Float)
fricas
Compiling function dene with type (PositiveInteger, PositiveInteger)
       -> Stream(Float)
fricas
Compiling function cfe with type (PositiveInteger, PositiveInteger)
       -> ContinuedFraction(Float)
fricas
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.

\label{eq23}0.1991482734 \<u> 7145577188(23)
Type: Expression(Float)
fricas
E1(x)== gamcfe(0,x)
Type: Void
fricas
E1(2.0)
Cannot compile map: gamcfe We will attempt to interpret the code.
fricas
Compiling function nume with type NonNegativeInteger -> Stream(Float
      )
fricas
Compiling function dene with type (NonNegativeInteger, Float) -> 
      Stream(Float)
fricas
Compiling function cfe with type (NonNegativeInteger, Float) -> 
      ContinuedFraction(Float)
fricas
Compiling function ccfe with type (NonNegativeInteger, Float) -> 
      Stream(Fraction(Float))

\label{eq24}0.0489005107 \<u> 0806111956(24)
Type: Fraction(Float)

Fraction Float? --Bill Page, Thu, 02 Feb 2006 11:23:08 -0600 reply
There seems to be a problem with Axiom's ContinuedFraction? domain. The type of the result is shown as Fraction Float but this is nonesense.
fricas
ff:Fraction Float
Fraction(Float) is not a valid type.

Something similar happens if the argument is Fraction Integer

fricas
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
Type: Void
fricas
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
Type: Void
fricas
cfe(a,x) == continuedFraction(0,nume(a),dene(a,x));
1 old definition(s) deleted for function or rule cfe
Type: Void
fricas
ccfe(a,x) == convergents cfe(a,x);
1 old definition(s) deleted for function or rule ccfe
Type: Void
fricas
ccfe(0,2::Float)
fricas
Compiling function nume with type NonNegativeInteger -> Stream(
      Integer)
fricas
Compiling function dene with type (NonNegativeInteger, Float) -> 
      Stream(Float)
fricas
Compiling function cfe with type (NonNegativeInteger, Float) -> 
      ContinuedFraction(Float)
fricas
Compiling function ccfe with type (NonNegativeInteger, Float) -> 
      Stream(Fraction(Float))

\label{eq25}\begin{array}{@{}l}
\displaystyle
\left[{0.0}, \:{0.3333333333 \<u> 3333333333}, \:{0.3571428571 \</u> 4285714286}, \: \right.
\
\
\displaystyle
\left.{0.3604651162 \<u> 7906976744}, \:{0.3611111111 \</u> 111111
1111}, \: \right.
\
\
\displaystyle
\left.{0.3612656467 \<u> 3157162726}, \:{0.3613083856 \</u> 869707
7301}, \: \right.
\
\
\displaystyle
\left.{0.3613215638 \<u> 624830248}, \:{0.3613259896 \</u> 0595485
929}, \: \right.
\
\
\displaystyle
\left.{0.3613275827 \_ 1362214001}, \: \ldots \right] 
(25)
Type: Stream(Fraction(Float))
fricas
ccfe(0,2::Fraction Integer)
fricas
Compiling function dene with type (NonNegativeInteger, Fraction(
      Integer)) -> Stream(Fraction(Integer))
fricas
Compiling function cfe with type (NonNegativeInteger, Fraction(
      Integer)) -> ContinuedFraction(Fraction(Integer))
fricas
Compiling function ccfe with type (NonNegativeInteger, Fraction(
      Integer)) -> Stream(Fraction(Fraction(Integer)))

\label{eq26}\left[ 0, \:{1 \over 3}, \:{5 \over{14}}, \:{{31}\over{86}}, \:{{13}\over{36}}, \:{{1039}\over{2876}}, \:{{5291}\over{1464
4}}, \:{{20221}\over{55964}}, \:{{193003}\over{534152}}, \:{{3
85207}\over{1066088}}, \: \ldots \right](26)
Type: Stream(Fraction(Fraction(Integer)))

Fraction Fraction Integer is also nonesense.




  Subject: (replying)   Be Bold !!
  ( 15 subscribers )  
Please rate this page: