Fix random polynomial bias

src/sage/crypto/lattice.py

@@ -112,10 +112,10 @@ def gen_lattice(type='modular', n=4, m=8, q=11, seed=None,
         [ 0 11  0  0  0  0  0  0]
         [ 0  0 11  0  0  0  0  0]
         [ 0  0  0 11  0  0  0  0]
-        [-2 -3 -3  4  1  0  0  0]
-        [ 4 -2 -3 -3  0  1  0  0]
-        [-3  4 -2 -3  0  0  1  0]
-        [-3 -3  4 -2  0  0  0  1]
+        [-3 -3 -2  4  1  0  0  0]
+        [ 4 -3 -3 -2  0  1  0  0]
+        [-2  4 -3 -3  0  0  1  0]
+        [-3 -2  4 -3  0  0  0  1]
 
     Ideal bases also work with polynomials::
 
@@ -125,10 +125,10 @@ def gen_lattice(type='modular', n=4, m=8, q=11, seed=None,
         [ 0 11  0  0  0  0  0  0]
         [ 0  0 11  0  0  0  0  0]
         [ 0  0  0 11  0  0  0  0]
-        [ 1  4 -3  3  1  0  0  0]
-        [ 3  1  4 -3  0  1  0  0]
-        [-3  3  1  4  0  0  1  0]
-        [ 4 -3  3  1  0  0  0  1]
+        [-3  4  1  4  1  0  0  0]
+        [ 4 -3  4  1  0  1  0  0]
+        [ 1  4 -3  4  0  0  1  0]
+        [ 4  1  4 -3  0  0  0  1]
 
     Cyclotomic bases with n=2^k are SWIFFT bases::
 
@@ -137,10 +137,10 @@ def gen_lattice(type='modular', n=4, m=8, q=11, seed=None,
         [ 0 11  0  0  0  0  0  0]
         [ 0  0 11  0  0  0  0  0]
         [ 0  0  0 11  0  0  0  0]
-        [-2 -3 -3  4  1  0  0  0]
-        [-4 -2 -3 -3  0  1  0  0]
-        [ 3 -4 -2 -3  0  0  1  0]
-        [ 3  3 -4 -2  0  0  0  1]
+        [-3 -3 -2  4  1  0  0  0]
+        [-4 -3 -3 -2  0  1  0  0]
+        [ 2 -4 -3 -3  0  0  1  0]
+        [ 3  2 -4 -3  0  0  0  1]
 
     Dual modular bases are related to Regev's famous public-key
     encryption [Reg2005]_::

src/sage/crypto/lwe.py

@@ -670,7 +670,7 @@ def __init__(self, ringlwe):
             sage: lwe = RingLWEConverter(RingLWE(16, 257, D, secret_dist='uniform'))
             sage: set_random_seed(1337)
             sage: lwe()
-            ((32, 216, 3, 125, 58, 197, 171, 43), ...)
+            ((171, 197, 58, 125, 3, 216, 32, 130), ...)
         """
         self.ringlwe = ringlwe
         self._i = 0
@@ -686,7 +686,7 @@ def __call__(self):
             sage: lwe = RingLWEConverter(RingLWE(16, 257, D, secret_dist='uniform'))
             sage: set_random_seed(1337)
             sage: lwe()
-            ((32, 216, 3, 125, 58, 197, 171, 43), ...)
+            ((171, 197, 58, 125, 3, 216, 32, 130), ...)
         """
         R_q = self.ringlwe.R_q
 

src/sage/misc/randstate.pyx

@@ -56,22 +56,22 @@ results of these random number generators reproducible. ::
 
     sage: set_random_seed(0)
     sage: print(rtest())
-    (303, -0.266166246380421, 1/6, (1,2), [ 0, 1, 1, 0, 0 ], 265625921, 79302, 0.2450652680687958)
+    (303, -0.266166246380421, 1/2*x^2 - 1/95*x - 1/2, (1,3), [ 1, 0, 0, 1, 1 ], 265625921, 5842, 0.9661911734708414)
     sage: set_random_seed(1)
     sage: print(rtest())
-    (978, 0.0557699430711638, -1/8*x^2 - 1/2*x + 1/2, (1,2,3), [ 1, 0, 0, 0, 1 ], 807447831, 23865, 0.6170498912488264)
+    (978, 0.0557699430711638, -3*x^2 - 1/12, (1,2), [ 0, 0, 1, 1, 0 ], 807447831, 29982, 0.8335077654199736)
     sage: set_random_seed(2)
     sage: print(rtest())
-    (207, -0.0141049486533456, 0, (1,3)(4,5), [ 1, 0, 1, 1, 1 ], 1642898426, 16190, 0.9343331114872127)
+    (207, -0.0141049486533456, 4*x^2 + 1/2, (1,2)(4,5), [ 1, 0, 0, 1, 1 ], 1642898426, 41662, 0.19982565117278328)
     sage: set_random_seed(0)
     sage: print(rtest())
-    (303, -0.266166246380421, 1/6, (1,2), [ 0, 1, 1, 0, 0 ], 265625921, 79302, 0.2450652680687958)
+    (303, -0.266166246380421, 1/2*x^2 - 1/95*x - 1/2, (1,3), [ 1, 0, 0, 1, 1 ], 265625921, 5842, 0.9661911734708414)
     sage: set_random_seed(1)
     sage: print(rtest())
-    (978, 0.0557699430711638, -1/8*x^2 - 1/2*x + 1/2, (1,2,3), [ 1, 0, 0, 0, 1 ], 807447831, 23865, 0.6170498912488264)
+    (978, 0.0557699430711638, -3*x^2 - 1/12, (1,2), [ 0, 0, 1, 1, 0 ], 807447831, 29982, 0.8335077654199736)
     sage: set_random_seed(2)
     sage: print(rtest())
-    (207, -0.0141049486533456, 0, (1,3)(4,5), [ 1, 0, 1, 1, 1 ], 1642898426, 16190, 0.9343331114872127)
+    (207, -0.0141049486533456, 4*x^2 + 1/2, (1,2)(4,5), [ 1, 0, 0, 1, 1 ], 1642898426, 41662, 0.19982565117278328)
 
 Once we've set the random number seed, we can check what seed was used.
 (This is not the current random number state; it does not change when
@@ -81,7 +81,7 @@ random numbers are generated.)  ::
     sage: initial_seed()
     12345
     sage: print(rtest())
-    (720, -0.612180244315804, 0, (1,3), [ 1, 0, 1, 1, 0 ], 1911581957, 65175, 0.8043027951758298)
+    (720, -0.612180244315804, x^2 - x, (1,2,3), [ 1, 0, 0, 0, 1 ], 1911581957, 27093, 0.9205331599518184)
     sage: initial_seed()
     12345
 
@@ -216,19 +216,19 @@ that you get without intervening ``with seed``. ::
 
     sage: set_random_seed(0)
     sage: r1 = rtest(); print(r1)
-    (303, -0.266166246380421, 1/6, (1,2), [ 0, 1, 1, 0, 0 ], 265625921, 79302, 0.2450652680687958)
+    (303, -0.266166246380421, 1/2*x^2 - 1/95*x - 1/2, (1,3), [ 1, 0, 0, 1, 1 ], 265625921, 5842, 0.9661911734708414)
     sage: r2 = rtest(); print(r2)
-    (443, 0.185001351421963, -2, (1,3), [ 0, 0, 1, 1, 0 ], 53231108, 8171, 0.28363811590618193)
+    (105, 0.642309615982449, -x^2 - x - 6, (1,3)(4,5), [ 1, 0, 0, 0, 1 ], 53231108, 77132, 0.001767155077382232)
 
 We get slightly different results with an intervening ``with seed``. ::
 
     sage: set_random_seed(0)
     sage: r1 == rtest()
     True
     sage: with seed(1): rtest()
-    (978, 0.0557699430711638, -1/8*x^2 - 1/2*x + 1/2, (1,2,3), [ 1, 0, 0, 0, 1 ], 807447831, 23865, 0.6170498912488264)
+    (978, 0.0557699430711638, -3*x^2 - 1/12, (1,2), [ 0, 0, 1, 1, 0 ], 807447831, 29982, 0.8335077654199736)
     sage: r2m = rtest(); r2m
-    (443, 0.185001351421963, -2, (1,3), [ 0, 0, 1, 1, 0 ], 53231108, 51295, 0.28363811590618193)
+    (105, 0.642309615982449, -x^2 - x - 6, (1,3)(4,5), [ 1, 0, 0, 0, 1 ], 53231108, 40267, 0.001767155077382232)
     sage: r2m == r2
     False
 
@@ -245,8 +245,8 @@ case, as we see in this example::
     sage: with seed(1):
     ....:     print(rtest())
     ....:     print(rtest())
-    (978, 0.0557699430711638, -1/8*x^2 - 1/2*x + 1/2, (1,2,3), [ 1, 0, 0, 0, 1 ], 807447831, 23865, 0.6170498912488264)
-    (181, 0.607995392046754, -x + 1/2, (2,3)(4,5), [ 1, 0, 0, 1, 1 ], 1010791326, 9693, 0.5691716786307407)
+    (978, 0.0557699430711638, -3*x^2 - 1/12, (1,2), [ 0, 0, 1, 1, 0 ], 807447831, 29982, 0.8335077654199736)
+    (138, -0.0404945051288503, 2*x - 24, (1,2,3), [ 1, 1, 0, 1, 1 ], 1010791326, 91360, 0.0033332230808060803)
     sage: r2m == rtest()
     True
 
@@ -258,7 +258,7 @@ NTL random numbers were generated inside the ``with seed``.
     True
     sage: with seed(1):
     ....:     rtest()
-    (978, 0.0557699430711638, -1/8*x^2 - 1/2*x + 1/2, (1,2,3), [ 1, 0, 0, 0, 1 ], 807447831, 23865, 0.6170498912488264)
+    (978, 0.0557699430711638, -3*x^2 - 1/12, (1,2), [ 0, 0, 1, 1, 0 ], 807447831, 29982, 0.8335077654199736)
     sage: r2m == rtest()
     True
 

src/sage/rings/polynomial/polynomial_ring.py

@@ -1320,7 +1320,7 @@
             sage: R.monomial(m.degree()) == m
             True
         """
-        return self({exponent:self.base_ring().one()})
+        return self({exponent: self.base_ring().one()})
 
     def krull_dimension(self):
         """
@@ -1362,23 +1362,25 @@
         """
         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
 
-        -  ``*args, **kwds`` - Passed on to the ``random_element`` method for
-           the base ring
+        -  ``monic`` -- boolean (optional); indicate whether the sampled
+           polynomial should be monic
+
+        -  ``*args, **kwds`` -- additional keyword parameters passed on to the
+           ``random_element`` method for the base ring
 
         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.degree()
             10
             sage: f.parent() is R
@@ -1388,21 +1390,30 @@
             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
+        among all polynomials with degree between them. If the base ring can be sampled uniformly,
+        then this method also samples uniformly::
 
-            sage: R.random_element(degree=(0, 8)).degree() in range(0, 9)
+            sage: R.random_element(degree=(0, 4)).degree() in range(0, 5)
             True
-            sage: found = [False]*9
+            sage: found = [False]*5
             sage: while not all(found):
-            ....:     found[R.random_element(degree=(0, 8)).degree()] = True
+            ....:     found[R.random_element(degree=(0, 4)).degree()] = True
 
         Note that the zero polynomial has degree `-1`, so if you want to
         consider it set the minimum degree to `-1`::
 
             sage: while R.random_element(degree=(-1,2), x=-1, y=1) != R.zero():
             ....:     pass
 
+        Monic polynomials are chosen among all monic polynomials with degree between the given
+        ``degree`` argument::
+
+            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])
@@ -1419,14 +1430,42 @@
 
             sage: R = PolynomialRing(GF(2), 'z')
             sage: for _ in range(100):
-            ....:     d = randint(-1,20)
+            ....:     d = randint(-1, 20)
             ....:     P = R.random_element(degree=d)
-            ....:     assert P.degree() == d, "problem with {} which has not degree {}".format(P,d)
+            ....:     assert P.degree() == d
+
+        In :issue:`37118`, ranges including integers below `-1` no longer raise
+        an error::
+
+            sage: R.random_element(degree=(-2, 3))  # random
+            z^3 + z^2 + 1
+
+        ::
+
+            sage: 0 in [R.random_element(degree=(-1, 2), monic=True) for _ in range(500)]
+            False
+
+        Testing error handling::
+
+            sage: R.random_element(degree=-5)
+            Traceback (most recent call last):
+            ...
+            ValueError: degree (=-5) must be at least -1
 
-            sage: R.random_element(degree=-2)
+            sage: R.random_element(degree=(-3, -2))
             Traceback (most recent call last):
             ...
-            ValueError: degree should be an integer greater or equal than -1
+            ValueError: maximum degree (=-2) must be at least -1
+
+        Testing uniformity::
+
+            sage: from collections import Counter
+            sage: R = GF(3)["x"]
+            sage: samples = [R.random_element(degree=(-1, 2)) for _ in range(27000)]    # long time
+            sage: assert all(750 <= f <= 1250 for f in Counter(samples).values())       # long time
+
+            sage: samples = [R.random_element(degree=(-1, 2), monic=True) for _ in range(13000)] # long time
+            sage: assert all(750 <= f <= 1250 for f in Counter(samples).values())       # long time
         """
         R = self.base_ring()
 
@@ -1435,11 +1474,15 @@
                 raise ValueError("degree argument must be an integer or a tuple of 2 integers (min_degree, max_degree)")
             if degree[0] > degree[1]:
                 raise ValueError("minimum degree must be less or equal than maximum degree")
+            if degree[1] < -1:
+                raise ValueError(f"maximum degree (={degree[1]}) must be at least -1")
         else:
-            degree = (degree,degree)
+            if degree < -1:
+                raise ValueError(f"degree (={degree}) must be at least -1")
+            degree = (degree, degree)
 
         if degree[0] <= -2:
-            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,
@@ -1450,24 +1493,44 @@
             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 and monic:
+            if degree[1] == -1:
+                raise ValueError("the maximum degree of monic polynomials needs to be at least 0")
+            if degree[1] == 0:
+                return self.one()
+            degree = (0, degree[1])
 
         # Pick random coefficients
-        p = self([R.random_element(*args, **kwds) for _ in range(d)])
+        end = degree[1]
+        if degree[0] == -1:
+            return self([R.random_element(*args, **kwds) for _ in range(end + 1)])
+
+        nonzero = False
+        coefs = [None] * (end + 1)
+
+        while not nonzero:
+            # Pick leading coefficients, if `monic` is set it's handle here.
+            if monic:
+                for i in range(degree[1] - degree[0] + 1):
+                    coefs[end - i] = R.random_element(*args, **kwds)
+                    if not nonzero and not coefs[end - i].is_zero():
+                        coefs[end - i] = R.one()
+                        nonzero = True
+            else:
+                # Fast path
+                for i in range(degree[1] - degree[0] + 1):
+                    coefs[end - i] = R.random_element(*args, **kwds)
+                    nonzero |= not coefs[end - i].is_zero()
 
-        # Add non-zero leading coefficient
-        p += R._random_nonzero_element(*args, **kwds) * self.gen() ** d
+        # Now we pick the remaining coefficients.
+        for i in range(degree[1] - degree[0] + 1, degree[1] + 1):
+            coefs[end - i] = R.random_element(*args, **kwds)
 
-        return p
+        return self(coefs)
 
     def _monics_degree(self, of_degree):
         """
@@ -1587,11 +1650,11 @@
 
         INPUT: Pass exactly one of:
 
-        -  ``max_degree`` - an int; the iterator will generate
+        -  ``max_degree`` -- an int; the iterator will generate
            all polynomials which have degree less than or equal to
            ``max_degree``
 
-        -  ``of_degree`` - an int; the iterator will generate
+        -  ``of_degree`` -- an int; the iterator will generate
            all polynomials which have degree ``of_degree``
 
         OUTPUT: an iterator
@@ -1653,11 +1716,11 @@
         INPUT: Pass exactly one of:
 
 
-        -  ``max_degree`` - an int; the iterator will generate
+        -  ``max_degree`` -- an int; the iterator will generate
            all monic polynomials which have degree less than or equal to
            ``max_degree``
 
-        -  ``of_degree`` - an int; the iterator will generate
+        -  ``of_degree`` -- an int; the iterator will generate
            all monic polynomials which have degree ``of_degree``
 
 
@@ -1734,7 +1797,7 @@
 
         INPUT:
 
-        - ``f`` - either a polynomial in ``self``, or a principal
+        - ``f`` -- either a polynomial in ``self``, or a principal
           ideal of ``self``.
         - further named arguments that are passed to the quotient constructor.
 
@@ -2414,8 +2477,8 @@
 #                 P += (F[i] * prod)
 #             return P
 
-        # using Neville's method for recursively generating the
-        # Lagrange interpolation polynomial
+# using Neville's method for recursively generating the
+# Lagrange interpolation polynomial
         elif algorithm == "neville":
             if previous_row is None:
                 previous_row = []