Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update the rank function of elliptic curves to use ellrank in pari #35626

Merged
merged 19 commits into from
Jun 11, 2023

Conversation

chriswuthrich
Copy link
Contributor

@chriswuthrich chriswuthrich commented May 8, 2023

📚 Description

Fixes #35621

Pari has now a very good in-built algorithm to determine the rank and a set of generators of the Mordell-Weil group of elliptic curves over the rational. It is, as far as I understand, a port of Denis Simon's script ellQ.gp. It performs a 2-descent and uses the Cassels Tate pairing to conclude the rank in many cases - when Sha has no 4-torsion elements.
The functionality to search rational points is also very quick.

📝 Checklist

  • [x] The title is concise, informative, and self-explanatory.
  • [x] The description explains in detail what this PR is about.
  • [x] I have linked a relevant issue or discussion.
  • [x] I have created tests covering the changes.
  • [x] I have updated the documentation accordingly.

⌛ Dependencies

@chriswuthrich chriswuthrich marked this pull request as ready for review May 9, 2023 07:36
@chriswuthrich
Copy link
Contributor Author

chriswuthrich commented May 9, 2023

This is ready for review.

I implemented that .rank() and .prove_BSD() still take mwrank by default, but they have the option to take algorithm=pari now. Instead .selmer_bound() and rank_bound() now take pari by default, but allow mwrank as an option. For rank_bound this is justified as the implementation in pari is better in the sense that it also uses the Cassels-Tate pairing and hence determines a smaller upper bound in some cases.
I made simon_two_descent deprecated. Only for elliptic curves over Q, it is still available over number fields.

Maybe it is better to take pari as the default in all these, but I do not have enough experience in comparing the two. Obvioulsy @JohnCremona is asked for comment on this here :)

However, this question may be decided later with a small change. My choice here already improves sage and it could be the optimal choice...

@JohnCremona JohnCremona self-requested a review May 9, 2023 09:45
@JohnCremona JohnCremona self-assigned this May 9, 2023
@JohnCremona
Copy link
Member

I will review this. When Bill Allombert was implementing all this in pari he corresponded frequently with me, and I know that he did a really excellent job. As a result of that, mwrank was improved a bit too (which Sage has alreay benefited from), for example in saturation of points.

Copy link
Member

@JohnCremona JohnCremona left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have made a few suggestions. By the way, the code changes here are on top of those in your other PR, so I have not re-reviewed those parts -- but it would have been a bit better to have made a fresh branch for this.

Anyway -- thanks a million for doing this! I have been hoping that someone would ever since pari's ellrank has been available. Also, I would be quite happy for 'pari' to replace 'mwrank' as default.


- ``E`` -- an elliptic curve

OUTPUT:
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Say that the output is a tuple of 5 elements, ...

# dim Sha[2] = dim Sel2 - rank E(Q) - dim tors
# rank_upper_bd = dim Sel_2 - dim tors - s
sha_upper_bd = rank_upper_bd - len(gens) + s
return len(gens), rank_upper_bd, s, sha_upper_bd, gens
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why return len(gens) instead of pari's 'lower', which you don't use? Can the lower rank bound be greater than the number of gens? But if pari increases the lower rank bound by 1 conditional on parity then I approve of using len(gens) instead.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That is my understanding of pari, although their documentation does not mention that the lower bound is conditional on the finiteness of the 2-primary part of the Tate-Shafarevich group.

In line 1953 of file ellrank.c the function ell2selmer, we have   if (odd(dim - nbpoints)) mwrank++;

I was planning to write to the pari-gp list asking about this (and whether they have plans to port the code for the number fields, too.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK -- let's just put a comment in that we don't use pari's lower bound as that is sometimes incremented by 1 based on a standard conjecture.

This function is deprecated as the functionality of
Simon's script for elliptic curves over the rationals
has been ported over to pari.
Use :meth:`.rank` with the keyword ``algorithm='pari;`` instead.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

sytax of this line looks dodgy, and I think the semicolon should be another single quote.

proof=None):
"""
proof=None,
pari_effort=0):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is 0 really the best default?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good question. It is the default in pari. I am amazed how many points are already found with effort 0. If I understand correctly, the algorithm calculates the 2-Selmer group first, then searches until either the effort bound limits it or enough points are found. In this way having a larger default value wouldn't harm the speed of the easy examples. However, maybe it is best to leave it to the user to do it with minimum effort first and then to increase. The way I implemented it the points found already in the first run are passed over to pari's second run.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good answer! Leave it as 0.

@@ -2038,14 +2061,16 @@ def rank(self, use_database=True, verbose=False,
sage: E.minimal_model().rank()
1

A large example where mwrank doesn't determine the result with certainty::
A large example where mwrank doesn't determine the result with certainty, but pari does::
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nice!

@@ -2218,13 +2283,18 @@ def gens(self, proof=None, **kwds):
:meth:`~gens_certain` method to find out afterwards
whether the generators were proved.

IMPLEMENTATION: Uses Cremona's mwrank C library.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's actually C++, but this has been here for years...

RuntimeError: generators could not be determined. So far we found []. Hint: increase pari_effort.
sage: E.gens(use_database=False, algorithm="pari",pari_effort=10)
[(-166136231668185267540804/2825630694251145858025 : 167661624456834335404812111469782006/150201095200135518108761470235125 : 1)]

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I like the last two examples -- are they fast enough or should they be tagged 'long time'?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The example -127^2 and effort 4 is done in 100ms on my old laptop. The last -157^2 with effort 10 takes seconds, so I put long time to that.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK

@@ -2333,7 +2426,26 @@ def _compute_gens(self, proof,
except RuntimeError:
pass
# end if (not_use_mwrank)
if algorithm == "mwrank_lib":
if algorithm == "pari":
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This block repeats lines 2210-2227 rather closely. Would it be possible to have a helper function used by both? The helper would not itself raise a runtime error, but that condition would be caught by the calling function.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do you insist? I am too lazy to put it into one for just two cases. I agree it might be better.

ep = self.pari_curve()
lower, upper, s, pts = ep.ellrank()
T = self.torsion_subgroup().invariants()
tor = sum(x%2==0 for x in T)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You could use self.two_torsion_rank(), provided just for this!

@@ -1,7 +1,7 @@
r"""
Îlu algorithm for elliptic-curve isogenies
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think that the funny square-root character in place of V is intentional, as that is what @yyyyx4 and friends call the improved Velu algorithm 'velusqrt'. So perhaps revert all these?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh, I didn't know. I saw all of these (find them ugly and possibly insulting to Monsieur Vélu qui n'est pas un élu). I am personally against it in headers, as it just looks like a bug. If they want to use it as a curiosity in some of the documentation I am not against, it, but I would still want most of them to be reverted. Sorry to be spoil fun.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@yyyyx4 How strongly do you feel about this? Every time I see it I think there is something wrong with my computer screen!

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Well, √élu is the name given to the algorithm by its authors (not me). Other names that have appeared in the literature are all variants of "square‑root Vélu", with varying choices of punctuation, capitalization, accentuation, abbreviation, and word order. I obviously prefer √élu, but most of the other variants are acceptable to me. Identifying this algorithm with Vélu's seems wrong, though.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for the comments. (Due to my bad choice of basing this issue on top of the other, the change of Vélu actually belongs to the dependency and not to this pull request. I will correct it there.)

I will of course correct this and I am sorry that I did this change without looking closer.
I will use "square-root Vélu". If there are other who prefer keeping it in the documentation and in the _repr_ then they should comment here (or better on the other pull request). I don't have a strong feeling and I am ok to revert to √élu.

@chriswuthrich
Copy link
Contributor Author

I still have not decided which algorithm to take as default. I think they are similar in speed but pari determines the rank more often.

Here some very naive timing on my old laptop:

curves #curves mwrank_lib pari
1<N<=1000 5113 2.9 s 2.94 s
8888<=N<=10000 7444 4.46 s 4.49 s
11111<=N<=12345 8201 5.4 s 4.85 s
111111<=N<=112345 7609 5.02 s 4.51 s
[1,B] with 1<=B<=1000 1000 40.2 s but failed in 31 cases 2min 17s
[3,B] with 1<=B<=1000 1000 40.7 s but failed in 35 cases 2min 21s
[A,1] with 1<=A<=1000 1000 59.7 s but failed in 67 cases 3min 23s but failed in 4 cases

Based on this I would suggest to change for pari as default. But I welcome any other opinion on this.

@JohnCremona
Copy link
Member

I think that I have re-reviewed and checked that the opints I raised have all been resolved, but I am not certain. GitHub's reviewing system is very sophisticated but also quite confusing. I did not see how it was possible for me as a reviewer to mark each point I had raised as "resolved".

About the default, let's leave it as it is for now. Changing it later would be easy.

@chriswuthrich
Copy link
Contributor Author

I agree with leaving the default as it is. I am not confident at all about my rough timings above. I will check with pari if they have experience in that respect.

Many thanks, John, for the helpful comments and the detailed review of the two tickets. Now that I am getting more used to this interface here, I am contactable for any reviews myself.

@vbraun
Copy link
Member

vbraun commented May 24, 2023

Merge conflict

@JohnCremona
Copy link
Member

Same as #35614 -- I'll give it +ve review when tests have passed

@chriswuthrich
Copy link
Contributor Author

Probably not the canonical solution, but I updated the dependant branch by resolving the merge conflicts with develop there, then I rebased this branch on the new version of the dependent branch. (Sorry if that is not good, but I think the history is cleaner this way.) [[ In the future I will not pile pull requests on top of each other, promised...]]

@github-actions
Copy link

github-actions bot commented Jun 4, 2023

Documentation preview for this PR (built with commit 0f8b28a) is ready! 🎉

@vbraun vbraun merged commit 86d2be0 into sagemath:develop Jun 11, 2023
@mkoeppe mkoeppe added this to the sage-10.1 milestone Jun 11, 2023
vbraun pushed a commit that referenced this pull request Jun 21, 2023
#35626 introduced a test failure on some machines as some tests return
random choices of generators of Mordell-Weil groups for elliptic curves.

This adds "random" to these tests and add some tests that do test if the
returned answer is correct. This is only possible when the rank is 1
(and the returned point must be in a finite set). For larger ranks at
best we can test if the number of generators is ok and if the points are
saturated.

URL: #35787
Reported by: Chris Wuthrich
Reviewer(s):
@chriswuthrich chriswuthrich deleted the pariell branch July 4, 2023 12:05
Comment on lines +2301 to +2303
sage: E = EllipticCurve([0,2429469980725060,0,275130703388172136833647756388,0])
sage: len(E.gens(algorithm="pari"))
14
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This test takes 50s, and brings testing the whole file (around 900 tests) from 22s in sagemath-10.0 to 73s in sagemath-10.1.

This doesn't seem ok even for --long.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This looks like an example where the old default (mwrank) is better than the newly implemented pari code.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Indeed, mwrank takes 13s while pari takes 51s. The default seems to be still mwrank.

In any case, I would suggest that this test be marked "# not tested, too long" or something (or be replaced by a smaller rank example which takes <1s).

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fully agree. It should stay in the documentation as a complicated example, but not be part of the tests. (Just back from holidays and fairly busy for a while) Please feel free to change this, I will do it but it will take a while.

vbraun pushed a commit to vbraun/sage that referenced this pull request Sep 10, 2023
…too long

    
A doctest introduced in sagemath#35626 takes ~50s on current hardware (for
comparision, doctesting the whole of
`src/sage/schemes/elliptic_curves/ell_rational_field.py` takes ~22s for
889 tests).

This is too long even for `--long`, but the example is still
interesting, so we just mark it `# not tested`.

### 📝 Checklist

- [x] The title is concise, informative, and self-explanatory.
- [x] The description explains in detail what this PR is about.
- [x] I have linked a relevant issue or discussion.


CC: @chriswuthrich (should be a trivial review 🙏)
    
URL: sagemath#36195
Reported by: Gonzalo Tornaría
Reviewer(s):
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Update the rank function of elliptic curves to use pari
8 participants