Opened 8 years ago

Closed 8 years ago

#15299 closed defect (fixed)

Incorrect results for analytic Sha due to low precision

Reported by: jdemeyer Owned by:
Priority: major Milestone: sage-5.13
Component: elliptic curves Keywords:
Cc: cremona, was Merged in: sage-5.13.beta4
Authors: Jeroen Demeyer Reviewers: Peter Bruin
Report Upstream: Fixed upstream, but not in a stable release. Work issues:
Branch: Commit:
Dependencies: #15337, #15402 Stopgaps:

Status badges

Description (last modified by pbruin)

See https://groups.google.com/forum/#!topic/sage-support/rYQ4rWyncG4

sage: E = EllipticCurve(QQ,[0, 0, 1, -79, 342]); E.sha().an(),E.sha().an_numerical()
(0, 1.00000000000000)
sage: E.lseries().deriv_at1(100*sqrt(E.conductor()) + 10)  # L1, error_bound
(1.94655218772191e-15, 1.82478252135394e-270)

This seems to be due to inappropriate use of Python floats instead of more precise real numbers. After the patch:

sage: E = EllipticCurve(QQ,[0, 0, 1, -79, 342]); E.sha().an(),E.sha().an_numerical()
(1.00000000000000, 1.00000000000000)
sage: E.lseries().deriv_at1(100*sqrt(E.conductor()) + 10)  # L1, error_bound
(-4.32787398660869448751904675450772492666840247314688171540527473331725818170217268435223462033366791557160872926179439894639315476270837428785657638252738603056742447337636326343956370276624493916496382120766160023620812331280787034239648552009947468067829864968026720015203778821069593806584e-277,
 1.82478252137476307223140369768561190028055347258560054363485475966241792307587640145132294203994875344783110100551912347495775160520204557245032474939095251969168953786545612090565728262067746413119194690260652692254781091147749697957445424152473292233020112755190503925812425294821095313979e-270)

While working on this, we found an upstream PARI bug: the precision for exponential_integal_1() was not as good as it could be.

Apply: 15299_lseries_prec.patch, 15299_reviewer.patch

Attachments (2)

15299_lseries_prec.patch (38.5 KB) - added by jdemeyer 8 years ago.
15299_reviewer.patch (1.5 KB) - added by pbruin 8 years ago.
replace page number by section number in Cohen reference

Download all attachments as: .zip

Change History (35)

comment:1 Changed 8 years ago by jdemeyer

  • Description modified (diff)

comment:2 follow-up: Changed 8 years ago by was

It gets worse!

sage: E = EllipticCurve('11a')
sage: E.lseries().deriv_at1()
0
sage: E.lseries().dokchitser().derivative(1)
0.308708533963172

I.e., it is wrong for all curves of rank 0 too... (This isn't what we wrote the code for. Ugh.)

comment:3 in reply to: ↑ 2 ; follow-up: Changed 8 years ago by jdemeyer

So this is plain wrong then (I don't know enough mathematics to judge this):

        if self.__E.root_number() == 1:
            return 0

comment:4 Changed 8 years ago by jdemeyer

I don't know about 11a1, but this at least fixes the original problem.

comment:5 in reply to: ↑ 3 ; follow-up: Changed 8 years ago by cremona

Replying to jdemeyer:

So this is plain wrong then (I don't know enough mathematics to judge this):

        if self.__E.root_number() == 1:
            return 0

The root number is the sign of the functional equation so is +1 for even analytic rank and -1 for odd. This function computes the first derivative. *In practice* this is something one would only want to do if the 0'th derivative was already known to be 0, in which case the code you quote would be OK since if the value is 0 and the order is even then the order is at least 2 so the first derivative is exactly 0. But of course this function then lies in wait for the user who decides they want the first derivative's value even when the value is nonzero (as for 11a1). The trouble is that (1) Formulas for the r'th derivative which are implemented are *only* valid under the assumption that all previous derivatives are 0; and of course (2) proving the earlier derivatives are exactly 0 is usually impossible with current theory.

Where does that leave this deriv_at1 function? At the very least it should come with a huge warning about all this. And it really should return 0 when the root number is +1 unless the user has made an explicit assumption (assume_order_of_vanishing_is_positive=True, say) and otherwise raise a NotImplemented? error (or attempt to prove that L(1)=0).

comment:6 Changed 8 years ago by jdemeyer

  • Status changed from new to needs_review

comment:7 Changed 8 years ago by jdemeyer

  • Description modified (diff)

comment:8 Changed 8 years ago by jdemeyer

  • Description modified (diff)

comment:9 Changed 8 years ago by cremona

...review in progress...and all works fine, including docbuilding. the changes look good to the human eye (this one at least). Oops, forgot the --long when testing..... and it's still all good. Thanks, Jeroen!

Last edited 8 years ago by cremona (previous) (diff)

comment:10 follow-up: Changed 8 years ago by jdemeyer

Corrected the error computation for at1(). I believe this is rigorous now:

        for n in xrange(1,k+1):
            term = (zpow * an[n])/n
            zpow *= z
            L += term
            # 8n+1 is the relative error in half-ulps to compute term.
            # For addition, multiplication, division, sqrt, this is
            # bounded by the number of operations. exp(x) multiplies the
            # relative error by abs(x) and adds 1 half-ulp. The relative
            # error for -2*pi/sqrtN is 3 half-ulps. Assuming that
            # 2*pi/sqrtN <= 2, the relative error for z is 7 half-ulps.
            # This implies a relative error of 8n-1 half-ulps for zpow.
            # Adding 2 for the computation of term gives:
            error += term.ulp()*(8*n+1) + L.ulp()

comment:11 in reply to: ↑ 10 ; follow-up: Changed 8 years ago by cremona

Replying to jdemeyer:

Corrected the error computation for at1(). I believe this is rigorous now:

        for n in xrange(1,k+1):
            term = (zpow * an[n])/n
            zpow *= z
            L += term
            # 8n+1 is the relative error in half-ulps to compute term.
            # For addition, multiplication, division, sqrt, this is
            # bounded by the number of operations. exp(x) multiplies the
            # relative error by abs(x) and adds 1 half-ulp. The relative
            # error for -2*pi/sqrtN is 3 half-ulps. Assuming that
            # 2*pi/sqrtN <= 2, the relative error for z is 7 half-ulps.
            # This implies a relative error of 8n-1 half-ulps for zpow.
            # Adding 2 for the computation of term gives:
            error += term.ulp()*(8*n+1) + L.ulp()

I can see where this is in the code -- can you say how it affects any doctest outputs? I am not a numerical analyst...

comment:12 in reply to: ↑ 11 Changed 8 years ago by jdemeyer

Replying to cremona:

can you say how it affects any doctest outputs?

I would only make the error results slightly larger than before. Perhaps more importantly, I believe the error computation is now rigorous, in the sense that one could prove that the actual error is bound by the error result.

comment:13 Changed 8 years ago by jdemeyer

If you care less about speed, I could also write a version using interval arithmetic, which will be simpler but slower.

comment:14 Changed 8 years ago by jdemeyer

Added tolerance to E.lseries().twist_values(1, -12, -4) doctest to account for doctest failure on ia64.

comment:15 Changed 8 years ago by jdemeyer

Passes tests on the buildbots now.

comment:16 in reply to: ↑ 5 ; follow-up: Changed 8 years ago by pbruin

Replying to cremona:

Where does that leave this deriv_at1 function? At the very least it should come with a huge warning about all this. And it really should return 0 when the root number is +1 unless the user has made an explicit assumption (assume_order_of_vanishing_is_positive=True, say) and otherwise raise a NotImplemented? error (or attempt to prove that L(1)=0).

The Lseries_ell class has a method L1_vanishes(), written by William Stein. According to the documentation, it is provably correct if the Manin constant is <= 2 for the optimal quotient in the isogeny class. This method is used in some places, but not currently in deriv_at1(), where we might use it to check the assumption L(E, 1) = 0. If that is too slow, a check=False flag could be added.

comment:17 follow-up: Changed 8 years ago by pbruin

Jeroen, do you have a particular reason for writing QQ() and R() instead of QQ(0) and R(0)? It saves one character, but is less readable.

comment:18 in reply to: ↑ 17 Changed 8 years ago by jdemeyer

Replying to pbruin:

Jeroen, do you have a particular reason for writing QQ() and R() instead of QQ(0) and R(0)? It saves one character, but is less readable.

No reason, it's just a habit to think in term of default constructors (I guess this is my C++ background). I just benchmarked QQ(), QQ(0) and QQ.zero() and the latter is the fastest, so perhaps we should use that.

Last edited 8 years ago by jdemeyer (previous) (diff)

comment:19 in reply to: ↑ 16 ; follow-up: Changed 8 years ago by jdemeyer

Replying to pbruin:

The Lseries_ell class has a method L1_vanishes(), written by William Stein. According to the documentation, it is provably correct if the Manin constant is <= 2 for the optimal quotient in the isogeny class. This method is used in some places, but not currently in deriv_at1(), where we might use it to check the assumption L(E, 1) = 0. If that is too slow, a check=False flag could be added.

Do you propose that this change should be made, or is it just an observation? Given that the function deriv_at1() is in practice only called when we know that L(E,1) = 0, I personally think the warning suffices.

comment:20 in reply to: ↑ 19 ; follow-up: Changed 8 years ago by pbruin

Replying to jdemeyer:

Do you propose that this change should be made, or is it just an observation? Given that the function deriv_at1() is in practice only called when we know that L(E,1) = 0, I personally think the warning suffices.

I agree, it was more an observation that we could in principle use L1_vanishes() here than a proposal to actually do it.

There is a formula for the r-th derivative which is valid when all lower derivatives vanish. As far as I know, only for the 0-th derivative is there a known way to prove that it vanishes by a numerical computation. For parity reasons (the root number), this means that if the order of vanishing is 0, 1 or 2, then we can prove this. If the order of vanishing is 3, then in general we don't know how to prove that it is not 1.

This means that if and when the formula mentioned above is implemented, we won't be able to verify the condition "all lower derivatives are 0" when r is at least 3. Hence we should probably not insist on verifying it when r = 1.

comment:21 follow-up: Changed 8 years ago by pbruin

Another question: is it necessary to compute the error bound to the same precision as the result, i.e. in RealField(prec)? It seems sufficient, and more efficient, to compute it in a lower-precision RealField or even just using Python floats.

comment:22 in reply to: ↑ 20 ; follow-up: Changed 8 years ago by cremona

Replying to pbruin:

Replying to jdemeyer:

Do you propose that this change should be made, or is it just an observation? Given that the function deriv_at1() is in practice only called when we know that L(E,1) = 0, I personally think the warning suffices.

I agree, it was more an observation that we could in principle use L1_vanishes() here than a proposal to actually do it.

There is a formula for the r-th derivative which is valid when all lower derivatives vanish. As far as I know, only for the 0-th derivative is there a known way to prove that it vanishes by a numerical computation. For parity reasons (the root number), this means that if the order of vanishing is 0, 1 or 2, then we can prove this. If the order of vanishing is 3, then in general we don't know how to prove that it is not 1.

You can go one step further thanks to Gross-Zagier: if the parity is odd and L'(1) looks zero then you can prove it, since if in fact L'(1)!=0 then the curve would have rank 1, but you can disprove that by finding three (or oeven only 2) independent points. See by talk http://homepages.warwick.ac.uk/staff/J.E.Cremona/papers/bsd50.pdf if you want to read more!

This means that if and when the formula mentioned above is implemented, we won't be able to verify the condition "all lower derivatives are 0" when r is at least 3. Hence we should probably not insist on verifying it when r = 1.

comment:23 in reply to: ↑ 22 Changed 8 years ago by pbruin

Replying to cremona:

You can go one step further thanks to Gross-Zagier: if the parity is odd and L'(1) looks zero then you can prove it, since if in fact L'(1)!=0 then the curve would have rank 1, but you can disprove that by finding three (or oeven only 2) independent points.

That is true (in fact I seem to remember learning this from the talk you linked to). However, it does require you to search for points; there seems to be no "analytic" way of proving that L'(1) = 0 by computing it to finite precision, like the L1_vanishes() function.

comment:24 in reply to: ↑ 21 Changed 8 years ago by jdemeyer

Replying to pbruin:

Another question: is it necessary to compute the error bound to the same precision as the result, i.e. in RealField(prec)? It seems sufficient, and more efficient, to compute it in a lower-precision RealField

Done. Needs #15337.

or even just using Python floats.

Not a good idea, as these have limited range and it's a lot harder (maybe even impossible) to control the rounding. One doctest has an error of 2.74997188336901e-449, which would be rounded to 0.0 as Python float.

comment:25 follow-up: Changed 8 years ago by pbruin

  • Dependencies set to #15337
  • Reviewers set to Peter Bruin
  • Status changed from needs_review to positive_review

This looks very good now. The error analysis appears to be completely rigorous for at1() and almost completely rigorous for deriv_at1(), the only source of non-rigorousness being due to the unknown error in the exponential integral function eint1() from PARI. This ticket is not the place to try to fix this, though.

Are the PARI developers aware of this precision issue? Should it be regarded it as a bug, or does PARI not strive for proven error bounds for functions such as eint1()?

comment:26 in reply to: ↑ 25 ; follow-up: Changed 8 years ago by jdemeyer

  • Milestone changed from sage-5.13 to sage-pending

Replying to pbruin:

Are the PARI developers aware of this precision issue?

No idea. I might report it.

does PARI not strive for proven error bounds for functions such as eint1()?

I don't think PARI does. However, the errors are quite large (over 30 bits can be wrong), so perhaps that's a bug indeed.

comment:27 in reply to: ↑ 26 Changed 8 years ago by jdemeyer

  • Description modified (diff)
  • Report Upstream changed from N/A to Fixed upstream, but not in a stable release.

Reported the precision issue, they fixed it. There is supposed to be an absolute error bound (not relative), but I don't know what the bound is...

comment:28 Changed 8 years ago by jdemeyer

  • Dependencies changed from #15337 to #15337, #15402

Moved part of the patch to #15402.

Changed 8 years ago by jdemeyer

comment:29 Changed 8 years ago by jdemeyer

  • Status changed from positive_review to needs_review

Changed error bounds because of #15402, needs review.

comment:30 Changed 8 years ago by pbruin

  • Status changed from needs_review to positive_review

Looks even better than before, the precision is now much better under control thanks to #15402, and the remaining "arbitrary" precision increase is clearly motivated.

As in #15402, just a trivial review patch to refer to a section instead of a page number in Cohen's book.

Changed 8 years ago by pbruin

replace page number by section number in Cohen reference

comment:31 Changed 8 years ago by pbruin

  • Description modified (diff)

comment:32 Changed 8 years ago by jdemeyer

  • Milestone changed from sage-pending to sage-5.13

comment:33 Changed 8 years ago by jdemeyer

  • Merged in set to sage-5.13.beta4
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.