Opened 9 years ago
Closed 9 years ago
#11767 closed defect (fixed)
elliptic_logarithm of high precision points often hangs forever
Reported by: | was | Owned by: | cremona |
---|---|---|---|
Priority: | major | Milestone: | sage-4.7.2 |
Component: | elliptic curves | Keywords: | ecc2011 |
Cc: | Merged in: | sage-4.7.2.alpha3 | |
Authors: | Paul Zimmermann | Reviewers: | John Cremona, Leif Leonhardy, William Stein |
Report Upstream: | N/A | Work issues: | |
Branch: | Commit: | ||
Dependencies: | Stopgaps: |
Description (last modified by )
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.
Attachments (4)
Change History (34)
comment:1 Changed 9 years ago by
comment:2 Changed 9 years ago by
- 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.
Changed 9 years ago by
comment:3 Changed 9 years ago by
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
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
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
another typo in that file:
...the the returned value z...
Paul
comment:7 Changed 9 years ago by
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
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
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
Replying to zimmerma:
John, if you change
eps
to2^(-prec2+1)
instead of2^(-prec2)
then the case wherer = 1 - 2^(-prec2)
will also terminate with the < condition. But the case wherer = 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
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
Replying to 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 usingprec2 = 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.
comment:13 follow-up: ↓ 14 Changed 9 years ago by
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
Replying to 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
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
- 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
Oh, the superfluous semicolons survived Paul's patch, but I don't mind.
comment:17 follow-ups: ↓ 18 ↓ 19 Changed 9 years ago by
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
Replying to zimmerma:
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
Replying to 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 wherer
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
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
comment:22 in reply to: ↑ 21 Changed 9 years ago by
Replying to leif:
Replying to 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).
comment:23 follow-up: ↓ 24 Changed 9 years ago by
- 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
Replying to zimmerma:
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
- 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:26 Changed 9 years ago by
- Keywords ecc2011 added
comment:27 Changed 9 years ago by
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
- Status changed from needs_review to positive_review
final positive review.
comment:29 Changed 9 years ago by
- Reviewers changed from John Cremona, Leif Leonhardy to John Cremona, Leif Leonhardy, William Stein
comment:30 Changed 9 years ago by
- Merged in set to sage-4.7.2.alpha3
- Resolution set to fixed
- Status changed from positive_review to closed
In case you provide a patch, you could also fix this:
(And remove the superfluous semicolon from
eps = R(1)>>(prec2);
twice.)