Opened 9 years ago
Closed 9 years ago
#15299 closed defect (fixed)
Incorrect results for analytic Sha due to low precision
Reported by: | Jeroen Demeyer | Owned by: | |
---|---|---|---|
Priority: | major | Milestone: | sage-5.13 |
Component: | elliptic curves | Keywords: | |
Cc: | John Cremona, William Stein | 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: |
Description (last modified by )
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.
Attachments (2)
Change History (35)
comment:1 Changed 9 years ago by
Description: | modified (diff) |
---|
comment:2 follow-up: 3 Changed 9 years ago by
comment:3 follow-up: 5 Changed 9 years ago by
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 9 years ago by
I don't know about 11a1
, but this at least fixes the original problem.
comment:5 follow-up: 16 Changed 9 years ago by
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 9 years ago by
Status: | new → needs_review |
---|
comment:7 Changed 9 years ago by
Description: | modified (diff) |
---|
comment:8 Changed 9 years ago by
Description: | modified (diff) |
---|
comment:9 Changed 9 years ago by
...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.....
comment:10 follow-up: 11 Changed 9 years ago by
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 follow-up: 12 Changed 9 years ago by
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 Changed 9 years ago by
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 9 years ago by
If you care less about speed, I could also write a version using interval arithmetic, which will be simpler but slower.
comment:14 Changed 9 years ago by
Added tolerance to E.lseries().twist_values(1, -12, -4)
doctest to account for doctest failure on ia64.
comment:16 follow-up: 19 Changed 9 years ago by
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: 18 Changed 9 years ago by
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 Changed 9 years ago by
Replying to pbruin:
Jeroen, do you have a particular reason for writing
QQ()
andR()
instead ofQQ(0)
andR(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.
comment:19 follow-up: 20 Changed 9 years ago by
Replying to pbruin:
The
Lseries_ell
class has a methodL1_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 inderiv_at1()
, where we might use it to check the assumption L(E, 1) = 0. If that is too slow, acheck=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 follow-up: 22 Changed 9 years ago by
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 thatL(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: 24 Changed 9 years ago by
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 follow-up: 23 Changed 9 years ago by
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 thatL(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 Changed 9 years ago by
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 Changed 9 years ago by
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-precisionRealField
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: 26 Changed 9 years ago by
Dependencies: | → #15337 |
---|---|
Reviewers: | → Peter Bruin |
Status: | needs_review → 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 follow-up: 27 Changed 9 years ago by
Milestone: | sage-5.13 → 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 Changed 9 years ago by
Description: | modified (diff) |
---|---|
Report Upstream: | N/A → 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 9 years ago by
Dependencies: | #15337 → #15337, #15402 |
---|
Moved part of the patch to #15402.
Changed 9 years ago by
Attachment: | 15299_lseries_prec.patch added |
---|
comment:29 Changed 9 years ago by
Status: | positive_review → needs_review |
---|
Changed error bounds because of #15402, needs review.
comment:30 Changed 9 years ago by
Status: | needs_review → positive_review |
---|
Changed 9 years ago by
Attachment: | 15299_reviewer.patch added |
---|
replace page number by section number in Cohen reference
comment:31 Changed 9 years ago by
Description: | modified (diff) |
---|
comment:32 Changed 9 years ago by
Milestone: | sage-pending → sage-5.13 |
---|
comment:33 Changed 9 years ago by
Merged in: | → sage-5.13.beta4 |
---|---|
Resolution: | → fixed |
Status: | positive_review → closed |
It gets worse!
I.e., it is wrong for all curves of rank 0 too... (This isn't what we wrote the code for. Ugh.)