Annotation of OpenXM_contrib/gmp/tests/rand/stat.c, Revision 1.1
1.1 ! maekawa 1: /* stat.c -- statistical tests of random number sequences. */
! 2:
! 3: /*
! 4: Copyright (C) 1999, 2000 Free Software Foundation, Inc.
! 5:
! 6: This file is part of the GNU MP Library.
! 7:
! 8: The GNU MP Library is free software; you can redistribute it and/or modify
! 9: it under the terms of the GNU Lesser General Public License as published by
! 10: the Free Software Foundation; either version 2.1 of the License, or (at your
! 11: option) any later version.
! 12:
! 13: The GNU MP Library is distributed in the hope that it will be useful, but
! 14: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
! 15: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
! 16: License for more details.
! 17:
! 18: You should have received a copy of the GNU Lesser General Public License
! 19: along with the GNU MP Library; see the file COPYING.LIB. If not, write to
! 20: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
! 21: MA 02111-1307, USA.
! 22: */
! 23:
! 24: /* Examples:
! 25:
! 26: $ gen 1000 | stat
! 27: Test 1000 real numbers.
! 28:
! 29: $ gen 30000 | stat -2 1000
! 30: Test 1000 real numbers 30 times and then test the 30 results in a
! 31: ``second level''.
! 32:
! 33: $ gen -f mpz_urandomb 1000 | stat -i 0xffffffff
! 34: Test 1000 integers 0 <= X <= 2^32-1.
! 35:
! 36: $ gen -f mpz_urandomb -z 34 1000 | stat -i 0x3ffffffff
! 37: Test 1000 integers 0 <= X <= 2^34-1.
! 38:
! 39: */
! 40:
! 41: #include <stdio.h>
! 42: #include <stdlib.h>
! 43: #include <unistd.h>
! 44: #include <math.h>
! 45: #include "gmp.h"
! 46: #include "gmpstat.h"
! 47:
! 48: #if !HAVE_DECL_OPTARG
! 49: extern char *optarg;
! 50: extern int optind, opterr;
! 51: #endif
! 52:
! 53: #define FVECSIZ (100000L)
! 54:
! 55: int g_debug = 0;
! 56:
! 57: static void
! 58: print_ks_results (mpf_t f_p, mpf_t f_p_prob,
! 59: mpf_t f_m, mpf_t f_m_prob,
! 60: FILE *fp)
! 61: {
! 62: double p, pp, m, mp;
! 63:
! 64: p = mpf_get_d (f_p);
! 65: m = mpf_get_d (f_m);
! 66: pp = mpf_get_d (f_p_prob);
! 67: mp = mpf_get_d (f_m_prob);
! 68:
! 69: fprintf (fp, "%.4f (%.0f%%)\t", p, pp * 100.0);
! 70: fprintf (fp, "%.4f (%.0f%%)\n", m, mp * 100.0);
! 71: }
! 72:
! 73: static void
! 74: print_x2_table (unsigned int v, FILE *fp)
! 75: {
! 76: double t[7];
! 77: int f;
! 78:
! 79:
! 80: fprintf (fp, "Chi-square table for v=%u\n", v);
! 81: fprintf (fp, "1%%\t5%%\t25%%\t50%%\t75%%\t95%%\t99%%\n");
! 82: x2_table (t, v);
! 83: for (f = 0; f < 7; f++)
! 84: fprintf (fp, "%.2f\t", t[f]);
! 85: fputs ("\n", fp);
! 86: }
! 87:
! 88:
! 89:
! 90: /* Pks () -- Distribution function for KS results with a big n (like 1000
! 91: or so): F(x) = 1 - pow(e, -2*x^2) [Knuth, vol 2, p.51]. */
! 92: /* gnuplot: plot [0:1] Pks(x), Pks(x) = 1-exp(-2*x**2) */
! 93:
! 94: static void
! 95: Pks (mpf_t p, mpf_t x)
! 96: {
! 97: double dt; /* temp double */
! 98:
! 99: mpf_set (p, x);
! 100: mpf_mul (p, p, p); /* p = x^2 */
! 101: mpf_mul_ui (p, p, 2); /* p = 2*x^2 */
! 102: mpf_neg (p, p); /* p = -2*x^2 */
! 103: /* No pow() in gmp. Use doubles. */
! 104: /* FIXME: Use exp()? */
! 105: dt = pow (M_E, mpf_get_d (p));
! 106: mpf_set_d (p, dt);
! 107: mpf_ui_sub (p, 1, p);
! 108: }
! 109:
! 110: /* f_freq() -- frequency test on real numbers 0<=f<1*/
! 111: static void
! 112: f_freq (const unsigned l1runs, const unsigned l2runs,
! 113: mpf_t fvec[], const unsigned long n)
! 114: {
! 115: unsigned f;
! 116: mpf_t f_p, f_p_prob;
! 117: mpf_t f_m, f_m_prob;
! 118: mpf_t *l1res; /* level 1 result array */
! 119:
! 120: mpf_init (f_p); mpf_init (f_m);
! 121: mpf_init (f_p_prob); mpf_init (f_m_prob);
! 122:
! 123:
! 124: /* Allocate space for 1st level results. */
! 125: l1res = (mpf_t *) malloc (l2runs * 2 * sizeof (mpf_t));
! 126: if (NULL == l1res)
! 127: {
! 128: fprintf (stderr, "stat: malloc failure\n");
! 129: exit (1);
! 130: }
! 131:
! 132: printf ("\nEquidistribution/Frequency test on real numbers (0<=X<1):\n");
! 133: printf ("\tKp\t\tKm\n");
! 134:
! 135: for (f = 0; f < l2runs; f++)
! 136: {
! 137: /* f_printvec (fvec, n); */
! 138: mpf_freqt (f_p, f_m, fvec + f * n, n);
! 139:
! 140: /* what's the probability of getting these results? */
! 141: ks_table (f_p_prob, f_p, n);
! 142: ks_table (f_m_prob, f_m, n);
! 143:
! 144: if (l1runs == 0)
! 145: {
! 146: /*printf ("%u:\t", f + 1);*/
! 147: print_ks_results (f_p, f_p_prob, f_m, f_m_prob, stdout);
! 148: }
! 149: else
! 150: {
! 151: /* save result */
! 152: mpf_init_set (l1res[f], f_p);
! 153: mpf_init_set (l1res[f + l2runs], f_m);
! 154: }
! 155: }
! 156:
! 157: /* Now, apply the KS test on the results from the 1st level rounds
! 158: with the distribution
! 159: F(x) = 1 - pow(e, -2*x^2) [Knuth, vol 2, p.51] */
! 160:
! 161: if (l1runs != 0)
! 162: {
! 163: /*printf ("-------------------------------------\n");*/
! 164:
! 165: /* The Kp's. */
! 166: ks (f_p, f_m, l1res, Pks, l2runs);
! 167: ks_table (f_p_prob, f_p, l2runs);
! 168: ks_table (f_m_prob, f_m, l2runs);
! 169: printf ("Kp:\t");
! 170: print_ks_results (f_p, f_p_prob, f_m, f_m_prob, stdout);
! 171:
! 172: /* The Km's. */
! 173: ks (f_p, f_m, l1res + l2runs, Pks, l2runs);
! 174: ks_table (f_p_prob, f_p, l2runs);
! 175: ks_table (f_m_prob, f_m, l2runs);
! 176: printf ("Km:\t");
! 177: print_ks_results (f_p, f_p_prob, f_m, f_m_prob, stdout);
! 178: }
! 179:
! 180: mpf_clear (f_p); mpf_clear (f_m);
! 181: mpf_clear (f_p_prob); mpf_clear (f_m_prob);
! 182: free (l1res);
! 183: }
! 184:
! 185: /* z_freq(l1runs, l2runs, zvec, n, max) -- frequency test on integers
! 186: 0<=z<=MAX */
! 187: static void
! 188: z_freq (const unsigned l1runs,
! 189: const unsigned l2runs,
! 190: mpz_t zvec[],
! 191: const unsigned long n,
! 192: unsigned int max)
! 193: {
! 194: mpf_t V; /* result */
! 195: double d_V; /* result as a double */
! 196:
! 197: mpf_init (V);
! 198:
! 199:
! 200: printf ("\nEquidistribution/Frequency test on integers (0<=X<=%u):\n", max);
! 201: print_x2_table (max, stdout);
! 202:
! 203: mpz_freqt (V, zvec, max, n);
! 204:
! 205: d_V = mpf_get_d (V);
! 206: printf ("V = %.2f (n = %lu)\n", d_V, n);
! 207:
! 208: mpf_clear (V);
! 209: }
! 210:
! 211: unsigned int stat_debug = 0;
! 212:
! 213: int
! 214: main (argc, argv)
! 215: int argc;
! 216: char *argv[];
! 217: {
! 218: const char usage[] =
! 219: "usage: stat [-d] [-2 runs] [-i max | -r max] [file]\n" \
! 220: " file filename\n" \
! 221: " -2 runs perform 2-level test with RUNS runs on 1st level\n" \
! 222: " -d increase debugging level\n" \
! 223: " -i max input is integers 0 <= Z <= MAX\n" \
! 224: " -r max input is real numbers 0 <= R < 1 and use MAX as\n" \
! 225: " maximum value when converting real numbers to integers\n" \
! 226: "";
! 227:
! 228: mpf_t fvec[FVECSIZ];
! 229: mpz_t zvec[FVECSIZ];
! 230: unsigned long int f, n, vecentries;
! 231: char *filen;
! 232: FILE *fp;
! 233: int c;
! 234: int omitoutput = 0;
! 235: int realinput = -1; /* 1: input is real numbers 0<=R<1;
! 236: 0: input is integers 0 <= Z <= MAX. */
! 237: long l1runs = 0, /* 1st level runs */
! 238: l2runs = 1; /* 2nd level runs */
! 239: mpf_t f_temp;
! 240: mpz_t z_imax; /* max value when converting between
! 241: real number and integer. */
! 242: mpf_t f_imax_plus1; /* f_imax + 1 stored in an mpf_t for
! 243: convenience */
! 244: mpf_t f_imax_minus1; /* f_imax - 1 stored in an mpf_t for
! 245: convenience */
! 246:
! 247:
! 248: mpf_init (f_temp);
! 249: mpz_init_set_ui (z_imax, 0x7fffffff);
! 250: mpf_init (f_imax_plus1);
! 251: mpf_init (f_imax_minus1);
! 252:
! 253: while ((c = getopt (argc, argv, "d2:i:r:")) != -1)
! 254: switch (c)
! 255: {
! 256: case '2':
! 257: l1runs = atol (optarg);
! 258: l2runs = -1; /* set later on */
! 259: break;
! 260: case 'd': /* increase debug level */
! 261: stat_debug++;
! 262: break;
! 263: case 'i':
! 264: if (1 == realinput)
! 265: {
! 266: fputs ("stat: options -i and -r are mutually exclusive\n", stderr);
! 267: exit (1);
! 268: }
! 269: if (mpz_set_str (z_imax, optarg, 0))
! 270: {
! 271: fprintf (stderr, "stat: bad max value %s\n", optarg);
! 272: exit (1);
! 273: }
! 274: realinput = 0;
! 275: break;
! 276: case 'r':
! 277: if (0 == realinput)
! 278: {
! 279: fputs ("stat: options -i and -r are mutually exclusive\n", stderr);
! 280: exit (1);
! 281: }
! 282: if (mpz_set_str (z_imax, optarg, 0))
! 283: {
! 284: fprintf (stderr, "stat: bad max value %s\n", optarg);
! 285: exit (1);
! 286: }
! 287: realinput = 1;
! 288: break;
! 289: case 'o':
! 290: omitoutput = atoi (optarg);
! 291: break;
! 292: case '?':
! 293: default:
! 294: fputs (usage, stderr);
! 295: exit (1);
! 296: }
! 297: argc -= optind;
! 298: argv += optind;
! 299:
! 300: if (argc < 1)
! 301: fp = stdin;
! 302: else
! 303: filen = argv[0];
! 304:
! 305: if (fp != stdin)
! 306: if (NULL == (fp = fopen (filen, "r")))
! 307: {
! 308: perror (filen);
! 309: exit (1);
! 310: }
! 311:
! 312: if (-1 == realinput)
! 313: realinput = 1; /* default is real numbers */
! 314:
! 315: /* read file and fill appropriate vec */
! 316: if (1 == realinput) /* real input */
! 317: {
! 318: for (f = 0; f < FVECSIZ ; f++)
! 319: {
! 320: mpf_init (fvec[f]);
! 321: if (!mpf_inp_str (fvec[f], fp, 10))
! 322: break;
! 323: }
! 324: }
! 325: else /* integer input */
! 326: {
! 327: for (f = 0; f < FVECSIZ ; f++)
! 328: {
! 329: mpz_init (zvec[f]);
! 330: if (!mpz_inp_str (zvec[f], fp, 10))
! 331: break;
! 332: }
! 333: }
! 334: vecentries = n = f; /* number of entries read */
! 335: fclose (fp);
! 336:
! 337: if (FVECSIZ == f)
! 338: fprintf (stderr, "stat: warning: discarding input due to lazy allocation "\
! 339: "of only %ld entries. sorry.\n", FVECSIZ);
! 340:
! 341: printf ("Got %lu numbers.\n", n);
! 342:
! 343: /* convert and fill the other vec */
! 344: /* since fvec[] contains 0<=f<1 and we want ivec[] to contain
! 345: 0<=z<=imax and we are truncating all fractions when
! 346: converting float to int, we have to add 1 to imax.*/
! 347: mpf_set_z (f_imax_plus1, z_imax);
! 348: mpf_add_ui (f_imax_plus1, f_imax_plus1, 1);
! 349: if (1 == realinput) /* fill zvec[] */
! 350: {
! 351: for (f = 0; f < n; f++)
! 352: {
! 353: mpf_mul (f_temp, fvec[f], f_imax_plus1);
! 354: mpz_init (zvec[f]);
! 355: mpz_set_f (zvec[f], f_temp); /* truncating fraction */
! 356: if (stat_debug > 1)
! 357: {
! 358: mpz_out_str (stderr, 10, zvec[f]);
! 359: fputs ("\n", stderr);
! 360: }
! 361: }
! 362: }
! 363: else /* integer input; fill fvec[] */
! 364: {
! 365: /* mpf_set_z (f_imax_minus1, z_imax);
! 366: mpf_sub_ui (f_imax_minus1, f_imax_minus1, 1);*/
! 367: for (f = 0; f < n; f++)
! 368: {
! 369: mpf_init (fvec[f]);
! 370: mpf_set_z (fvec[f], zvec[f]);
! 371: mpf_div (fvec[f], fvec[f], f_imax_plus1);
! 372: if (stat_debug > 1)
! 373: {
! 374: mpf_out_str (stderr, 10, 0, fvec[f]);
! 375: fputs ("\n", stderr);
! 376: }
! 377: }
! 378: }
! 379:
! 380: /* 2 levels? */
! 381: if (1 != l2runs)
! 382: {
! 383: l2runs = n / l1runs;
! 384: printf ("Doing %ld second level rounds "\
! 385: "with %ld entries in each round", l2runs, l1runs);
! 386: if (n % l1runs)
! 387: printf (" (discarding %ld entr%s)", n % l1runs,
! 388: n % l1runs == 1 ? "y" : "ies");
! 389: puts (".");
! 390: n = l1runs;
! 391: }
! 392:
! 393: #ifndef DONT_FFREQ
! 394: f_freq (l1runs, l2runs, fvec, n);
! 395: #endif
! 396: #ifdef DO_ZFREQ
! 397: z_freq (l1runs, l2runs, zvec, n, mpz_get_ui (z_imax));
! 398: #endif
! 399:
! 400: mpf_clear (f_temp); mpz_clear (z_imax);
! 401: mpf_clear (f_imax_plus1);
! 402: mpf_clear (f_imax_minus1);
! 403: for (f = 0; f < vecentries; f++)
! 404: {
! 405: mpf_clear (fvec[f]);
! 406: mpz_clear (zvec[f]);
! 407: }
! 408:
! 409: return 0;
! 410: }
! 411:
! 412:
! 413:
! 414:
! 415:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to stat.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / tests / rand