upiterbarg / mpmath

Automatically exported from code.google.com/p/mpmath
Other
0 stars 0 forks source link

problem with ellipk / hyp2f1 #201

Closed GoogleCodeExporter closed 9 years ago

GoogleCodeExporter commented 9 years ago
What steps will reproduce the problem?
1. ellipk(eps^2) for eps <> 0, 0.5 gives incorrect results for
   extended precision
2. ellipk(eps^2) is identical to pi/2*hyp2f1(0.5,0.5,1,my_epsi*my_epsi)
   for all 58 values of eps from 0,0.0001,....,0.99999 I've tested
3. so I assume you use hyp2f1 to calculate ellipk 

What is the expected output? What do you see instead?
I've obtained the expected output by two different means: 1.) using Maxima 
and 2.) an Excel-Spreadsheet with the xNumbers-AddIn implementing the 
AGM-Method to calculate K(k). Both outputs agree to within the first 60 of 64 
digits for all tested values of eps, whereas the mpmath output differs after 
the first 14(27) digits in the worst(best) case of big(small) eps.

Examples:
eps = 0.3
"correct" K(k=0.09)=
1,608048619930512801267207222238687157112176728802652558496349254
"incorrect" mpmath ellipk(0.09) = pi/2*hyp2f1(0.5,0.5,1,0.09)=
1,608048619930512799813156244396991240851380691072008068964509435

eps = 0.9
"correct" K(k=0.81)=
2,280549138422770204613751944555530438743237966278793336928341064
"incorrect" mpmath ellipk(0.81) = pi/2*hyp2f1(0.5,0.5,1,0.81)=
2,280549138422770332454780412866039777284362348618189389837282963

eps = 0.5
"correct" K(k=0.25)=
1,685750354812596042871203657799076989500800894141089044119948298
"correct" mpmath ellipk(0.25) = pi/2*hyp2f1(0.5,0.5,1,0.25)=
1,685750354812596042871203657799076989500800894141089044119948298

What version of the product are you using? On what operating system?
Py 2.6.6, mpmath 0.16, sympy 0.6.7, gmpy 1.13 with WinXP 32

Regards

Original issue reported on code.google.com by dr...@web.de on 26 Oct 2010 at 9:38

GoogleCodeExporter commented 9 years ago
Your inputs only have double precision (where 0 and 0.5 happen to be exact).

You need to convert decimal literals to full-precision numbers with mpf:

>>> print ellipk(mpf('0.09'))
1.608048619930512801267207222238687157112176728802652558496349254

Or implicitly:

>>> print ellipk('0.09')
1.608048619930512801267207222238687157112176728802652558496349254

Original comment by fredrik....@gmail.com on 26 Oct 2010 at 10:05

GoogleCodeExporter commented 9 years ago
oh thank you; my code was (multiple lines in this comment form only!):
<---snip--->
for my_epsi in epsi:
    epsi_sq=my_epsi*my_epsi
    my_str = str(my_epsi) + ';x'+ str(ellipk(epsi_sq)) + ';x' + 
             str(pi/2*hyp2f1(0.5,0.5,1.0,epsi_sq)) + '\n'
<---snip--->

I changed that to :
<---snip--->
for my_epsi in epsi:
    mpf_epsi_sq=mpf(my_epsi)*mpf(my_epsi)
    my_str = str(my_epsi) + ';x'+ str(ellipk(mpf_epsi_sq)) + ';x' + 
             str(pi/mpf(2.0)*hyp2f1(mpf(0.5),mpf(0.5),mpf(1.0),mpf_epsi_sq)) + '\n'
<---snip--->
and that did not help!? I then realised that mpf() only works on strings... so 
after I changed the definiton of the list epsi from
epsi=[0,0.00001,0.00002,0.00005,...,0.99999]
to
epsi=['0.0','0.00001','0.00002','0.00005',...,'0.99999']
all is fine now... newbie error... sorry... but it might be a good idea to 
mention this behaviour of mpf() a bit ealier than in the penultimate paragraph 
of the manual page...

Regards

Original comment by dr...@web.de on 26 Oct 2010 at 3:11

GoogleCodeExporter commented 9 years ago
No worry :)

I think the documentation is adequate. People will run into this "bug" no 
matter what the documentation says. It's just something you have to live with 
in the Python language (or C, or most other programming languages).

Original comment by fredrik....@gmail.com on 30 Oct 2010 at 9:52