2009/6/23 Robert Bradshaw <rober...@math.washington.edu>: > > On Jun 21, 2009, at 12:54 PM, John Cremona wrote: > >> 2009/6/21 William Stein <wst...@gmail.com>: >>> >>> On Sun, Jun 21, 2009 at 8:38 PM, gsw<georgswe...@googlemail.com> >>> wrote: >>>> >>>> Hi John, >>>> >>>> On 21 Jun., 17:47, John Cremona <john.crem...@gmail.com> wrote: >>>>> This should be of interest to anyone who has ever had to manage >>>>> precision issues between Sage and pari real and complex numbers >>>>> (e.g. >>>>> Alex Ghitza). Others can move on. >>>>> >>>>> In the conversion of a pari complex number back to Sage (in >>>>> sage/libs/pari/gen_py.py in the function python(z)), the >>>>> precision of >>>>> the ComplexField in which the returned value should lie is taken >>>>> to be >>>>> the *minimum* of the precisions of the real and imaginary parts >>>>> of the >>>>> pari complex. Here I am assuming that both real and imaginary >>>>> parts >>>>> are inexact, otherwise there is no issue; and note that a pari gen >>>>> object (type t_COMPLEX) can have different types, and also >>>>> different >>>>> precisions, for the real and imaginary parts. And we *are* >>>>> converting >>>>> properly between pari's precision measured in words and Sage's >>>>> measured in bits. >>>>> >>>>> I propose to change this to the *maximum* of the two precisions. >>>>> >>>>> I have been experimenting with the wrapping of the pari functions >>>>> ellwp0(), which essentially takes ((w1,w2),z) where (w1,w2) is a >>>>> pair >>>>> of complex numbers defining a lattice L in C and z represents and >>>>> element of C/L, and returns the pair (P(z),P'(z)) where P() is the >>>>> Weierstrass P function. In some cases I find that the real and >>>>> imaginary parts returned by pari's function can have very different >>>>> precisions. >>>>> >>>>> Example (using code now developed, so this will not work in the >>>>> released version of Sage): >>>>> >>>>> sage: E = EllipticCurve([1,1,1,-8,6]) >>>>> sage: P = E(0,2) >>>>> sage: L = E.period_lattice() >>>>> sage: z = L.elliptic_logarithm(P, prec=129) >>>>> >>>>> Using the current min of the precisions (the exact value should >>>>> be (0,2)): >>>>> sage: L.elliptic_exponential(z) >>>>> (-3.8805107275644938151041666666176877354e-11 - >>>>> 2.7369110631344083416479093423631174439e-48*I, >>>>> 2.0000000000194025536378224690755208333 + >>>>> 4.1053665947016125124718640135446761659e-48*I) >>>>> >>>>> Using the proposed max: >>>>> sage: L.elliptic_exponential(z) >>>>> (8.8162076311671563097655240291668425836e-39 - >>>>> 2.7369110631344083416479093423631174439e-48*I, >>>>> 2.0000000000000000000000000000000000000 + >>>>> 4.1053665947016125124718640135446761659e-48*I) >>>>> >>>>> I swear that the only change between those two is changing line >>>>> 166 of >>>>> sage/libs/pari/gen_py.py from >>>>> prec = >>>>> min(gen.prec_words_to_bits(xprec),gen.prec_words_to_bits(yprec)) >>>>> to >>>>> prec = >>>>> max(gen.prec_words_to_bits(xprec),gen.prec_words_to_bits(yprec)) >>>>> >>>>> Note that in this example I started with 129 bits precision; at >>>>> 128 bits we get >>>>> >>>>> sage: z = L.elliptic_logarithm(P, prec=128) >>>>> sage: L.elliptic_exponential(z) >>>>> (8.6692708373143703712694319620140618739e-38, >>>>> 1.9999999999999999999999999999999999999) >>>>> (using min) >>>>> and the identical >>>>> sage: L.elliptic_exponential(z) >>>>> (8.6692708373143703712694319620140618739e-38, >>>>> 1.9999999999999999999999999999999999999) >>>>> (using max) >>>>> which are the same. Here 128 bits is an exact number of words, >>>>> while >>>>> 129 forces pari to use an extra word most of which is garbage on >>>>> input, and somehow the effect is to catastrophically reduce the >>>>> precision of the output! >>>>> >>>>> I just checked (using testall) that this change has no effect at >>>>> all >>>>> on any doctests, so I propose to include it in my upcoming patch >>>>> (which implements the complex elliptic exponential for elliptic >>>>> curves, by a simple call to the already wrapped pari.ellwp0(), but >>>>> which has taken ages to sort out the precision properly. >>>>> >>>>> John >>>> >>>> let me try to rephrase this, I'm not sure whether I did >>>> understand you >>>> correctly: >>>> >>>> a) On certain inputs a certain pari function returns complex values, >>>> whose real parts and imaginary parts it reports as being quite >>>> different (in terms of precision). >> >> True. >> >>>> >>>> b) Closer inspection yields that the part with the lower precision >>>> reported actually holds more information than the reported precision >>>> would allow for, especially it is OK to regard this part to be of >>>> the >>>> (greater) precision of the other part. >>>> >> >> Hmm. In the case in question the imaginary part should have been 0 >> exactly if the inputs were known to infinite precision. But the value >> returned was approximately 0 with a precision of only 3 words while >> the real part had precision 7 words. >> >>>> c) If so, I would consider this behaviour a bug in pari, to be fixed >>>> in pari. >> >> I was planning to report this to pari, but it is not easy to give a >> reproducible example --precision is handled differently in pari's >> library mode than in gp, so it would mean writing a C program. I'm >> not sure how the pari people would like getting a bug report which >> relied on them installing Sage to verify! It does not help that the >> pari documentation is less than clear about what precision should be >> expected. (There is a prec parameter to ellwp0 but that parameter is >> frequently ignored in cases where the precision of the other >> parameters can be used to determine the best output precision). >> >>>> >>>> d) In the meantime, you propose as a workaround to alter the >>>> conversion behaviour from pari to Sage for complex numbers in >>>> general. >>>> >>>> If I got right --- I vote +1 for your proposal. >>>> >> >> Thanks, but I'll take William's suggestion on board too. >> >>> >>> If he got it right above, then I would suggest instead making a >>> function in gen.pyx that does what you want and using it directly. >>> You could also add a non-default flag to the conversion function so >>> that when that option is set the conversion function uses the max >>> instead of the min. >> > > > I'm wary of automatically promoting low precision numbers to high > ones, so I think the flag sounds like a good option too. This also > seems strange: > > sage: pari(ComplexField(500).gen()/3) > 0.E-154 + 0.333333333333333*I
Here's what I ended up doing in the patch (#6386, needs review ;)): I left alone the conversion of pari t_COMPLEX to Sage. Instead I use this line C(t.real().python(),t.imag().python()) instead of C(t.python()) where C is a ComplexField of the precision I started with. For best results of course one should use a bit precision in Sage which corresponds to an exact number of words, but I do not want to make that a rule! John > > - Robert > > > > > --~--~---------~--~----~------------~-------~--~----~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send 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 -~----------~----~----~----~------~----~------~--~---
