#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] and
> [attachment:trac_11767c.patch] (in this order) 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
1. [attachment:trac_11767b.patch]
1. [attachment:trac_11767c.patch]
to the Sage library.
--
Comment(by leif):
FWIW, 11767b+c also passes doctests here.
[[BR]]
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.)
--
Ticket URL: <http://trac.sagemath.org/sage_trac/ticket/11767#comment:25>
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.
