home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: OtherApp / OtherApp.zip / entrand.zip / ent / randtest.c < prev    next >
C/C++ Source or Header  |  1998-10-19  |  4KB  |  180 lines
  1. /*
  2.  
  3.      Apply various randomness tests to a stream of bytes
  4.  
  5.           by John Walker  --  September 1996
  6.                http://www.fourmilab.ch/
  7.  
  8. */
  9.  
  10. #include <math.h>
  11.  
  12. #define FALSE 0
  13. #define TRUE  1
  14.  
  15. #define log2of10 3.32192809488736234787
  16.  
  17. static int binary = FALSE;       /* Treat input as a bitstream */
  18.  
  19. static long ccount[256],       /* Bins to count occurrences of values */
  20.         totalc = 0;        /* Total bytes counted */
  21. static double prob[256];       /* Probabilities per bin for entropy */
  22.  
  23. /*  LOG2  --  Calculate log to the base 2  */
  24.  
  25. static double log2(double x)
  26. {
  27.     return log2of10 * log10(x);
  28. }
  29.  
  30. #define MONTEN    6              /* Bytes used as Monte Carlo
  31.                      co-ordinates.    This should be no more
  32.                      bits than the mantissa of your
  33.                                          "double" floating point type. */
  34.  
  35. static int mp, sccfirst;
  36. static unsigned int monte[MONTEN];
  37. static long inmont, mcount;
  38. static double a, cexp, incirc, montex, montey, montepi,
  39.           scc, sccun, sccu0, scclast, scct1, scct2, scct3,
  40.           ent, chisq, datasum;
  41.  
  42. /*  RT_INIT  --  Initialise random test counters.  */
  43.  
  44. void rt_init(int binmode)
  45. {
  46.     int i;
  47.  
  48.     binary = binmode;           /* Set binary / byte mode */
  49.  
  50.     /* Initialise for calculations */
  51.  
  52.     ent = 0.0;               /* Clear entropy accumulator */
  53.     chisq = 0.0;           /* Clear Chi-Square */
  54.     datasum = 0.0;           /* Clear sum of bytes for arithmetic mean */
  55.  
  56.     mp = 0;               /* Reset Monte Carlo accumulator pointer */
  57.     mcount = 0;            /* Clear Monte Carlo tries */
  58.     inmont = 0;            /* Clear Monte Carlo inside count */
  59.     incirc = 65535.0 * 65535.0;/* In-circle distance for Monte Carlo */
  60.  
  61.     sccfirst = TRUE;           /* Mark first time for serial correlation */
  62.     scct1 = scct2 = scct3 = 0.0; /* Clear serial correlation terms */
  63.  
  64.     incirc = pow(pow(256.0, (double) (MONTEN / 2)) - 1, 2.0);
  65.  
  66.     for (i = 0; i < 256; i++) {
  67.     ccount[i] = 0;
  68.     }
  69.     totalc = 0;
  70. }
  71.  
  72. /*  RT_ADD  --    Add one or more bytes to accumulation.    */
  73.  
  74. void rt_add(void *buf, int bufl)
  75. {
  76.     unsigned char *bp = buf;
  77.     int oc, c, bean;
  78.  
  79.     while (bean = 0, (bufl-- > 0)) {
  80.        oc = *bp++;
  81.  
  82.        do {
  83.       if (binary) {
  84.          c = !!(oc & 0x80);
  85.       } else {
  86.          c = oc;
  87.       }
  88.       ccount[c]++;          /* Update counter for this bin */
  89.       totalc++;
  90.  
  91.       /* Update inside / outside circle counts for Monte Carlo
  92.          computation of PI */
  93.  
  94.       if (bean == 0) {
  95.           monte[mp++] = oc;       /* Save character for Monte Carlo */
  96.           if (mp >= MONTEN) {     /* Calculate every MONTEN character */
  97.          int mj;
  98.  
  99.          mp = 0;
  100.          mcount++;
  101.          montex = montey = 0;
  102.          for (mj = 0; mj < MONTEN / 2; mj++) {
  103.             montex = (montex * 256.0) + monte[mj];
  104.             montey = (montey * 256.0) + monte[(MONTEN / 2) + mj];
  105.          }
  106.          if ((montex * montex + montey *  montey) <= incirc) {
  107.             inmont++;
  108.          }
  109.           }
  110.       }
  111.  
  112.       /* Update calculation of serial correlation coefficient */
  113.  
  114.       sccun = c;
  115.       if (sccfirst) {
  116.          sccfirst = FALSE;
  117.          scclast = 0;
  118.          sccu0 = sccun;
  119.       } else {
  120.          scct1 = scct1 + scclast * sccun;
  121.       }
  122.       scct2 = scct2 + sccun;
  123.       scct3 = scct3 + (sccun * sccun);
  124.       scclast = sccun;
  125.       oc <<= 1;
  126.        } while (binary && (++bean < 8));
  127.     }
  128. }
  129.  
  130. /*  RT_END  --    Complete calculation and return results.  */
  131.  
  132. void rt_end(double *r_ent, double *r_chisq, double *r_mean,
  133.         double *r_montepicalc, double *r_scc)
  134. {
  135.     int i;
  136.  
  137.     /* Complete calculation of serial correlation coefficient */
  138.  
  139.     scct1 = scct1 + scclast * sccu0;
  140.     scct2 = scct2 * scct2;
  141.     scc = totalc * scct3 - scct2;
  142.     if (scc == 0.0) {
  143.        scc = -100000;
  144.     } else {
  145.        scc = (totalc * scct1 - scct2) / scc;
  146.     }
  147.  
  148.     /* Scan bins and calculate probability for each bin and
  149.        Chi-Square distribution */
  150.  
  151.     cexp = totalc / (binary ? 2.0 : 256.0);  /* Expected count per bin */
  152.     for (i = 0; i < (binary ? 2 : 256); i++) {
  153.        prob[i] = (double) ccount[i] / totalc;
  154.        a = ccount[i] - cexp;
  155.        chisq = chisq + (a * a) / cexp;
  156.        datasum += ((double) i) * ccount[i];
  157.     }
  158.  
  159.     /* Calculate entropy */
  160.  
  161.     for (i = 0; i < (binary ? 2 : 256); i++) {
  162.        if (prob[i] > 0.0) {
  163.       ent += prob[i] * log2(1 / prob[i]);
  164.        }
  165.     }
  166.  
  167.     /* Calculate Monte Carlo value for PI from percentage of hits
  168.        within the circle */
  169.  
  170.     montepi = 4.0 * (((double) inmont) / mcount);
  171.  
  172.     /* Return results through arguments */
  173.  
  174.     *r_ent = ent;
  175.     *r_chisq = chisq;
  176.     *r_mean = datasum / totalc;
  177.     *r_montepicalc = montepi;
  178.     *r_scc = scc;
  179. }
  180.