Opened 10 years ago
Closed 9 years ago
#12509 closed defect (fixed)
computation of height of point on elliptic curve over Q(sqrt(5)) is WRONG
Reported by: | was | Owned by: | was |
---|---|---|---|
Priority: | critical | Milestone: | sage-5.10 |
Component: | elliptic curves | Keywords: | heights |
Cc: | Merged in: | sage-5.10.beta3 | |
Authors: | John Cremona | Reviewers: | Peter Bruin, Chris Wuthrich |
Report Upstream: | N/A | Work issues: | |
Branch: | Commit: | ||
Dependencies: | Stopgaps: | #12692 |
Description (last modified by )
There are evidently many examples in which computing P.height()
, for P
a point on an elliptic curve over Q(sqrt(5)) yields a completely wrong answer. This is very serious, since it is a blatantly wrong mathematical answer -- raising NotImplementedError? would be much better! Here's an example that Ashwath Rabindranath (Princeton) found, where Sage and Magma do not agree. According to BSD, Sha has order 1 using the Magma answer, and a crazy order with the Sage answer.
sage: K.<a> = NumberField(x^2-x-1) sage: v = [0, a + 1, 1, 28665*a - 46382, 2797026*a - 4525688] sage: E = EllipticCurve(v) sage: E == E.global_minimal_model() True sage: F.a_invariants() (0, a + 1, 1, 28665*a - 46382, 2797026*a - 4525688) sage: P = E([72*a - 509/5, -682/25*a - 434/25]) sage: P.height() 1.35648516097058 sage: Q = magma(E)(magma([P[0], P[1]])) sage: Q (1/5*(360*a - 509) : 1/25*(-682*a - 434) : 1) sage: Q.Height() 1.38877711688727252538242306
Apply: trac12509-heights.patch, trac12509-heights2.patch, trac12509-heights3.patch
Attachments (3)
Change History (30)
comment:1 Changed 10 years ago by
comment:2 Changed 10 years ago by
A little hg blame sage/schemes/elliptic_curves/ell_point.py
and hg history sage/schemes/elliptic_curves/ell_point.py
shows that this code comes from #8496 , merged in sage-4.4.alpha1. The next change to this file happened with #8827 in a later release. I have verified that the example in the ticket gives the same result on sage-4.4, so the error has been there all along.
Note that:
sage: [P.height(precision=p) for p in [53..70]] [1.35648516097058, 1.37988660372465, 1.392904136135677, 1.388358180620918, 1.390651822960756, 1.3901262440498197, 1.3886249896939375, 1.3887137996432686, 1.38889123047486600, 1.38868605331097548, 1.38879145662338488, 1.38876788728230994, 1.388772046900726933, 1.388778979290526601, 1.388773433409450769, 1.3887772462291215145, 1.3887758164353761832, 1.3887774411997285045]
so the computation is probably not really wrong. It's probably just a matter of precision estimates being off.
comment:3 Changed 10 years ago by
This is also pretty damning:
sage: (3*P).height()/P.height() 9.19316602595430 sage: (2*P).height()/P.height() 4.09522243754782
My guess is that the local contributions from the finite primes are fine (each is rational times log(p)), it's just the (two) real places. One could compare with the real local heights as computed by Magma to verify this.
comment:4 Changed 10 years ago by
- Stopgaps set to #12692
comment:5 Changed 9 years ago by
This is a completely TRIVIAL thing, and I cannot understand why I did not spot this earlier. In the example, the height is computed as the sum of archimedean and non-arch components. The latter is log(25)/2 and is returned as such *using the python builtin log function*! This is where the precision is lost, Compare
sage: P.height() 1.35648516097058 sage: P.archimedian_local_height() + P.nonarchimedian_local_height() 1/2*log(25) - 0.2206607955468278522013468194382 sage: RR(P.archimedian_local_height() + P.nonarchimedian_local_height()) 1.38877711688727
and you wonder how it is possible to get this wrong:
sage: P.archimedian_local_height(prec=53) + P.nonarchimedian_local_height(prec=53) 1.35648516097
I will make a patch now.
comment:6 Changed 9 years ago by
I was a little hasty. There's no problem in the non-arch local height being an exact log:
sage: P.nonarchimedian_local_height() 1/2*log(25) sage: P.nonarchimedian_local_height(prec=53) 1.60943791243410 sage: RR(P.nonarchimedian_local_height()) 1.60943791243410
is all correct. The problem is with this:
sage: P.archimedian_local_height() -0.2206607955468278522013468194382 sage: P.archimedian_local_height(prec=53) -0.252952751464
where the first is correct as this shows:
sage: P.archimedian_local_height()+ P.nonarchimedian_local_height(prec=53) 1.38877711688727
but the second is not.
comment:7 Changed 9 years ago by
It seems to be a numerical stability problem. With the same example, and using the first embedding of K, we have, first to default 53-bit precision:
sage: b2,b4,b6,b8 (1.527864045, -128195.888575, -25017379.5417, -4118102250.56) sage: x -146.29844719 sage: t -0.00683534254264 sage: 1 - (b4 * t**2) - (2*b6 * t**3) - (b8 * t**4) 3.5527136788e-15
while at 100-bit precision:
sage: bb2,bb4,bb6,bb8 (1.5278640450004206071816526625, -128195.88857503147164756896321, -2.5017379541668653550777130780e7, -4.1181022505609323099570853156e9) sage: xx -146.29844718999242907073025207 sage: tt -0.0068353425426405016233236537484 sage: 1 - (bb4 * tt**2) - (2*bb6 * tt**3) - (bb8 * tt**4) 3.5142315895588997690642550530e-15
This means that whenever we evaluate one of the four polynomials (which coefficients are essentially the b-invariants) at some point we may lose precision badly.
One way out would be to use a "working precision" higher than the requested precision. Another would be to use algebraic values instead of approximate ones for t,z,w in the algorithm, which is likely to be slow, but I will try it.
A better solution for the longer term would be to implement Mestre's AGM algorithm instead BUT this only exists at present for real embeddings (I have worked on extending to the complex case, but that is still in progress).
comment:8 Changed 9 years ago by
Patch ready for review. Basically we double the working precision while computing each archidean component of the height, and this suffices. I had to add prec as an alias for precision for CDF (as it already was for RDF), and a couple of doctest outputs had the last digit changed by 1. There are doctests to show that the problem is fixed.
comment:9 Changed 9 years ago by
- Component changed from number theory to elliptic curves
- Description modified (diff)
- Keywords heights added
- Status changed from new to needs_review
comment:10 Changed 9 years ago by
- Reviewers set to Peter Bruin
comment:11 Changed 9 years ago by
The original precision problem seems to be solved by the patch, but there is a remaining problem if the given place v
has lower precision than prec
(the answer ends with a string of zeroes):
sage: K.<a> = NumberField(x^2-x-1) sage: v = [0, a + 1, 1, 28665*a - 46382, 2797026*a - 4525688] sage: E = EllipticCurve(v) sage: P = E([72*a - 509/5, -682/25*a - 434/25]) sage: P.archimedian_local_height(v=K.places()[1],prec=1000) 2.81105581458853745263599796447519791747641640992814370122945754557165333149900887291328934930695388487942633879651281942220091469088330120430293759305481969374444866348067978423230783846520353108644485473632812500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
The precision of v
shouldn't influence the result, of course.
Also, the default value of prec
used when applying Silverman's Theorem 4.2 is currently the precision of the refined place vv
, but it suffices to take the same default precision to which the result is returned, namely the precision of v
.
I wrote a patch that fixes both issues; I am now going to test it.
comment:12 Changed 9 years ago by
Thanks: your second point is that while we still work in double the precision, the number of terms needed can be computed from the original precision? That makes sense.
On the second point, it is certainly possible that the user asks for two different things by using a low-precision embedding but then asking for a higher precision without refining the embedding. I would fix this by first replacing v (if necessary) by a refinement to precision "prec". Which may be what you patch does, so I should look at it.... and indeed that is what you do.
The second patch looks good to me, hope the tests pass!
comment:13 follow-up: ↓ 14 Changed 9 years ago by
The tests in sage/schemes/elliptic_curves/ all pass for me, now testing the whole Sage library.
If this also succeeds, do my additional changes still allow me to give a positive review, or do we need a new reviewer?
comment:14 in reply to: ↑ 13 Changed 9 years ago by
Replying to pbruin:
The tests in sage/schemes/elliptic_curves/ all pass for me, now testing the whole Sage library.
If this also succeeds, do my additional changes still allow me to give a positive review, or do we need a new reviewer?
Maybe, maybe not. William?
comment:15 Changed 9 years ago by
OK, all tests passed.
Changed 9 years ago by
comment:16 Changed 9 years ago by
trac12509-heights2.patch now also makes the documentation slightly more accurate
comment:17 Changed 9 years ago by
- Description modified (diff)
comment:18 Changed 9 years ago by
I did a small rebase of the first patch for 5.9.beta5; both patches still need to be applied. We need a reviewer.
comment:19 follow-up: ↓ 20 Changed 9 years ago by
I'm keeping myself on the reviewers list for the moment, but only for John's patch. For the additional changes I made, someone else should probably review the whole ticket.
Should I rebase my patch too?
comment:20 in reply to: ↑ 19 Changed 9 years ago by
Replying to pbruin:
I'm keeping myself on the reviewers list for the moment, but only for John's patch. For the additional changes I made, someone else should probably review the whole ticket.
Should I rebase my patch too?
It's not necessary, your patch still applies OK. (The only bit of mine which did not apply was the alias of prec=precision in sage.rings.complex_double.pyx)
comment:21 Changed 9 years ago by
- Status changed from needs_review to positive_review
I checked the two patches. They apply cleanly to 5.9 and the tests pass. As far as I understand they do the right thing.
comment:22 Changed 9 years ago by
- Reviewers changed from Peter Bruin to Peter Bruin, Chris Wuthrich
comment:23 Changed 9 years ago by
- Description modified (diff)
comment:24 Changed 9 years ago by
- Status changed from positive_review to needs_work
There is a documentation problem:
dochtml.log:[plane_cur] /mazur/release/merger/sage-5.10.beta2/local/lib/python2.7/site-packages/sage/schemes/elliptic_curves/ell_point.py:docstring of sage.schemes.elliptic_curves.ell_point.EllipticCurvePoint_number_field.archimedian_local_height:49: WARNING: Literal block expected; none found.
comment:25 Changed 9 years ago by
- Description modified (diff)
- Status changed from needs_work to needs_review
comment:26 Changed 9 years ago by
- Status changed from needs_review to positive_review
comment:27 Changed 9 years ago by
- Merged in set to sage-5.10.beta3
- Resolution set to fixed
- Status changed from positive_review to closed
I cannot remember who implemented this, but have looked at the code recently and it is using the standard Silverman algorithm.
I think it is likely that there is already a ticket for this.