Annotation of OpenXM_contrib/gmp/tests/rand/stat.c, Revision 1.1.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