#include <stdio.h> 
#include <stdlib.h> 
#include <math.h>
//#include "gaussian.h"
 
// Returns the next pseudorandom, Gaussian ("normally") distributed double 
// value with mean 0.0 and standard deviation 1.0 from this random 
// number generator's sequence.

// This uses the polar method of G. E. P. Box, M. E. Muller, and G. Marsaglia,
// as described by Donald E. Knuth in The Art of Computer Programming, Volume 2:
// Seminumerical Algorithms, section 3.4.1, subsection C, algorithm P. 

// Note that it generates two independent values at the cost of only one call 
// to log() and one call to sqrt().

unsigned char haveNext = 0;
double savedGaussian;

void setGaussianSeed(unsigned int seed)
{ srand(seed);
}
 
double nextGaussian() 
{ if (haveNext) 
  { haveNext = 0;
    return savedGaussian;
  } 

  else 
  { //printf("find two new gaussians\n");
    double v1 = 0.0;
    double v2 = 0.0;
    double s = 0.0;
    while (s >= 1.0 || s == 0.0)
    { v1 = (2.0 * rand())/RAND_MAX - 1.0;   // between -1.0 and 1.0
      v2 = (2.0 * rand())/RAND_MAX - 1.0;   // between -1.0 and 1.0
      //printf("%f %f %f\n", v1, v2, s);
      s = v1 * v1 + v2 * v2;
    }
    double multiplier = sqrt(-2.0 * log(s)/s);
    savedGaussian = v2 * multiplier;
    haveNext = 1;
    return v1 * multiplier;
  }
}

int roundToInt(double number)
{ if (number >= 0)  return (int)(number + 0.5); 
  return (int)(number - 0.5);
}




int test_main()
{ int i;
  int positiveBins = 7;
  int binCount = 2*positiveBins+1;
  int bins[binCount];
  for (i=0; i<binCount; i++)
  { bins[i] = 0; 
  }

  setGaussianSeed(400);

  for (i=0; i<1000000; i++)
  { //printf("%f\n",nextGaussian());
    double r = nextGaussian();
    int idx = roundToInt(r) + positiveBins;
    if (idx >=0 && idx <binCount)
    { //printf("%d\n", idx); 
      bins[idx]++;
    }
  }

  for (i=0; i<binCount; i++)
  { printf("%3d) %d\n", i-positiveBins, bins[i]); 
  }

  printf("%d, %d, %d, %d, %d\n", roundToInt(1.2), roundToInt(1.6), roundToInt(-2.2), roundToInt(-2.8), roundToInt(-0.45));
}