Opened 9 years ago

Closed 9 years ago

#15483 closed defect (fixed)

Simon 2-descent gives RuntimeError for an elliptic curve over a quadratic field

Reported by: Peter Bruin Owned by:
Priority: major Milestone: sage-6.2
Component: elliptic curves Keywords: simon_two_descent
Cc: John Cremona Merged in:
Authors: Peter Bruin Reviewers: Jeroen Demeyer
Report Upstream: Workaround found; Bug reported upstream. Work issues:
Branch: b662aba (Commits, GitHub, GitLab) Commit: b662aba98ad4122a82fc135cf9b5dbf58de2ba0e
Dependencies: #11005 Stopgaps:

Status badges

Description (last modified by Peter Bruin)

[See #15608 for a list of open simon_two_descent tickets]

The following tries to do 2-descent on an elliptic curve over a quadratic field:

K.<s>=QuadraticField(229)
c4=2173-235*(1-s)/2
c6=-124369+15988*(1-s)/2
E=EllipticCurve([-c4/48,-c6/864])
E.simon_two_descent()

It fails with a RuntimeError:

Traceback (most recent call last):
...
RuntimeError: 
  ***   at top-level: ans=bnfellrank(K,[0,0,0,
  ***                     ^--------------------
  ***   in function bnfellrank: ...eqtheta,rnfeq,bbnf];rang=
  ***   bnfell2descent_gen(b
  ***   ^--------------------
  ***   in function bnfell2descent_gen: ...riv,r=nfsqrt(nf,norm(zc))
  ***   [1];if(DEBUGLEVEL_el
  ***   ^--------------------
  ***   array index (1) out of allowed range [none].
An error occurred while running Simon's 2-descent program

The problem is apparently that the square root of norm(zc) (which is indeed a square) is not computed due to a precision bug; see comment:7 below.

Change History (23)

comment:1 Changed 9 years ago by John Cremona

There are many similar bug reports with separate trac tickets. I have no idea if this is a new bug or not. Anyway, the version Sage has of DS's scripts are not the most recent one so there is no point in reporting this back to him unless you check that it happens with the scripts on his web page and you email him without mentioning Sage, otherwise he will just assume that the problem lies with Sage (which is possible).

Several of us spend hours or even days in 2011 sage Days (more than one) working on this sort of thing and gave up. We were trying to do more, namely rewrite the interface so that it can use a gp2c-compiled version of the scripts. One day I'll have another go at that.

Meanwhile, this entry is to make sure that no-one wastes any time at all looking into the versions of the scripts which Sage currently includes; you should test against the current versions.

comment:2 in reply to:  1 Changed 9 years ago by Peter Bruin

Description: modified (diff)

Replying to cremona:

There are many similar bug reports with separate trac tickets. I have no idea if this is a new bug or not. Anyway, the version Sage has of DS's scripts are not the most recent one so there is no point in reporting this back to him unless you check that it happens with the scripts on his web page and you email him without mentioning Sage, otherwise he will just assume that the problem lies with Sage (which is possible).

I searched for other tickets related to Simon's code and did not find reports of this particular problem. We should first upgrade to the latest version of his script and then see which of the various problems listed at #11041 persist. I added this ticket to that list for easier reference.

comment:3 Changed 9 years ago by wuthrich

Report Upstream: N/ANot yet reported upstream; Will do shortly.

Lorenzini reported me a similar bug and I had a look. I now believe that this is a bug in Simon's script as I can reproduce it outside sage. I used the version of his file dated 06/04/2011, which is the newest I found on the web, while sage has one dated the 25/03/2009. So besides that bug, the file needs to be updated anyway. I tried both in sage's gp (2.5.5) and the current alpha (2.6.1).

I send an email to Denis and see if and how fast we can hope for a solution of this. Otherwise, I may at least update the script, later.

Lorenzini's example is a bit simpler:

R.<t> = QQ[]
L.<g> = NumberField(t^3 - 9*t^2 + 13*t - 4)
E1 = EllipticCurve(L,[1-g*(g-1),-g^2*(g-1),-g^2*(g-1),0,0])
E1.rank()

comment:4 Changed 9 years ago by John Cremona

Description: modified (diff)

comment:5 Changed 9 years ago by For batch modifications

Milestone: sage-6.1sage-6.2

comment:6 Changed 9 years ago by Peter Bruin

Dependencies: #11005

The problem still exists after #11005, which upgrades Simon's script to the latest version.

comment:7 Changed 9 years ago by Peter Bruin

Description: modified (diff)

It seems that the bug is caused by a precision problem in subst, called by nfrealsign in ell.gp. The following GP input demonstrates this (using the latest version of ell.gp, see #11005):

\r src/ext/pari/simon/ell.gp
b = -1554544300737274875964190134520312870631312460283689944298138572669148295776039072867720281361776956435252620954745928376624817557704277432961924925312*y + 23524523971732905757341977352314040726186200302188191824300117738073539522011689544444863977622786771332621915440577829842674416407299864303146477224320
a = Mod(b, y^2 - 229);
\p 38                    \\ default precision
K = bnfinit(y^2 - 229);
nfrealsign(K, a, 1)      \\ result: 1
nfrealsign(K, a, 2)      \\ result: -1  (wrong)
subst(b, y, K.roots[2])  \\ result: -7.695704335233296721 E112, incorrect to this precision
\p 500                   \\ increase precision
K = bnfinit(y^2 - 229);
nfrealsign(K, a, 1)      \\ result: 1
nfrealsign(K, a, 2)      \\ result: 1
subst(b, y, K.roots[2])  \\ result: 55550556624985845118007242443189926820306719407.796355955086894963616120700422308457186069825115395206111444495243228840597412527828358943327267422898011087182852030228685210492253493260963313348652457717717703991690402543456279919008587884729307897541432624377818611055431792706380900043638904108307170236335925883131494247446074384269328376967381902948861789570537169455258630011202296209679102466188417

The value of K.roots does not seem to be the problem; the roots are calculated correctly to the requested precision.

comment:8 Changed 9 years ago by Peter Bruin

The basic problem seems to be that PARI is too optimistic about the resulting precision when adding real numbers:

? b = -1554544300737274875964190134520312870631312460283689944298138572669148295776039072867720281361776956435252620954745928376624817557704277432961924925312*y + 23524523971732905757341977352314040726186200302188191824300117738073539522011689544444863977622786771332621915440577829842674416407299864303146477224320
? K = bnfinit(y^2 - 229);
? u=polcoeff(b,1)*K.roots[2]
%167 = -2.3524523971732905757341977352314040726 E151
? v=polcoeff(b,0)*1.
%168 = 2.3524523971732905757341977352314040726 E151
? u+v
%169 = -7.695704335233296721 E112
? precision(u)
%170 = 38
? precision(v)
%171 = 38
? precision(u+v)
%172 = 19

comment:9 Changed 9 years ago by Peter Bruin

Report Upstream: Not yet reported upstream; Will do shortly.Reported upstream. No feedback yet.

E-mailed pari-dev and Denis Simon about the precision issue, also suggesting a quick fix for ell.gp.

comment:10 Changed 9 years ago by Peter Bruin

The quick fix is to increase the bound in the precision check in nfrealsign from 10 digits to 38 digits (128 bits). This at least makes the error in my example and that reported by Lorenzini disappear.

Last edited 9 years ago by Peter Bruin (previous) (diff)

comment:11 Changed 9 years ago by Peter Bruin

Branch: u/pbruin/15483-two_descent_precision
Commit: d34a77ab7abe4072349e10e0154a598da94b7c31
Status: newneeds_review

comment:12 Changed 9 years ago by Peter Bruin

Status: needs_reviewneeds_work

This needs a doctest...

comment:13 Changed 9 years ago by git

Commit: d34a77ab7abe4072349e10e0154a598da94b7c31b662aba98ad4122a82fc135cf9b5dbf58de2ba0e

Branch pushed to git repo; I updated commit sha1. New commits:

b662abaadd/fix doctests for simon_two_descent

comment:14 Changed 9 years ago by Peter Bruin

Status: needs_workneeds_review

comment:15 Changed 9 years ago by Peter Bruin

Authors: Peter Bruin
Report Upstream: Reported upstream. No feedback yet.Workaround found; Bug reported upstream.

comment:16 Changed 9 years ago by Jeroen Demeyer

I quite like the suggestion by Karim Belabas to use nfsign() together with nfnewprec. I think that's a lot better than the current patch, which is really just a hack.

Maybe you should ask Denis Simon what he thinks (perhaps with a proof-of-concept patch)? Ideally, he fixes his script "upstream".

comment:17 in reply to:  16 Changed 9 years ago by Peter Bruin

Replying to jdemeyer:

I quite like the suggestion by Karim Belabas to use nfsign() together with nfnewprec. I think that's a lot better than the current patch, which is really just a hack.

Maybe you should ask Denis Simon what he thinks (perhaps with a proof-of-concept patch)? Ideally, he fixes his script "upstream".

I completely agree with all of this and hope that we will soon hear Denis's opinion (I e-mailed him and pari-dev simultaneously). The current patch is indeed fundamentally a hack and will hopefully be superseded by a better solution.

comment:18 Changed 9 years ago by John Cremona

Sorry that I have not had time to join this discussion -- I am very interested in the outcome.

It would be good to get Denis to adapt his script. He has done this in the past several times. I have found that it does not work for someone else to revise his script and send him the result; he want to do it himself and remain sole author. But a proof-of-concept patch might help him.

comment:19 in reply to:  18 Changed 9 years ago by Jeroen Demeyer

Replying to cremona:

It would be good to get Denis to adapt his script. He has done this in the past several times. I have found that it does not work for someone else to revise his script and send him the result; he want to do it himself and remain sole author. But a proof-of-concept patch might help him.

Good to know. I also needed to make some changes to his scripts for the PARI-2.7 port (#15767).

comment:20 Changed 9 years ago by Jeroen Demeyer

Did you get any answer from Denis Simon? If not, we should try to implement the nfsign() solution ourselves.

comment:21 in reply to:  20 Changed 9 years ago by Peter Bruin

Replying to jdemeyer:

Did you get any answer from Denis Simon?

Unfortunately not (yet).

If not, we should try to implement the nfsign() solution ourselves.

If you want to do it, go ahead, but to be honest, it seems a bit pointless to me at the moment.

Although I agreed above that the increase from 10 to 38 in the precision check was just a hack, there is actually a reason for the 38. (The following assumes a 64-bit word length.) For any bound n < 38 in this test, a 1-bit error in the result already has catastrophic consequences: a number that should be negative is taken to be a positive number with 64 bits of precision, due to the precision being only 1 bit, and that bit being affected by a round-off error). For n >= 38, the round-off errors have to accumulate to at least 65 bits of precision for the problem to occur (wrongly trusting the positivity of a 128-bit number). Of course that can happen, but it is extremely unlikely and I would be surprised if you could find any example.

comment:22 Changed 9 years ago by Jeroen Demeyer

Reviewers: Jeroen Demeyer
Status: needs_reviewpositive_review

comment:23 Changed 9 years ago by Volker Braun

Branch: u/pbruin/15483-two_descent_precisionb662aba98ad4122a82fc135cf9b5dbf58de2ba0e
Resolution: fixed
Status: positive_reviewclosed
Note: See TracTickets for help on using tickets.