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
Cheers and thanks everyone for sharing how to get the best out of the
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