Opened 10 years ago

Closed 9 years ago

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

Reported by: Owned by: was was critical sage-5.10 elliptic curves heights sage-5.10.beta3 John Cremona Peter Bruin, Chris Wuthrich N/A #12692

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
```

### comment:1 Changed 10 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 10 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 10 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 10 years ago by jen

• Stopgaps set to #12692

### comment:5 Changed 10 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 10 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 10 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 10 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 10 years ago by cremona

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

### comment:10 Changed 10 years ago by pbruin

• Reviewers set to Peter Bruin

### comment:11 Changed 10 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 10 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: ↓ 14 Changed 10 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 10 years ago by cremona

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

OK, all tests passed.

### comment:16 Changed 10 years ago by pbruin

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

### Changed 9 years ago by cremona

rebased to 5.9beta5

### comment:17 Changed 9 years ago by cremona

• Description modified (diff)

### comment:18 Changed 9 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: ↓ 20 Changed 9 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 9 years ago by cremona

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

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

### comment:23 Changed 9 years ago by jdemeyer

• Description modified (diff)

### comment:24 Changed 9 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 9 years ago by pbruin

fix documentation syntax

### comment:25 Changed 9 years ago by pbruin

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

### comment:26 Changed 9 years ago by jdemeyer

• Status changed from needs_review to positive_review

### comment:27 Changed 9 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.