Re: Numerical accuracy of Hypergeometric2F1

• To: mathgroup at smc.vnet.net
• Subject: [mg55803] Re: Numerical accuracy of Hypergeometric2F1
• From: "Christos Argyropoulos M.D." <chrisarg at fuse.net>
• Date: Wed, 6 Apr 2005 03:12:02 -0400 (EDT)
• Sender: owner-wri-mathgroup at wolfram.com

```Hello,
Thanx everyone for pointing out the subtleties of making Mathematica
understand precision of input (I always thought that SetPrecision threads
over its first input). However I still think that the asymptotic expansion
(or whatever numerical method is used internally) used to obtain the
numerical value of Hypergeometric2F1 for large (floating point, arbitrary
precision) values of one of the parameters needs to be re-examined by WRI.
I am attaching an example where the parameters to hypergeometric2F1 were
arbitrary precision numbers (using SetPrecision before passing them to the
function).
When rational numbers or integers are used to represent input, Mathematica
gets it right (I believe that in those cases it symbolically expands the
hypergeometric in a polynomial and then evaluates the polynomial).
For the example attached  (i.e. Hypergeometric2F1[80+1/2,1,3/2,90/10]) and
floating point representations for [90+1/2, 90/10] one has to increase
precision to 25 digits to get an accurate result (and relying on machine
arithmetic is off the mark). By comparison, Forrey's implementation of the
Hypergeometric2F1 ( R. C. Forrey, Computing the Hypergeometric Function,
Journal of Computational Physics 137, 79-100 (1997). ) which uses machine
floating point arithmetic gives the following answer (also obtainable by
Mathematica for integer/rational input) :

0.9860651611217461E+89

Using a smaller first argument (i.e. 10) results in no "significant"
difference in the answers provided by Mathematica irrespective of the
encoding of the parameters (rational, machine arithmetic floating point,
arbitrary precision numbers). My interpretation of these results is that the
algorithm that provides the numerical answer needs some
tweaking/optimization

arbirtrary-precision point arithmetic with Mathematica

Christos Argyropoulos MD PhD

Mathematica Code
--------------------------------------------------------------------------------------------------------------------
In[1]:=
\!\(fpol[x_, k_] :=
Module[{}, coeff = Table[1, {k + 1}];
For[i = 2, i \[LessEqual] k + 1, \(i++\),
coeff[\([i]\)] =
coeff[\([i -
1]\)]\ \(\(-k\) + i - 1\)\/\(i - k - 1/2\)]; \
{Hypergeometric2F1[k + 1/2, 1, 3/2,
x], \(1\/\(2  k -
1\)\) \(\[Sum]\+\(i = 1\)\%\(k + 1\)\((\((1 - x)\)\^\(-i\)\
\
coeff[\([i]\)])\)\)}]\[IndentingNewLine]
fhorn[x_, k_] :=
Module[{}, coeff = Table[1, {k + 1}];
For[i = 2, i \[LessEqual] k + 1, \(i++\),
coeff[\([i]\)] =
coeff[\([i -
1]\)]\ \(\(-k\) + i - 1\)\/\(i - k - 1/2\)]; \
{Hypergeometric2F1[k + 1/2, 1, 3/2, x], \(1\/\(2  k - 1\)\)
Module[{y = 1/\((1 - x)\), res = 1}, res = 0;
For[i = k + 1, i \[GreaterEqual] \ 1\ , \(i--\),
res = \((res + coeff[\([i]\)])\)\ y\ ]; res]}]\)
In[3]:=
x=9/10;
k=90.;
workprec=16;
printprec=16;
In[7]:=

Print["Workprec = ",workprec];
fpol[x,k]//SetPrecision[#,printprec]&
fpol[N[x],k]//SetPrecision[#,16]&
fpol[SetPrecision[x,workprec],k]//SetPrecision[#,printprec]&
fpol[SetPrecision[x,workprec],SetPrecision[k,workprec]]//
SetPrecision[#,printprec]&
fpol[x,Rationalize[k]]//SetPrecision[#,printprec]&
workprec=20;
Print["Workprec = ",workprec];
fpol[x,k]//SetPrecision[#,printprec]&
fpol[N[x],k]//SetPrecision[#,16]&
fpol[SetPrecision[x,workprec],k]//SetPrecision[#,printprec]&
fpol[SetPrecision[x,workprec],SetPrecision[k,workprec]]//
SetPrecision[#,printprec]&
fpol[x,Rationalize[k]]//SetPrecision[#,printprec]&
workprec=25;
Print["Workprec = ",workprec];
fpol[x,k]//SetPrecision[#,printprec]&
fpol[N[x],k]//SetPrecision[#,16]&
fpol[SetPrecision[x,workprec],k]//SetPrecision[#,printprec]&
fpol[SetPrecision[x,workprec],SetPrecision[k,workprec]]//
SetPrecision[#,printprec]&
fpol[x,Rationalize[k]]//SetPrecision[#,printprec]&
From In[7]:=
Workprec = \[InvisibleSpace]16
Out[8]=
\!\({7.20338647551431814577488138862587`15.9546*^95,
9.860651611217270453862251879131`15.9546*^88}\)
Out[9]=
\!\({7.20338647551431814577488138862587`15.9546*^95,
9.860651611217469754211158880883`15.9546*^88}\)
Out[10]=
\!\({7.20338647551431814577488138862587`15.9546*^95,
9.860651611217466927255855235469`15.9546*^88}\)
Out[11]=
\!\({9.990653929442517915012215795621`15.9546*^88,
9.860651611217461201585792023911`15.9546*^88}\)
Out[12]=
\!\({9.860651611217264799951644588301`15.9546*^88,
9.860651611217264799951644588301`15.9546*^88}\)
From In[7]:=
Workprec = \[InvisibleSpace]20
Out[15]=
\!\({7.20338647551431814577488138862587`15.9546*^95,
9.860651611217270453862251879131`15.9546*^88}\)
Out[16]=
\!\({7.20338647551431814577488138862587`15.9546*^95,
9.860651611217469754211158880883`15.9546*^88}\)
Out[17]=
\!\({7.20338647551431814577488138862587`15.9546*^95,
9.860651611217270453862251879131`15.9546*^88}\)
Out[18]=
\!\({9.8606020395819286190853862`15.9546*^88,
9.860651611217264267879523610443`15.9546*^88}\)
Out[19]=
\!\({9.860651611217264799951644588301`15.9546*^88,
9.860651611217264799951644588301`15.9546*^88}\)
From In[7]:=
Workprec = \[InvisibleSpace]25
Out[22]=
\!\({7.20338647551431814577488138862587`15.9546*^95,
9.860651611217270453862251879131`15.9546*^88}\)
Out[23]=
\!\({7.20338647551431814577488138862587`15.9546*^95,
9.860651611217469754211158880883`15.9546*^88}\)
Out[24]=
\!\({7.20338647551431814577488138862587`15.9546*^95,
9.860651611217270453862251879131`15.9546*^88}\)
Out[25]=
\!\({9.8606516120571455699381674`15.9546*^88,
9.8606516112172642678795376`15.9546*^88}\)
Out[26]=
\!\({9.860651611217264799951644588301`15.9546*^88,
9.860651611217264799951644588301`15.9546*^88}\)
In[27]:=
workprec=25;
x=9/10;
k=10.;
Print["Workprec = ",workprec];
fpol[x,k]//SetPrecision[#,printprec]&
fpol[N[x],k]//SetPrecision[#,16]&
fpol[SetPrecision[x,workprec],k]//SetPrecision[#,printprec]&
fpol[SetPrecision[x,workprec],SetPrecision[k,workprec]]//
SetPrecision[#,printprec]&
fpol[x,Rationalize[k]]//SetPrecision[#,printprec]&
From In[27]:=
Workprec = \[InvisibleSpace]25
Out[31]=
\!\({2.991232093254027843475341796875`15.9546*^9,
2.991232093254022121429443359375`15.9546*^9}\)
Out[32]=
\!\({2.991232093254027843475341796875`15.9546*^9,
2.9912320932540283203125`15.9546*^9}\)
Out[33]=
\!\({2.991232093254027843475341796875`15.9546*^9,
2.991232093254022121429443359375`15.9546*^9}\)
Out[34]=
\!\({2.9912320932540215202753902439`15.9546*^9,
2.9912320932540215202753902444`15.9546*^9}\)
Out[35]=
\!\({2.99123209325402164459228515625`15.9546*^9,
2.99123209325402164459228515625`15.9546*^9}\)
------------------------------------------------------------------------------------------------

Fortran Implementation
source for the hypergeometric function :
http://cfa-www.harvard.edu/~rcf/codes/hyp.f
test file source code:

c********************************************************************
c
c   program :  hyptest

c***********************************************************************
c

real*8 a,b,c,z,w,re,im
a=1.d0
b=90.5d0
c=1.5d0
z=0.9d0
call hyp(z,a,b,c,re,im)
write(*,900)re
900   format (g22.16)
1000  stop
end

compile as:
g77  hyp.f  hyptest.f

```

• Prev by Date: Re: Re: Mathematica bug in handling trigonometric functions? (and more)
• Next by Date: ploting functions slowly
• Previous by thread: Re: Re: Re[2]: Re: Numerical accuracy of Hypergeometric2F1
• Next by thread: Re: Numerical accuracy of Hypergeometric2F1