support calling PARI's qfbcornacchia() from BinaryQF

src/sage/quadratic_forms/binary_qf.py

@@ -51,7 +51,7 @@
 # ****************************************************************************
 from functools import total_ordering
 
-from sage.libs.pari.all import pari_gen
+from sage.libs.pari.all import pari_gen, pari
 from sage.rings.all import ZZ, is_fundamental_discriminant
 from sage.arith.all import gcd
 from sage.structure.sage_object import SageObject
@@ -1540,7 +1540,7 @@ def small_prime_value(self, Bmax=1000):
                 raise ValueError("Unable to find a prime value of %s" % self)
             B += 10
 
-    def solve_integer(self, n):
+    def solve_integer(self, n, *, algorithm="general"):
         r"""
         Solve `Q(x, y) = n` in integers `x` and `y` where `Q` is this
         quadratic form.
@@ -1549,19 +1549,40 @@ def solve_integer(self, n):
 
         - ``n`` -- a positive integer
 
+        - ``algorithm`` -- ``"general"`` (default) or ``"cornacchia"``
+
+        To use the Cornacchia algorithm, the quadratic form must have
+        `a=1` and `b=0` and `c>0`, and ``n`` must be a prime or four
+        times a prime (but this is not checked).
+
         OUTPUT:
 
         A tuple `(x, y)` of integers satisfying `Q(x, y) = n`, or ``None``
         if no solution exists.
 
-        ALGORITHM: :pari:`qfbsolve`
+        ALGORITHM: :pari:`qfbsolve` or :pari:`qfbcornacchia`
 
         EXAMPLES::
 
             sage: Q = BinaryQF([1, 0, 419])
             sage: Q.solve_integer(773187972)
             (4919, 1337)
 
+        If `Q` is of the form `[1,0,c]` as above and `n` is a prime or
+        four times a prime, Cornacchia's algorithm can be used, which is
+        typically much faster than the general method::
+
+            sage: Q = BinaryQF([1, 0, 12345])
+            sage: n = 2^99 + 5273
+            sage: Q.solve_integer(n)
+            (-67446480057659, 7139620553488)
+            sage: Q.solve_integer(n, algorithm='cornacchia')
+            (67446480057659, 7139620553488)
+            sage: timeit('Q.solve_integer(n)')                          # not tested
+            125 loops, best of 3: 3.13 ms per loop
+            sage: timeit('Q.solve_integer(n, algorithm="cornacchia")')  # not tested
+            625 loops, best of 3: 18.6 μs per loop
+
         ::
 
             sage: Qs = BinaryQF_reduced_representatives(-23, primitive_only=True)
@@ -1584,6 +1605,20 @@ def solve_integer(self, n):
             sage: xy is None or Q(*xy) == n
             True
 
+        ...also when using the ``"cornacchia"`` algorithm::
+
+            sage: abc = [1,0,randrange(1,10^3)]
+            sage: Q = BinaryQF(abc)
+            sage: n = random_prime(10^9)
+            sage: if randrange(2):
+            ....:     n *= 4
+            sage: xy1 = Q.solve_integer(n, algorithm='cornacchia')
+            sage: xy1 is None or Q(*xy1) == n
+            True
+            sage: xy2 = Q.solve_integer(n)
+            sage: (xy1 is None) == (xy2 is None)
+            True
+
         Test for square discriminants specifically (:trac:`33026`)::
 
             sage: n = randrange(-10^3, 10^3)
@@ -1625,6 +1660,13 @@ def solve_integer(self, n):
             Q = self.matrix_action_right(M)
             assert not Q._c
 
+            if not Q._b:
+                # at this point, Q = a*x^2
+                if Q._a.divides(n) and (n // Q._a).is_square():
+                    x = (n // Q._a).isqrt()
+                    return tuple(row[0] * x for row in M.rows())
+                return None
+
             # at this point, Q = a*x^2 + b*x*y
             if not n:
                 return tuple(M.columns()[1])
@@ -1636,6 +1678,15 @@ def solve_integer(self, n):
 
             return None
 
+        if algorithm == 'cornacchia':
+            if not (self._a.is_one() and self._b.is_zero() and self._c > 0):
+                raise ValueError("Cornacchia's algorithm requires a=1 and b=0 and c>0")
+            sol = pari.qfbcornacchia(self._c, n)
+            return tuple(map(ZZ, sol)) if sol else None
+
+        if algorithm != 'general':
+            raise NotImplementedError(f'algorithm {algorithm!r} is unknown')
+
         flag = 2  # single solution, possibly imprimitive
         sol = self.__pari__().qfbsolve(n, flag)
         return tuple(map(ZZ, sol)) if sol else None