Opened 12 years ago
Last modified 5 years ago
#8829 closed enhancement
Saturation for MWgroups of elliptic curves over number fields. — at Version 36
Reported by:  robertwb  Owned by:  cremona 

Priority:  major  Milestone:  sage8.1 
Component:  elliptic curves  Keywords:  saturation 
Cc:  pbruin, kedlaya  Merged in:  
Authors:  Robert Bradshaw, John Cremona  Reviewers:  
Report Upstream:  N/A  Work issues:  
Branch:  u/cremona/8829 (Commits, GitHub, GitLab)  Commit:  8686677fa5cdf4359fd6f6b8d8e25925f6893a4c 
Dependencies:  #8828  Stopgaps: 
Description (last modified by )
Implementation of saturation of points on elliptic curves over number fields. Original patch by Robert Bradshaw in 2010, refactored and made into a git branch by John Cremona in 2015.
I also implemented the simple case of E.gens() for E(K) when E/Q and [K:Q] = 2.
Change History (37)
comment:1 Changed 12 years ago by
 Status changed from new to needs_review
comment:2 followup: ↓ 3 Changed 12 years ago by
I have had a quick look and will go through this in more detail later (after #8828 is completed, probably). I spent a long time on my C++ implementation of this (over QQ but the algorithm is general) so am quite familiar with the details.
Here are two references you should give: [1] S. Siksek "Infinite descent on elliptic curves", Rocky Mountain J of M, Vol 25 No. 4 (1995), 15011538. [2] M. Prickett, "Saturation of MordellWeil groups of elliptic curves over number fields", U of Nottingham PhD thesis (2004), http://etheses.nottingham.ac.uk/52/.
Martin Prickett implemented this in Magma, but the code was very slow and hard to read so it never got incorporated into Magma releases.
Incidentally, it was for this that I implemented group structure for curves over GF(q) in the first place! In my C++ implementation I cache a lot of the information of this group structure so that when you do psaturation for larger and larger p, the structures are already there. A good example is to take one of those curves of very high rank: I think I once successfully psaturated the rank 24 curve at all p < 10^6
(the bound was totally out of reach, something like 10^100
).
Another point which might be useful over number fields: it suffices to use degree one primes to reduce modulo.
Changed 12 years ago by
comment:3 in reply to: ↑ 2 Changed 12 years ago by
Replying to cremona:
I have had a quick look and will go through this in more detail later (after #8828 is completed, probably). I spent a long time on my C++ implementation of this (over QQ but the algorithm is general) so am quite familiar with the details.
Here are two references you should give: [1] S. Siksek "Infinite descent on elliptic curves", Rocky Mountain J of M, Vol 25 No. 4 (1995), 15011538. [2] M. Prickett, "Saturation of MordellWeil groups of elliptic curves over number fields", U of Nottingham PhD thesis (2004), http://etheses.nottingham.ac.uk/52/.
Ah, those look like good references to read too :).
Martin Prickett implemented this in Magma, but the code was very slow and hard to read so it never got incorporated into Magma releases.
Incidentally, it was for this that I implemented group structure for curves over GF(q) in the first place! In my C++ implementation I cache a lot of the information of this group structure so that when you do psaturation for larger and larger p, the structures are already there.
The way I do it is consider many p at once, and for each curve over GF(q) I see which primes in my set it could help with, though this won't scale as far. I'm sure there's still lots of room for improvement.
A good example is to take one of those curves of very high rank: I think I once successfully psaturated the rank 24 curve at all p <
10^6
(the bound was totally out of reach, something like10^100
).
That reminds meI was wondering if there's any way to go from min(h(P)) to a bound on the regulator for rank > 1.
Another point which might be useful over number fields: it suffices to use degree one primes to reduce modulo.
Good point. Those get pretty rare for large degree number fields though, right?
comment:4 Changed 12 years ago by
You might also like to look at my C++ code which is in eclib, in src/qcurves. I can point to the right files if it is not clear. In case you wonder, "TLSS" stands for "TateLichtenbaumSamir_Siksek" since I use the TL map when the ptorsion in E(GF(q)) is not cyclic and Samir's original method when it is. Samir only used reduction modulo primes where p exactly divided the order, and in particular for which the reduction had cyclic ppart. But Martin and I discovered that this can fail when there is a pisogeny. Here, fail means in the sense that there can exist points which are not multiples of p in E(QQ) but which map to zero in E(GF(q))/p for all q.
In MP's thesis he proves that this cannot happen if you use all q, or all but a finite number, or all but a finite number of degree 1 primes, .... some of these results we then found had been proved elsewhere (3 or 4 times, independently, within 3 or 4 years!). But it can happen if you leave out the q for which the quotient has noncyclic ppart.
comment:5 Changed 12 years ago by
 Reviewers set to John Cremona
 Status changed from needs_review to needs_work
Patch applies fine to 4.4.1 and tests pass.
This functionality is badly needed!
We now have heights for points on curves defined over number_fields but no associated regulator function. I suggest that the function regulator_of_points() be moved up from ell_rational_field to ell_number_field. This tcan then be called instead of the code in lines 424432 [line numbers are from the patched file, not the patch].
Line 439 uses a function self.height_function() which does not exist. This block needs to be replaced by something sensible. If one has a lower bound on the height of nontorsion points, then a bound on the index can be deduced from standard geometry of numbers estimates. To get such a lower bound, see papers of Cremona & Siksek (over Q) and Thongjunthug (over number fields) for an algorithm which would need to be implemented. (Not hard over Q, not much harder for totally real fields, quite a lot worse over fields with a complex embedding). Until this is done, I don't think this saturation function can allow maxprime==0.
In the rank one code: when large primes p (say, over 20!) are being tested then the division_points code will be expensive since it involves constructing the multiplicationbyp map. I would recommend using a reduction strategy here just as in the general case. To check psaturation just find primes q such that #E(Fq) is divisible by p and then see if the reduction of P mod q is a multiple of p. This will almost always prove saturation very quickly. If it fails for several (say 5) q then try to divide P by p; else use more q, and so on. There is one theoretical drawback here: this strategy might fail if there is a rational pisogeny. Over Q, we know which p this might happen for and I would first test for the existence of isogenies of primes degree, and use the division test (as here) for any primes that show up. Over number fields that's harder to deal with, but again we can fall back on the division test to rpove that P cannot be divided by p.
The function list_of_multiples(P,n) duplicated the generic function multiples() which I wrote for just this sort of purpose!
I don't like the loop through all linear combinations for small
primes. Even with p=2 there are curves with 24 independent points out
there and 2^24
divisions is not nice to contemplate. If you want this
short cut, do it based on the size of p^r
.
The main code with reduction etc looks fine to me (but I did not check the logic rigorously).
The gens function for E(K) when E is defined over Q and [K:Q]=2 looks fine. For a more general case we could always try using simon_two_descent (followed by saturation). Trying such an examples led me to:
sage: K.<a> = NumberField(x^22) sage: E = EllipticCurve([a,0]) sage: P = E(0,0) sage: P.has_finite_order() True sage: P.order() 2 sage: P.height() 0 sage: E.saturation([P], verbose=True, max_prime=5) ## infinite loop
This is caused as follows: The height matrix is [0] with det=0 but reg / min(heights) is NaN so reg / min(heights) < 1e6 is False!. This will need fixing. At the very least, I would discard any points of finite order before doing anything else with them. Then min(heights) will never be 0.
Most of the above is easy to deal with. The hard part is computing a suitable max_prime form a lower height bound on points. I suggest that for now you make it compulsory to have a positive max_prime and add a TODO.
comment:6 followup: ↓ 7 Changed 12 years ago by
Thank you for all your input. self.height_function
comes from #8828, though as you suggest we could make max_prime mandatory for now (and for rank > 1 once that goes in). That's a good point about large primes in the rank one case. I found the loop through all linear combinations to be much faster in practice for small primes, but the hard coded p == 2
case was left by accident, I meant to cap that on p^r
as I did the others.
I probably won't fix and polish this up before finishing my thesis, but at the latest we should be able to get it done during the workshop at MSRI next month.
comment:7 in reply to: ↑ 6 Changed 12 years ago by
Replying to robertwb:
Thank you for all your input.
self.height_function
comes from #8828, though as you suggest we could make max_prime mandatory for now (and for rank > 1 once that goes in). That's a good point about large primes in the rank one case. I found the loop through all linear combinations to be much faster in practice for small primes, but the hard codedp == 2
case was left by accident, I meant to cap that onp^r
as I did the others.I probably won't fix and polish this up before finishing my thesis, but at the latest we should be able to get it done during the workshop at MSRI next month.
OK  looking forward to it! I'll take a look at #8828 by then as well.
comment:8 Changed 12 years ago by
Since rwb is now busy at Google, and I want this functionality, I am now implementing the changes I suggested above!
comment:9 Changed 12 years ago by
I made a separate ticket for the regulator functions. See #9372.
comment:10 Changed 10 years ago by
How is this going John? It seems to have been awhile....
comment:11 followup: ↓ 12 Changed 9 years ago by
See #12509: until we can fix the height computation, saturation cannot be carried out properly. It's still on my todo list though.
comment:12 in reply to: ↑ 11 Changed 9 years ago by
comment:13 Changed 9 years ago by
 Milestone changed from sage5.11 to sage5.12
comment:14 Changed 8 years ago by
 Dependencies set to #8828
comment:15 Changed 8 years ago by
 Summary changed from Saturation for curves over number fields. to Saturation for MWgroups of elliptic curves over number fields.
comment:16 Changed 8 years ago by
 Milestone changed from sage6.1 to sage6.2
comment:17 Changed 8 years ago by
 Cc pbruin added
comment:18 Changed 8 years ago by
 Milestone changed from sage6.2 to sage6.3
comment:19 Changed 8 years ago by
 Milestone changed from sage6.3 to sage6.4
comment:20 Changed 7 years ago by
I do not know why this was left drifting, but I really need it myself now so will look at it again, rebase on 6.8 and see what we can do. But I only have one day before a week off, so...
comment:21 Changed 7 years ago by
 Branch set to u/cremona/8829
 Commit set to 0e1e35f624edb087d3fb1ddc21226fec7acfafad
 Description modified (diff)
 Keywords saturation added
comment:22 Changed 7 years ago by
Current branch works but more doctests and testing are needed; so not ready for review yet.
I did a lot of rewriting of the main saturation routine, separating off psaturation and also allowing saturation to be done at just one prime. This is a useful special case, since if you take the images of some saturated points under a pisogeny the images may not be psaturated but will still be saturated at all other primes.
The code for computing E(K) when K is quadratic and E is a base change has been completely rewritten and will now work in many more cases (not just when the coefficients of E are in Q).
comment:23 Changed 7 years ago by
 Commit changed from 0e1e35f624edb087d3fb1ddc21226fec7acfafad to a85f40629df89da829e3a38d9e0443576bced01d
Branch pushed to git repo; I updated commit sha1. New commits:
a85f406  #8829: finished doctests for ec/nf saturation

comment:24 Changed 7 years ago by
The latest commit involves more than adding more doctests to the new functions, since bugs were revealed which led to a rewrite of the sieving code for the two cases where the prank of the reduction is 1 or 2; the former uses discrete log in the reduction, the latter uses Weil pairing and discrte log in the multiplicative group. In the sieve I restrict to primes of degree 1. It is a Theorem (see http://eprints.nottingham.ac.uk/10052/) that this will suffice to prove psaturation, provided that one does use reductions with prank 2 and not just those of prank 1 as originally suggested by Siksek in https://ore.exeter.ac.uk/repository/handle/10871/8323 .
I will mark this as ready for review so the bots get to work on it, and of course humans are welcome to look at the new code, but I will now test it thoroughly on the LMFDB curves and report back.
comment:25 Changed 7 years ago by
 Reviewers John Cremona deleted
 Status changed from needs_work to needs_review
comment:26 followup: ↓ 27 Changed 7 years ago by
Hello,
1) indent correctly the INPUT and OUTPUT fields:
INPUT:  first thing goes one there (note the shift of 2 characters)
2) use the new syntax for raise:
raise MyError("is rich")
Point 1 may be the source of the doc build failure found by the bot:
OSError: [plane_cur] /home/patchbot/sagepatchbot/local/lib/python2.7/sitepackages/sage/schemes/elliptic_curves/ell_number_field.py:docstring of sage.schemes.elliptic_curves.ell_number_field.EllipticCurve_number_field.saturation:9: WARNING: Bullet list ends without a blank line; unexpected unindent.
comment:27 in reply to: ↑ 26 Changed 7 years ago by
Replying to chapoton:
Hello,
1) indent correctly the INPUT and OUTPUT fields:
INPUT:  first thing goes one there (note the shift of 2 characters)2) use the new syntax for raise:
raise MyError("is rich")Point 1 may be the source of the doc build failure found by the bot:
OSError: [plane_cur] /home/patchbot/sagepatchbot/local/lib/python2.7/sitepackages/sage/schemes/elliptic_curves/ell_number_field.py:docstring of sage.schemes.elliptic_curves.ell_number_field.EllipticCurve_number_field.saturation:9: WARNING: Bullet list ends without a blank line; unexpected unindent.
Thanks, I will fix those.
comment:28 Changed 7 years ago by
 Commit changed from a85f40629df89da829e3a38d9e0443576bced01d to 4f511b89fd199d070bad0d1054efa148d765fdf6
Branch pushed to git repo; I updated commit sha1. New commits:
4f511b8  fixup doc errors and raise() syntax

comment:29 Changed 7 years ago by
 Status changed from needs_review to needs_work
many failing doctests, see bot report
comment:30 Changed 7 years ago by
Apologies, it was a mistake to set this to needs_review prematurely. Next time I do, I will mean it.
comment:31 Changed 7 years ago by
 Commit changed from 4f511b89fd199d070bad0d1054efa148d765fdf6 to 2f20f7373a13e2267bc5ae3d2a7cae93278a76c4
comment:32 Changed 7 years ago by
Progress report: I am currently running the psaturation (for single primes) on lots of LMFDB curves and all is well so far. This is almost always for very small p (mainly 2 and 3) though, since I am starting with some saturated points on one curve (provided by Magma) and using pisogenies to map to other curves in the isogeny class. Higher degree isogenies are not so common.
I did start to veryify that the points from Magma were fully saturated, but ran into problems computing the saturation index, using (line 3717) the lower bound on the height of all nontorsion points  previously implemented and merged i n6.3 (see #8828). For example, I had a curve where the value of 5 in that line was insufficient *and led to an infinite loop in the call to min()*, while 10 worked fine, but now I have a curve where I have not yet found a value which gives anything. For the record I will give that example here:
K.<phi> = NumberField(x^2x1) # Q(sqrt(5)) E = EllipticCurve([phi + 1, phi + 1, 1, 20*phi  39, 196*phi + 237]) H = E.height_function() H.min(.1,10,verbose=True) # does not appear to terminate
Strictly this is about the code merged in #8828, but it will need fixing here before we can let this (useful!) function out into the world.
comment:33 Changed 7 years ago by
 Commit changed from 2f20f7373a13e2267bc5ae3d2a7cae93278a76c4 to 8686677fa5cdf4359fd6f6b8d8e25925f6893a4c
Branch pushed to git repo; I updated commit sha1. New commits:
7dbb186  Merge branch 'develop' of trac.sagemath.org:sage into isogs

c66d49f  Merge branch 'develop' of trac.sagemath.org:sage into isogs

ce2472b  Merge branch 'isogs' into sat+isogs

8c80b6f  Merge branch 'saturate' into sat+isogs

8686677  8829: two bug fixes, more tests

comment:34 Changed 7 years ago by
The latest branch I just pushed has some merges in it which were not intended but I hope that will not cause any problems  as well as merging 6.9.beta6 I also merged by branch 'isogs' which has been merged into beta7.
One bug fix addresses the previous comment  after rereading my own 2006 paper I found that the original implementer from #8828 had missed one point (when mu is halved one must increment n_max in order to guarantee termination). A small additional improvement in the same place (the method min_gr() in height.py) now gives a small improvement in the bound, which is why one doctest there has been changed.
The second bug was to do with mutability of lists giving unwanted side effects, and is commented at the point in the source which has changed.
It is likely that users who call the saturation() method will also want to lll_reduce() the output but I have not made that automatic.
I will set this to needs_review once my own full test has completed.
comment:35 Changed 7 years ago by
 Status changed from needs_work to needs_review
comment:36 Changed 7 years ago by
 Description modified (diff)
 Status changed from needs_review to needs_work
After further testing (on many thousands of curves but only psaturating for small p) I saw that it was bad to use discrete_log_lambda() for the dlog in the multiplcative group (in the rarer case where the prank of the reduction is 2 and the Weil pairing is used), both unnecessary and less efficient than simply w.log(zeta). One additional commit coming up...
Some dependance on #8828.