gh-40317: Implement partial integer factorization using flint
    
As in the title. This is more efficient than trial division, which is
the only currently implemented partial integer factorization algorithm.

Side note: it may be possible to reuse `limit=` argument and set
`bits=limit.bit_length()`, but this *kind of* break backwards
compatibility. (or maybe not?)

### 📝 Checklist

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

- [x] The title is concise and informative.
- [x] The description explains in detail what this PR is about.
- [ ] I have linked a relevant issue or discussion.
- [x] 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: #40317
Reported by: user202729
Reviewer(s): Sahil Jain, Travis Scrimshaw, user202729

Author: Release Manager

src/sage/arith/misc.py

@@ -2560,7 +2560,7 @@ def trial_division(n, bound=None):
         return ZZ(n).trial_division(bound)
 
 
-def factor(n, proof=None, int_=False, algorithm='pari', verbose=0, **kwds):
+def factor(n, proof=None, int_=False, algorithm=None, verbose=0, **kwds):
     """
     Return the factorization of ``n``.  The result depends on the
     type of ``n``.
@@ -2734,6 +2734,11 @@ def factor(n, proof=None, int_=False, algorithm='pari', verbose=0, **kwds):
 
         sage: len(factor(2^2203-1,proof=false))
         1
+
+    Test ``limit``::
+
+        sage: factor(2990, limit=10)
+        2 * 5 * 299
     """
     try:
         m = n.factor

src/sage/matrix/matrix_integer_dense.pyx

@@ -2978,6 +2978,10 @@ cdef class Matrix_integer_dense(Matrix_dense):
                           precision=precision)
 
             R = A.to_matrix(self.new_matrix())
+
+        else:
+            raise ValueError("unknown algorithm")
+
         return R
 
     def LLL(self, delta=None, eta=None, algorithm='fpLLL:wrapper', fp=None, prec=0, early_red=False, use_givens=False, use_siegel=False, transformation=False, **kwds):

src/sage/misc/sageinspect.py

@@ -1448,7 +1448,7 @@ def sage_getargspec(obj):
                     annotations={})
         sage: sage_getargspec(factor)
         FullArgSpec(args=['n', 'proof', 'int_', 'algorithm', 'verbose'],
-                    varargs=None, varkw='kwds', defaults=(None, False, 'pari', 0),
+                    varargs=None, varkw='kwds', defaults=(None, False, None, 0),
                     kwonlyargs=[], kwonlydefaults=None, annotations={})
 
     In the case of a class or a class instance, the :class:`FullArgSpec` of the
@@ -1808,7 +1808,7 @@ def sage_signature(obj):
         sage: sage_signature(identity_matrix)                                          # needs sage.modules
         <Signature (ring, n=0, sparse=False)>
         sage: sage_signature(factor)
-        <Signature (n, proof=None, int_=False, algorithm='pari', verbose=0, **kwds)>
+        <Signature (n, proof=None, int_=False, algorithm=None, verbose=0, **kwds)>
 
     In the case of a class or a class instance, the :class:`Signature` of the
     ``__new__``, ``__init__`` or ``__call__`` method is returned::
@@ -2671,8 +2671,8 @@ def __internal_tests():
 
     A cython function with default arguments (one of which is a string)::
 
-        sage: sage_getdef(sage.rings.integer.Integer.factor, obj_name='factor')
-        "factor(algorithm='pari', proof=None, limit=None, int_=False, verbose=0)"
+        sage: sage_getdef(sage.rings.integer.Integer.binomial, obj_name='binomial')
+        "binomial(m, algorithm='gmp')"
 
     This used to be problematic, but was fixed in :issue:`10094`::
 

src/sage/rings/factorint_flint.pyx

@@ -24,7 +24,7 @@ from sage.libs.flint.fmpz_factor_sage cimport *
 from sage.rings.integer cimport Integer
 
 
-def factor_using_flint(Integer n):
+def factor_using_flint(Integer n, unsigned bits=0):
     r"""
     Factor the nonzero integer ``n`` using FLINT.
 
@@ -35,6 +35,7 @@ def factor_using_flint(Integer n):
     INPUT:
 
     - ``n`` -- a nonzero sage Integer; the number to factor
+    - ``bits`` -- if nonzero, passed to ``fmpz_factor_smooth``
 
     OUTPUT:
 
@@ -89,7 +90,10 @@ def factor_using_flint(Integer n):
     fmpz_factor_init(factors)
 
     sig_on()
-    fmpz_factor(factors, p)
+    if bits:
+        fmpz_factor_smooth(factors, p, bits, 0)  # TODO make proved=* customizable
+    else:
+        fmpz_factor(factors, p)
     sig_off()
 
     pairs = fmpz_factor_to_pairlist(factors)

src/sage/rings/integer.pyx

@@ -3912,8 +3912,8 @@ cdef class Integer(sage.structure.element.EuclideanDomainElement):
             sig_off()
             return x
 
-    def factor(self, algorithm='pari', proof=None, limit=None, int_=False,
-               verbose=0):
+    def factor(self, algorithm=None, proof=None, limit=None, int_=False,
+               verbose=0, *, flint_bits=None):
         """
         Return the prime factorization of this integer as a
         formal Factorization object.
@@ -3932,7 +3932,7 @@ cdef class Integer(sage.structure.element.EuclideanDomainElement):
           - ``'magma'`` -- use the MAGMA computer algebra system (requires
             an installation of MAGMA)
 
-          - ``'qsieve'`` -- use Bill Hart's quadratic sieve code;
+          - ``'qsieve'`` -- use ``qsieve_factor`` in the FLINT library;
             WARNING: this may not work as expected, see qsieve? for
             more information
 
@@ -3946,6 +3946,10 @@ cdef class Integer(sage.structure.element.EuclideanDomainElement):
           given it must fit in a ``signed int``, and the factorization is done
           using trial division and primes up to limit
 
+        - ``flint_bits`` -- integer or ``None`` (default: ``None``); if specified,
+          perform only a partial factorization, primes at most ``2^flint_bits``
+          have a high probability of being detected
+
         OUTPUT: a Factorization object containing the prime factors and
         their multiplicities
 
@@ -3994,6 +3998,13 @@ cdef class Integer(sage.structure.element.EuclideanDomainElement):
             sage: n.factor(algorithm='flint')                                           # needs sage.libs.flint
             2 * 3 * 11 * 13 * 41 * 73 * 22650083 * 1424602265462161
 
+        Example with ``flint_bits``. The small prime factor is found since it is
+        much smaller than `2^{50}`, but not the large ones::
+
+            sage: n = next_prime(2^256) * next_prime(2^257) * next_prime(2^40)
+            sage: n.factor(algorithm='flint', flint_bits=50)                            # needs sage.libs.flint
+            1099511627791 * 2681...6291
+
         We factor using a quadratic sieve algorithm::
 
             sage: # needs sage.libs.pari
@@ -4021,14 +4032,11 @@ cdef class Integer(sage.structure.element.EuclideanDomainElement):
             sage: n.factor(algorithm='foobar')
             Traceback (most recent call last):
             ...
-            ValueError: Algorithm is not known
+            ValueError: algorithm is not known
         """
         from sage.structure.factorization import Factorization
         from sage.structure.factorization_integer import IntegerFactorization
 
-        if algorithm not in ['pari', 'flint', 'kash', 'magma', 'qsieve', 'ecm']:
-            raise ValueError("Algorithm is not known")
-
         cdef Integer n, p, unit
 
         if mpz_sgn(self.value) == 0:
@@ -4039,18 +4047,27 @@ cdef class Integer(sage.structure.element.EuclideanDomainElement):
             unit = one
         else:
             n = PY_NEW(Integer)
-            unit = PY_NEW(Integer)
             mpz_neg(n.value, self.value)
-            mpz_set_si(unit.value, -1)
-
-        if mpz_cmpabs_ui(n.value, 1) == 0:
-            return IntegerFactorization([], unit=unit, unsafe=True,
-                                        sort=False, simplify=False)
+            unit = smallInteger(-1)
 
         if limit is not None:
+            if algorithm is not None:
+                raise ValueError('trial division will always be used when limit is provided')
             from sage.rings.factorint import factor_trial_division
             return factor_trial_division(self, limit)
 
+        if algorithm is None:
+            algorithm = 'pari'
+        elif algorithm not in ['pari', 'flint', 'kash', 'magma', 'qsieve', 'ecm']:
+            raise ValueError("algorithm is not known")
+
+        if algorithm != 'flint' and flint_bits is not None:
+            raise ValueError("cannot specify flint_bits when algorithm is not flint")
+
+        if n.is_one():
+            return IntegerFactorization([], unit=unit, unsafe=True,
+                                        sort=False, simplify=False)
+
         if mpz_fits_slong_p(n.value):
             global n_factor_to_list
             if n_factor_to_list is None:
@@ -4080,7 +4097,9 @@ cdef class Integer(sage.structure.element.EuclideanDomainElement):
                                         sort=False, simplify=False)
         elif algorithm == 'flint':
             from sage.rings.factorint_flint import factor_using_flint
-            F = factor_using_flint(n)
+            if flint_bits is not None and flint_bits <= 0:
+                raise ValueError('flint_bits must be positive or None')
+            F = factor_using_flint(n, flint_bits or 0)
             F.sort()
             return IntegerFactorization(F, unit=unit, unsafe=True,
                                         sort=False, simplify=False)
@@ -6698,13 +6717,13 @@ cdef class Integer(sage.structure.element.EuclideanDomainElement):
         if not mpz_sgn(n.value):
             mpz_set_ui(t.value, 0)
             mpz_abs(g.value, self.value)
-            mpz_set_si(s.value, 1 if mpz_sgn(self.value) >= 0 else -1)
+            s = smallInteger(1 if mpz_sgn(self.value) >= 0 else -1)
             return g, s, t
 
         if not mpz_sgn(self.value):
             mpz_set_ui(s.value, 0)
             mpz_abs(g.value, n.value)
-            mpz_set_si(t.value, 1 if mpz_sgn(n.value) >= 0 else -1)
+            t = smallInteger(1 if mpz_sgn(n.value) >= 0 else -1)
             return g, s, t
 
         # both n and self are nonzero, so we need to do a division and