Trac #29935: implicitly fuzz RNG-dependent doctests with a random random seed

Our documentation and tests include a variety of
examples and tests involving randomness.

Those depend on a randomness seed, and up to Sage 9.1
they always used the same one: 0. Thus every run
of `sage -t` or `make [p]test[all][long]` would run
those "random" doctests deterministically.
In many cases, the output of those tests even relied
on that. As a result, random examples and tests
were actually testing that they were run non-randomly!

This is reminiscent of
an [https://xkcd.com/221/ xkcd comic illustration of random number
generators].

Based on these considerations and the
related [https://groups.google.com/d/msg/sage-
devel/c4UbKSdt3Aw/UQAo1iYoAAAJ sage-devel discussion],
we propose to:

- allow specifying a randomness when running tests (#29962),
- adapt tests involving randomness, making sure they
  test mathematical functionality independent of
  the randomness seed used to run them (see roadmap),
- default to a random randomness seed when none is specified
  (present ticket),

thus making those tests more robust and more useful,
by becoming more likely to reveal bugs in a variety
of cases including corner cases.

The first step (see #29962, merged in Sage 9.2)
adds a `--random-seed` flag to `sage -t`, allowing:
{{{
sage -t --long --random-seed=9876543210 src/sage/
}}}

Still, when no randomness seed is specified,
the default seed 0 is used. This means most testers
test with the same randomness seed, making "random"
doctests still mostly deterministic in practice.

Here is a way to pick a random randomness seed
and run tests with it (can be used to work on
the tickets in the roadmap):
{{{
$ randseed() {
    sage -c "import sage.misc.randstate as randstate; \
             randstate.set_random_seed(); \
             print(randstate.initial_seed())";
    }
$ SEED=$(randseed)
$ DIR=src/sage
$ echo "$ sage -t --long --random-seed=${SEED} ${DIR}" \
  && ./sage -t --long --random-seed=${SEED} ${DIR}
}}}

Once examples and tests involving randomness
have been adapted, the present ticket puts
the final touch by making it so that
when running tests with no seed specified,
a random one will be used:
{{{
sage -t src/sage/all.py
Running doctests with ID 2020-06-23-23-19-03-8003eea5.
...
sage -t --warn-long 89.5 --random-
seed=273987373473433732996760115183658447263 src/sage/all.py
    [16 tests, 0.73 s]
----------------------------------------------------------------------
All tests passed!
----------------------------------------------------------------------
Total time for all tests: 0.8 seconds
    cpu time: 0.7 seconds
    cumulative wall time: 0.7 seconds
}}}

Being displayed in the output, the seed used
can be used again if needed:
{{{
sage -t --warn-long 89.5 --random-
seed=273987373473433732996760115183658447263 src/sage/all.py
}}}
allowing to explore any problematic case
revealed by running the tests with that seed.

Roadmap:

- Allow fuzzing: #29962
- Make all parts of sage ready for default fuzzing:
  - #29945: coding
  - #29963: geometry
  - #29964: libs
  - #29965: graphs
  - #32107: groups
  - #29967: interfaces
  - #29968: algebras
  - #29969: misc
  - #29970: arith
  - #29971: categories
  - #29972: stats
  - #29973: sets
  - #29974: combinat
  - #29975: numerical, probability
  - #29976: matrix
  - #29977: modular
  - #29978: modules
  - #29979: rings
  - #29980: crypto
  - #29981: documentation
  - #29982: dynamics
  - #29983: finance
  - #29984: symbolic
  - #29985: schemes
  - #29986: plot
  - #32188: various missed tests along the way
- #32216: Update the developers guide for implicitly fuzzing doctests
- Finally make fuzzing default with this ticket.

Follow-up:

- #32544: Meta-ticket: Fix unstable doctests detected after #29935
- Remove `set_random_seed` in doctests where possible.

Errors discovered by this ticket:

- #29936: Hyperbolic geometry bugs revealed by fuzzing
- #29945: Failing doctest in `src/sage/coding/linear_code.py`
  (just a wrong doctest, will be fixed here)
- #29954: Unstable plotting
- #29956: Bug in `KlyachkoBundle_class.random_deformation`
- #29957: Bug in `ContinuedFraction` rounding
- #29958: Too many strong articulation points
- #29961: Random symbolic expression is completely unstable
- #30045: Bug in Reed-Solomon encoder and error-erasure decoder
- Index error with random derangement (fixed in #29974)
- #31890: simplify_hypergeometric is unstable
- #31891: `ZeroDivisonError` when creating polynomial system
- #31892: Conic parametrization broken
- #32075: Polynomial generic power trunk broken
- #32083: Various errors with polybori including segmentation fault
- #32084 `_nth_root_naive` fails for integer mod
- #32085: Errors when computing norms of padic elements
- #32086: apply_homography unstable for continue fraction
- #32095: DiFUB algorithm fails on some random graph
- #32108: Fix random tree on one or less vertices
- #32109: Fix 0/0 in ore function field
- #32111: Unstable minimal polynomial for element of 2-adic Eisenstein
Extension Field in pi defined by x^4 - 2*a
- #32117: Random relative number field checks only irreducibility over
QQ
- #32118: AlgebraicForm checks invariance with random matrix that can be
the identity
- #32124: SL2Z.random_element unstable
- #32125: random Ore polynomials do not respect minimum degree bound
- #32126: padic QpLC.random_element is broken
- #32127: gosper_iterator of continued fractions is unstable
- #32129: sage_input is unreliable for elements of ComplexField
- #32131: Cut width of graph with one edge incorrect
- #32132: Wrong gyration orbit length
- #32138: is_groebner fails over fraction fields
- #32141: Unstable doctest involving permutation groups
- #32169: Bug in edge disjoint spanning trees
- #32185: Failing weak order assertion on random symbolic expression
- #32186: Random bounded tolerance graph
- #32187: permutation group generated by list perms in L of degree n
incorrect when compared to GAP
- #32657: `plot_vector_field` unstable

URL: https://trac.sagemath.org/29935
Reported by: dimpase
Ticket author(s): Jonathan Kliem
Reviewer(s): Michael Orlitzky

Author: Release Manager

src/bin/sage

@@ -466,7 +466,7 @@ usage_advanced() {
     echo "                         labeled \"# optional\" or labeled"
     echo "                         \"# optional tag\" for any of the tags given."
     echo "     --randorder[=seed] -- randomize order of tests"
-    echo "     --random-seed[=seed]  -- random seed for fuzzing doctests"
+    echo "     --random-seed[=seed]  -- random seed (integer) for fuzzing doctests"
     echo "     --new            -- only test files modified since last commit"
     echo "     --initial        -- only show the first failure per block"
     echo "     --debug          -- drop into PDB after an unexpected error"

src/bin/sage-runtests

@@ -56,7 +56,7 @@ if __name__ == "__main__":
              'if "build" is listed, will also run tests specific to Sage\'s build/packaging system; '
              'if set to "all", then all tests will be run')
     parser.add_argument("--randorder", type=int, metavar="SEED", help="randomize order of tests")
-    parser.add_argument("--random-seed", dest="random_seed", type=int, metavar="SEED", help="random seed for fuzzing doctests")
+    parser.add_argument("--random-seed", dest="random_seed", type=int, metavar="SEED", help="random seed (integer) for fuzzing doctests")
     parser.add_argument("--global-iterations", "--global_iterations", type=int, default=0, help="repeat the whole testing process this many times")
     parser.add_argument("--file-iterations", "--file_iterations", type=int, default=0, help="repeat each file this many times, stopping on the first failure")
     parser.add_argument("--environment", type=str, default="sage.repl.ipython_kernel.all_jupyter", help="name of a module that provides the global environment for tests")

src/doc/en/developer/doctesting.rst

@@ -803,7 +803,33 @@ You can also pass in an explicit amount of time::
 Finally, you can disable any warnings about long tests with
 ``--warn-long 0``.
 
-Doctests may start from a random seed::
+Doctests start from a random seed::
+
+    [kliem@sage sage-9.2]$ sage -t src/sage/doctest/tests/random_seed.rst
+    Running doctests with ID 2020-06-23-23-22-59-49f37a55.
+    ...
+    Doctesting 1 file.
+    sage -t --warn-long 89.5 --random-seed=112986622569797306072457879734474628454 src/sage/doctest/tests/random_seed.rst
+    **********************************************************************
+    File "src/sage/doctest/tests/random_seed.rst", line 3, in sage.doctest.tests.random_seed
+    Failed example:
+        randint(5, 10)
+    Expected:
+        9
+    Got:
+        8
+    **********************************************************************
+    1 item had failures:
+       1 of   2 in sage.doctest.tests.random_seed
+        [1 test, 1 failure, 0.00 s]
+    ----------------------------------------------------------------------
+    sage -t --warn-long 89.5 --random-seed=112986622569797306072457879734474628454 src/sage/doctest/tests/random_seed.rst  # 1 doctest failed
+    ----------------------------------------------------------------------
+    Total time for all tests: 0.0 seconds
+        cpu time: 0.0 seconds
+        cumulative wall time: 0.0 seconds
+
+This seed can be set explicitly to reproduce possible failures::
 
     [kliem@sage sage-9.2]$ sage -t --warn-long 89.5 --random-seed=112986622569797306072457879734474628454 src/sage/doctest/tests/random_seed.rst
     Running doctests with ID 2020-06-23-23-24-28-14a52269.

src/sage/arith/misc.py

@@ -5667,12 +5667,25 @@ def sort_complex_numbers_for_display(nums):
         ....:                     RDF.random_element()))
         sage: shuffle(nums)
         sage: nums = sort_c(nums)
-        sage: nums[:3]
-        [0.0, 1.0, 2.0]
-        sage: for i in range(3, len(nums)-1):
-        ....:     assert nums[i].real() <= nums[i+1].real() + 1e-10
-        ....:     if abs(nums[i].real() - nums[i+1].real()) < 1e-10:
-        ....:         assert nums[i].imag() <= nums[i+1].imag() + 1e-10
+        sage: for i in range(len(nums)):
+        ....:     if nums[i].imag():
+        ....:         first_non_real = i
+        ....:         break
+        ....: else:
+        ....:     first_non_real = len(nums)
+        sage: assert first_non_real >= 3
+        sage: for i in range(first_non_real - 1):
+        ....:     assert nums[i].real() <= nums[i + 1].real()
+
+        sage: def truncate(n):
+        ....:     if n.real() < 1e-10:
+        ....:         return 0
+        ....:     else:
+        ....:         return n.real().n(digits=9)
+        sage: for i in range(first_non_real, len(nums)-1):
+        ....:     assert truncate(nums[i]) <= truncate(nums[i + 1])
+        ....:     if truncate(nums[i]) == truncate(nums[i + 1]):
+        ....:         assert nums[i].imag() <= nums[i+1].imag()
     """
     if not nums:
         return nums

src/sage/crypto/util.py

@@ -334,6 +334,8 @@ def carmichael_lambda(n):
         ....:     L = coprime(n)
         ....:     return list(map(power_mod, L, [k]*len(L), [n]*len(L)))
         sage: def my_carmichael(n):
+        ....:     if n == 1:
+        ....:         return 1
         ....:     for k in range(1, n):
         ....:         L = znpower(n, k)
         ....:         ones = [1] * len(L)

src/sage/doctest/control.py

@@ -28,6 +28,7 @@
 import re
 import types
 import sage.misc.flatten
+import sage.misc.randstate as randstate
 from sage.structure.sage_object import SageObject
 from sage.env import DOT_SAGE, SAGE_LIB, SAGE_SRC, SAGE_VENV, SAGE_EXTCODE
 from sage.misc.temporary_file import tmp_dir
@@ -461,7 +462,8 @@ def __init__(self, options, args):
         self._init_warn_long()
 
         if self.options.random_seed is None:
-            self.options.random_seed = 0
+            randstate.set_random_seed()
+            self.options.random_seed = randstate.initial_seed()
 
     def __del__(self):
         if getattr(self, 'logfile', None) is not None:

src/sage/functions/exp_integral.py

@@ -1497,7 +1497,7 @@ def exponential_integral_1(x, n=0):
         ....:     n = 2^ZZ.random_element(14)
         ....:     x = exponential_integral_1(a, n)
         ....:     y = exponential_integral_1(S(a), n)
-        ....:     c = RDF(2 * max(1.0, y[0]))
+        ....:     c = RDF(4 * max(1.0, y[0]))
         ....:     for i in range(n):
         ....:         e = float(abs(S(x[i]) - y[i]) << prec)
         ....:         if e >= c:

src/sage/functions/orthogonal_polys.py

@@ -2083,7 +2083,7 @@ class Func_ultraspherical(GinacFunction):
         32*t^3 - 12*t
         sage: _ = var('x')
         sage: for N in range(100):
-        ....:     n = ZZ.random_element().abs() + 5
+        ....:     n = ZZ.random_element(5, 5001)
         ....:     a = QQ.random_element().abs() + 5
         ....:     assert ((n+1)*ultraspherical(n+1,a,x) - 2*x*(n+a)*ultraspherical(n,a,x) + (n+2*a-1)*ultraspherical(n-1,a,x)).expand().is_zero()
         sage: ultraspherical(5,9/10,3.1416)

src/sage/graphs/generic_graph.py

@@ -6331,8 +6331,11 @@ def edge_disjoint_spanning_trees(self, k, root=None, solver=None, verbose=0):
         By Edmond's theorem, a graph which is `k`-connected always has `k`
         edge-disjoint arborescences, regardless of the root we pick::
 
-            sage: g = digraphs.RandomDirectedGNP(28, .3) # reduced from 30 to 28, cf. #9584
+            sage: g = digraphs.RandomDirectedGNP(11, .3)  # reduced from 30 to 11, cf. #32169
             sage: k = Integer(g.edge_connectivity())
+            sage: while not k:
+            ....:     g = digraphs.RandomDirectedGNP(11, .3)
+            ....:     k = Integer(g.edge_connectivity())
             sage: arborescences = g.edge_disjoint_spanning_trees(k)  # long time (up to 15s on sage.math, 2011)
             sage: all(a.is_directed_acyclic() for a in arborescences)  # long time
             True
@@ -6341,7 +6344,9 @@ def edge_disjoint_spanning_trees(self, k, root=None, solver=None, verbose=0):
 
         In the undirected case, we can only ensure half of it::
 
-            sage: g = graphs.RandomGNP(30, .3)
+            sage: g = graphs.RandomGNP(14, .3)  # reduced from 30 to 14, see #32169
+            sage: while not g.is_biconnected():
+            ....:     g = graphs.RandomGNP(14, .3)
             sage: k = Integer(g.edge_connectivity()) // 2
             sage: trees = g.edge_disjoint_spanning_trees(k)
             sage: all(t.is_tree() for t in trees)

src/sage/graphs/tutte_polynomial.py

@@ -544,10 +544,12 @@ def tutte_polynomial(G, edge_selector=None, cache=None):
         + 105*x^2*y^2 + 65*x*y^3 + 35*y^4 + 180*x^3 + 240*x^2*y + 171*x*y^2
         + 75*y^3 + 120*x^2 + 168*x*y + 84*y^2 + 36*x + 36*y
 
-    The Tutte polynomial of `G` evaluated at (1,1) is the number of
+    The Tutte polynomial of a connected graph `G` evaluated at (1,1) is the number of
     spanning trees of `G`::
 
         sage: G = graphs.RandomGNP(10,0.6)
+        sage: while not G.is_connected():
+        ....:     G = graphs.RandomGNP(10,0.6)
         sage: G.tutte_polynomial()(1,1) == G.spanning_trees_count()
         True
 

src/sage/groups/perm_gps/permgroup_morphism.py

@@ -117,11 +117,10 @@ def image(self, J):
             sage: G = L.galois_group()
             sage: D4 = DihedralGroup(4)
             sage: h = D4.isomorphism_to(G)
-            sage: h.image(D4)
-            Subgroup generated by [(1,2)(3,4)(5,7)(6,8), (1,6,4,7)(2,5,3,8)] of (Galois group 8T4 ([4]2) with order 8 of x^8 + 4*x^7 + 12*x^6 + 22*x^5 + 23*x^4 + 14*x^3 + 28*x^2 + 24*x + 16)
-            sage: r, s = D4.gens()
-            sage: h.image(r)
-            (1,6,4,7)(2,5,3,8)
+            sage: h.image(D4).is_isomorphic(G)
+            True
+            sage: all(h.image(g) in G for g in D4.gens())
+            True
         """
         H = self.codomain()
         if J in self.domain():

src/sage/matrix/matrix0.pyx

@@ -4830,7 +4830,10 @@ cdef class Matrix(sage.structure.element.Matrix):
             sage: B.multiplicative_order()
             1
 
-            sage: E = MatrixSpace(GF(11^2,'e'),5).random_element()
+            sage: M = MatrixSpace(GF(11^2,'e'),5)
+            sage: E = M.random_element()
+            sage: while E.det() == 0:
+            ....:     E = M.random_element()
             sage: (E^E.multiplicative_order()).is_one()
             True
 

src/sage/misc/misc.py

@@ -803,6 +803,8 @@ def __rmul__(self, left):
         EXAMPLES::
 
             sage: A = random_matrix(ZZ, 4)
+            sage: while A.rank() != 4:
+            ....:     A = random_matrix(ZZ, 4)
             sage: B = random_matrix(ZZ, 4)
             sage: temp = A * BackslashOperator()
             sage: temp.left is A

src/sage/modular/modform/numerical.py

@@ -444,7 +444,7 @@ def systems_of_eigenvalues(self, bound):
 
         EXAMPLES::
 
-            sage: numerical_eigenforms(61).systems_of_eigenvalues(10)  # rel tol 1e-12
+            sage: numerical_eigenforms(61).systems_of_eigenvalues(10)  # rel tol 1e-11
             [
             [-1.4811943040920152, 0.8060634335253695, 3.1563251746586642, 0.6751308705666477],
             [-1.0, -2.0000000000000027, -3.000000000000003, 1.0000000000000044],
@@ -471,7 +471,7 @@ def systems_of_abs(self, bound):
 
         EXAMPLES::
 
-            sage: numerical_eigenforms(61).systems_of_abs(10)  # rel tol 1e-12
+            sage: numerical_eigenforms(61).systems_of_abs(10)  # rel tol 1e-11
             [
             [0.3111078174659775, 2.903211925911551, 2.525427560843529, 3.214319743377552],
             [1.0, 2.0000000000000027, 3.000000000000003, 1.0000000000000044],

src/sage/stats/hmm/chmm.pyx

@@ -903,7 +903,7 @@ cdef class GaussianHiddenMarkovModel(HiddenMarkovModel):
             sage: m = hmm.GaussianHiddenMarkovModel([[.1,.9],[.5,.5]], [(1,.5), (-1,3)], [.1,.9])
             sage: v = m.sample(10)
             sage: l = stats.TimeSeries([m.baum_welch(v,max_iter=1)[0] for _ in range(len(v))])
-            sage: all(l[i] <= l[i+1] for i in range(9))
+            sage: all(l[i] <= l[i+1] + 0.0001 for i in range(9))
             True
             sage: l  # random
             [-20.1167, -17.7611, -16.9814, -16.9364, -16.9314, -16.9309, -16.9309, -16.9309, -16.9309, -16.9309]

src/sage/tests/book_stein_ent.py

@@ -267,8 +267,9 @@
 ....:         if g != 1 and g != n:
 ....:             return g
 sage: n=32295194023343; e=29468811804857; d=11127763319273
-sage: crack_given_decrypt(n, e*d - 1)
-737531
+sage: p = crack_given_decrypt(n, e*d - 1)
+sage: p in (737531, n/737531)  # could be other prime divisor
+True
 sage: factor(n)
 737531 * 43788253
 sage: e = 22601762315966221465875845336488389513

src/sage/tests/books/computational-mathematics-with-sagemath/linalg_doctest.py

@@ -269,6 +269,8 @@
   sage: A = random_matrix(R,2,3); A                    # random
   [      3*x^2 + x       x^2 + 2*x       2*x^2 + 2]
   [    x^2 + x + 2 2*x^2 + 4*x + 3   x^2 + 4*x + 3]
+  sage: while A.rank() < 2:
+  ....:     A = random_matrix(R,2,3)
 
 Sage example in ./linalg.tex, line 1830::
 

src/sage/tests/books/computational-mathematics-with-sagemath/sol/graphique_doctest.py

@@ -42,7 +42,7 @@
   sage: g = plot([c*e^(-1/x) for c in srange(-8, 8, 0.4)], (x, -3, 3))
   sage: y = var('y')
   sage: g += plot_vector_field((x^2, y), (x,-3,3), (y,-5,5))
-  sage: g.show()
+  sage: g.show()  # not tested, known bug, see :trac:`32657`
 
 Sage example in ./sol/graphique.tex, line 124::