Fix random polynomial bias

src/sage/rings/polynomial/polynomial_ring.py

@@ -1360,7 +1360,7 @@ 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.
 
@@ -1370,6 +1370,10 @@ def random_element(self, degree=(-1,2), *args, **kwds):
            or a tuple of minimum and maximum degrees. By default set to
            ``(-1,2)``.
 
+        -  ``monic`` - optional boolean to indicate whether the sample
+           polynomial should be monic or not. If degree includes -1, it is still
+           possible for the algorithm to return ``0``, which is not monic.
+
         -  ``*args, **kwds`` - Passed on to the ``random_element`` method for
            the base ring
 
@@ -1386,8 +1390,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 +1405,16 @@ 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::
+
+            sage: R.random_element(5, monic=True).is_monic()
+            True
+
+        Note that if the degree range includes `-1`, then the resulting polynomial might not be monic::
+
+            sage: all(R.random_element(degree=(-1, 1), monic=True).is_monic() for _ in range(10^3))
+            False
+
         TESTS::
 
             sage: R.random_element(degree=[5])
@@ -1439,6 +1453,10 @@ def random_element(self, degree=(-1,2), *args, **kwds):
         if degree[0] <= -2:
             raise ValueError("degree should be an integer greater or equal than -1")
 
+        # Actually, it's probably more useful to let the user check this.
+        # if monic and degree[0] == -1:
+        #     raise ValueError("degree cannot include -1 when monic is set")
+
         # If the coefficient range only contains 0, then
         # * if the degree range includes -1, return the zero polynomial,
         # * otherwise raise a value error
@@ -1448,24 +1466,64 @@ 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 == (-1, -1):
+            return 0
+
+        if degree[0] == -1:
+            allow_zero = True
+            degree = (0, degree[1])
+        else:
+            allow_zero = False
+
+        while True:
+            # Pick random coefficients
+            coefs = []
+            nonzero = False
+
+            for _ in range(degree[1] - degree[0] + 1):
+                c = R.random_element(*args, **kwds)
+                coefs.append(c)
+                if c != 0 and not nonzero:
+                    coefs[-1] = 1
+                    nonzero = True
+
+            if not (allow_zero or nonzero):
+                continue
 
-        # If degree is -1, return the 0 polynomial
-        if d == -1:
-            return self.zero()
+            for _ in range(degree[0]):
+                c = R.random_element(*args, **kwds)
+                coefs.append(c)
 
-        # If degree is 0, return a random constant term
-        if d == 0:
-            return self(R._random_nonzero_element(*args, **kwds))
+            coefs.reverse()
+            return self(coefs)
 
-        # Pick random coefficients
-        p = self([R.random_element(*args, **kwds) for _ in range(d)])
+    def random_monic_element(self, degree=(0, 2), *args, **kwargs):
+        r"""
+        Return a random monic polynomial of given degree or with 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)``. If the bounds include ``-1``, it is still possible for the
+           algorithm to return ``0``, which is not monic.
 
-        # Add non-zero leading coefficient
-        p += R._random_nonzero_element(*args, **kwds) * self.gen() ** d
+        -  ``*args, **kwds`` - Passed on to the ``random_element`` method for
+           the base ring
+
+        EXAMPLES::
 
-        return p
+            sage: R.<x> = GF(3)[]
+            sage: R.random_element()  # random
+            x^2 + 11*x + 5
+            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
+        """
+        return self.random_element(degree=degree, monic=True, *args, **kwargs)
 
     def _monics_degree(self, of_degree):
         """