On Sat, Aug 15, 2009 at 4:40 PM, Minh Nguyen<nguyenmi...@gmail.com> wrote:
>
> Hi folks,
>
> I noticed the following thread from the Maxima mailing list.
>
> --
> Regards
> Minh Van Nguyen
>
>
> ---------- Forwarded message ----------
> From: Richard Fateman <fate...@cs.berkeley.edu>
> Date: Sun, Aug 16, 2009 at 12:35 AM
> Subject: Re: [Maxima] mpmath + sage + hypergeometric numerics
> To: Barton Willis <will...@unk.edu>
> Cc: Maxima List <max...@math.utexas.edu>
>
>
> I tried Mathematica 6.0 on
> HypergeometricPFQ[{112, 51/10}, {-9/10}, -99999/100000]
>
> and got the (presumably exact ??) result
>
> -183933169187931829794351104521526197029774379000000000000000000000000\
> 0000000000000000000000000000000000000000000000000000000000000000000000\
> 0000000000000000000000000000000000000000000000000000000000000000000000\
> 0000000000000000000000000000000000000000000000000000000000000000000000\
> 0000000000000000000000000000000000000000000000000000000000000000000000\
> 0000000000000000000000000000000000000000000000000000000000000000000000\
> 0000000000000000000000000000000000000000000000000000000000000000000000\
> 0000000000000000000000000000000000000000000000000000000000000000000000\
> 0000000000000000000000000000000000000000000000/
> 11324984934739307220163426028747102350310899924999124109622447068264\
> 2089509049347831693103349943262869210969949172830357221037131460936617\
> 2457317140303388728526568056839394088427274999192970896030049564669168\
> 5199955047630663915447045353398112210710118309005643847603222953120252\
> 7472676579674432971253271831088126529101868279771770153695910346445313\
> 9991441759761159325993368148026744708619920408656674968525340473179327\
> 4543394221143555750764130509074977665133986904900063950461891643545329\
> 6240971237497707332066495140483266031003755059081400674335749107079469\
> 7608945169169804679596517702998625671016983228364785324615691195240034\
> 1
> which evaluates to the Mathematica result via N[%,100].
>
> In the documentation for HypergeometricPFQ there are illustrations where the
> internal $MaxExtraPrecision is insufficient to obtain the requested
> precision for this function,
> so apparently Mathematica fails to do adequate numerical evaluation
> using N[] without resetting some global system parameter.
>
> I don't know how to run your code or mpmath on exact inputs.
>
> I have, in the past, used python-wrapped code from Allegro CL with
> almost no effort.
> This was for GMP.
> I found, however, that the python wrappers were too unsophisticated for me,
> and
> implemented only the "Functional" style in GMP. That is, A:=A*B+C could not
> access the combined multiply-add-replace capability. I don't know if
> that is still the case.
>
>
>
> Barton Willis wrote:
>
> For Maxima 5.19, I contributed code for numerical and symbolic
> evaluation of the hypergeometric functions. The other day, I started
> looking into how Sage evaluates these functions (but I haven't tried to
> install Sage). One way seems to be mpmath
> (http://code.google.com/p/mpmath/).
> The NSF supported the development of mpmath:
>
> http://groups.google.com/group/mpmath/browse_thread/thread/960b02f0a2ed5dfa
>
>
> Let's try mpmath on a hypergeometric function; first let's try 50 digits
> and
> second 100 digits. The signs of the two values aren't the same:
>
>
>
> mp.dps = 50
> hyper(['112.0', '5.1'],['-0.9'],'-0.99999')
>
>
> mpf('0.00000000000073150973254264463451967741826767251189694769140430165')
>
>
>
> mp.dps = 100
> hyper(['112.0', '5.1'],['-0.9'],'-0.99999')
>
>
> mpf('-0.00000000000000000000000162413610479<junk>')
>
> Wolfram Alpha (input: hypergeometricPQ[{112,5.1}, {-0.9}, -0.99999]) gives
>
> -1.62413610479708629611767409391575176600337612526370... x 10^-24
>
> And my (non-NSF) code gives
>
> (%i4) hypergeometric([112.0b0, 5.1b0],[-0.9b0],-0.99999b0);
> (%o4) - 1.6331655633849675761598634419620944688583571143885b-24
>
> I'm not sure why my value differs from mma by about 9.0e-27, but
> likely it is due to differences in the precisions of the inputs. At
> least my code gets the sign correct. Grump^infinity.
>
> Barton
This is indeed a bug in mpmath, and I'm grateful that it was brought
to my attention. After fixing it, I get:
>>> from mpmath import *
>>> mp.dps = 200
>>> a,b,c,z = map(mpmathify, ['112', '5.1', '-0.9', '-0.99999'])
>>> for dps in [5, 15, 30, 50, 100]:
... mp.dps = dps
... print hyp2f1(a,b,c,z)
...
-1.6241e-24
-1.62413610479709e-24
-1.62413610479708629611767409392e-24
-1.6241361047970862961176740939157517660033761252637e-24
-0.000000000000000000000001624136104797086296117674093915751766003376125263700886074038894673955627186461247032427607279195267
This seems to agree to the full relative precision with Mathematica
and Wolfram Alpha.
I have not committed the bugfix yet, as I want to test it some more first.
Regarding exact/approximate parameters: Note that I precomputed each
argument above to a much higher precision than the target precision.
It will also work if you pass numbers with an exact binary
representation. Generally, mpmath assumes that inputs are exact binary
floating-point numbers, and the user has to ensure that they are
accurate enough if the function is very sensitive to perturbations (as
is the case here).
If you pass decimal strings as parameters, they will be converted to
binary using the calling precision, which is a convenient feature but
only gives fully accurate answers if not much extra working precision
is required. So for example (with the bugfix applied):
>>> for dps in [5, 15, 30, 50, 80]:
... mp.dps = dps
... print hyp2f1(112,'5.1','-0.9','-0.99999')
...
-8.4964e-13
-1.49991019565633e-22
-1.62413610479695451985358601598e-24
-1.6241361047970862961176740939157535519047737042344e-24
-0.0000000000000000000000016241361047970862961176740939157517660033761252637008860740388960827834203160782
Here the signs are correct, and presumably each output is correct for
the binary input parameters that result from perturbing the inputs
(but not correct for the given decimal numbers). The reason why Maxima
gives a nearly, but not quite, correct result could be the same one.
This behavior is rather by design than a bug, although it can be
surprising and inconvenient. Possible solutions are:
1. Support exact rational parameters. Mpmath generally supports these
for the parameters of hypergeometric functions (they can be specified
as integer pairs). The case of hyp2f1 near the unit circle,
unfortunately, is one of very few algorithms that don't support exact
rational parameters yet. I should be able to fix this relatively soon.
2. Interval arithmetic (to handle rational as well as nonrational
parameters). Much of the hypergeometric code in mpmath is generic
enough that it should work with interval arithmetic with some tweaks.
The algorithm used above worked with intervals by changing 4 lines,
and then it returned something like [-234982034.239784,
+234679263.394857] without the bugfix, correctly indicating that
catastrophic cancellation had occurred. It then gave a narrow,
presumably correct, interval with the precision increased. It might be
a while before I enable support for intervals, however.
Fredrik
--~--~---------~--~----~------------~-------~--~----~
To post to this group, send an email to sage-devel@googlegroups.com
To unsubscribe from this group, send an email to
sage-devel-unsubscr...@googlegroups.com
For more options, visit this group at http://groups.google.com/group/sage-devel
URLs: http://www.sagemath.org
-~----------~----~----~----~------~----~------~--~---
