elliptic_logarithm of high precision points often hangs forever

[williamstein]

I am doing a project to compute Chow-Heegner points with Darmon, Rotger, et al., and am using Sage's elliptic_logarithm function with high precision points as input. Unfortunately, it completely hangs on many input points over number fields. Here is an example below, where computing to precision 500 works fine, but precision 600 hangs Sage forever:

sage: R.<x> = QQ[]
sage: K.<a> = NumberField(x^2 + x + 5)
sage: E = EllipticCurve(K, [0,0,1,-3,-5])
sage: P = E([0,a])
sage: L = P.curve().period_lattice(K.embeddings(ComplexField(600))[0])
sage: time L.elliptic_logarithm(P, prec=500)

-0.842248166487739393375018008381693990800588864069506187033873183845246233548058477561706400464057832396643843146464236956684557207157300006542470428494 - 0.571366031453267388121279381354098224265947866751130917440598461117775339240176310729173301979590106474259885638797913383502735083088736326391919063211*I
Time: CPU 0.08 s, Wall: 0.09 s

BUT:

sage: L.elliptic_logarithm(P, prec=600)
HANGS FOREVER

Hitting control-c and using the debugger suggests that the termination condition is impossible. You end up in this line if (r.abs()-1).abs() < eps: break with actually having (r.abs()-1).abs() == eps, so the strict inequality isnt' satisfied, and maybe for some reason it can't be???

ipdb> l
   1357             r = C(((xP-e3)/(xP-e2)).sqrt())
   1358             if r.real()<0: r=-r
   1359             t = -C(wP)/(2*r*(xP-e2))
   1360             eps = R(1)>>(prec2);
   1361             while True:        
-> 1362                 s = b*r+a
   1363                 a, b = (a+b)/2, (a*b).sqrt()
   1364                 if (a+b).abs() < (a-b).abs():  b=-b
   1365                 r = (a*(r+1)/s).sqrt()
   1366                 if (r.abs()-1).abs() < eps: break
   1367                 if r.real()<0: r=-r

ipdb> print eps
5.80771375621750318328344999898952221581714435905885826948966319082059051450573607513973642929306226703333487164506974570166210185861027247933383445853709821724799129294394972880718063974478259360634645856449058721344643691320490207325897321618392592839150233975392885056947380044108965624364302537017698344822203678107744035231793749824965395444912652822986530e-362
ipdb> print r    
1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
ipdb> print (r.abs()-1).abs()
5.80771375621750318328344999898952221581714435905885826948966319082059051450573607513973642929306226703333487164506974570166210185861027247933383445853709821724799129294394972880718063974478259360634645856449058721344643691320490207325897321618392592839150233975392885056947380044108965624364302537017698344822203678107744035231793749824965395444912652822986530e-362
ipdb> print eps
5.80771375621750318328344999898952221581714435905885826948966319082059051450573607513973642929306226703333487164506974570166210185861027247933383445853709821724799129294394972880718063974478259360634645856449058721344643691320490207325897321618392592839150233975392885056947380044108965624364302537017698344822203678107744035231793749824965395444912652822986530e-362
ipdb> print (r.abs()-1).abs() < eps
False
ipdb> print (r.abs()-1).abs() <= eps
True

Changing < eps to <= eps in two spots in the relevant file seems to fix the problem for me.

---

Apply

1. attachment: trac_11767b.patch
2. attachment: trac_11767c.patch to the Sage library.

Component: elliptic curves

Keywords: ecc2011

Author: Paul Zimmermann

Reviewer: John Cremona, Leif Leonhardy, William Stein

Merged: sage-4.7.2.alpha3

Issue created by migration from https://trac.sagemath.org/ticket/11767

[83660e46-0051-498b-a8c1-f7a7bd232b5a]

comment:1

In case you provide a patch, you could also fix this:

        # embeddings uses only real arithmetic i nthe iteraion, and is

(And remove the superfluous semicolon from eps = R(1)>>(prec2); twice.)

[williamstein]

comment:2

I do not claim to understand the < versus <= issue, since I do not know the algorithm in question. But I can attach a patch with this fix, which fixes the problem, has an example to illustrate that the problem is fixed, and fixes the typos leif mentions above. Perhaps this will be easy for John Cremona to referee.

[williamstein]

Attachment: trac_11767.patch.gz

[83660e46-0051-498b-a8c1-f7a7bd232b5a]

comment:3

FWIW, using

        R = RealField(prec2+1)
        C = ComplexField(prec2+1)

instead also works for (and seems a bit cleaner to) me, though of course a bit more expensive.

John will certainly know better.

[JohnCremona]

comment:4

I cannot fix this now and maybe someone else should anyway, since my numerical analysis expertise is clearly lacking -- apologies. The underlying algorithm here is very good indeed (double exponential) but the stopping condition needs to be handled properly.

In the function elliptic_logarithm, first of all the working precision is doubled (prec2 is 2*prec)! Overkill perhaps but with this algorithm the number of correct digits doubles with each step anyway. The stopping condition (line 1366) is that the relative error is < eps where eps is 2^-prec2. In the loop, a and b are sequences of complex numbers which converge to a common limit and we ant to stop when a/b is close to 1. Similarly on line 1400 where a,b are real.

I think that the right thing to do would be to replace eps by 2^-(prec2-1) or something similar. Also, replacing "while True" by something else so that infinite looping is impossible?

I'm glad this function is useful, by the way -- there is no other software which can take complex elliptic logs!

[williamstein]

comment:5

I'm glad this function is useful, by the way -- there is no other software which can take complex elliptic logs!

Thanks -- it's absolutely critical for my research right now!

[zimmermann6]

comment:6

another typo in that file:

...the the returned value z...

Paul

[zimmermann6]

this patch should be applied after the previous one

[zimmermann6]

comment:7

Attachment: trac_11767a.patch.gz

I have added a 2nd patch (to be applied after the first one) which solves the typo, and modifies the termination condition.

As far as I understand, the termination is when the (complex) variable r becomes 1. We compute the complex norm of r, subtract one from it, and compare to eps which is 2^(-prec2).

If r is near one, typically r.abs() will be 1, or one of the floating-point numbers near from 1, which are 1-2^(-prec2) and 1+2^(1-prec2). The subtraction of 1 is exact, and the only case where the comparison r.abs()-1 < eps can be true is when r=1.

What I've done is that I've changed the termination condition to also check whether the relative difference r.abs()-1 does not change from one step to the next one. Since r.abs()-1 converges to one, this should catch cases where it stays at say 1-2^(-prec2) due to rounding errors. However it will not catch cases where it oscillates between two different values near from 1.

Paul

[JohnCremona]

comment:8

Paul, thanks. I trust your knowledge of binary floating point issues as being much better than my own. Are you saying that even with your patch the termination condition may never hold and hence we loop forever?

Bearing in mind that I already doubled the precision (prec2 = 2*prec) should we perhaps change 1-2^(-prec2) to 1-2^(-prec2+1)? If that would mean that termination is then guaranteed?

[zimmermann6]

comment:9

John, if you change eps to 2^(-prec2+1) instead of 2^(-prec2) then the case where r = 1 - 2^(-prec2) will also terminate with the < condition. But the case where r = 1 + 2^(-prec2+1) will still loop.

If the target precision is prec and you compute internally with 2*prec bits, then eps could be 2^(-prec), no?

Paul

[83660e46-0051-498b-a8c1-f7a7bd232b5a]

comment:10

Replying to @zimmermann6:

John, if you change eps to 2^(-prec2+1) instead of 2^(-prec2) then the case where r = 1 - 2^(-prec2) will also terminate with the < condition. But the case where r = 1 + 2^(-prec2+1) will still loop.

Well, either r.abs() converges to one, or it doesn't. So if you compute with sufficiently high precision, especially not choosing epsilon to be the minimal positive value of the floating point type used, the termination condition will get true after some iterations.

[zimmermann6]

comment:11

ok, after looking again, it seems that using prec2 = 2*prec is overkill. Usually, what we do in MPFR is that we use as internal precision prec + log2(prec) + a few bits. To avoid computing a log(), I suggest using prec2 = prec + 40 which should be large enough. With this change, William's example with prec=600 runs ok. John, if you agree with that, I can make a new patch.

Paul

[JohnCremona]

comment:12

Replying to @zimmermann6:

ok, after looking again, it seems that using prec2 = 2*prec is overkill. Usually, what we do in MPFR is that we use as internal precision prec + log2(prec) + a few bits. To avoid computing a log(), I suggest using prec2 = prec + 40 which should be large enough. With this change, William's example with prec=600 runs ok. John, if you agree with that, I can make a new patch.

Paul

I agree, that is fine. I put in 2*prec for testing and it worked (well, for several months!) and the algorithm was fast enough that I never bothered to work out a more sensible increase.

After today (Thursday 15/9/11) I'll be away for 8 days so please don't wait for my approval to give the new patch a positive review assuming that the tests pass, including the new doctest with William's example.

[zimmermann6]

attach only that patch

[zimmermann6]

comment:13

Attachment: trac_11767b.patch.gz

John, thank you for your feedback. I've now uploaded a new patch trac_11767b.patch which should be applied solely. William, please could you review that ticket?

Paul

[JohnCremona]

comment:14

Replying to @zimmermann6:

John, thank you for your feedback. I've now uploaded a new patch trac_11767b.patch which should be applied solely. William, please could you review that ticket?

Paul

The patch looks fine to me though I have not tested it. With a Positive Review the testbot will test it thoroughly anyway, but I think that at least one human who is not the author should test it too.

[83660e46-0051-498b-a8c1-f7a7bd232b5a]

comment:15

Well, I think even I can review this. ;-)

By the way, I guess the double "the" should have read "then the", i.e., an "n" was missing rather than a word duplicated, but it's ok as is, too.

So I think we have a positive review.

Of course William could again stress-test this, but at least there should no longer be endless loops, and the algorithm converges fast (with an unmodified termination condition, except that epsilon is now larger for high precision inputs), such that I don't expect very different or much less precise results w.r.t. the previous version either. Note that epsilon will now even be smaller for lower precision inputs.

[83660e46-0051-498b-a8c1-f7a7bd232b5a]

Reviewer: John Cremona, Leif Leonhardy

[83660e46-0051-498b-a8c1-f7a7bd232b5a]

Description changed:

--- 
+++ 
@@ -50,4 +50,9 @@
 True

-Changing < eps to <= eps in two spots in the relevant file seems to fix the problem for me. +Changing < eps to <= eps in two spots in the relevant file seems to fix the problem for me. + +--- + +Apply only attachment: trac_11767b.patch to the Sage library. +

[83660e46-0051-498b-a8c1-f7a7bd232b5a]

Author: Paul Zimmermann

[83660e46-0051-498b-a8c1-f7a7bd232b5a]

comment:16

Oh, the superfluous semicolons survived Paul's patch, but I don't mind.

[zimmermann6]

comment:17

Leif,

So I think we have a positive review.

you didn't hit "positive review" thus you are waiting for William's review?

Note that William's example seems to exercise only the case where r is real in the code. Maybe one should check that there is at least one example where r is complex with large precision, or add one if it is not the case.

Paul

[83660e46-0051-498b-a8c1-f7a7bd232b5a]

comment:18

Replying to @zimmermann6:

Leif,

So I think we have a positive review.

you didn't hit "positive review" thus you are waiting for William's review?

Yes, it's 4 am over there...

[JohnCremona]

comment:19

Replying to @zimmermann6:

Leif,

So I think we have a positive review.

you didn't hit "positive review" thus you are waiting for William's review?

Note that William's example seems to exercise only the case where r is real in the code. Maybe one should check that there is at least one example where r is complex with large precision, or add one if it is not the case.

That is a very good observation. The real case algorithm has been known for a very long time, but the complex case (using a generalised complex AGM) is still in preprint form! I will come up with an example.

Paul

[JohnCremona]

comment:20

Try this:

sage: K.<a> = QuadraticField(-5)
sage: E = EllipticCurve([1,1,a,a,0])
sage: P = E(0,0)
sage: P.order()
+Infinity
sage: L = P.curve().period_lattice(K.embeddings(ComplexField())[0])
sage: time L.elliptic_logarithm(P, prec=500)
1.17058357737548897849026170185581196033579563441850967539191867385734983296504066660506637438866628981886518901958717288150400849746892393771983141354 - 1.13513899565966043682474529757126359416758251309237866586896869548539516543734207347695898664875799307727928332953834601460994992792519799260968053875*I
Time: CPU 0.35 s, Wall: 0.36 s
sage: L = P.curve().period_lattice(K.embeddings(ComplexField(1000))[0])
sage: time L.elliptic_logarithm(P, prec=1000)

The last line hangs with vanilla 4.7.1. prec=900 is still OK. Note that it is irrelevant to specify the precision when constructing the period lattice, since the only use made of the embedding at that stage is to determine which emebedding is required; the elliptic log computation is where the floating point work is done.

[83660e46-0051-498b-a8c1-f7a7bd232b5a]

comment:21

Replying to @JohnCremona:

Try this:

Both work for me with Paul's patch.

[JohnCremona]

comment:22

Replying to @nexttime:

Replying to @JohnCremona:

Try this:

Both work for me with Paul's patch.

Thanks for testing. Paul, could you add the new example to the tests? You can use exactly the code I provided though the line

sage: L = P.curve().period_lattice(K.embeddings(ComplexField(1000))[0])

can be omitted for the reasons I explained (the lattice only needs to be defined once, with no precision specified).

[zimmermann6]

Attachment: trac_11767c.patch.gz

apply this patch after trac_11767b.patch

[zimmermann6]

Description changed:

--- 
+++ 
@@ -54,5 +54,6 @@

 ---

-Apply only [attachment: trac_11767b.patch](https://github.com/sagemath/sage-prod/files/10653607/trac_11767b.patch.gz) to the Sage library.
+Apply only [attachment: trac_11767b.patch](https://github.com/sagemath/sage-prod/files/10653607/trac_11767b.patch.gz) and
+[attachment: trac_11767c.patch](https://github.com/sagemath/sage-prod/files/10653608/trac_11767c.patch.gz) (in this order) to the Sage library.

[zimmermann6]

comment:23

I knew I should not have opened my mouth... Attached is a second patch. Don't ask me to make a unified patch, I prefer to spend my time on other tickets.

Paul

[JohnCremona]

comment:24

Replying to @zimmermann6:

I knew I should not have opened my mouth... Attached is a second patch. Don't ask me to make a unified patch, I prefer to spend my time on other tickets.

Paul

Thank you!

[83660e46-0051-498b-a8c1-f7a7bd232b5a]

comment:25

FWIW, 11767b+c also passes doctests here.

P.S.:

There are various ways to make one patch out of two. E.g., take a fresh branch, import the first with --no-commit, and the second with -f(orce), then commit the cumulative changes. A diff is even easier, just compare the tip to the changeset two below the tip. (If you export two changesets at once, they will end up in one file, but effectively be two patches which are simply concatenated, but you can in turn re-import such a patch as one changeset into a fresh branch and re-export it as a true single patch.)

[83660e46-0051-498b-a8c1-f7a7bd232b5a]

Description changed:

--- 
+++ 
@@ -54,6 +54,8 @@

 ---

-Apply only [attachment: trac_11767b.patch](https://github.com/sagemath/sage-prod/files/10653607/trac_11767b.patch.gz) and
-[attachment: trac_11767c.patch](https://github.com/sagemath/sage-prod/files/10653608/trac_11767c.patch.gz) (in this order) to the Sage library.
+Apply
+1. [attachment: trac_11767b.patch](https://github.com/sagemath/sage-prod/files/10653607/trac_11767b.patch.gz)
+2. [attachment: trac_11767c.patch](https://github.com/sagemath/sage-prod/files/10653608/trac_11767c.patch.gz)
+to the Sage library.

[zimmermann6]

Changed keywords from none to ecc2011

[83660e46-0051-498b-a8c1-f7a7bd232b5a]

comment:27

Can someone (William?) give this a final positive review?

If this happens within a few hours, I can still merge it into Sage 4.7.2.alpha3.

[williamstein]

comment:28

final positive review.

[83660e46-0051-498b-a8c1-f7a7bd232b5a]

Changed reviewer from John Cremona, Leif Leonhardy to John Cremona, Leif Leonhardy, William Stein

[83660e46-0051-498b-a8c1-f7a7bd232b5a]

Merged: sage-4.7.2.alpha3