#11767: elliptic_logarithm of high precision points often hangs forever
--------------------------------------------+-------------------------------
Reporter: was | Owner: cremona
Type: defect | Status: needs_review
Priority: major | Milestone: sage-4.7.2
Component: elliptic curves | Keywords:
Work_issues: | Upstream: N/A
Reviewer: John Cremona, Leif Leonhardy | Author: Paul Zimmermann
Merged: | Dependencies:
--------------------------------------------+-------------------------------
Changes (by newvalueoldvalue):
* reviewer: => John Cremona, Leif Leonhardy
* author: => Paul Zimmermann
Old description:
> 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.
New description:
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 only [attachment:trac_11767b.patch] to the Sage library.
--
Comment:
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.
--
Ticket URL: <http://trac.sagemath.org/sage_trac/ticket/11767#comment:15>
Sage <http://www.sagemath.org>
Sage: Creating a Viable Open Source Alternative to Magma, Maple, Mathematica,
and MATLAB
--
You received this message because you are subscribed to the Google Groups
"sage-trac" group.
To post to this group, send email to [email protected].
To unsubscribe from this group, send email to
[email protected].
For more options, visit this group at
http://groups.google.com/group/sage-trac?hl=en.
