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 Gam(0,1.1)
Error signalled from user code in function Gam: a too large, ITMAX too small Gam(1,1)

\label{eq1}0.3678794411_7144232159(1)
Type: Float
fricas
Gam(1,1.1)

\label{eq2}0.3328710836_9807955329(2)
Type: Float
fricas
Gam(5,10)

\label{eq3}0.7020645138_4706574415(3)
Type: Float
fricas
Gam(5,11)

\label{eq4}0.3625104156_5228203538(4)
Type: Float
fricas
Gam(7,0)

\label{eq5}720.0000000000_0000000000_0008(5)
Type: Float

Let's try that again with more precision:

fricas
digits(100)

\label{eq6}20(6)
fricas
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)

\label{eq7}0.3678794411_7144232159_5523770161_4608674458_1113103176_7834
507836_8016974614_9574489980_3357147274_3459196437(7)
Type: Float
fricas
Gam(1,1.1)

\label{eq8}0.3328710836_9807955328_8846906431_3155216124_7952156921_2491
793331_3867507470_8541284431_1612617072_7005478519(8)
Type: Float
fricas
Gam(5,10)

\label{eq9}0.7020645138_4706574414_6387196628_3546367191_6532623256_0684
622278_6705870550_5584357048_3474646670_2985365058(9)
Type: Float
fricas
Gam(5,11)

\label{eq10}0.3625104156_5228203538_0753904311_4079803866_4530925132_7036
797697_4190490374_2658968752_0305953551_1648548436(10)
Type: Float
fricas
Gam(7,0)

\label{eq11}719.9999999999_9999999999_9999999999_9999999999_9999999999_99
99999999_9999999999_9999999999_9999999999_9999999999_9996(11)
Type: Float

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.01 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 14 JUL 2013 10:47:12 AM):
; /var/aw/var/LatexWiki/SFX.NRLIB/SFX.fasl written ; compilation finished in 0:00:00.042 ------------------------------------------------------------------------ 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{eq12}\Gamma \left({a , \: x}\right)(12)
Type: Expression(Integer)
fricas
Gamma(0,1)$SFX

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

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

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

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

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

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

\label{eq19}720.0000000000001(19)
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
      AXIOM-XL compiler and options 
-O -Fasy -Fao -Flsp -laxiom -Mno-ALDOR_W_WillObsolete -DAxiom -Y $AXIOM/algebra -I $AXIOM/algebra
      Use the system command )set compiler args to change these 
      options.
#1 (Warning) Could not use archive file `libaxiom.al'.
#2 (Warning) Could not use archive file `libaxiom.al'.
"/usr/local/aldor/linux/1.1.0/include/axiom.as", line 4: 
import from AxiomLib;
............^
[L4 C13] #3 (Error) No meaning for identifier `AxiomLib'.
"/usr/local/aldor/linux/1.1.0/include/axiom.as", line 15: import { true: %, false: % } from Boolean; ..................................^ [L15 C35] #4 (Error) No meaning for identifier `Boolean'.
"/usr/local/aldor/linux/1.1.0/include/axiom.as", line 17: string: Literal -> %; ........................^.......^ [L17 C25] #5 (Error) No meaning for identifier `Literal'. [L17 C33] #6 (Error) There are no suitable meanings for the operator `->'.
"/usr/local/aldor/linux/1.1.0/include/axiom.as", line 18: } from String; .......^ [L18 C8] #8 (Error) No meaning for identifier `String'.
"/var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/sfx2.as", line 8: Gamma:(DoubleFloat,DoubleFloat)->DoubleFloat ...........^.......................^.^ [L8 C12] #9 (Error) (After Macro Expansion) No meaning for identifier `DoubleFloat'. Expanded expression was: DoubleFloat [L8 C36] #11 (Error) (After Macro Expansion) There are no suitable meanings for the operator `->'. Expanded expression was: -> [L8 C38] #10 (Error) (After Macro Expansion) No meaning for identifier `DoubleFloat'. Expanded expression was: DoubleFloat [L8 C38] #13 (Fatal Error) (After Macro Expansion) Too many errors (use `-M emax=n' or `-M no-emax' to change the limit). Expanded expression was: ()
The )library system command was not called after compilation.

fricas
Gamma(a,x)

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

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

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

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

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

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

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

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

\label{eq28}\Gamma \left({7, \: 0}\right)(28)
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{eq29}719.9999999999_9999999999_9999999999_9999999999_9999999999_99
99999999_9999999999_9999999999_9999999999_9999999999_9996(29)
Type: Float
fricas
Gam(7,0.1)

\label{eq30}719.9999999869_1035963050_9717349012_8710490499_2980729057_24
68299528_7297035908_3400303838_7582842912_4885996(30)
Type: Float
fricas
Gam(7,0.2)

\label{eq31}719.9999984646_1597708521_8246906824_6256588795_3666381154_46
81930285_0299928361_6981388256_1798599100_3380769(31)
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{eq32}120(32)
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 35 exposed and 33 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{eq33}0.1991482734_7145577188(33)
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{eq34}0.0489005107_0806111956(34)
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{eq35}\begin{array}{@{}l}
\displaystyle
\left[{0.0}, \:{0.3333333333_3333333333}, \:{0.3571428571_428
5714286}, \: \right.
\
\
\displaystyle
\left.{0.3604651162_7906976744}, \:{0.3611111111_1111111111}, \: \right.
\
\
\displaystyle
\left.{0.3612656467_3157162726}, \:{0.3613083856_8697077301}, \: \right.
\
\
\displaystyle
\left.{0.3613215638_624830248}, \:{0.3613259896_0595485929}, \: \right.
\
\
\displaystyle
\left.{0.3613275827_1362214001}, \:...\right] 
(35)
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{eq36}\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}}, \:...\right](36)
Type: Stream(Fraction(Fraction(Integer)))

Fraction Fraction Integer is also nonesense.




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