#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:
--------------------------------------------+-------------------------------
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.
>
> ----
>
> Apply only [attachment:trac_11767b.patch] to the Sage library.
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] and
[attachment:trac_11767c.patch] (in this order) to the Sage library.
--
Comment(by 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
--
Ticket URL: <http://trac.sagemath.org/sage_trac/ticket/11767#comment:23>
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.
