/* 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;
}
