Fix random polynomial bias

src/sage/rings/polynomial/polynomial_ring.py

@@ -1360,23 +1360,32 @@ def ngens(self):
         """
         return 1
 
-    def random_element(self, degree=(-1,2), *args, **kwds):
+    def random_element(self, degree=(-1, 2), monic=False, *args, **kwds):
         r"""
-        Return a random polynomial of given degree or with given degree bounds.
+        Return a random polynomial of given degree (bounds).
 
         INPUT:
 
-        -  ``degree`` - optional integer for fixing the degree
-           or a tuple of minimum and maximum degrees. By default set to
-           ``(-1,2)``.
+        -  ``degree`` - (default: ``(-1, 2)``) integer for fixing the degree or
+           a tuple of minimum and maximum degrees
+
+        -  ``monic`` - optional boolean to indicate whether the sampled
+           polynomial should be monic, or not. If this is set, `0` is not a
+           possible output of the method
 
         -  ``*args, **kwds`` - Passed on to the ``random_element`` method for
            the base ring
 
+        .. SEEALSO::
+
+            :meth:`random_monic_element`
+
         EXAMPLES::
 
             sage: R.<x> = ZZ[]
-            sage: f = R.random_element(10, 5, 10)
+            sage: f = R.random_element(10, x=5, y=10)
+            sage: f  # random
+            5*x^10 + 6*x^9 + 5*x^8 + 8*x^7 + 8*x^6 + 5*x^5 + 7*x^4 + 8*x^3 + 6*x^2 + 9*x + 6
             sage: f.degree()
             10
             sage: f.parent() is R
@@ -1386,8 +1395,8 @@ def random_element(self, degree=(-1,2), *args, **kwds):
             sage: R.random_element(6).degree()
             6
 
-        If a tuple of two integers is given for the ``degree`` argument, a degree
-        is first uniformly chosen, then a polynomial of that degree is given::
+        If a tuple of two integers is given for the ``degree`` argument, a polynomial is chosen
+        uniformly among all polynomials with degree between them::
 
             sage: R.random_element(degree=(0, 8)).degree() in range(0, 9)
             True
@@ -1401,6 +1410,21 @@ def random_element(self, degree=(-1,2), *args, **kwds):
             sage: while R.random_element(degree=(-1,2), x=-1, y=1) != R.zero():
             ....:     pass
 
+        It is possible to sample a monic polynomial, but it is preferable to
+        use :meth:`random_monic_element`::
+
+            sage: R.random_element(degree=(-1, 5), monic=True).is_monic()
+            True
+            sage: R.random_monic_element(degree=(-1, 5)).is_monic()
+            True
+
+        Note that if the degree range includes `-1` and ``monic`` is set, it is silently ignored::
+
+            sage: all(R.random_element(degree=(-1, 1), monic=True).is_monic() for _ in range(10^3))
+            True
+            sage: all(R.random_element(degree=(0, 1), monic=True).is_monic() for _ in range(10^3))
+            True
+
         TESTS::
 
             sage: R.random_element(degree=[5])
@@ -1421,10 +1445,17 @@ def random_element(self, degree=(-1,2), *args, **kwds):
             ....:     P = R.random_element(degree=d)
             ....:     assert P.degree() == d, "problem with {} which has not degree {}".format(P,d)
 
-            sage: R.random_element(degree=-2)
-            Traceback (most recent call last):
-            ...
-            ValueError: degree should be an integer greater or equal than -1
+        In :issue:`37118`, ranges including integers below `-1` no longer error::
+
+            sage: R.random_element(degree=(-2, 3))  # random
+            1
+
+        ::
+
+            sage: 0 in [R.random_element(degree=(-1, 2), monic=True) for _ in range(500)]
+            False
+            sage: 0 in [R.random_monic_element(degree=(-1, 2)) for _ in range(500)]
+            False
         """
         R = self.base_ring()
 
@@ -1437,7 +1468,9 @@ def random_element(self, degree=(-1,2), *args, **kwds):
             degree = (degree,degree)
 
         if degree[0] <= -2:
-            raise ValueError("degree should be an integer greater or equal than -1")
+            # This error has been removed in issue #37118.
+            # raise ValueError("degree should be an integer greater or equal than -1")
+            degree = (-1, degree[1])
 
         # If the coefficient range only contains 0, then
         # * if the degree range includes -1, return the zero polynomial,
@@ -1448,24 +1481,101 @@ def random_element(self, degree=(-1,2), *args, **kwds):
             else:
                 raise ValueError("No polynomial of degree >= 0 has all coefficients zero")
 
-        # Pick a random degree
-        d = randint(degree[0], degree[1])
-
-        # If degree is -1, return the 0 polynomial
-        if d == -1:
+        if degree == (-1, -1):
             return self.zero()
 
-        # If degree is 0, return a random constant term
-        if d == 0:
-            return self(R._random_nonzero_element(*args, **kwds))
+        # If `monic` is set, zero should be ignored
+        if degree[0] == -1:
+            degree = (0, degree[1])
+            has_zero = not monic
+        else:
+            has_zero = False
+
+        while True:
+            # Pick random coefficients
+            coefs = [None] * (degree[1] + 1)
+            nonzero = False
+
+            # Pick leading coefficients, if `monic` is set it's handle here.
+            if monic:
+                for i in range(degree[1] - degree[0] + 1):
+                    coefs[i] = R.random_element(*args, **kwds)
+                    if not nonzero and not coefs[i].is_zero():
+                        coefs[i] = R.one()
+                        nonzero = True
+
+                # Leading terms must be nonzero
+                if not nonzero:
+                    continue
+            else:
+                # Fast path
+                for i in range(degree[1] - degree[0] + 1):
+                    coefs[i] = R.random_element(*args, **kwds)
+                    nonzero |= not coefs[i].is_zero()
+
+                if not nonzero and not has_zero:
+                    continue
 
-        # Pick random coefficients
-        p = self([R.random_element(*args, **kwds) for _ in range(d)])
+            # Now we pick the remaining coefficients. Zeros still should be
+            # tracked to handle `has_zero`.
+            for i in range(degree[1] - degree[0] + 1, degree[1] + 1):
+                coefs[i] = R.random_element(*args, **kwds)
+                nonzero |= not coefs[i].is_zero()
 
-        # Add non-zero leading coefficient
-        p += R._random_nonzero_element(*args, **kwds) * self.gen() ** d
+            # If we don't want zero (not has_zero), but coefs is zero (not
+            # nonzero), then reject
+            if not has_zero and not nonzero:
+                continue
 
-        return p
+            coefs.reverse()
+            return self(coefs)
+
+    def random_monic_element(self, degree=(0, 2), *args, **kwargs):
+        r"""
+        Return a random monic polynomial of given degree (bounds).
+
+        Calls :meth:`random_element` with ``monic=True``.
+
+        INPUT:
+
+        -  ``degree`` - optional integer for fixing the degree
+           or a tuple of minimum and maximum degrees. By default set to
+           ``(0, 2)``
+
+        - ``zero`` - boolean, if set the algorithm may return `0`, even though
+          it is not a monic polynomial. By default set to ``False``
+
+        -  ``*args, **kwds`` - Passed on to the ``random_element`` method for
+           the base ring
+
+        .. SEEALSO::
+
+            :meth:`random_element`
+
+        EXAMPLES::
+
+            sage: R.<x> = GF(3)[]
+            sage: R.random_monic_element()  # random
+            x^2 + x + 1
+            sage: all(R.random_monic_element().is_monic() for _ in range(100))
+            True
+            sage: all(R.random_monic_element(degree=(-1, 1)).is_monic() for _ in range(100))
+            True
+
+        TESTS:
+
+        There are 7 monic polynomials of degree not exceeding 2 over GF(2). We
+        apply the chi-square test::
+
+            sage: from collections import Counter
+            sage: from scipy.stats import chisquare
+            sage: N = 10^5
+            sage: R.<x> = GF(2)[]
+            sage: cnts = Counter(R.random_monic_element() for _ in range(N))
+            sage: chisquare(list(cnts.values()), [N / 7] * 7).pvalue < 0.1
+            False
+        """
+        return self.random_element(degree=degree, monic=True, *args, **kwargs)
 
     def _monics_degree(self, of_degree):
         """