From octave-sources-request at bevo dot che dot wisc dot edu  Thu Apr  9 09:28:49 1998
Subject: RE: sanity test for gammai
From: (Ted Harding) <Ted dot Harding at nessie dot mcc dot ac dot uk>
To: Jim Van Zandt <jrv at mbunix dot mitre dot org>
Cc: jrv at vanzandt dot mv dot com, octave-sources@bevo.che.wisc.edu
Date: Thu, 09 Apr 1998 15:25:50 +0100 (BST)

On 09-Apr-98 Jim Van Zandt wrote:
> Octave 2.0.5 finds:
> [snip]
> which are correct.  (Unfortunately, as has been pointed out, 2.0.11
> results are not correct.)  

Recompiled with John Eaton's patches, they are the right way round now!

> I would like to propose the following for
> test/octave.test/arith/gammai-1.m:
> 
>     # reference: Abramowitz & Stegun section
>     # 6.5, where P(a,x) is the Octave
>     # function gammai(a,x)
>     x = [sqrt(pi)*erf(1), gammai(.5,1)*gamma(.5), gammainc(1,.5)*gamma(.5)];
>     v = [1.49364826562485, 1.49364826562485, 1.49364826562485];
>     all (abs (x - v) < sqrt(eps))
> 
> Unfortunately, octave 2.0.5 chokes on this (infinite loop or
> something).  I have not been able to determine why.

2.0.11 does not give "infinite loop" or other nasties. However, in 2.0.11 the
computed value of gammai(x,a) is simply the integral of (t^(a-1))*exp(-t), not
(as defined in A&S) of (t^(a-1))*exp(-t)/gamma(a) and also contrary to what
"help gammai" says (see below). With the appropriate change, the above test
becomes

>     x = [sqrt(pi)*erf(1), gammai(.5,1), gammainc(1,.5)];
>     v = [1.49364826562485, 1.49364826562485, 1.49364826562485];
>     all (abs (x - v) < sqrt(eps))

and this succeeds (result=1).

octave:130> gammai(6,100)
ans = 120.00   [ = 5! ].

One can live with either interpretation of gammai, but it should be made
definite which it is supposed to be.

For instance, Bender & Orszag, "Advanced Mathematical Methods for Scientists
and Engineers, McGraw-Hill define it as the integral; Jeffreys & Jeffreys also
do not divide by gamma(a); while Abramowitz & Stegun adopt integral divided by
gamma(a), as do Press, Flannery Teukolsky & Vetterling "Numerical recipes in C".
Whittaker and Watson also do not divide by gamma(a) and even do not normalise
erf(x) so that erf(inf)=sqrt(pi)/2.

Although A&S tends to be what people look up, the fact that conventions for
special functions tend to differ from author to author means that one should
always check which definition is used, and also verify what definition is used
in software.

Octave-2.0.11 gives

octave:130> gammai(6,100)
ans = 120.00   [ = 5! ].


In octave-1.1.1 I get

octave:6> gammai(6,100)
ans = 1

so there has been a change along the line (octave-2.0.8 seg-faults on gammai so
I can't test it there ... )

Best wishes,
Ted.


--------------------------------------------------------------------
E-Mail: (Ted Harding) <Ted dot Harding at nessie dot mcc dot ac dot uk>
Date: 09-Apr-98                                       Time: 15:25:50
--------------------------------------------------------------------