/* Measure the quality of random numbers for different seed values.

   This file exists because I found that the rand() that comes with
   MingW appears to be absymal.  Basically the first output of rand()
   is not very random at all with respect to seed value such that
   consecutive seeds will likely yield similar random numbers, and
   subsequent calls yields only 15 bits of random numbers.

   Discovery of this poor rand() is why I have the rand_patch.h hack.
                                                                      */
#include<math.h>
#include<stdio.h>
#include<stdlib.h>
#include<time.h>

#define SEED_COUNT 1000
#define RAND_COUNT 5

static int rand_values[RAND_COUNT][SEED_COUNT];

static int CmpDouble(const void *a, const void *b)
{
   double x = *(double*)a;
   double y = *(double*)b;
   return x < y ? -1 : x > y ? 1 : 0;
}

static void OutputStatistics(double *values, int count, double expected_range)
{
   double mean, stddev, d;
   int i;

   qsort(values, count, sizeof(double), CmpDouble);
   mean = *values;
   for(i = 1; i < count; i++)
      mean += values[i];
   mean /= count;

   stddev = 0;
   for(i = 0; i < count; i++)
   {
      d = values[i] - mean;
      stddev += d * d;
   }
   stddev = sqrt(stddev / count);

   printf("%.5f  range=%g, min=%g, max=%g, mean=%g, "
          "median=%g, stddev=%g\n",
          stddev / expected_range,
          values[count - 1] - values[0],
          values[0],
          values[count - 1],
          mean,
          values[count / 2],
          stddev);
}

int main(void)
{
   int s, r, t;
   double table[SEED_COUNT];

   /* Collect random values. */
   t = time(NULL);
   for(s = 0; s < SEED_COUNT; s++)
   {
      srand(t + s);
      for(r = 0; r < RAND_COUNT; r++)
         rand_values[r][s] = rand();
   }

   /* Output statistics about random values according to their position.
      That is, rand[0] are stats for the first return value from rand()
      for some number of seeds, rand[1] is the second, and so on.        */
   for(r = 0; r < RAND_COUNT; r++)
   {
      for(s = 0; s < SEED_COUNT; s++)
         table[s] = (double)rand_values[r][s] / RAND_MAX;
      printf("rand[%d]:    ", r);
      OutputStatistics(table, SEED_COUNT, 1.0);

      for(s = 0; s < SEED_COUNT; s++)
         table[s] = (double)((rand_values[r][s] >> 16) & 0xffff);
      printf("rand[%d].hi: ", r);
      OutputStatistics(table, SEED_COUNT, 0xffff);

      for(s = 0; s < SEED_COUNT; s++)
         table[s] = (double)(rand_values[r][s] & 0xffff);
      printf("rand[%d].lo: ", r);
      OutputStatistics(table, SEED_COUNT, 0xffff);

      putchar('\n');
   }

   printf("RAND_MAX = %d = 0x%x\n", RAND_MAX, RAND_MAX);
   return 0;
}