Opened 9 years ago

Closed 9 years ago

elliptic_logarithm of high precision points often hangs forever

Reported by: Owned by: was cremona major sage-4.7.2 elliptic curves ecc2011 sage-4.7.2.alpha3 Paul Zimmermann John Cremona, Leif Leonhardy, William Stein N/A

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

to the Sage library.

comment:1 Changed 9 years ago by leif

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

comment:2 Changed 9 years ago by was

• Status changed from new to needs_review

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.

comment:3 Changed 9 years ago by leif

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.

comment:4 Changed 9 years ago by cremona

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!

comment:5 Changed 9 years ago by was

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!

comment:6 Changed 9 years ago by zimmerma

another typo in that file:

```...the the returned value z...
```

Paul

Changed 9 years ago by zimmerma

this patch should be applied after the previous one

comment:7 Changed 9 years ago by zimmerma

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

comment:8 Changed 9 years ago by cremona

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?

comment:9 follow-up: ↓ 10 Changed 9 years ago by zimmerma

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

comment:10 in reply to: ↑ 9 Changed 9 years ago by leif

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.

comment:11 follow-up: ↓ 12 Changed 9 years ago by zimmerma

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

comment:12 in reply to: ↑ 11 Changed 9 years ago by cremona

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.

Changed 9 years ago by zimmerma

attach only that patch

comment:13 follow-up: ↓ 14 Changed 9 years ago by zimmerma

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

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

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.

comment:15 Changed 9 years ago by leif

• Authors set to Paul Zimmermann
• Description modified (diff)
• Reviewers set to John Cremona, Leif Leonhardy

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.

comment:16 Changed 9 years ago by leif

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

comment:17 follow-ups: ↓ 18 ↓ 19 Changed 9 years ago by zimmerma

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

comment:18 in reply to: ↑ 17 Changed 9 years ago by leif

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

comment:19 in reply to: ↑ 17 Changed 9 years ago by cremona

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

comment:20 follow-up: ↓ 21 Changed 9 years ago by cremona

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.

comment:21 in reply to: ↑ 20 ; follow-up: ↓ 22 Changed 9 years ago by leif

Try this:

Both work for me with Paul's patch.

comment:22 in reply to: ↑ 21 Changed 9 years ago by cremona

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

Changed 9 years ago by zimmerma

apply this patch after trac_11767b.patch

comment:23 follow-up: ↓ 24 Changed 9 years ago by zimmerma

• Description modified (diff)

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

comment:24 in reply to: ↑ 23 Changed 9 years ago by cremona

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!

comment:25 Changed 9 years ago by leif

• Description modified (diff)

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

comment:27 Changed 9 years ago by leif

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.

comment:28 Changed 9 years ago by was

• Status changed from needs_review to positive_review

final positive review.

comment:29 Changed 9 years ago by leif

• Reviewers changed from John Cremona, Leif Leonhardy to John Cremona, Leif Leonhardy, William Stein

comment:30 Changed 9 years ago by leif

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