gh-40535: Implement construction of hyperelliptic curves from defining equation
    
<!-- ^ Please provide a concise and informative title. -->
<!-- ^ Don't put issue numbers in the title, do this in the PR
description below. -->
<!-- ^ For example, instead of "Fixes #12345" use "Introduce new method
to calculate 1 + 2". -->
<!-- v Describe your changes below in detail. -->
<!-- v Why is this change required? What problem does it solve? -->
<!-- v If this PR resolves an open issue, please link to it here. For
example, "Fixes #12345". -->



### 📝 Checklist

<!-- Put an `x` in all the boxes that apply. -->

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

### ⌛ Dependencies

<!-- List all open PRs that this PR logically depends on. For example,
-->
<!-- - #12345: short description why this is a dependency -->
<!-- - #34567: ... -->
    
URL: #40535
Reported by: user202729
Reviewer(s): Travis Scrimshaw

Author: Release Manager

src/sage/schemes/elliptic_curves/constructor.py

@@ -161,7 +161,7 @@ class EllipticCurveFactory(UniqueFactory):
 
         sage: R.<x,y> = GF(5)[]
         sage: EllipticCurve(x^3 + x^2 + 2 - y^2 - y*x)
-        Elliptic Curve defined by y^2 + x*y  = x^3 + x^2 + 2 over Finite Field of size 5
+        Elliptic Curve defined by y^2 + x*y = x^3 + x^2 + 2 over Finite Field of size 5
 
     We can also create elliptic curves by giving a smooth plane cubic with a rational point::
 
@@ -579,41 +579,19 @@ def coefficients_from_Weierstrass_polynomial(f):
         sage: R.<w,z> = QQ[]
         sage: coefficients_from_Weierstrass_polynomial(-w^2 + z^3 + 1)
         [0, 0, 0, 0, 1]
+        sage: R.<u,v> = GF(13)[]
+        sage: EllipticCurve(u^2 + 2*v*u + 3*u - (v^3 + 4*v^2 + 5*v + 6))  # indirect doctest
+        Elliptic Curve defined by y^2 + 2*x*y + 3*y = x^3 + 4*x^2 + 5*x + 6 over Finite Field of size 13
     """
-    R = f.parent()
-    cubic_variables = [ x for x in R.gens() if f.degree(x) == 3 ]
-    quadratic_variables = [ y for y in R.gens() if f.degree(y) == 2 ]
-    try:
-        x = cubic_variables[0]
-        y = quadratic_variables[0]
-    except IndexError:
+    from sage.schemes.hyperelliptic_curves.constructor import _parse_multivariate_defining_equation
+    f, h = _parse_multivariate_defining_equation(f)
+    # OUTPUT: tuple (f, h), each of them given as a list of coefficients.
+    if len(f) != 4 or len(h) > 2:
         raise ValueError('polynomial is not in long Weierstrass form')
-
-    a1 = a2 = a3 = a4 = a6 = 0
-    x3 = y2 = None
-    for coeff, mon in f:
-        if mon == x**3:
-            x3 = coeff
-        elif mon == x**2:
-            a2 = coeff
-        elif mon == x:
-            a4 = coeff
-        elif mon == 1:
-            a6 = coeff
-        elif mon == y**2:
-            y2 = -coeff
-        elif mon == x*y:
-            a1 = -coeff
-        elif mon == y:
-            a3 = -coeff
-        else:
-            raise ValueError('polynomial is not in long Weierstrass form')
-
-    if x3 != y2:
+    if not f[3].is_one():
         raise ValueError('the coefficient of x^3 and -y^2 must be the same')
-    elif x3 != 1:
-        a1, a2, a3, a4, a6 = a1/x3, a2/x3, a3/x3, a4/x3, a6/x3
-    return [a1, a2, a3, a4, a6]
+    h += [0] * (2 - len(h))
+    return [h[1], f[2], h[0], f[1], f[0]]
 
 
 def EllipticCurve_from_c4c6(c4, c6):
@@ -625,7 +603,7 @@ def EllipticCurve_from_c4c6(c4, c6):
 
         sage: E = EllipticCurve_from_c4c6(17, -2005)
         sage: E
-        Elliptic Curve defined by y^2  = x^3 - 17/48*x + 2005/864 over Rational Field
+        Elliptic Curve defined by y^2 = x^3 - 17/48*x + 2005/864 over Rational Field
         sage: E.c_invariants()
         (17, -2005)
     """

src/sage/schemes/hyperelliptic_curves/constructor.py

@@ -29,7 +29,69 @@
 from sage.structure.dynamic_class import dynamic_class
 
 
-def HyperellipticCurve(f, h=0, names=None, PP=None, check_squarefree=True):
+def _parse_multivariate_defining_equation(g):
+    """
+    Parse a defining equation for a hyperelliptic curve.
+    The input `g` should have the form `g(x, y) = y^2 + h(x) y - f(x)`,
+    or a constant multiple of that.
+
+    OUTPUT: tuple (f, h), each of them given as a list of coefficients.
+
+    TESTS::
+
+        sage: from sage.schemes.hyperelliptic_curves.constructor import _parse_multivariate_defining_equation
+        sage: R.<x,y> = QQ[]
+        sage: _parse_multivariate_defining_equation(y^2 + 3*x^2*y - (x^5 + x + 1))
+        ([1, 1, 0, 0, 0, 1], [0, 0, 3])
+        sage: _parse_multivariate_defining_equation(2*y^2 + 3*x^2*y - (x^5 + x + 1))
+        ([1/2, 1/2, 0, 0, 0, 1/2], [0, 0, 3/2])
+
+    The variable names are arbitrary::
+
+        sage: S.<z,t> = GF(13)[]
+        sage: _parse_multivariate_defining_equation(2*t^2 + 3*z^2*t - (z^5 + z + 1))
+        ([7, 7, 0, 0, 0, 7], [0, 0, 8])
+    """
+    from sage.rings.polynomial.multi_polynomial import MPolynomial
+    if not isinstance(g, MPolynomial):
+        raise ValueError("must be a multivariate polynomial")
+
+    variables = g.variables()
+    if len(variables) != 2:
+        raise ValueError("must be a polynomial in two variables")
+
+    y, x = sorted(variables, key=g.degree)
+    if g.degree(y) != 2:
+        raise ValueError("must be a polynomial of degree 2 in a variable")
+
+    f = []
+    h = []
+    for k, v in g:
+        dx = v.degree(x)
+        dy = v.degree(y)
+        if dy == 2:
+            if dx != 0:
+                raise ValueError(f"cannot have a term y*x^{dx}")
+            y2 = k
+        elif dy == 1:
+            while len(h) <= dx:
+                h.append(0)
+            h[dx] = k
+        else:
+            assert dy == 0
+            while len(f) <= dx:
+                f.append(0)
+            f[dx] = -k
+
+    if not y2.is_one():
+        y2_inv = y2.inverse_of_unit()
+        f = [c * y2_inv for c in f]
+        h = [c * y2_inv for c in h]
+
+    return f, h
+
+
+def HyperellipticCurve(f, h=None, names=None, PP=None, check_squarefree=True):
     r"""
     Return the hyperelliptic curve `y^2 + h y = f`, for
     univariate polynomials `h` and `f`. If `h`
@@ -76,6 +138,12 @@ def HyperellipticCurve(f, h=0, names=None, PP=None, check_squarefree=True):
         Hyperelliptic Curve over Finite Field in a of size 3^2
          defined by y^2 + (x + a)*y = x^3 + x + 2
 
+    Construct from defining polynomial::
+
+        sage: R.<x,y> = QQ[]
+        sage: HyperellipticCurve(y^2 + 3*x^2*y - (x^5 + x + 1))
+        Hyperelliptic Curve over Rational Field defined by y^2 + 3*x^2*y = x^5 + x + 1
+
     Characteristic two::
 
         sage: # needs sage.rings.finite_rings
@@ -200,6 +268,17 @@ def HyperellipticCurve(f, h=0, names=None, PP=None, check_squarefree=True):
     """
     # F is the discriminant; use this for the type check
     # rather than f and h, one of which might be constant.
+    if h is None:
+        from sage.rings.polynomial.multi_polynomial import MPolynomial
+        if isinstance(f, MPolynomial) and len(f.parent().gens()) == 2:
+            from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
+            from sage.structure.element import get_coercion_model
+            P = PolynomialRing(f.base_ring(), 'x')
+            f, h = _parse_multivariate_defining_equation(f)
+            f, h = P(f), P(h)
+        else:
+            h = 0
+
     F = h**2 + 4 * f
     if not isinstance(F, Polynomial):
         raise TypeError(f"arguments f = {f} and h = {h} must be polynomials")