patch from Charles Karney:

> Here is the patch which repeats the chi-squared tests in the case of
> suspiciously high values.  The probability of a false positive is now
> 1/10^9.  I also got rid of the mpfr_printf's.


git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@9121 280ebfd0-de03-0410-8827-d642c229c3f4

4 files changed, 287 insertions, 164 deletions

diff --git a/tests/chi-squared-tests.txt b/tests/chi-squared-tests.txt
index 8832473c6..fa64ed9fb 100644
--- a/tests/chi-squared-tests.txt
+++ b/tests/chi-squared-tests.txt
@@ -1,61 +1,94 @@
-[About the new t[ne]random_chisq.c tests by Charles Karney]
+About the trandom_deviate.c and t[ne]random_chisq.c
 
-Here is a patch which:
+trandom_deviate tests the functions in random_deviate.c, which provides
+support routines for mpfr_[ne]random.  The tests
 
-* Expands the tests carried out by trandom_deviate.c beyond the
-  simplistic goal of improving the code coverage.  Specifically, it now
-  checks
+  * improve the code coverage of the tests by constructing random
+    deviates with special properties: (a) the initial part of the
+    fraction matches that of another random deviate; (b) with the
+    leading bit in various positions.
 
   * that the values provided by mpfr_random_deviate_value at a given
-    precision does not change if additional random words are added to
-    the deviate;
+    precision do not change if additional random words are added to the
+    deviate;
 
   * that the values provided by mpfr_random_deviate_value at a given
     precision but different rounding modes are consistent.
 
-* Provides chi-squared tests for mpfr_[ne]random.  Two type of tests are
-  performed
+t[ne]random_chisq.c provide chi-squared tests for mpfr_[ne]random.  Note
+that there's a fair amount of code duplicated in these files, but with
+enough differences to make putting the code into a separate file a
+little messy.
+
+Two types of tests are performed
 
   * using equal width bins at moderate precision (prec = 64) treating
     the deviates are continuous quantities;
 
-  * assigning a bin to each possible random deviates (within a range) at
-    low precisions (prec = 2, 3, 4).
-
-  By default, these tests are run with 100000 samples.  This value is
-  high enough to catch most glaring problems with the algorithms and is
-  small enough that program completes in a under a second.  More subtle
-  errors in the implementation would be uncovered by running the test
-  with an argument of 10^6 or more.
-
-  These programs abort with status 1 if any of the chi-squared values
-  are too large or too small.  But by the nature of the tests, extreme
-  values are expected occasionally, and I've set the the "too large" and
-  "too small" thresholds to be hit with probability 1/1000.  Since there
-  are 8 such tests in each program (the low+high test on 4 separate
-  chi-squared runs), each test will fail about 8 times out of 1000.  I'm
-  not sure how these tests are best integrated into MPFR's testing
-  framework.  If the automatic tests are run with fixed seeds then it
-  will be easy to ensure that they don't fail.
-
-  The chief virtue of the chi-squared test will be to catch nearly all
-  mistakes made by a maintainer trying to improve the implementations in
-  some way.  In my experience, any such mistake leads to enormous
-  chi-squared values.
-
-  It might also be sensible to eliminate the tests for low values of
-  chi-squared.  Low values are only possible if the deviates are too
-  close to ideally distributed and that is only possible if the deviates
-  are correlated in some way.  This in turn is only possible if gmp's
-  random numbers are correlated and I assume that you have evidence (and
-  tests) that this is not so.
-
-  I've left t[ne]random as is, even though t[ne]random_chisq are
-  performing the same tests in a more comprehensive fashion.  There's
-  some utility in having a set of really simple tests.
-
-At this point, I think you have a comprehensive set of tests for
-mpre_[ne]random.  So I don't have any plans to add more tests.  Let me
-know if have any questions.
-
-  --Charles
+  * assigning a bin to each possible value for the random deviates
+    (within a range) at low precisions (prec = 2, 3, 4).
+
+Only abnormally large values of chi-squared are flagged, since
+abnormally low values can only result from
+
+  * the expected variation in chi-squared results;
+
+  * a correlation between the random deviates; but this can only happen
+    if gmp's random numbers are correlated (in a rather pernicious way);
+    the tests on gmp's random number generator should catch this.
+
+The chi-squared tests have to accomplish contradictory goals:
+
+  * reliably detect when there's a problem with the routines under test
+    (the probability of a false negative is low);
+
+  * rarely report a problem when there isn't one (the probability of a
+    false positive is low);
+
+  * complete in about 1 second.
+
+Some attributes of the testing come to our aid:
+
+  * the tests are executed many times on many different platforms;
+
+  * typical coding errors in the routines for normal and exponential
+    sampling lead to distributions which are far from the correct one;
+
+  * carrying out a longer (by a factor of N) test on a bad routine
+    multiplies chi-squared by about N, a result which has a vanishingly
+    small probability for a good routine.
+
+The strategy for testing is
+
+  * use 100000 samples (with this number, the test of mpfr_nrandom takes
+    about 0.3 sec); compute Q, the probability that chi-square exceeds
+    the value obtained.  If Q > 0.01, the test passes; if Q < 1.0e-9,
+    the test fails; otherwise:
+
+  * repeat the test with 10 times the number of samples and compute Q
+    again.  If Q > 0.001, the test passes; if Q < 1.0e-9, the test
+    fails; otherwise:
+
+  * repeat the test with 10 times the number of samples and compute Q
+    again.  If Q > 0.0001, the test passes; if Q < 1.0e-9, the test
+    fails; otherwise:
+
+  * the test fails.
+
+If the mpfr_[ne]random is OK, then the probability that the test fails
+(a false positive) is about 1.0e-9.  If there's a bad error in
+mpfr_[ne]random, failure will probably occur at one of the stages.  It
+there's an error in mpfr_[ne]random which results in a slight deviation
+from the true distribution (which I don't think is likely given the
+nature of the algorithms), then there a reasonable chance that the
+problem will be discovered "eventually".  Because "eventually" is not
+very satisfactory, it is important that the tests be run with a large
+number (e.g., 10^9) of samples whenever changes are made to
+mpfr_[ne]random.
+
+I've left the t[ne]random tests as is, even though t[ne]random_chisq are
+performing the same tests in a more comprehensive fashion.  There's some
+utility in having a set of really simple tests.
+
+  --Charles Karney <charles.karney@sri.com>
+
diff --git a/tests/terandom_chisq.c b/tests/terandom_chisq.c
index 16bc1bb18..9fbcda44d 100644
--- a/tests/terandom_chisq.c
+++ b/tests/terandom_chisq.c
@@ -35,6 +35,39 @@ exponential_cumulative (mpfr_t z, mpfr_t x, mpfr_rnd_t rnd)
   mpfr_neg (z, z, rnd);
 }
 
+/* Given nu and chisqp, compute probability that chisq > chisqp.  This uses,
+ * A&S 26.4.16,
+ *
+ * Q(nu,chisqp) =
+ *     erfc( (3/2)*sqrt(nu) * ( cbrt(chisqp/nu) - 1 + 2/(9*nu) ) ) / 2
+ *
+ * which is valid for nu > 30.  This is the basis for the formula in Knuth,
+ * TAOCP, Vol 2, 3.3.1, Table 1.  It more accurate than the similar formula,
+ * DLMF 8.11.10. */
+static void
+chisq_prob (mpfr_t q, long nu, mpfr_t chisqp)
+{
+  mpfr_t t;
+  mpfr_rnd_t rnd;
+
+  rnd = MPFR_RNDN;  /* This uses an approx formula.  Might as well use RNDN. */
+  mpfr_init2 (t, mpfr_get_prec (q));
+
+  mpfr_div_si (q, chisqp, nu, rnd); /* chisqp/nu */
+  mpfr_cbrt (q, q, rnd);            /* (chisqp/nu)^(1/3) */
+  mpfr_sub_ui (q, q, 1, rnd);       /* (chisqp/nu)^(1/3) - 1 */
+  mpfr_set_ui (t, 2, rnd);
+  mpfr_div_si (t, t, 9*nu, rnd); /* 2/(9*nu) */
+  mpfr_add (q, q, t, rnd);       /* (chisqp/nu)^(1/3) - 1 + 2/(9*nu) */
+  mpfr_sqrt_ui (t, nu, rnd);     /* sqrt(nu) */
+  mpfr_mul_d (t, t, 1.5, rnd);   /* (3/2)*sqrt(nu) */
+  mpfr_mul (q, q, t, rnd);       /* arg to erfc */
+  mpfr_erfc (q, q, rnd);         /* erfc(...) */
+  mpfr_div_ui (q, q, 2, rnd);    /* erfc(...)/2 */
+
+  mpfr_clear (t);
+}
+
 /* The continuous chi-squared test on with a set of bins of equal width.
  *
  * A single precision is picked for sampling and the chi-squared calculation.
@@ -49,7 +82,7 @@ exponential_cumulative (mpfr_t z, mpfr_t x, mpfr_rnd_t rnd)
  *
  * The testing of low precision exponential deviates is done by
  * test_erandom_chisq_disc. */
-static void
+static double
 test_erandom_chisq_cont (long num, mpfr_prec_t prec, int nu,
                          double xmin, double xmax, int verbose)
 {
@@ -58,7 +91,7 @@ test_erandom_chisq_cont (long num, mpfr_prec_t prec, int nu,
   int i, inexact;
   long k;
   mpfr_rnd_t rnd, rndd;
-  double chilim, xp, sqrt2dof;
+  double Q, chisq;
 
   rnd = MPFR_RNDN;              /* For chi-squared calculation */
   rndd = MPFR_RNDD;             /* For sampling and figuring the bins */
@@ -84,7 +117,8 @@ test_erandom_chisq_cont (long num, mpfr_prec_t prec, int nu,
       if (inexact == 0)
         {
           /* one call in the loop pretended to return an exact number! */
-          printf ("Error: mpfr_erandom() returns a zero ternary value.\n");
+          fprintf (stderr,
+                   "Error: mpfr_erandom() returns a zero ternary value.\n");
           exit (1);
         }
       if (mpfr_signbit (x))
@@ -122,32 +156,19 @@ test_erandom_chisq_cont (long num, mpfr_prec_t prec, int nu,
       mpfr_swap (pa, pb);       /* i.e., pa = pb */
     }
 
+  chisq = mpfr_get_d (t, rnd);
+  chisq_prob (t, nu, t);
+  Q = mpfr_get_d (t, rnd);
   if (verbose)
-    mpfr_printf ("num = %ld, discrete (prec = %ld) bins in [%.2f, %.2f], "
-                 "nu = %d: chi2 = %.2Rf\n", num, prec, xmin, xmax, nu, t);
-
-  /* See Knuth, TAOCP, Vol 2, 3.3.1, Table 1, approx formula for nu > 30; xp =
-   * +/- 3.1 gives 0.1% and 99.9% percentile points. */
-  mpfr_set_si (pa, 2*nu, rnd);
-  mpfr_sqrt (pa, pa, rnd);
-  sqrt2dof = mpfr_get_d (pa, rnd);
-  xp = 3.1;
-  chilim = nu - sqrt2dof * xp  +   2 * xp * xp / 3;
-  if (mpfr_cmp_d (t, chilim) < 0)
-    {
-      mpfr_printf ("chi2 = %.2Rf is less than %.2f"
-                   " should only happen 1 in 1000 times\n", t, chilim);
-      exit (1);
-    }
-  chilim = nu + sqrt2dof * xp  +   2 * xp * xp / 3;
-  if (mpfr_cmp_d (t, chilim) > 0)
     {
-      mpfr_printf ("chi2 = %.2Rf is greater than %.2f"
-                   " should only happen 1 in 1000 times\n", t, chilim);
-      exit (1);
+      printf ("num = %ld, equal bins in [%.2f, %.2f], nu = %d: chisq = %.2f\n",
+              num, xmin, xmax, nu, chisq);
+      if (Q < 0.05)
+        printf ("    WARNING: probability (less than 5%%) = %.2e\n", Q);
     }
 
   mpfr_clears (x, a, b, dx, z, pa, pb, ps, t, (mpfr_ptr) 0);
+  return Q;
 }
 
 /* Return a sequential number for a positive low-precision x.  x is altered by
@@ -185,8 +206,8 @@ sequential (mpfr_t x)
  * The sampling is with MPFR_RNDN.  This is the rounding mode which elicits the
  * most information.  trandom_deviate includes checks on the consistency of the
  * results extracted from a random_deviate with other rounding modes.  */
-static void
-test_erandom_chisq_disc (long num, mpfr_prec_t wprec, mpfr_prec_t prec,
+static double
+test_erandom_chisq_disc (long num, mpfr_prec_t wprec, int prec,
                          double xmin, double xmax, int verbose)
 {
   mpfr_t x, v, pa, pb, z, t;
@@ -194,7 +215,7 @@ test_erandom_chisq_disc (long num, mpfr_prec_t wprec, mpfr_prec_t prec,
   int i, inexact, nu;
   long *counts;
   long k, seqmin, seqmax, seq;
-  double chilim, xp, sqrt2dof;
+  double Q, chisq;
 
   rnd = MPFR_RNDN;
   mpfr_init2 (x, prec);
@@ -261,33 +282,54 @@ test_erandom_chisq_disc (long num, mpfr_prec_t wprec, mpfr_prec_t prec,
       mpfr_add (t, t, z, rnd);
       mpfr_swap (pa, pb);       /* i.e., pa = pb */
     }
+
+  chisq = mpfr_get_d (t, rnd);
+  chisq_prob (t, nu, t);
+  Q = mpfr_get_d (t, rnd);
   if (verbose)
-    mpfr_printf ("num = %ld, discrete (prec = %ld) bins in [%.6f, %.2f], "
-                 "nu = %d: chi2 = %.2Rf\n", num, prec, xmin, xmax, nu, t);
-
-  /* See Knuth, TAOCP, Vol 2, 3.3.1, Table 1, approx formula for nu > 30; xp =
-   * +/- 3.1 gives 0.1% and 99.9% percentile points. */
-  mpfr_set_si (pa, 2*nu, rnd);
-  mpfr_sqrt (pa, pa, rnd);
-  sqrt2dof = mpfr_get_d (pa, rnd);
-  xp = 3.1;
-  chilim = nu - sqrt2dof * xp  +   2 * xp * xp / 3;
-  if (mpfr_cmp_d (t, chilim) < 0)
     {
-      mpfr_printf ("chi2 = %.2Rf is less than %.2f"
-                   " should only happen 1 in 1000 times\n", t, chilim);
-      exit (1);
-    }
-  chilim = nu + sqrt2dof * xp  +   2 * xp * xp / 3;
-  if (mpfr_cmp_d (t, chilim) > 0)
-    {
-      mpfr_printf ("chi2 = %.2Rf is greater than %.2f"
-                   " should only happen 1 in 1000 times\n", t, chilim);
-      exit (1);
+      printf ("num = %ld, discrete (prec = %d) bins in [%.6f, %.2f], "
+              "nu = %d: chisq = %.2f\n", num, prec, xmin, xmax, nu, chisq);
+      if (Q < 0.05)
+        printf ("    WARNING: probability (less than 5%%) = %.2e\n", Q);
     }
 
   free (counts);
   mpfr_clears (x, v, pa, pb, z, t, (mpfr_ptr) 0);
+  return Q;
+}
+
+static void
+run_chisq ( double (*f)(long, mpfr_prec_t, int, double, double, int),
+            long num, mpfr_prec_t prec, int bin,
+            double xmin, double xmax, int verbose)
+{
+  double Q, Qcum, Qbad, Qthresh;
+  int i;
+  Qcum = 1;
+  Qbad = 1.e-9;
+  Qthresh = 0.01;
+  for (i = 0; i < 3; ++i)
+    {
+      Q = (*f)(num, prec, bin, xmin, xmax, verbose);
+      Qcum *= Q;
+      if (Q > Qthresh)
+        return;
+      else if (Q < Qbad)
+        {
+          fprintf (stderr, "Error: mpfr_erandom chi-squared failure "
+                   "(prob = %.2e)\n", Q);
+          exit (1);
+        }
+      num *= 10;
+      Qthresh /= 10;
+    }
+  if (Qcum < Qbad)              /* Presumably this is true */
+    {
+      fprintf (stderr, "Error: mpfr_erandom combined chi-squared failure "
+               "(prob = %.2e)\n", Qcum);
+      exit (1);
+    }
 }
 
 int
@@ -307,10 +349,10 @@ main (int argc, char *argv[])
         nbtests = a;
     }
 
-  test_erandom_chisq_cont (nbtests, 64, 60, 0, 7, verbose);
-  test_erandom_chisq_disc (nbtests, 64, 2, 0.002, 6, verbose);
-  test_erandom_chisq_disc (nbtests, 64, 3, 0.02, 7, verbose);
-  test_erandom_chisq_disc (nbtests, 64, 4, 0.04, 8, verbose);
+  run_chisq (test_erandom_chisq_cont, nbtests, 64, 60, 0, 7, verbose);
+  run_chisq (test_erandom_chisq_disc, nbtests, 64, 2, 0.002, 6, verbose);
+  run_chisq (test_erandom_chisq_disc, nbtests, 64, 3, 0.02, 7, verbose);
+  run_chisq (test_erandom_chisq_disc, nbtests, 64, 4, 0.04, 8, verbose);
 
   tests_end_mpfr ();
   return 0;
diff --git a/tests/tnrandom_chisq.c b/tests/tnrandom_chisq.c
index 89740761f..49c02e40f 100644
--- a/tests/tnrandom_chisq.c
+++ b/tests/tnrandom_chisq.c
@@ -36,6 +36,39 @@ normal_cumulative (mpfr_t z, mpfr_t x, mpfr_rnd_t rnd)
   mpfr_div_ui (z, z, 2, rnd);
 }
 
+/* Given nu and chisqp, compute probability that chisq > chisqp.  This uses,
+ * A&S 26.4.16,
+ *
+ * Q(nu,chisqp) =
+ *     erfc( (3/2)*sqrt(nu) * ( cbrt(chisqp/nu) - 1 + 2/(9*nu) ) ) / 2
+ *
+ * which is valid for nu > 30.  This is the basis for the formula in Knuth,
+ * TAOCP, Vol 2, 3.3.1, Table 1.  It more accurate than the similar formula,
+ * DLMF 8.11.10. */
+static void
+chisq_prob (mpfr_t q, long nu, mpfr_t chisqp)
+{
+  mpfr_t t;
+  mpfr_rnd_t rnd;
+
+  rnd = MPFR_RNDN;  /* This uses an approx formula.  Might as well use RNDN. */
+  mpfr_init2 (t, mpfr_get_prec (q));
+
+  mpfr_div_si (q, chisqp, nu, rnd); /* chisqp/nu */
+  mpfr_cbrt (q, q, rnd);            /* (chisqp/nu)^(1/3) */
+  mpfr_sub_ui (q, q, 1, rnd);       /* (chisqp/nu)^(1/3) - 1 */
+  mpfr_set_ui (t, 2, rnd);
+  mpfr_div_si (t, t, 9*nu, rnd); /* 2/(9*nu) */
+  mpfr_add (q, q, t, rnd);       /* (chisqp/nu)^(1/3) - 1 + 2/(9*nu) */
+  mpfr_sqrt_ui (t, nu, rnd);     /* sqrt(nu) */
+  mpfr_mul_d (t, t, 1.5, rnd);   /* (3/2)*sqrt(nu) */
+  mpfr_mul (q, q, t, rnd);       /* arg to erfc */
+  mpfr_erfc (q, q, rnd);         /* erfc(...) */
+  mpfr_div_ui (q, q, 2, rnd);    /* erfc(...)/2 */
+
+  mpfr_clear (t);
+}
+
 /* The continuous chi-squared test on with a set of bins of equal width.
  *
  * A single precision is picked for sampling and the chi-squared calculation.
@@ -50,7 +83,7 @@ normal_cumulative (mpfr_t z, mpfr_t x, mpfr_rnd_t rnd)
  *
  * The testing of low precision normal deviates is done by
  * test_nrandom_chisq_disc. */
-static void
+static double
 test_nrandom_chisq_cont (long num, mpfr_prec_t prec, int nu,
                          double xmin, double xmax, int verbose)
 {
@@ -59,7 +92,7 @@ test_nrandom_chisq_cont (long num, mpfr_prec_t prec, int nu,
   int i, inexact;
   long k;
   mpfr_rnd_t rnd, rndd;
-  double chilim, xp, sqrt2dof;
+  double Q, chisq;
 
   rnd = MPFR_RNDN;              /* For chi-squared calculation */
   rndd = MPFR_RNDD;             /* For sampling and figuring the bins */
@@ -85,7 +118,8 @@ test_nrandom_chisq_cont (long num, mpfr_prec_t prec, int nu,
       if (inexact == 0)
         {
           /* one call in the loop pretended to return an exact number! */
-          printf ("Error: mpfr_nrandom() returns a zero ternary value.\n");
+          fprintf (stderr,
+                   "Error: mpfr_nrandom() returns a zero ternary value.\n");
           exit (1);
         }
       mpfr_sub (x, x, a, rndd);
@@ -118,32 +152,19 @@ test_nrandom_chisq_cont (long num, mpfr_prec_t prec, int nu,
       mpfr_swap (pa, pb);       /* i.e., pa = pb */
     }
 
+  chisq = mpfr_get_d (t, rnd);
+  chisq_prob (t, nu, t);
+  Q = mpfr_get_d (t, rnd);
   if (verbose)
-    mpfr_printf ("num = %ld, discrete (prec = %ld) bins in [%.2f, %.2f], "
-                 "nu = %d: chi2 = %.2Rf\n", num, prec, xmin, xmax, nu, t);
-
-  /* See Knuth, TAOCP, Vol 2, 3.3.1, Table 1, approx formula for nu > 30; xp =
-   * +/- 3.1 gives 0.1% and 99.9% percentile points. */
-  mpfr_set_si (pa, 2*nu, rnd);
-  mpfr_sqrt (pa, pa, rnd);
-  sqrt2dof = mpfr_get_d (pa, rnd);
-  xp = 3.1;
-  chilim = nu - sqrt2dof * xp  +   2 * xp * xp / 3;
-  if (mpfr_cmp_d (t, chilim) < 0)
-    {
-      mpfr_printf ("chi2 = %.2Rf is less than %.2f"
-                   " should only happen 1 in 1000 times\n", t, chilim);
-      exit (1);
-    }
-  chilim = nu + sqrt2dof * xp  +   2 * xp * xp / 3;
-  if (mpfr_cmp_d (t, chilim) > 0)
     {
-      mpfr_printf ("chi2 = %.2Rf is greater than %.2f"
-                   " should only happen 1 in 1000 times\n", t, chilim);
-      exit (1);
+      printf ("num = %ld, equal bins in [%.2f, %.2f], nu = %d: chisq = %.2f\n",
+              num, xmin, xmax, nu, chisq);
+      if (Q < 0.05)
+        printf ("    WARNING: probability (less than 5%%) = %.2e\n", Q);
     }
 
   mpfr_clears (x, a, b, dx, z, pa, pb, ps, t, (mpfr_ptr) 0);
+  return Q;
 }
 
 /* Return a sequential number for a positive low-precision x.  x is altered by
@@ -184,8 +205,8 @@ sequential (mpfr_t x)
  * MPFR_RNDN.  This is the rounding mode which elicits the most information.
  * trandom_deviate includes checks on the consistency of the results extracted
  * from a random_deviate with other rounding modes.  */
-static void
-test_nrandom_chisq_disc (long num, mpfr_prec_t wprec, mpfr_prec_t prec,
+static double
+test_nrandom_chisq_disc (long num, mpfr_prec_t wprec, int prec,
                          double xmin, double xmax, int verbose)
 {
   mpfr_t x, v, pa, pb, z, t;
@@ -193,7 +214,7 @@ test_nrandom_chisq_disc (long num, mpfr_prec_t wprec, mpfr_prec_t prec,
   int i, inexact, nu;
   long *counts;
   long k, seqmin, seqmax, seq;
-  double chilim, xp, sqrt2dof;
+  double Q, chisq;
 
   rnd = MPFR_RNDN;
   mpfr_init2 (x, prec);
@@ -262,33 +283,54 @@ test_nrandom_chisq_disc (long num, mpfr_prec_t wprec, mpfr_prec_t prec,
       mpfr_add (t, t, z, rnd);
       mpfr_swap (pa, pb);       /* i.e., pa = pb */
     }
+
+  chisq = mpfr_get_d (t, rnd);
+  chisq_prob (t, nu, t);
+  Q = mpfr_get_d (t, rnd);
   if (verbose)
-    mpfr_printf ("num = %ld, discrete (prec = %ld) bins in [%.6f, %.2f], "
-                 "nu = %d: chi2 = %.2Rf\n", num, prec, xmin, xmax, nu, t);
-
-  /* See Knuth, TAOCP, Vol 2, 3.3.1, Table 1, approx formula for nu > 30; xp =
-   * +/- 3.1 gives 0.1% and 99.9% percentile points. */
-  mpfr_set_si (pa, 2*nu, rnd);
-  mpfr_sqrt (pa, pa, rnd);
-  sqrt2dof = mpfr_get_d (pa, rnd);
-  xp = 3.1;
-  chilim = nu - sqrt2dof * xp  +   2 * xp * xp / 3;
-  if (mpfr_cmp_d (t, chilim) < 0)
     {
-      mpfr_printf ("chi2 = %.2Rf is less than %.2f"
-                   " should only happen 1 in 1000 times\n", t, chilim);
-      exit (1);
-    }
-  chilim = nu + sqrt2dof * xp  +   2 * xp * xp / 3;
-  if (mpfr_cmp_d (t, chilim) > 0)
-    {
-      mpfr_printf ("chi2 = %.2Rf is greater than %.2f"
-                   " should only happen 1 in 1000 times\n", t, chilim);
-      exit (1);
+      printf ("num = %ld, discrete (prec = %d) bins in [%.6f, %.2f], "
+              "nu = %d: chisq = %.2f\n", num, prec, xmin, xmax, nu, chisq);
+      if (Q < 0.05)
+        printf ("    WARNING: probability (less than 5%%) = %.2e\n", Q);
     }
 
   free (counts);
   mpfr_clears (x, v, pa, pb, z, t, (mpfr_ptr) 0);
+  return Q;
+}
+
+static void
+run_chisq ( double (*f)(long, mpfr_prec_t, int, double, double, int),
+            long num, mpfr_prec_t prec, int bin,
+            double xmin, double xmax, int verbose)
+{
+  double Q, Qcum, Qbad, Qthresh;
+  int i;
+  Qcum = 1;
+  Qbad = 1.e-9;
+  Qthresh = 0.01;
+  for (i = 0; i < 3; ++i)
+    {
+      Q = (*f)(num, prec, bin, xmin, xmax, verbose);
+      Qcum *= Q;
+      if (Q > Qthresh)
+        return;
+      else if (Q < Qbad)
+        {
+          fprintf (stderr, "Error: mpfr_nrandom chi-squared failure "
+                   "(prob = %.2e)\n", Q);
+          exit (1);
+        }
+      num *= 10;
+      Qthresh /= 10;
+    }
+  if (Qcum < Qbad)              /* Presumably this is true */
+    {
+      fprintf (stderr, "Error: mpfr_nrandom combined chi-squared failure "
+               "(prob = %.2e)\n", Qcum);
+      exit (1);
+    }
 }
 
 int
@@ -308,10 +350,10 @@ main (int argc, char *argv[])
         nbtests = a;
     }
 
-  test_nrandom_chisq_cont (nbtests, 64, 60, -4, 4, verbose);
-  test_nrandom_chisq_disc (nbtests, 64, 2, 0.0005, 3, verbose);
-  test_nrandom_chisq_disc (nbtests, 64, 3, 0.002, 4, verbose);
-  test_nrandom_chisq_disc (nbtests, 64, 4, 0.004, 4, verbose);
+  run_chisq (test_nrandom_chisq_cont, nbtests, 64, 60, -4, 4, verbose);
+  run_chisq (test_nrandom_chisq_disc, nbtests, 64, 2, 0.0005, 3, verbose);
+  run_chisq (test_nrandom_chisq_disc, nbtests, 64, 3, 0.002, 4, verbose);
+  run_chisq (test_nrandom_chisq_disc, nbtests, 64, 4, 0.004, 4, verbose);
 
   tests_end_mpfr ();
   return 0;
diff --git a/tests/trandom_deviate.c b/tests/trandom_deviate.c
index 04647adaa..668b53fa9 100644
--- a/tests/trandom_deviate.c
+++ b/tests/trandom_deviate.c
@@ -73,7 +73,8 @@ test_compare (long nbtests, int verbose)
           t2 = mpfr_random_deviate_less (u, v, RANDS);
           if (t1 != t2)
             {
-              printf ("Error: mpfr_random_deviate_less() inconsistent.\n");
+              fprintf (stderr,
+                       "Error: mpfr_random_deviate_less() inconsistent.\n");
               exit (1);
             }
           if (t1)
@@ -84,13 +85,14 @@ test_compare (long nbtests, int verbose)
       t1 = mpfr_random_deviate_less (u, u, RANDS);
       if (t1)
         {
-          printf ("Error: mpfr_random_deviate_less() gives u < u.\n");
+          fprintf (stderr, "Error: mpfr_random_deviate_less() gives u < u.\n");
           exit (1);
         }
       t1 = mpfr_random_deviate_tstbit (u, 0, RANDS);
       if (t1)
         {
-          printf ("Error: mpfr_random_deviate_tstbit() says 1 bit is on.\n");
+          fprintf (stderr,
+                   "Error: mpfr_random_deviate_tstbit() says 1 bit is on.\n");
           exit (1);
         }
     }
@@ -150,7 +152,7 @@ test_value_consistency (long nbtests)
               inexa2 == inexb2 && mpfr_equal_p (a2, b2) &&
               inexa3 == inexb3 && mpfr_equal_p (a3, b3) ) )
         {
-          printf ("Error: random_deviate values are inconsistent.\n");
+          fprintf (stderr, "Error: random_deviate values are inconsistent.\n");
           exit (1);
         }
     }
@@ -190,9 +192,10 @@ test_value_round (long nbtests, mpfr_prec_t prec)
               inexd     < 0 && inexu     > 0 &&
               inexz * s < 0 && inexa * s > 0 ) )
         {
-          printf ("n %d t %d d %d u %d z %d a %d s %d\n",
+          fprintf (stderr, "n %d t %d d %d u %d z %d a %d s %d\n",
                   inexn, inext, inexd, inexu, inexz, inexa, s);
-          printf ("Error: random_deviate has wrong values for inexact.\n");
+          fprintf (stderr,
+                   "Error: random_deviate has wrong values for inexact.\n");
           exit (1);
         }
       if (inexn < 0)
@@ -206,11 +209,13 @@ test_value_round (long nbtests, mpfr_prec_t prec)
               mpfr_equal_p(xz, SAME_SIGN(inexn, inexz) ? xn : t) &&
               mpfr_equal_p(xa, SAME_SIGN(inexn, inexa) ? xn : t) ) )
         {
-          printf ("n %d t %d d %d u %d z %d a %d s %d\n",
+          fprintf (stderr, "n %d t %d d %d u %d z %d a %d s %d\n",
                   inexn, inext, inexd, inexu, inexz, inexa, s);
-          mpfr_printf ("n %.4Rf t %.4Rf d %.4Rf u %.4Rf z %.4Rf a %.4Rf\n",
-                       xn, t, xd, xu, xz, xa);
-          printf ("Error: random_deviate rounds inconsistently.\n");
+          fprintf (stderr, "n %.4Rf t %.4Rf d %.4Rf u %.4Rf z %.4Rf a %.4Rf\n",
+                   mpfr_get_d (xn, MPFR_RNDN), mpfr_get_d (t, MPFR_RNDN),
+                   mpfr_get_d (xd, MPFR_RNDN), mpfr_get_d (xu, MPFR_RNDN),
+                   mpfr_get_d (xz, MPFR_RNDN), mpfr_get_d (xa, MPFR_RNDN));
+          fprintf (stderr, "Error: random_deviate rounds inconsistently.\n");
           exit (1);
         }
     }
@@ -267,7 +272,8 @@ test_value (long nbtests, mpfr_prec_t prec, mpfr_rnd_t rnd,
         }
       if (exact)
         {
-          printf ("Error: random_deviate() returns a zero ternary value.\n");
+          fprintf (stderr,
+                   "Error: random_deviate() returns a zero ternary value.\n");
           exit (1);
         }
       mpfr_random_deviate_reset (u);