 /*
 Code is in the public domain.
 Can you seriously envisage ANYONE paying to use this mess? I can't.
	- Peter Thompson 14 Feb 92
 */

#include <stdio.h>
#include "random.h"

main(int argc, char *argv[])
{
  RandState *randbase;
  long i;
  long j;
  long *counts;
  long N;
  long R;
  ULONG seed;
  double sum = 0;
  double chisqr = 0;

  seed = (ULONG) atol(argv[1]);
  N = atol(argv[2]);
  R = atol(argv[3]);
  printf("Seed:%li N:%li R:%li\n", seed, N, R);
  counts = (long *) calloc(R,sizeof(long int));
  if ((long *) NULL == counts)
    exit(21);
  randbase = (RandState *) malloc(sizeof(RandState));
  if ((RandState *) NULL == randbase)
    exit(21);
  rand32init(randbase, seed);
  for (i = N; (i>0); i--) {
    ULONG result;
    result = rnda32(randbase);
 /* Debug stuff
    for (j = 0; j < 64; j++)
      printf("%8.8x%c",(ULONG)randbase->rs_State[j],(((j+1)&7)==0)? '\n':' ');
    printf("%li %x\n",randbase->rs_Index,2*result);	
 */
    counts[(long)(R*(result/(float) MAXRNDNUM))]++;
  }
  for (i = 0, sum = 0; i < R; i++) {
    sum += counts[i]*counts[i];
    if (i < 128)
      printf("%5.5li%c",counts[i],(((i+1)&7)==0)? '\n':' ');
  }
  chisqr = R*(sum/(double)N)-N;
  if (((chisqr-R)*(chisqr-R)) > (4*R))
    printf("\nFailed: %f\n",chisqr);
  else
    printf("\nPassed: %f\n",chisqr);
}
