computation of height of point on elliptic curve over Q(sqrt(5)) is WRONG

[was]

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​

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

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

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

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

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

[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).

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

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

[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!

[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?

[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?

[pbruin]

OK, all tests passed.

[pbruin]

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

[cremona]

rebased to 5.9beta5

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

[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?

[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)

[wuthrich]

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.

[jdemeyer]

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.

[pbruin]

fix documentation syntax