help links subscribe changes refresh edit

	SandBox
SandBox Functions ←	↑	→ SandBox Integration

SandBox Gamma

Edit detail for SandBox Gamma revision 1 of 2

	1 2
Editor: Time: 2007/11/18 18:03:02 GMT-8
Note: Fraction Float?

changed:
-
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)
\begin{axiom}
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
\end{axiom}

\begin{axiom}
Gam(0,1)
Gam(0,1.1)
Gam(1,1)
Gam(1,1.1)
Gam(5,10)
Gam(5,11)
Gam(7,0)
\end{axiom}

Let's try that again with more precision:
\begin{axiom}
digits(100)
Gam(0,1)
Gam(0,1.1)
Gam(1,1)
Gam(1,1.1)
Gam(5,10)
Gam(5,11)
Gam(7,0)
\end{axiom}

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

\begin{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
\end{spad}

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.

\begin{axiom}
Gamma(a,x)
Gamma(0,1)
Gamma(0,2)
Gamma(1,1)
Gamma(1,1.1)
Gamma(5,10)
Gamma(5,11)
Gamma(7,0)
\end{axiom}

Exactly the same code can also be compiled using Aldor (Axiom library
compiler version 2):
\begin{aldor}[sfx2]
#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
\end{aldor}

\begin{axiom}
Gamma(a,x)
Gamma(0,1)
Gamma(0,2)
Gamma(1,1)
Gamma(1,1.1)
Gamma(5,10)
Gamma(5,11)
Gamma(1,0)
Gamma(7,0)
\end{axiom}

From greg Sun Jan 29 18:12:10 -0600 2006
From: greg
Date: Sun, 29 Jan 2006 18:12:10 -0600
Subject: 
Message-ID: <20060129181210-0600@wiki.axiom-developer.org>

This is false (result is 720)::

  Gam(7,0)
              0.0
                                             Type: Float

From BillPage Sun Jan 29 19:42:30 -0600 2006
From: Bill Page
Date: Sun, 29 Jan 2006 19:42:30 -0600
Subject: Better?
Message-ID: <20060129194230-0600@wiki.axiom-developer.org>

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.

\begin{axiom}
Gam(7,0)
Gam(7,0.1)
Gam(7,0.2)
\end{axiom}

From greg Mon Jan 30 09:58:09 -0600 2006
From: greg
Date: Mon, 30 Jan 2006 09:58:09 -0600
Subject: yigal's code
Message-ID: <20060130095809-0600@wiki.axiom-developer.org>

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
\begin{axiom}
)set functions compile on
\end{axiom}

\begin{axiom}
j:=120
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)
E1(x)== gamcfe(0,x)
E1(2.0)
\end{axiom}


From BillPage Thu Feb 2 11:23:08 -0600 2006
From: Bill Page
Date: Thu, 02 Feb 2006 11:23:08 -0600
Subject: Fraction Float?
Message-ID: <20060202112308-0600@wiki.axiom-developer.org>

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.
\begin{axiom}
ff:Fraction Float
\end{axiom}

Something similar happens if the argument is Fraction Integer
\begin{axiom}
nume(a) == cons(1,[((a-i)*i) for i in 1..]);
dene(a,x) == [(x+2*i+1-a) for i in 0..];
cfe(a,x) == continuedFraction(0,nume(a),dene(a,x));
ccfe(a,x) == convergents cfe(a,x);
ccfe(0,2::Float)
ccfe(0,2::Fraction Integer)
\end{axiom}

'Fraction Fraction Integer' is also nonesense.

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)

axiom

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

axiom

Gam(0,1)

axiom

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.36787944117144232159$

(1)

Type: Float

axiom

Gam(1,1.1)

$\label{eq2}0.33287108369807955329$

(2)

Type: Float

axiom

Gam(5,10)

$\label{eq3}0.70206451384706574415$

(3)

Type: Float

axiom

Gam(5,11)

$\label{eq4}0.36251041565228203538$

(4)

Type: Float

axiom

Gam(7,0)

$\label{eq5}720.000000000000000000000008$

(5)

Type: Float

Let's try that again with more precision:

axiom

digits(100)

$\label{eq6}20$

(6)

Type: PositiveInteger?

axiom

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}\begin{array}{@{}l}\displaystyle 0.36787944117144232159552377016146086744581113103176 \ 78345078368016974614957448998033571472743459196437$

(7)

Type: Float

axiom

Gam(1,1.1)

$\label{eq8}\begin{array}{@{}l}\displaystyle 0.33287108369807955328884690643131552161247952156921 \ 24917933313867507470854128443116126170727005478519$

(8)

Type: Float

axiom

Gam(5,10)

$\label{eq9}\begin{array}{@{}l}\displaystyle 0.70206451384706574414638719662835463671916532623256 \ 06846222786705870550558435704834746466702985365058$

(9)

Type: Float

axiom

Gam(5,11)

$\label{eq10}\begin{array}{@{}l}\displaystyle 0.36251041565228203538075390431140798038664530925132 \ 70367976974190490374265896875203059535511648548436$

(10)

Type: Float

axiom

Gam(7,0)

$\label{eq11}\begin{array}{@{}l}\displaystyle 719.9999999999999999999999999999999999999999999999999 \ 9999999999999999999999999999999999999999999999999999 \ 996$

(11)

Type: Float

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

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

spad

   Compiling FriCAS source code from file 
      /var/zope2/var/LatexWiki/6599990542159375921-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
****** comp fails at level 3 with expression: ******
error in function Gamma 

(SEQ
 (IF (< |x| 0)
     (|error| #1="Invalid arguments")
     (IF (< |a| 0)
         (|error| #1#)
         |noBranch|))
 (IF (= |x| 0)
     (|return| 1 (|Gamma| |a|))
     |noBranch|)
 (LET EPS
   | << |
   ((|elt| (|DoubleFloat|) **)
    (|::| ((|elt| (|Float|) |float|) 10 0 10) (|DoubleFloat|))
    (- ((|elt| (|DoubleFloat|) |digits|))))
   | >> |)
 (LET (|:| |b| (|DoubleFloat|))
   (- (+ |x| 1) |a|))
 (LET (|:| |c| (|DoubleFloat|))
   (/ 1 (|::| ((|elt| (|Float|) |float|) 1 -323 10) (|DoubleFloat|))))
 (LET (|:| |d| (|DoubleFloat|))
   (/ 1 |b|))
 (LET (|:| |h| (|DoubleFloat|))
   |d|)
 (LET (|:| |i| (|DoubleFloat|))
   1)
 (REPEAT
  (SEQ
   (LET (|:| |an| (|DoubleFloat|))
     (- (* |i| (- |i| |a|))))
   (LET |b|
     (+ |b| (|::| ((|elt| (|Float|) |float|) 2 0 10) (|DoubleFloat|))))
   (LET |d|
     (+ (* |an| |d|) |b|))
   (SEQ
    (LET #2=#:G725
      (< (|abs| |d|)
         (|::| ((|elt| (|Float|) |float|) 1 -323 10) (|DoubleFloat|))))
    (|exit| 1
     (IF #2#
         (LET |d|
           (|::| ((|elt| (|Float|) |float|) 1 -323 10) (|DoubleFloat|)))
         |noBranch|)))
   (LET |c|
     (+ |b| (/ |an| |c|)))
   (SEQ
    (LET #3=#:G726
      (< (|abs| |c|)
         (|::| ((|elt| (|Float|) |float|) 1 -323 10) (|DoubleFloat|))))
    (|exit| 1
     (IF #3#
         (LET |c|
           (|::| ((|elt| (|Float|) |float|) 1 -323 10) (|DoubleFloat|)))
         |noBranch|)))
   (LET |d|
     (/ 1 |d|))
   (LET (|:| |del| (|DoubleFloat|))
     (* |d| |c|))
   (LET |h|
     (* |h| |del|))
   (SEQ
    (LET #4=#:G727
      (< (|::| ((|elt| (|Float|) |float|) 100 0 10) (|DoubleFloat|)) |i|))
    (|exit| 1
     (IF #4#
         (|leave| 1 |$NoValue|)
         (SEQ
          (LET #5=#:G728
            (< (|abs| (- |del| 1)) EPS))
          (|exit| 1
           (IF #5#
               (|leave| 1 |$NoValue|)
               |noBranch|))))))
   (|exit| 1
    (LET |i|
      (+ |i| 1)))))
 (SEQ
  (LET #6=#:G729
    (< (|::| ((|elt| (|Float|) |float|) 100 0 10) (|DoubleFloat|)) |i|))
  (|exit| 1
   (IF #6#
       (|error| "a too large, ITMAX too small")
       |noBranch|)))
 (|exit| 1 (|return| 1 (* (|exp| (+ (- |x|) (* (|log| |x|) |a|))) |h|))))
****** level 3  ******
$x:= ((elt (DoubleFloat) **) (:: ((elt (Float) float) 10 (Zero) 10) (DoubleFloat)) (- ((elt (DoubleFloat) digits))))
$m:= $EmptyMode
$f:=
((((|x| # . #1=#) (|a| # #) (|x| . #1#) (* #) ...)))

   >> Apparent user error:
   Cannot coerce (call (ELT $ 12) 10 (call (XLAM ignore 0)) 10) 
      of mode (Float) 
      to mode (DoubleFloat)

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.

axiom

Gamma(a,x)

$\label{eq12}\Gamma \left({a , \: x}\right)$

(12)

Type: Expression(Integer)

axiom

Gamma(0,1)

$\label{eq13}\Gamma \left({0, \: 1}\right)$

(13)

Type: Expression(Integer)

axiom

Gamma(0,2)

$\label{eq14}\Gamma \left({0, \: 2}\right)$

(14)

Type: Expression(Integer)

axiom

Gamma(1,1)

$\label{eq15}\Gamma \left({1, \: 1}\right)$

(15)

Type: Expression(Integer)

axiom

Gamma(1,1.1)

$\label{eq16}\Gamma \left({{1.0}, \:{1.1}}\right)$

(16)

Type: Expression(Float)

axiom

Gamma(5,10)

$\label{eq17}\Gamma \left({5, \:{10}}\right)$

(17)

Type: Expression(Integer)

axiom

Gamma(5,11)

$\label{eq18}\Gamma \left({5, \:{11}}\right)$

(18)

Type: Expression(Integer)

axiom

Gamma(7,0)

$\label{eq19}\Gamma \left({7, \: 0}\right)$

(19)

Type: Expression(Integer)

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

aldor

   Compiling FriCAS source code from file 
      /var/zope2/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.
"/var/zope2/var/LatexWiki/sfx2.as", line 16: 
      EPS := (10.0::DoubleFloat**(-digits()$DoubleFloat))$DoubleFloat  -- Relative accuracy
...................................^.....................^
[L16 C36] #3 (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()
[L16 C58] #2 (Error) (After Macro Expansion) There are no suitable meanings for the operator `**$DoubleFloat'.
Expanded expression was: **$DoubleFloat

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

   The )library system command was not called after compilation.

axiom

Gamma(a,x)

$\label{eq20}\Gamma \left({a , \: x}\right)$

(20)

Type: Expression(Integer)

axiom

Gamma(0,1)

$\label{eq21}\Gamma \left({0, \: 1}\right)$

(21)

Type: Expression(Integer)

axiom

Gamma(0,2)

$\label{eq22}\Gamma \left({0, \: 2}\right)$

(22)

Type: Expression(Integer)

axiom

Gamma(1,1)

$\label{eq23}\Gamma \left({1, \: 1}\right)$

(23)

Type: Expression(Integer)

axiom

Gamma(1,1.1)

$\label{eq24}\Gamma \left({{1.0}, \:{1.1}}\right)$

(24)

Type: Expression(Float)

axiom

Gamma(5,10)

$\label{eq25}\Gamma \left({5, \:{10}}\right)$

(25)

Type: Expression(Integer)

axiom

Gamma(5,11)

$\label{eq26}\Gamma \left({5, \:{11}}\right)$

(26)

Type: Expression(Integer)

axiom

Gamma(1,0)

$\label{eq27}\Gamma \left({1, \: 0}\right)$

(27)

Type: Expression(Integer)

axiom

Gamma(7,0)

$\label{eq28}\Gamma \left({7, \: 0}\right)$

(28)

Type: Expression(Integer)

... --greg, Sun, 29 Jan 2006 18:12:10 -0600 reply

This is false (result is 720):

  Gam(7,0)
              0.0
                                             Type: Float

Better? --Bill Page, Sun, 29 Jan 2006 19:42:30 -0600 reply

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

to compute $\Gamma$ since it converges more rapidly.

axiom

Gam(7,0)

$\label{eq29}\begin{array}{@{}l}\displaystyle 719.9999999999999999999999999999999999999999999999999 \ 9999999999999999999999999999999999999999999999999999 \ 996$

(29)

Type: Float

axiom

Gam(7,0.1)

$\label{eq30}\begin{array}{@{}l}\displaystyle 719.9999999869103596305097173490895137595484268324346 \ 065775193165312727417661992245691022941558764196$

(30)

Type: Float

axiom

Gam(7,0.2)

$\label{eq31}\begin{array}{@{}l}\displaystyle 719.9999984646159770852182469157013222705579469380749 \ 722296525136047137980713842586005969219440451807$

(31)

Type: Float

yigal's code --greg, Mon, 30 Jan 2006 09:58:09 -0600 reply

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

axiom

)set functions compile on

axiom

j:=120

$\label{eq32}120$

(32)

Type: PositiveInteger?

axiom

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

Type: Void

axiom

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

Type: Void

axiom

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

Type: Void

axiom

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

Type: Void

axiom

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

Type: Void

axiom

gamcfe(2,3)

axiom

Compiling function nume with type PositiveInteger -> Stream(Float)

axiom

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

axiom

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

axiom

Compiling function ccfe with type (PositiveInteger,PositiveInteger)
       -> Stream(Fraction(Float)) 
   There are 35 exposed and 24 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}\begin{array}{@{}l}\displaystyle 0.19914827347145577191736966260024710652679836875369 \ 28622705109104242426709207982006162169769465333782$

(33)

Type: Expression(Float)

axiom

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

Type: Void

axiom

E1(2.0)

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

axiom

Compiling function nume with type NonNegativeInteger -> Stream(Float
      )

axiom

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

axiom

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

axiom

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

$\label{eq34}\begin{array}{@{}l}\displaystyle 0.04890051070806111956723982691434728982121544510421 \ 327725184171637798809149832755994923592819658821724$

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

axiom

ff:Fraction Float

   Fraction(Float) is not a valid type.

Something similar happens if the argument is Fraction Integer

axiom

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

axiom

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

axiom

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

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

Type: Void

axiom

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

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

Type: Void

axiom

ccfe(0,2::Float)

axiom

Compiling function nume with type NonNegativeInteger -> Stream(
      Integer)

axiom

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

axiom

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

axiom

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

$\label{eq35}\begin{array}{@{}l} \displaystyle \left[{0.0}, \: \right. \ \ \displaystyle \left.{ \begin{array}{@{}l}\displaystyle 0.33333333333333333333333333333333333333333333333333 \ 33333333333333333333333333333333333333333333333333$

(35)

Type: Stream(Fraction(Float))

axiom

ccfe(0,2::Fraction Integer)

axiom

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

axiom

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

axiom

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.