stat.c
上传用户:qaz666999
上传日期:2022-08-06
资源大小:2570k
文件大小:10k
源码类别:

数学计算

开发平台:

Unix_Linux

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