Add hyperelliptic curves using the smooth model

[GiacomoPope]

Overview

This PR includes a whole new module into schemes which implements hyperelliptic curves over the smooth model, and focuses on implementing decent arithmetic for Jac(H) for all reasonable cases, something which is not done properly in the current implementation due to limitations of the curve model used.

The idea is that hyperelliptic_curves_sm will be available with the old model for some time, until someone decides we should depreciate the old model in favour of this new implementation.

Motivation

The current implementation of hyperelliptic curves in SageMath uses the projective plane model. Although this works nicely enough for imaginary curves with only one point at infinity, it is not descriptive enough for the real models. My gut feeling is the projective plane model was used as in early Sage days this was easier to do and allowed hyperelliptic curves to be supported at all. Mathematically, I believe it makes much more sense to use a weighted projective model to have this new description but it requires more tools which we only have thanks to the work of others in sage since it was released.

This PR is a total rewrite of the hyperelliptic curve classes to instead use the smooth model for hyperelliptic curves which facilitates implementing arithmetic of Jacobians of hyperelliptic curves for (almost) all cases. In particular, now if one tries to perform arithmetic on Jac(C) one either gets the correct value or a well-handled error instead of just returning something which is wrong.

The hope is that with this new model for the curves, more recent research in the area of hyperelliptic curves and their jacobians can be more easily integrated. As a personal example, being able to compute arithmetic in Jac(H) for the real model unconditionally would help with the implementation of genus two isogenies which are "in vogue" right now in the cryptography world.

Content of PR

This PR looks WAY bigger than it is, as it duplicates ALL code from hyperelliptic_curves and then has modifications to allow this to work in the new curve model. I believe the best thing to do in the long run is to depreciate the old hyperelliptic curve impl and have this replace it, but this will take years(?) so I think we'll just have to have both side by side for a while.

This PR still needs a lot of work including:

• more tests
• better comments and docstrings
• potentially refactoring of which files belong where

but the code has sat at this stage for a long time without much more progress from the three of us, so I think the best thing to do is open up the PR and try and get some extra attention on getting this code ready to be included.

One benefit is that all this code is "new" in that it should not conflict with anyone else's work, so we can get this code into a good position and keep adding to it with more interesting maths once the base layer is in place

Example of better arithmetic

For example, this code fixed the following issue:

#32024

sage: R.<x> = QQ[]
sage: f = 144*x^6 - 240*x^5 + 148*x^4 + 16*x^3 - 16*x^2 - 4*x + 1
sage: H = HyperellipticCurveSmoothModel(f)
sage: J = Jacobian(H)
sage: P = J(H(0,1))-J(H(0,-1))
sage: (5*P).is_zero()
False
sage: 
sage: H = HyperellipticCurve(f)
sage: J = Jacobian(H)
sage: P = J(H(0,1))-J(H(0,-1))
sage: (5*P).is_zero()
True

and is also related to the issues #37109 #37093 #37101 #37626

Help!

This code needs more work, and lots of time spent on the review which I don't think is a reasonable thing to do for one person. So if this area interests you any help would be appreciated (either in the review of some of the code or some extra commits which tidy up aspects)

[github-actions]

Documentation preview for this PR (built with commit 130cbd2; changes) is ready! 🎉
This preview will update shortly after each push to this PR.

[GiacomoPope]

This branch:

sage: R.<x> = GF(7)[]
....: H = HyperellipticCurve(3*x^6 + 2*x^2 + 1)
....: J = Jacobian(H)
....: JK = J(GF(7))
....: e = JK.an_element()
....: isinstance(e, JK.category().element_class)
False
sage: type(e)
<class 'sage.schemes.hyperelliptic_curves.jacobian_morphism.MumfordDivisorClassFieldInert'>
sage: JK.category().element_class
<class 'sage.categories.homsets.HomsetsOf.element_class'>

Main:

sage: R.<x> = GF(7)[]
....: H = HyperellipticCurve(3*x^6 + 2*x^2 + 1)
....: J = Jacobian(H)
....: JK = J(GF(7))
....: e = JK.an_element()
....: isinstance(e, JK.category().element_class)
False
sage: type(e)
<class 'sage.schemes.hyperelliptic_curves.jacobian_morphism.JacobianMorphism_divisor_class_field'>
sage: JK.category().element_class
<class 'sage.categories.homsets.HomsetsOf.element_class'>

[grhkm21]

I committed a few small changes. They address the following bugs:

#!/usr/bin/env sage
"""
Test script to demonstrate fixes for PR #39161

Issues fixed:
1. Pickle/unpickle equality (loads(dumps(P)) == P should be True)
2. J.is_parent_of(J.zero()) should return True
3. lift_u(0) should not raise NameError on split model
"""

def test_pickle_ramified():
    """Test pickle/unpickle on ramified model"""
    x = polygen(GF(5))
    H = HyperellipticCurve(x^5 + 3*x + 1)
    J = H.jacobian()
    P = J.random_element()
    result = loads(dumps(P)) == P
    print(f"[{'PASS' if result else 'FAIL'}] Pickle ramified: loads(dumps(P)) == P -> {result}")
    return result

def test_pickle_split():
    """Test pickle/unpickle on split model"""
    x = polygen(GF(5))
    H = HyperellipticCurve(x^6 + x + 1)
    J = H.jacobian()
    P = J.random_element()
    result = loads(dumps(P)) == P
    print(f"[{'PASS' if result else 'FAIL'}] Pickle split: loads(dumps(P)) == P -> {result}")
    return result

def test_pickle_inert():
    """Test pickle/unpickle on inert model"""
    x = polygen(GF(5))
    H = HyperellipticCurve(2*x^6 + 1)
    J = H.jacobian()
    P = J.random_element()
    result = loads(dumps(P)) == P
    print(f"[{'PASS' if result else 'FAIL'}] Pickle inert: loads(dumps(P)) == P -> {result}")
    return result

def test_is_parent_of():
    """Test J.is_parent_of(J.zero())"""
    x = polygen(GF(19))
    H = HyperellipticCurve(x^5 + x, x^2 + 1)
    J = H.jacobian()
    result = J.is_parent_of(J.zero())
    print(f"[{'PASS' if result else 'FAIL'}] is_parent_of: J.is_parent_of(J.zero()) -> {result}")
    return result

def test_lift_u_zero():
    """Test lift_u(0) on split model - should not raise NameError"""
    x = polygen(GF(7))
    H = HyperellipticCurve(x^6 + x + 1)
    J = H.jacobian()
    try:
        result = J.lift_u(0)
        print(f"[PASS] lift_u(0): returned {result}")
        return True
    except NameError as e:
        print(f"[FAIL] lift_u(0): NameError - {e}")
        return False
    except Exception as e:
        print(f"[FAIL] lift_u(0): {type(e).__name__} - {e}")
        return False

def run_all_tests():
    """Run all tests and return overall result"""
    print("=" * 60)
    print("Testing hyperelliptic Jacobian fixes")
    print("=" * 60)

    results = []
    results.append(test_pickle_ramified())
    results.append(test_pickle_split())
    results.append(test_pickle_inert())
    results.append(test_is_parent_of())
    results.append(test_lift_u_zero())

    print("=" * 60)
    passed = sum(results)
    total = len(results)
    print(f"Results: {passed}/{total} tests passed")
    print("=" * 60)

    return all(results)

success = run_all_tests()
import sys
sys.exit(0 if success else 1)

This gave 0/5 tests passed before the commit, and 5/5 now

[vincentmacri]

I've opened GiacomoPope#2 with some minor fixes. Otherwise, code structure looks good. All that remains I think is for me to review the correctness of the balanced divisors implementation. Is that pretty much all in jacobian_morphism?

Also, don't forget to move the citation to the bibliography (I pointed this out above).

[vincentmacri]

The diffs are too long to easily read in my browser. Is it alright if I submit the rest of my suggestions via a PR to your branch?

[vincentmacri]

I've opened another PR with suggestions and a question, but otherwise I think this is almost ready to merge.

Just a few last things:

1. There are several TODO notes in the code. Did you want to discuss any of them now and address them in this PR, or do you want to leave that for future work?
2. There is a merge conflict that needs to be resolved.
3. After everything else is done, since you've basically refactored all the hyperelliptic code, would you mind running our stricter linting autofix settings on it to save us from future code cleanup? I can do another PR if you prefer, but then my name will be all over the git blame on stuff that is mostly your work. Just run these two commands:
   i. To run our stricter linting rules and autofix violations: ruff check --preview --fix src/sage/schemes/hyperelliptic_curves
   ii. To fix some rules that are probably going to be added to our stricter list of linting rules soon anyway: ruff check --preview --fix --select RET src/sage/schemes/hyperelliptic_curves

[GiacomoPope]

0. I merged your PR
1. I think the TODO can be left for later
2. I have sorted this
   3.a. I get failures, mainly from imports from lazy_import, however there are two failures which we could fix, one I did myself but the other is in a scary function and I don't know if something will break if I change it

PLW0642 Reassigned `self` variable in instance method
   --> src/sage/schemes/hyperelliptic_curves/hyperelliptic_padic_field.py:734:17
    |
732 |                 A = PolynomialRing(RationalField(), "x")
733 |                 f = A(self.hyperelliptic_polynomials()[0])
734 |                 self = HyperellipticCurve(f).change_ring(K)
    |                 ^^^^
735 |                 xP = P[0]
736 |                 xPv = xP.valuation()
    |
help: Consider using a different variable name

3.b ruff check --preview --fix --select RET src/sage/schemes/hyperelliptic_curves I fixed all errors from this.

[vincentmacri]

Alright, as long as you fixed as much as you could that's fine. The linter is pretty reliable at not breaking things, but circular import issues are one thing it struggles with.

Positive review from me then! Thank you for all this work, and for bearing through this very long code review! I think this is the second-longest code review we've had on the Sage GitHub by number of comments!

I do have one last question, just so there is something written down for future reference:

For the inert case, what did you reference when writing the implementation? Most of the papers I've seen on the subject ignore the inert case.

[GiacomoPope]

I honestly can't remember where I learnt about the inert case. Is it not in the morales paper or Galbraith book? Maybe I figured it out with help from Magma and Sabrina?

[GiacomoPope]

Thank you so much for taking on the gargantuan task of reviewing this PR. I know it was a huge amount of work but I hope in the long run it really helps with Sage

[vincentmacri]

I honestly can't remember where I learnt about the inert case. Is it not in the morales paper or Galbraith book? Maybe I figured it out with help from Magma and Sabrina?

Galbraith's book discusses it a bit but in much less detail than the other cases. Anyway, that's fine! Just wanted to make sure I wasn't missing some resource.

[yyyyx4]

Big thanks to all those who contributed to this monumental effort!

[GiacomoPope]

Yeah, thanks everyone! Really happy about this PR and excited to build on it