problem with ellipk / hyp2f1
GoogleCodeExporter opened this issue · 3 comments
GoogleCodeExporter commented
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
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
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
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
- Changed state: WontFix