diff options
Diffstat (limited to 'gmp-6.3.0/tests/rand/stat.c')
| -rw-r--r-- | gmp-6.3.0/tests/rand/stat.c | 406 |
1 files changed, 406 insertions, 0 deletions
diff --git a/gmp-6.3.0/tests/rand/stat.c b/gmp-6.3.0/tests/rand/stat.c new file mode 100644 index 0000000..fa06f4c --- /dev/null +++ b/gmp-6.3.0/tests/rand/stat.c @@ -0,0 +1,406 @@ +/* stat.c -- statistical tests of random number sequences. */ + +/* +Copyright 1999, 2000 Free Software Foundation, Inc. + +This file is part of the GNU MP Library test suite. + +The GNU MP Library test suite 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 3 of the License, +or (at your option) any later version. + +The GNU MP Library test suite 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 +the GNU MP Library test suite. If not, see https://www.gnu.org/licenses/. */ + +/* Examples: + + $ gen 1000 | stat +Test 1000 real numbers. + + $ gen 30000 | stat -2 1000 +Test 1000 real numbers 30 times and then test the 30 results in a +``second level''. + + $ gen -f mpz_urandomb 1000 | stat -i 0xffffffff +Test 1000 integers 0 <= X <= 2^32-1. + + $ gen -f mpz_urandomb -z 34 1000 | stat -i 0x3ffffffff +Test 1000 integers 0 <= X <= 2^34-1. + +*/ + +#include <stdio.h> +#include <stdlib.h> +#include <unistd.h> +#include <math.h> +#include "gmpstat.h" + +#if !HAVE_DECL_OPTARG +extern char *optarg; +extern int optind, opterr; +#endif + +#define FVECSIZ (100000L) + +int g_debug = 0; + +static void +print_ks_results (mpf_t f_p, mpf_t f_p_prob, + mpf_t f_m, mpf_t f_m_prob, + FILE *fp) +{ + double p, pp, m, mp; + + p = mpf_get_d (f_p); + m = mpf_get_d (f_m); + pp = mpf_get_d (f_p_prob); + mp = mpf_get_d (f_m_prob); + + fprintf (fp, "%.4f (%.0f%%)\t", p, pp * 100.0); + fprintf (fp, "%.4f (%.0f%%)\n", m, mp * 100.0); +} + +static void +print_x2_table (unsigned int v, FILE *fp) +{ + double t[7]; + int f; + + + fprintf (fp, "Chi-square table for v=%u\n", v); + fprintf (fp, "1%%\t5%%\t25%%\t50%%\t75%%\t95%%\t99%%\n"); + x2_table (t, v); + for (f = 0; f < 7; f++) + fprintf (fp, "%.2f\t", t[f]); + fputs ("\n", fp); +} + + + +/* Pks () -- Distribution function for KS results with a big n (like 1000 + or so): F(x) = 1 - pow(e, -2*x^2) [Knuth, vol 2, p.51]. */ +/* gnuplot: plot [0:1] Pks(x), Pks(x) = 1-exp(-2*x**2) */ + +static void +Pks (mpf_t p, mpf_t x) +{ + double dt; /* temp double */ + + mpf_set (p, x); + mpf_mul (p, p, p); /* p = x^2 */ + mpf_mul_ui (p, p, 2); /* p = 2*x^2 */ + mpf_neg (p, p); /* p = -2*x^2 */ + /* No pow() in gmp. Use doubles. */ + /* FIXME: Use exp()? */ + dt = pow (M_E, mpf_get_d (p)); + mpf_set_d (p, dt); + mpf_ui_sub (p, 1, p); +} + +/* f_freq() -- frequency test on real numbers 0<=f<1*/ +static void +f_freq (const unsigned l1runs, const unsigned l2runs, + mpf_t fvec[], const unsigned long n) +{ + unsigned f; + mpf_t f_p, f_p_prob; + mpf_t f_m, f_m_prob; + mpf_t *l1res; /* level 1 result array */ + + mpf_init (f_p); mpf_init (f_m); + mpf_init (f_p_prob); mpf_init (f_m_prob); + + + /* Allocate space for 1st level results. */ + l1res = (mpf_t *) malloc (l2runs * 2 * sizeof (mpf_t)); + if (NULL == l1res) + { + fprintf (stderr, "stat: malloc failure\n"); + exit (1); + } + + printf ("\nEquidistribution/Frequency test on real numbers (0<=X<1):\n"); + printf ("\tKp\t\tKm\n"); + + for (f = 0; f < l2runs; f++) + { + /* f_printvec (fvec, n); */ + mpf_freqt (f_p, f_m, fvec + f * n, n); + + /* what's the probability of getting these results? */ + ks_table (f_p_prob, f_p, n); + ks_table (f_m_prob, f_m, n); + + if (l1runs == 0) + { + /*printf ("%u:\t", f + 1);*/ + print_ks_results (f_p, f_p_prob, f_m, f_m_prob, stdout); + } + else + { + /* save result */ + mpf_init_set (l1res[f], f_p); + mpf_init_set (l1res[f + l2runs], f_m); + } + } + + /* Now, apply the KS test on the results from the 1st level rounds + with the distribution + F(x) = 1 - pow(e, -2*x^2) [Knuth, vol 2, p.51] */ + + if (l1runs != 0) + { + /*printf ("-------------------------------------\n");*/ + + /* The Kp's. */ + ks (f_p, f_m, l1res, Pks, l2runs); + ks_table (f_p_prob, f_p, l2runs); + ks_table (f_m_prob, f_m, l2runs); + printf ("Kp:\t"); + print_ks_results (f_p, f_p_prob, f_m, f_m_prob, stdout); + + /* The Km's. */ + ks (f_p, f_m, l1res + l2runs, Pks, l2runs); + ks_table (f_p_prob, f_p, l2runs); + ks_table (f_m_prob, f_m, l2runs); + printf ("Km:\t"); + print_ks_results (f_p, f_p_prob, f_m, f_m_prob, stdout); + } + + mpf_clear (f_p); mpf_clear (f_m); + mpf_clear (f_p_prob); mpf_clear (f_m_prob); + free (l1res); +} + +/* z_freq(l1runs, l2runs, zvec, n, max) -- frequency test on integers + 0<=z<=MAX */ +static void +z_freq (const unsigned l1runs, + const unsigned l2runs, + mpz_t zvec[], + const unsigned long n, + unsigned int max) +{ + mpf_t V; /* result */ + double d_V; /* result as a double */ + + mpf_init (V); + + + printf ("\nEquidistribution/Frequency test on integers (0<=X<=%u):\n", max); + print_x2_table (max, stdout); + + mpz_freqt (V, zvec, max, n); + + d_V = mpf_get_d (V); + printf ("V = %.2f (n = %lu)\n", d_V, n); + + mpf_clear (V); +} + +unsigned int stat_debug = 0; + +int +main (argc, argv) + int argc; + char *argv[]; +{ + const char usage[] = + "usage: stat [-d] [-2 runs] [-i max | -r max] [file]\n" \ + " file filename\n" \ + " -2 runs perform 2-level test with RUNS runs on 1st level\n" \ + " -d increase debugging level\n" \ + " -i max input is integers 0 <= Z <= MAX\n" \ + " -r max input is real numbers 0 <= R < 1 and use MAX as\n" \ + " maximum value when converting real numbers to integers\n" \ + ""; + + mpf_t fvec[FVECSIZ]; + mpz_t zvec[FVECSIZ]; + unsigned long int f, n, vecentries; + char *filen; + FILE *fp; + int c; + int omitoutput = 0; + int realinput = -1; /* 1: input is real numbers 0<=R<1; + 0: input is integers 0 <= Z <= MAX. */ + long l1runs = 0, /* 1st level runs */ + l2runs = 1; /* 2nd level runs */ + mpf_t f_temp; + mpz_t z_imax; /* max value when converting between + real number and integer. */ + mpf_t f_imax_plus1; /* f_imax + 1 stored in an mpf_t for + convenience */ + mpf_t f_imax_minus1; /* f_imax - 1 stored in an mpf_t for + convenience */ + + + mpf_init (f_temp); + mpz_init_set_ui (z_imax, 0x7fffffff); + mpf_init (f_imax_plus1); + mpf_init (f_imax_minus1); + + while ((c = getopt (argc, argv, "d2:i:r:")) != -1) + switch (c) + { + case '2': + l1runs = atol (optarg); + l2runs = -1; /* set later on */ + break; + case 'd': /* increase debug level */ + stat_debug++; + break; + case 'i': + if (1 == realinput) + { + fputs ("stat: options -i and -r are mutually exclusive\n", stderr); + exit (1); + } + if (mpz_set_str (z_imax, optarg, 0)) + { + fprintf (stderr, "stat: bad max value %s\n", optarg); + exit (1); + } + realinput = 0; + break; + case 'r': + if (0 == realinput) + { + fputs ("stat: options -i and -r are mutually exclusive\n", stderr); + exit (1); + } + if (mpz_set_str (z_imax, optarg, 0)) + { + fprintf (stderr, "stat: bad max value %s\n", optarg); + exit (1); + } + realinput = 1; + break; + case 'o': + omitoutput = atoi (optarg); + break; + case '?': + default: + fputs (usage, stderr); + exit (1); + } + argc -= optind; + argv += optind; + + if (argc < 1) + fp = stdin; + else + filen = argv[0]; + + if (fp != stdin) + if (NULL == (fp = fopen (filen, "r"))) + { + perror (filen); + exit (1); + } + + if (-1 == realinput) + realinput = 1; /* default is real numbers */ + + /* read file and fill appropriate vec */ + if (1 == realinput) /* real input */ + { + for (f = 0; f < FVECSIZ ; f++) + { + mpf_init (fvec[f]); + if (!mpf_inp_str (fvec[f], fp, 10)) + break; + } + } + else /* integer input */ + { + for (f = 0; f < FVECSIZ ; f++) + { + mpz_init (zvec[f]); + if (!mpz_inp_str (zvec[f], fp, 10)) + break; + } + } + vecentries = n = f; /* number of entries read */ + fclose (fp); + + if (FVECSIZ == f) + fprintf (stderr, "stat: warning: discarding input due to lazy allocation "\ + "of only %ld entries. sorry.\n", FVECSIZ); + + printf ("Got %lu numbers.\n", n); + + /* convert and fill the other vec */ + /* since fvec[] contains 0<=f<1 and we want ivec[] to contain + 0<=z<=imax and we are truncating all fractions when + converting float to int, we have to add 1 to imax.*/ + mpf_set_z (f_imax_plus1, z_imax); + mpf_add_ui (f_imax_plus1, f_imax_plus1, 1); + if (1 == realinput) /* fill zvec[] */ + { + for (f = 0; f < n; f++) + { + mpf_mul (f_temp, fvec[f], f_imax_plus1); + mpz_init (zvec[f]); + mpz_set_f (zvec[f], f_temp); /* truncating fraction */ + if (stat_debug > 1) + { + mpz_out_str (stderr, 10, zvec[f]); + fputs ("\n", stderr); + } + } + } + else /* integer input; fill fvec[] */ + { + /* mpf_set_z (f_imax_minus1, z_imax); + mpf_sub_ui (f_imax_minus1, f_imax_minus1, 1);*/ + for (f = 0; f < n; f++) + { + mpf_init (fvec[f]); + mpf_set_z (fvec[f], zvec[f]); + mpf_div (fvec[f], fvec[f], f_imax_plus1); + if (stat_debug > 1) + { + mpf_out_str (stderr, 10, 0, fvec[f]); + fputs ("\n", stderr); + } + } + } + + /* 2 levels? */ + if (1 != l2runs) + { + l2runs = n / l1runs; + printf ("Doing %ld second level rounds "\ + "with %ld entries in each round", l2runs, l1runs); + if (n % l1runs) + printf (" (discarding %ld entr%s)", n % l1runs, + n % l1runs == 1 ? "y" : "ies"); + puts ("."); + n = l1runs; + } + +#ifndef DONT_FFREQ + f_freq (l1runs, l2runs, fvec, n); +#endif +#ifdef DO_ZFREQ + z_freq (l1runs, l2runs, zvec, n, mpz_get_ui (z_imax)); +#endif + + mpf_clear (f_temp); mpz_clear (z_imax); + mpf_clear (f_imax_plus1); + mpf_clear (f_imax_minus1); + for (f = 0; f < vecentries; f++) + { + mpf_clear (fvec[f]); + mpz_clear (zvec[f]); + } + + return 0; +} |