/*  "prng.scm" schlep compiles to prng.c for comparison with Scheme. */


#include "math.h"
#include "stdio.h"
#include "stdlib.h"

int integer_length(n)
     int n;
{
  switch (n) {
  case 0:
    return 0;
  case 1:
    return 1;
  case 2:
  case 3:
    return 2;
  case 4:
  case 5:
  case 6:
  case 7:
    return 3;
  default:
    return 4+(integer_length((n)>>4));
  }
}


int random_chunk(sta)
     unsigned char *sta;
{
  if (0 < (sta[0x102]))
    {
      sta[0x102] = 0;
      printf("random state called reentrantly\n");
      exit(1);
    }
  sta[0x102] = 1;
  {
    int idx = 0xff&(1+(sta[0x100]));
    int xtm = sta[idx];
    int idy = 0xff&((sta[0x101])+(xtm));
    sta[0x100] = idx;
    sta[0x101] = idy;
    {
      int ytm = sta[idy];
      sta[idy] = xtm;
      sta[idx] = ytm;
      {
	int ans = sta[0xff&((ytm)+(xtm))];
	sta[0x102] = 0;
	return ans;
      }
    }
  }
}


int random_random(modu, state)
     int modu;
     unsigned char *state;
{
  int bitlen = integer_length(-1+(modu));
  goto L_rnd;
L_rnd:
  {
    int bln = bitlen;
    int rbs = 0;
    while (!((bln)<=7)) {
      {
	bln = -8+(bln);
	rbs = ((rbs)<<8)+(random_chunk(state));
      }
    }
    if (0 < (bln))
      rbs = ((rbs)<<(bln))^(random_chunk(state));
    if ((rbs)<(modu))
      return rbs;
    else goto L_rnd;
  }
}


void seed_random_state(sta, seed)
     unsigned char *sta;
     unsigned char *seed;
{
  {
    int idx = 0xff;
    while (!(0 > (idx))) {
      sta[idx] = idx;
      {
	idx = -1+(idx);
      }
    }
  }
  {
    int i = 0;
    int j = 0;
    int seed_len = sizeof(seed)-1;
    int k = 0;
    while (!((i)>=0x100)) {
      {
	int swp = sta[i];
	k = 0xff&((k)+((((unsigned char*)(seed))[j]))+(swp));
	sta[i] = sta[k];
	sta[k] = swp;
      }
      {
	i = 1+(i);
	j = (1+(j))%(seed_len);
      }
    }
  }
}


double square_diff(x, y)
     int x;
     int y;
{
  int z = (x)-(y);
  return (z)*(z);
}


void prng(samples, modu, sta)
     int samples;
     int modu;
     unsigned char *sta;
{
  unsigned *sra = malloc((samples) * (sizeof (long)));
  {
    int cnt = -1+(samples);
    int num = random_random(modu, sta);
    int sum = 0;
    while (!(0 > (cnt))) {
      sra[cnt] = num;
      {
	int T_num = random_random(modu, sta);
	cnt = -1+(cnt);
	sum = (sum)+(num);
	num = T_num;
      }
    }
    sum = (sum)+(num);
    {
      double mean = (sum)/(samples);
      {
	int cnt = -1+(samples);
	int var2 = 0;
	while (!(0 > (cnt))) {
	  {
	    int T_cnt = -1+(cnt);
	    var2 = (square_diff(mean, sra[cnt]))+(var2);
	    cnt = T_cnt;
	  }
	}
	printf("%d / %d = %g +/- %g\n", sum, samples, mean, sqrt((var2)/(samples)));
      }
    }
  }
}


int main(argc, argv)
     int argc;
     char **argv;
{
  int n = (argc)>1
    ?atoi(argv[1])
    :0x3e8;
  unsigned char *sta = calloc(3+0x100, (sizeof (char)));
  seed_random_state(sta,
	 "http://swissnet.ai.mit.edu/~jaffer/SLIB.html");
  prng(n, 0x3e7, sta);
  return 0;
}