Opened 6 years ago

Closed 6 years ago

#15483 closed defect (fixed)

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

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

Description (last modified by pbruin)

[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 follow-up: Changed 6 years ago by 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 6 years ago by pbruin

  • 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 6 years ago by wuthrich

  • Report Upstream changed from N/A to Not 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 6 years ago by cremona

  • Description modified (diff)

comment:5 Changed 6 years ago by vbraun_spam

  • Milestone changed from sage-6.1 to sage-6.2

comment:6 Changed 6 years ago by pbruin

  • Dependencies set to #11005

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

comment:7 Changed 6 years ago by pbruin

  • 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 6 years ago by pbruin

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 6 years ago by pbruin

  • Report Upstream changed from Not yet reported upstream; Will do shortly. to 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 6 years ago by pbruin

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 6 years ago by pbruin (previous) (diff)

comment:11 Changed 6 years ago by pbruin

  • Branch set to u/pbruin/15483-two_descent_precision
  • Commit set to d34a77ab7abe4072349e10e0154a598da94b7c31
  • Status changed from new to needs_review

comment:12 Changed 6 years ago by pbruin

  • Status changed from needs_review to needs_work

This needs a doctest...

comment:13 Changed 6 years ago by git

  • Commit changed from d34a77ab7abe4072349e10e0154a598da94b7c31 to b662aba98ad4122a82fc135cf9b5dbf58de2ba0e

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

b662abaadd/fix doctests for simon_two_descent

comment:14 Changed 6 years ago by pbruin

  • Status changed from needs_work to needs_review

comment:15 Changed 6 years ago by pbruin

  • Authors set to Peter Bruin
  • Report Upstream changed from Reported upstream. No feedback yet. to Workaround found; Bug reported upstream.

comment:16 follow-up: Changed 6 years ago by 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".

comment:17 in reply to: ↑ 16 Changed 6 years ago by pbruin

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 follow-up: Changed 6 years ago by 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 6 years ago by jdemeyer

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 follow-up: Changed 6 years ago by jdemeyer

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 6 years ago by pbruin

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 6 years ago by jdemeyer

  • Reviewers set to Jeroen Demeyer
  • Status changed from needs_review to positive_review

comment:23 Changed 6 years ago by vbraun

  • Branch changed from u/pbruin/15483-two_descent_precision to b662aba98ad4122a82fc135cf9b5dbf58de2ba0e
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.