summaryrefslogtreecommitdiff
path: root/gsl-1.9/rng/test.c
diff options
context:
space:
mode:
Diffstat (limited to 'gsl-1.9/rng/test.c')
-rw-r--r--gsl-1.9/rng/test.c602
1 files changed, 602 insertions, 0 deletions
diff --git a/gsl-1.9/rng/test.c b/gsl-1.9/rng/test.c
new file mode 100644
index 0000000..8b267a2
--- /dev/null
+++ b/gsl-1.9/rng/test.c
@@ -0,0 +1,602 @@
+/* rng/test.c
+ *
+ * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or (at
+ * your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#include <config.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_test.h>
+#include <gsl/gsl_ieee_utils.h>
+
+void rng_test (const gsl_rng_type * T, unsigned long int seed, unsigned int n,
+ unsigned long int result);
+void rng_float_test (const gsl_rng_type * T);
+void generic_rng_test (const gsl_rng_type * T);
+void rng_state_test (const gsl_rng_type * T);
+void rng_parallel_state_test (const gsl_rng_type * T);
+void rng_read_write_test (const gsl_rng_type * T);
+int rng_max_test (gsl_rng * r, unsigned long int *kmax, unsigned long int ran_max) ;
+int rng_min_test (gsl_rng * r, unsigned long int *kmin, unsigned long int ran_min, unsigned long int ran_max) ;
+int rng_sum_test (gsl_rng * r, double *sigma);
+int rng_bin_test (gsl_rng * r, double *sigma);
+
+#define N 10000
+#define N2 200000
+
+int
+main (void)
+{
+ const gsl_rng_type ** rngs = gsl_rng_types_setup(); /* get all rng types */
+
+ const gsl_rng_type ** r ;
+
+ gsl_ieee_env_setup ();
+
+ gsl_rng_env_setup ();
+
+ /* specific tests of known results for 10000 iterations with seed = 1 */
+
+ rng_test (gsl_rng_rand, 1, 10000, 1910041713);
+ rng_test (gsl_rng_randu, 1, 10000, 1623524161);
+ rng_test (gsl_rng_cmrg, 1, 10000, 719452880);
+ rng_test (gsl_rng_minstd, 1, 10000, 1043618065);
+ rng_test (gsl_rng_mrg, 1, 10000, 2064828650);
+ rng_test (gsl_rng_taus, 1, 10000, 2733957125UL);
+ rng_test (gsl_rng_taus2, 1, 10000, 2733957125UL);
+ rng_test (gsl_rng_taus113, 1, 1000, 1925420673UL);
+ rng_test (gsl_rng_transputer, 1, 10000, 1244127297UL);
+ rng_test (gsl_rng_vax, 1, 10000, 3051034865UL);
+
+ /* Borosh13 test value from PARI: (1812433253^10000)%(2^32) */
+ rng_test (gsl_rng_borosh13, 1, 10000, 2513433025UL);
+
+ /* Fishman18 test value from PARI: (62089911^10000)%(2^31-1) */
+ rng_test (gsl_rng_fishman18, 1, 10000, 330402013UL);
+
+ /* Fishman2x test value from PARI:
+ ((48271^10000)%(2^31-1) - (40692^10000)%(2^31-249))%(2^31-1) */
+ rng_test (gsl_rng_fishman2x, 1, 10000, 540133597UL);
+
+ /* Knuthran2 test value from PARI:
+ { xn1=1; xn2=1; for (n=1,10000,
+ xn = (271828183*xn1 - 314159269*xn2)%(2^31-1);
+ xn2=xn1; xn1=xn; print(xn); ) } */
+ rng_test (gsl_rng_knuthran2, 1, 10000, 1084477620UL);
+
+ /* Knuthran test value taken from p188 in Knuth Vol 2. 3rd Ed */
+ rng_test (gsl_rng_knuthran, 310952, 1009 * 2009 + 1, 461390032);
+
+ /* Knuthran improved test value from Knuth's source */
+ rng_test (gsl_rng_knuthran2002, 310952, 1, 708622036);
+ rng_test (gsl_rng_knuthran2002, 310952, 2, 1005450560);
+ rng_test (gsl_rng_knuthran2002, 310952, 100 * 2009 + 1, 995235265);
+ rng_test (gsl_rng_knuthran2002, 310952, 1009 * 2009 + 1, 704987132);
+
+ /* Lecuyer21 test value from PARI: (40692^10000)%(2^31-249) */
+ rng_test (gsl_rng_lecuyer21, 1, 10000, 2006618587UL);
+
+ /* Waterman14 test value from PARI: (1566083941^10000)%(2^32) */
+ rng_test (gsl_rng_waterman14, 1, 10000, 3776680385UL);
+
+ /* specific tests of known results for 10000 iterations with seed = 6 */
+
+ /* Coveyou test value from PARI:
+ x=6; for(n=1,10000,x=(x*(x+1))%(2^32);print(x);) */
+
+ rng_test (gsl_rng_coveyou, 6, 10000, 1416754246UL);
+
+ /* Fishman20 test value from PARI: (6*48271^10000)%(2^31-1) */
+ rng_test (gsl_rng_fishman20, 6, 10000, 248127575UL);
+
+ /* FIXME: the ranlux tests below were made by running the fortran code and
+ getting the expected value from that. An analytic calculation
+ would be preferable. */
+
+ rng_test (gsl_rng_ranlux, 314159265, 10000, 12077992);
+ rng_test (gsl_rng_ranlux389, 314159265, 10000, 165942);
+
+ rng_test (gsl_rng_ranlxs0, 1, 10000, 11904320);
+ /* 0.709552764892578125 * ldexp(1.0,24) */
+
+ rng_test (gsl_rng_ranlxs1, 1, 10000, 8734328);
+ /* 0.520606517791748047 * ldexp(1.0,24) */
+
+ rng_test (gsl_rng_ranlxs2, 1, 10000, 6843140);
+ /* 0.407882928848266602 * ldexp(1.0,24) */
+
+ rng_test (gsl_rng_ranlxd1, 1, 10000, 1998227290UL);
+ /* 0.465248546261094020 * ldexp(1.0,32) */
+
+ rng_test (gsl_rng_ranlxd2, 1, 10000, 3949287736UL);
+ /* 0.919515205581550532 * ldexp(1.0,32) */
+
+ /* FIXME: the tests below were made by running the original code in
+ the ../random directory and getting the expected value from
+ that. An analytic calculation would be preferable. */
+
+ rng_test (gsl_rng_slatec, 1, 10000, 45776);
+ rng_test (gsl_rng_uni, 1, 10000, 9214);
+ rng_test (gsl_rng_uni32, 1, 10000, 1155229825);
+ rng_test (gsl_rng_zuf, 1, 10000, 3970);
+
+ /* The tests below were made by running the original code and
+ getting the expected value from that. An analytic calculation
+ would be preferable. */
+
+ rng_test (gsl_rng_r250, 1, 10000, 1100653588);
+ rng_test (gsl_rng_mt19937, 4357, 1000, 1186927261);
+ rng_test (gsl_rng_mt19937_1999, 4357, 1000, 1030650439);
+ rng_test (gsl_rng_mt19937_1998, 4357, 1000, 1309179303);
+ rng_test (gsl_rng_tt800, 0, 10000, 2856609219UL);
+
+ rng_test (gsl_rng_ran0, 0, 10000, 1115320064);
+ rng_test (gsl_rng_ran1, 0, 10000, 1491066076);
+ rng_test (gsl_rng_ran2, 0, 10000, 1701364455);
+ rng_test (gsl_rng_ran3, 0, 10000, 186340785);
+
+ rng_test (gsl_rng_ranmar, 1, 10000, 14428370);
+
+ rng_test (gsl_rng_rand48, 0, 10000, 0xDE095043UL);
+ rng_test (gsl_rng_rand48, 1, 10000, 0xEDA54977UL);
+
+ rng_test (gsl_rng_random_glibc2, 0, 10000, 1908609430);
+ rng_test (gsl_rng_random8_glibc2, 0, 10000, 1910041713);
+ rng_test (gsl_rng_random32_glibc2, 0, 10000, 1587395585);
+ rng_test (gsl_rng_random64_glibc2, 0, 10000, 52848624);
+ rng_test (gsl_rng_random128_glibc2, 0, 10000, 1908609430);
+ rng_test (gsl_rng_random256_glibc2, 0, 10000, 179943260);
+
+ rng_test (gsl_rng_random_bsd, 0, 10000, 1457025928);
+ rng_test (gsl_rng_random8_bsd, 0, 10000, 1910041713);
+ rng_test (gsl_rng_random32_bsd, 0, 10000, 1663114331);
+ rng_test (gsl_rng_random64_bsd, 0, 10000, 864469165);
+ rng_test (gsl_rng_random128_bsd, 0, 10000, 1457025928);
+ rng_test (gsl_rng_random256_bsd, 0, 10000, 1216357476);
+
+ rng_test (gsl_rng_random_libc5, 0, 10000, 428084942);
+ rng_test (gsl_rng_random8_libc5, 0, 10000, 1910041713);
+ rng_test (gsl_rng_random32_libc5, 0, 10000, 1967452027);
+ rng_test (gsl_rng_random64_libc5, 0, 10000, 2106639801);
+ rng_test (gsl_rng_random128_libc5, 0, 10000, 428084942);
+ rng_test (gsl_rng_random256_libc5, 0, 10000, 116367984);
+
+ rng_test (gsl_rng_ranf, 0, 10000, 2152890433UL);
+ rng_test (gsl_rng_ranf, 2, 10000, 339327233);
+
+ /* Test constant relationship between int and double functions */
+
+ for (r = rngs ; *r != 0; r++)
+ rng_float_test (*r);
+
+ /* Test save/restore functions */
+
+ for (r = rngs ; *r != 0; r++)
+ rng_state_test (*r);
+
+ for (r = rngs ; *r != 0; r++)
+ rng_parallel_state_test (*r);
+
+ for (r = rngs ; *r != 0; r++)
+ rng_read_write_test (*r);
+
+ /* generic statistical tests (these are just to make sure that we
+ don't get any crazy results back from the generator, i.e. they
+ aren't a test of the algorithm, just the implementation) */
+
+ for (r = rngs ; *r != 0; r++)
+ generic_rng_test (*r);
+
+ exit (gsl_test_summary ());
+}
+
+
+void
+rng_test (const gsl_rng_type * T, unsigned long int seed, unsigned int n,
+ unsigned long int result)
+{
+ gsl_rng *r = gsl_rng_alloc (T);
+ unsigned int i;
+ unsigned long int k = 0;
+ int status;
+
+ if (seed != 0)
+ {
+ gsl_rng_set (r, seed);
+ }
+
+ for (i = 0; i < n; i++)
+ {
+ k = gsl_rng_get (r);
+ }
+
+ status = (k != result);
+ gsl_test (status, "%s, %u steps (%u observed vs %u expected)",
+ gsl_rng_name (r), n, k, result);
+
+ gsl_rng_free (r);
+}
+
+void
+rng_float_test (const gsl_rng_type * T)
+{
+ gsl_rng *ri = gsl_rng_alloc (T);
+ gsl_rng *rf = gsl_rng_alloc (T);
+
+ double u, c ;
+ unsigned int i;
+ unsigned long int k = 0;
+ int status = 0 ;
+
+ do
+ {
+ k = gsl_rng_get (ri);
+ u = gsl_rng_get (rf);
+ }
+ while (k == 0) ;
+
+ c = k / u ;
+ for (i = 0; i < N2; i++)
+ {
+ k = gsl_rng_get (ri);
+ u = gsl_rng_get (rf);
+ if (c*k != u)
+ {
+ status = 1 ;
+ break ;
+ }
+ }
+
+ gsl_test (status, "%s, ratio of int to double (%g observed vs %g expected)",
+ gsl_rng_name (ri), c, k/u);
+
+ gsl_rng_free (ri);
+ gsl_rng_free (rf);
+}
+
+
+void
+rng_state_test (const gsl_rng_type * T)
+{
+ unsigned long int test_a[N], test_b[N];
+
+ int i;
+
+ gsl_rng *r = gsl_rng_alloc (T);
+ gsl_rng *r_save = gsl_rng_alloc (T);
+
+ for (i = 0; i < N; ++i)
+ {
+ gsl_rng_get (r); /* throw away N iterations */
+ }
+
+ gsl_rng_memcpy (r_save, r); /* save the intermediate state */
+
+ for (i = 0; i < N; ++i)
+ {
+ test_a[i] = gsl_rng_get (r);
+ }
+
+ gsl_rng_memcpy (r, r_save); /* restore the intermediate state */
+ gsl_rng_free (r_save);
+
+ for (i = 0; i < N; ++i)
+ {
+ test_b[i] = gsl_rng_get (r);
+ }
+
+ {
+ int status = 0;
+ for (i = 0; i < N; ++i)
+ {
+ status |= (test_b[i] != test_a[i]);
+ }
+ gsl_test (status, "%s, random number state consistency",
+ gsl_rng_name (r));
+ }
+
+ gsl_rng_free (r);
+}
+
+
+void
+rng_parallel_state_test (const gsl_rng_type * T)
+{
+ unsigned long int test_a[N], test_b[N];
+ unsigned long int test_c[N], test_d[N];
+ double test_e[N], test_f[N];
+
+ int i;
+
+ gsl_rng *r1 = gsl_rng_alloc (T);
+ gsl_rng *r2 = gsl_rng_alloc (T);
+
+ for (i = 0; i < N; ++i)
+ {
+ gsl_rng_get (r1); /* throw away N iterations */
+ }
+
+ gsl_rng_memcpy (r2, r1); /* save the intermediate state */
+
+ for (i = 0; i < N; ++i)
+ {
+ /* check that there is no hidden state intermixed between r1 and r2 */
+ test_a[i] = gsl_rng_get (r1);
+ test_b[i] = gsl_rng_get (r2);
+ test_c[i] = gsl_rng_uniform_int (r1, 1234);
+ test_d[i] = gsl_rng_uniform_int (r2, 1234);
+ test_e[i] = gsl_rng_uniform (r1);
+ test_f[i] = gsl_rng_uniform (r2);
+ }
+
+ {
+ int status = 0;
+ for (i = 0; i < N; ++i)
+ {
+ status |= (test_b[i] != test_a[i]);
+ status |= (test_c[i] != test_d[i]);
+ status |= (test_e[i] != test_f[i]);
+ }
+ gsl_test (status, "%s, parallel random number state consistency",
+ gsl_rng_name (r1));
+ }
+
+ gsl_rng_free (r1);
+ gsl_rng_free (r2);
+
+}
+
+void
+rng_read_write_test (const gsl_rng_type * T)
+{
+ unsigned long int test_a[N], test_b[N];
+
+ int i;
+
+ gsl_rng *r = gsl_rng_alloc (T);
+
+ for (i = 0; i < N; ++i)
+ {
+ gsl_rng_get (r); /* throw away N iterations */
+ }
+
+ { /* save the state to a binary file */
+ FILE *f = fopen("test.dat", "wb");
+ gsl_rng_fwrite(f, r);
+ fclose(f);
+ }
+
+ for (i = 0; i < N; ++i)
+ {
+ test_a[i] = gsl_rng_get (r);
+ }
+
+ { /* read the state from a binary file */
+ FILE *f = fopen("test.dat", "rb");
+ gsl_rng_fread(f, r);
+ fclose(f);
+ }
+
+ for (i = 0; i < N; ++i)
+ {
+ test_b[i] = gsl_rng_get (r);
+ }
+
+ {
+ int status = 0;
+ for (i = 0; i < N; ++i)
+ {
+ status |= (test_b[i] != test_a[i]);
+ }
+ gsl_test (status, "%s, random number generator read and write",
+ gsl_rng_name (r));
+ }
+
+ gsl_rng_free (r);
+}
+
+void
+generic_rng_test (const gsl_rng_type * T)
+{
+ gsl_rng *r = gsl_rng_alloc (T);
+ const char *name = gsl_rng_name (r);
+ unsigned long int kmax = 0, kmin = 1000;
+ double sigma = 0;
+ const unsigned long int ran_max = gsl_rng_max (r);
+ const unsigned long int ran_min = gsl_rng_min (r);
+
+ int status = rng_max_test (r, &kmax, ran_max);
+
+ gsl_test (status,
+ "%s, observed vs theoretical maximum (%lu vs %lu)",
+ name, kmax, ran_max);
+
+ status = rng_min_test (r, &kmin, ran_min, ran_max);
+
+ gsl_test (status,
+ "%s, observed vs theoretical minimum (%lu vs %lu)",
+ name, kmin, ran_min);
+
+ status = rng_sum_test (r, &sigma);
+
+ gsl_test (status,
+ "%s, sum test within acceptable sigma (observed %.2g sigma)",
+ name, sigma);
+
+ status = rng_bin_test (r, &sigma);
+
+ gsl_test (status,
+ "%s, bin test within acceptable chisq (observed %.2g sigma)",
+ name, sigma);
+
+ gsl_rng_set (r, 1); /* set seed to 1 */
+ status = rng_max_test (r, &kmax, ran_max);
+
+ gsl_rng_set (r, 1); /* set seed to 1 */
+ status |= rng_min_test (r, &kmin, ran_min, ran_max);
+
+ gsl_rng_set (r, 1); /* set seed to 1 */
+ status |= rng_sum_test (r, &sigma);
+
+ gsl_test (status, "%s, maximum and sum tests for seed=1", name);
+
+ gsl_rng_set (r, 12345); /* set seed to a "typical" value */
+ status = rng_max_test (r, &kmax, ran_max);
+
+ gsl_rng_set (r, 12345); /* set seed to a "typical" value */
+ status |= rng_min_test (r, &kmin, ran_min, ran_max);
+
+ gsl_rng_set (r, 12345); /* set seed to a "typical" value */
+ status |= rng_sum_test (r, &sigma);
+
+ gsl_test (status, "%s, maximum and sum tests for non-default seeds", name);
+
+ gsl_rng_free (r);
+}
+
+int
+rng_max_test (gsl_rng * r, unsigned long int *kmax, unsigned long int ran_max)
+{
+ unsigned long int actual_uncovered;
+ double expect_uncovered;
+ int status;
+ unsigned long int max = 0;
+ int i;
+
+ for (i = 0; i < N2; ++i)
+ {
+ unsigned long int k = gsl_rng_get (r);
+ if (k > max)
+ max = k;
+ }
+
+ *kmax = max;
+
+ actual_uncovered = ran_max - max;
+ expect_uncovered = (double) ran_max / (double) N2;
+
+ status = (max > ran_max) || (actual_uncovered > 7 * expect_uncovered) ;
+
+ return status;
+}
+
+int
+rng_min_test (gsl_rng * r, unsigned long int *kmin,
+ unsigned long int ran_min, unsigned long int ran_max)
+{
+ unsigned long int actual_uncovered;
+ double expect_uncovered;
+ int status;
+ unsigned long int min = 1000000000UL;
+ int i;
+
+ for (i = 0; i < N2; ++i)
+ {
+ unsigned long int k = gsl_rng_get (r);
+ if (k < min)
+ min = k;
+ }
+
+ *kmin = min;
+
+ actual_uncovered = min - ran_min;
+ expect_uncovered = (double) ran_max / (double) N2;
+
+ status = (min < ran_min) || (actual_uncovered > 7 * expect_uncovered);
+
+ return status;
+}
+
+int
+rng_sum_test (gsl_rng * r, double *sigma)
+{
+ double sum = 0;
+ int i, status;
+
+ for (i = 0; i < N2; ++i)
+ {
+ double x = gsl_rng_uniform (r) - 0.5;
+ sum += x;
+ }
+
+ sum /= N2;
+
+ /* expect the average to have a variance of 1/(12 n) */
+
+ *sigma = sum * sqrt (12.0 * N2);
+
+ /* more than 3 sigma is an error (increased to 3.1 to avoid false alarms) */
+
+ status = (fabs (*sigma) > 3.1 || fabs(*sigma) < 0.003);
+
+ if (status) {
+ fprintf(stderr,"sum=%g, sigma=%g\n",sum,*sigma);
+ }
+
+ return status;
+}
+
+#define BINS 17
+#define EXTRA 10
+int
+rng_bin_test (gsl_rng * r, double *sigma)
+{
+ int count[BINS+EXTRA];
+ double chisq = 0;
+ int i, status;
+
+ for (i = 0; i < BINS+EXTRA; i++)
+ count[i] = 0 ;
+
+
+ for (i = 0; i < N2; i++)
+ {
+ int j = gsl_rng_uniform_int (r, BINS);
+ count[j]++ ;
+ }
+
+ chisq = 0 ;
+ for (i = 0; i < BINS; i++)
+ {
+ double x = (double)N2/(double)BINS ;
+ double d = (count[i] - x) ;
+ chisq += (d*d) / x;
+ }
+
+ *sigma = sqrt(chisq/BINS) ;
+
+ /* more than 3 sigma is an error */
+
+ status = (fabs (*sigma) > 3 || fabs(*sigma) < 0.003);
+
+ for (i = BINS; i < BINS+EXTRA; i++)
+ {
+ if (count[i] != 0)
+ {
+ status = 1 ;
+ gsl_test (status,
+ "%s, wrote outside range in bin test "
+ "(%d observed vs %d expected)",
+ gsl_rng_name(r), i, BINS - 1);
+ }
+ }
+
+ return status;
+}
+
generated by cgit v1.2.3 (git 2.25.1) at 2024-05-01 00:00:42 +0000