Opened 8 years ago

Closed 7 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 pbruin)

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)

trac12509-heights2.patch (2.9 KB) - added by pbruin 7 years ago.
trac12509-heights.patch (9.0 KB) - added by cremona 7 years ago.
rebased to 5.9beta5
trac12509-heights3.patch (607 bytes) - added by pbruin 7 years ago.
fix documentation syntax

Download all attachments as: .zip

Change History (30)

comment:1 Changed 8 years ago by cremona

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.

comment:2 Changed 8 years ago by nbruin

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 8 years ago by cremona

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 8 years ago by jen

  • Stopgaps set to #12692

comment:5 Changed 8 years ago by cremona

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 8 years ago by cremona

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 8 years ago by cremona

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 7 years ago by cremona

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 7 years ago by cremona

  • Authors set to John Cremona
  • Component changed from number theory to elliptic curves
  • Description modified (diff)
  • Keywords heights added
  • Status changed from new to needs_review

comment:10 Changed 7 years ago by pbruin

  • Reviewers set to Peter Bruin

comment:11 Changed 7 years ago by pbruin

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 7 years ago by cremona

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: Changed 7 years ago by 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?

comment:14 in reply to: ↑ 13 Changed 7 years ago by cremona

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 7 years ago by pbruin

OK, all tests passed.

Changed 7 years ago by pbruin

comment:16 Changed 7 years ago by pbruin

trac12509-heights2.patch now also makes the documentation slightly more accurate

Changed 7 years ago by cremona

rebased to 5.9beta5

comment:17 Changed 7 years ago by cremona

  • Description modified (diff)

comment:18 Changed 7 years ago by cremona

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: Changed 7 years ago by 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?

comment:20 in reply to: ↑ 19 Changed 7 years ago by cremona

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 7 years ago by wuthrich

  • 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 7 years ago by jdemeyer

  • Reviewers changed from Peter Bruin to Peter Bruin, Chris Wuthrich

comment:23 Changed 7 years ago by jdemeyer

  • Description modified (diff)

comment:24 Changed 7 years ago by jdemeyer

  • 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.

Changed 7 years ago by pbruin

fix documentation syntax

comment:25 Changed 7 years ago by pbruin

  • Description modified (diff)
  • Status changed from needs_work to needs_review

comment:26 Changed 7 years ago by jdemeyer

  • Status changed from needs_review to positive_review

comment:27 Changed 7 years ago by jdemeyer

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