Sage: Ticket #12509: computation of height of point on elliptic curve over Q(sqrt(5)) is WRONG
https://trac.sagemath.org/ticket/12509
<p>
There are evidently many examples in which computing <code>P.height()</code>, for <code>P</code> 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 <a class="missing wiki">NotImplementedError?</a> 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.
</p>
<pre class="wiki">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
</pre><p>
Apply: <a class="attachment" href="https://trac.sagemath.org/attachment/ticket/12509/trac12509-heights.patch" title="Attachment 'trac12509-heights.patch' in Ticket #12509">trac12509-heights.patch</a><a class="trac-rawlink" href="https://trac.sagemath.org/raw-attachment/ticket/12509/trac12509-heights.patch" title="Download"></a>, <a class="attachment" href="https://trac.sagemath.org/attachment/ticket/12509/trac12509-heights2.patch" title="Attachment 'trac12509-heights2.patch' in Ticket #12509">trac12509-heights2.patch</a><a class="trac-rawlink" href="https://trac.sagemath.org/raw-attachment/ticket/12509/trac12509-heights2.patch" title="Download"></a>, <a class="attachment" href="https://trac.sagemath.org/attachment/ticket/12509/trac12509-heights3.patch" title="Attachment 'trac12509-heights3.patch' in Ticket #12509">trac12509-heights3.patch</a><a class="trac-rawlink" href="https://trac.sagemath.org/raw-attachment/ticket/12509/trac12509-heights3.patch" title="Download"></a>
</p>
en-usSagehttps://trac.sagemath.org/chrome/site/logo_sagemath_trac.png
https://trac.sagemath.org/ticket/12509
Trac 1.1.6cremonaWed, 15 Feb 2012 09:24:51 GMT
https://trac.sagemath.org/ticket/12509#comment:1
https://trac.sagemath.org/ticket/12509#comment:1
<p>
I cannot remember who implemented this, but have looked at the code recently and it is using the standard Silverman algorithm.
</p>
<p>
I think it is likely that there is already a ticket for this.
</p>
TicketnbruinThu, 16 Feb 2012 06:52:01 GMT
https://trac.sagemath.org/ticket/12509#comment:2
https://trac.sagemath.org/ticket/12509#comment:2
<p>
A little <code>hg blame sage/schemes/elliptic_curves/ell_point.py</code> and <code>hg history sage/schemes/elliptic_curves/ell_point.py</code> shows that this code comes from <a class="closed ticket" href="https://trac.sagemath.org/ticket/8496" title="enhancement: Implement canonical heights for elliptic curves over number fields (closed: fixed)">#8496</a> , merged in sage-4.4.alpha1. The next change to this file happened with <a class="closed ticket" href="https://trac.sagemath.org/ticket/8827" title="defect: Cache heights of points on elliptic curves (closed: fixed)">#8827</a> 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:
</p>
<pre class="wiki">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]
</pre><p>
so the computation is probably not really wrong. It's probably just a matter of precision estimates being off.
</p>
TicketcremonaThu, 16 Feb 2012 12:02:04 GMT
https://trac.sagemath.org/ticket/12509#comment:3
https://trac.sagemath.org/ticket/12509#comment:3
<p>
This is also pretty damning:
</p>
<pre class="wiki">sage: (3*P).height()/P.height()
9.19316602595430
sage: (2*P).height()/P.height()
4.09522243754782
</pre><p>
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.
</p>
TicketjenMon, 19 Mar 2012 00:48:06 GMTstopgaps set
https://trac.sagemath.org/ticket/12509#comment:4
https://trac.sagemath.org/ticket/12509#comment:4
<ul>
<li><strong>stopgaps</strong>
set to <em>#12692</em>
</li>
</ul>
TicketcremonaTue, 08 Jan 2013 10:23:07 GMT
https://trac.sagemath.org/ticket/12509#comment:5
https://trac.sagemath.org/ticket/12509#comment:5
<p>
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
</p>
<pre class="wiki">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
</pre><p>
and you wonder how it is possible to get this wrong:
</p>
<pre class="wiki">sage: P.archimedian_local_height(prec=53) + P.nonarchimedian_local_height(prec=53)
1.35648516097
</pre><p>
I will make a patch now.
</p>
TicketcremonaTue, 08 Jan 2013 10:30:55 GMT
https://trac.sagemath.org/ticket/12509#comment:6
https://trac.sagemath.org/ticket/12509#comment:6
<p>
I was a little hasty. There's no problem in the non-arch local height being an exact log:
</p>
<pre class="wiki">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
</pre><p>
is all correct. The problem is with this:
</p>
<pre class="wiki">sage: P.archimedian_local_height()
-0.2206607955468278522013468194382
sage: P.archimedian_local_height(prec=53)
-0.252952751464
</pre><p>
where the first is correct as this shows:
</p>
<pre class="wiki">sage: P.archimedian_local_height()+ P.nonarchimedian_local_height(prec=53)
1.38877711688727
</pre><p>
but the second is not.
</p>
TicketcremonaTue, 08 Jan 2013 12:46:21 GMT
https://trac.sagemath.org/ticket/12509#comment:7
https://trac.sagemath.org/ticket/12509#comment:7
<p>
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:
</p>
<pre class="wiki">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
</pre><p>
while at 100-bit precision:
</p>
<pre class="wiki">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
</pre><p>
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.
</p>
<p>
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.
</p>
<p>
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).
</p>
TicketcremonaWed, 09 Jan 2013 10:56:36 GMT
https://trac.sagemath.org/ticket/12509#comment:8
https://trac.sagemath.org/ticket/12509#comment:8
<p>
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.
</p>
TicketcremonaWed, 09 Jan 2013 10:58:10 GMTstatus, component, description changed; keywords, author set
https://trac.sagemath.org/ticket/12509#comment:9
https://trac.sagemath.org/ticket/12509#comment:9
<ul>
<li><strong>status</strong>
changed from <em>new</em> to <em>needs_review</em>
</li>
<li><strong>keywords</strong>
<em>heights</em> added
</li>
<li><strong>component</strong>
changed from <em>number theory</em> to <em>elliptic curves</em>
</li>
<li><strong>description</strong>
modified (<a href="/ticket/12509?action=diff&version=9">diff</a>)
</li>
<li><strong>author</strong>
set to <em>John Cremona</em>
</li>
</ul>
TicketpbruinThu, 10 Jan 2013 15:20:29 GMTreviewer set
https://trac.sagemath.org/ticket/12509#comment:10
https://trac.sagemath.org/ticket/12509#comment:10
<ul>
<li><strong>reviewer</strong>
set to <em>Peter Bruin</em>
</li>
</ul>
TicketpbruinFri, 11 Jan 2013 12:54:16 GMT
https://trac.sagemath.org/ticket/12509#comment:11
https://trac.sagemath.org/ticket/12509#comment:11
<p>
The original precision problem seems to be solved by the patch, but there is a remaining problem if the given place <code>v</code> has lower precision than <code>prec</code> (the answer ends with a string of zeroes):
</p>
<pre class="wiki">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
</pre><p>
The precision of <code>v</code> shouldn't influence the result, of course.
</p>
<p>
Also, the default value of <code>prec</code> used when applying Silverman's Theorem 4.2 is currently the precision of the refined place <code>vv</code>, but it suffices to take the same default precision to which the result is returned, namely the precision of <code>v</code>.
</p>
<p>
I wrote a patch that fixes both issues; I am now going to test it.
</p>
TicketcremonaFri, 11 Jan 2013 13:28:47 GMT
https://trac.sagemath.org/ticket/12509#comment:12
https://trac.sagemath.org/ticket/12509#comment:12
<p>
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.
</p>
<p>
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.
</p>
<p>
The second patch looks good to me, hope the tests pass!
</p>
TicketpbruinFri, 11 Jan 2013 13:39:02 GMT
https://trac.sagemath.org/ticket/12509#comment:13
https://trac.sagemath.org/ticket/12509#comment:13
<p>
The tests in sage/schemes/elliptic_curves/ all pass for me, now testing the whole Sage library.
</p>
<p>
If this also succeeds, do my additional changes still allow me to give a positive review, or do we need a new reviewer?
</p>
TicketcremonaFri, 11 Jan 2013 14:10:03 GMT
https://trac.sagemath.org/ticket/12509#comment:14
https://trac.sagemath.org/ticket/12509#comment:14
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/12509#comment:13" title="Comment 13">pbruin</a>:
</p>
<blockquote class="citation">
<p>
The tests in sage/schemes/elliptic_curves/ all pass for me, now testing the whole Sage library.
</p>
<p>
If this also succeeds, do my additional changes still allow me to give a positive review, or do we need a new reviewer?
</p>
</blockquote>
<p>
Maybe, maybe not. William?
</p>
TicketpbruinFri, 11 Jan 2013 19:22:32 GMT
https://trac.sagemath.org/ticket/12509#comment:15
https://trac.sagemath.org/ticket/12509#comment:15
<p>
OK, all tests passed.
</p>
TicketpbruinTue, 15 Jan 2013 00:10:55 GMTattachment set
https://trac.sagemath.org/ticket/12509
https://trac.sagemath.org/ticket/12509
<ul>
<li><strong>attachment</strong>
set to <em>trac12509-heights2.patch</em>
</li>
</ul>
TicketpbruinTue, 15 Jan 2013 00:13:57 GMT
https://trac.sagemath.org/ticket/12509#comment:16
https://trac.sagemath.org/ticket/12509#comment:16
<p>
trac12509-heights2.patch now also makes the documentation slightly more accurate
</p>
TicketcremonaMon, 22 Apr 2013 13:42:44 GMTattachment set
https://trac.sagemath.org/ticket/12509
https://trac.sagemath.org/ticket/12509
<ul>
<li><strong>attachment</strong>
set to <em>trac12509-heights.patch</em>
</li>
</ul>
<p>
rebased to 5.9beta5
</p>
TicketcremonaMon, 22 Apr 2013 13:43:18 GMTdescription changed
https://trac.sagemath.org/ticket/12509#comment:17
https://trac.sagemath.org/ticket/12509#comment:17
<ul>
<li><strong>description</strong>
modified (<a href="/ticket/12509?action=diff&version=17">diff</a>)
</li>
</ul>
TicketcremonaMon, 22 Apr 2013 13:47:46 GMT
https://trac.sagemath.org/ticket/12509#comment:18
https://trac.sagemath.org/ticket/12509#comment:18
<p>
I did a small rebase of the first patch for 5.9.beta5; both patches still need to be applied. We need a reviewer.
</p>
TicketpbruinMon, 22 Apr 2013 14:05:14 GMT
https://trac.sagemath.org/ticket/12509#comment:19
https://trac.sagemath.org/ticket/12509#comment:19
<p>
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.
</p>
<p>
Should I rebase my patch too?
</p>
TicketcremonaMon, 22 Apr 2013 15:39:35 GMT
https://trac.sagemath.org/ticket/12509#comment:20
https://trac.sagemath.org/ticket/12509#comment:20
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/12509#comment:19" title="Comment 19">pbruin</a>:
</p>
<blockquote class="citation">
<p>
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.
</p>
<p>
Should I rebase my patch too?
</p>
</blockquote>
<p>
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)
</p>
TicketwuthrichSat, 04 May 2013 19:22:29 GMTstatus changed
https://trac.sagemath.org/ticket/12509#comment:21
https://trac.sagemath.org/ticket/12509#comment:21
<ul>
<li><strong>status</strong>
changed from <em>needs_review</em> to <em>positive_review</em>
</li>
</ul>
<p>
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.
</p>
TicketjdemeyerSat, 04 May 2013 21:03:21 GMTreviewer changed
https://trac.sagemath.org/ticket/12509#comment:22
https://trac.sagemath.org/ticket/12509#comment:22
<ul>
<li><strong>reviewer</strong>
changed from <em>Peter Bruin</em> to <em>Peter Bruin, Chris Wuthrich</em>
</li>
</ul>
TicketjdemeyerTue, 07 May 2013 09:56:50 GMTdescription changed
https://trac.sagemath.org/ticket/12509#comment:23
https://trac.sagemath.org/ticket/12509#comment:23
<ul>
<li><strong>description</strong>
modified (<a href="/ticket/12509?action=diff&version=23">diff</a>)
</li>
</ul>
TicketjdemeyerTue, 07 May 2013 12:27:43 GMTstatus changed
https://trac.sagemath.org/ticket/12509#comment:24
https://trac.sagemath.org/ticket/12509#comment:24
<ul>
<li><strong>status</strong>
changed from <em>positive_review</em> to <em>needs_work</em>
</li>
</ul>
<p>
There is a documentation problem:
</p>
<pre class="wiki">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.
</pre>
TicketpbruinFri, 10 May 2013 07:17:05 GMTattachment set
https://trac.sagemath.org/ticket/12509
https://trac.sagemath.org/ticket/12509
<ul>
<li><strong>attachment</strong>
set to <em>trac12509-heights3.patch</em>
</li>
</ul>
<p>
fix documentation syntax
</p>
TicketpbruinFri, 10 May 2013 07:17:57 GMTstatus, description changed
https://trac.sagemath.org/ticket/12509#comment:25
https://trac.sagemath.org/ticket/12509#comment:25
<ul>
<li><strong>status</strong>
changed from <em>needs_work</em> to <em>needs_review</em>
</li>
<li><strong>description</strong>
modified (<a href="/ticket/12509?action=diff&version=25">diff</a>)
</li>
</ul>
TicketjdemeyerFri, 10 May 2013 08:28:13 GMTstatus changed
https://trac.sagemath.org/ticket/12509#comment:26
https://trac.sagemath.org/ticket/12509#comment:26
<ul>
<li><strong>status</strong>
changed from <em>needs_review</em> to <em>positive_review</em>
</li>
</ul>
TicketjdemeyerMon, 13 May 2013 13:25:55 GMTstatus changed; resolution, merged set
https://trac.sagemath.org/ticket/12509#comment:27
https://trac.sagemath.org/ticket/12509#comment:27
<ul>
<li><strong>status</strong>
changed from <em>positive_review</em> to <em>closed</em>
</li>
<li><strong>resolution</strong>
set to <em>fixed</em>
</li>
<li><strong>merged</strong>
set to <em>sage-5.10.beta3</em>
</li>
</ul>
Ticket