#include "randomizer.h"


/* --------------------------------------------------------------------------------------------- *
 *                                                                                               *
 *    Here are some random number generators                                                     *
 *                                                                                               *
 * --------------------------------------------------------------------------------------------- */


/* --------------------------------------------------------------------------------------------- */


float flat_random()   /* flat random distribution in interval [0,1[ */
{
  int r ;
  r = random() ;
  return (float)( (double)r / ((double)RAND_MAX+1) ) ;
}


/* --------------------------------------------------------------------------------------------- */


float exp_random(float mean)   /* exponential random distribution with average <mean>*/
{
  float r ;
  r = flat_random() ;
  if (r>=1)  return 9.999E+27 ;  else  return -mean*log(1-r) ;
}


/* --------------------------------------------------------------------------------------------- */


point spher_random()   /* 3-dim points uniformly distributed on a unitary sphere */
{
  point P ;
  float cos_theta, sin_theta, phi ;
  cos_theta = 2*(flat_random()-0.5) ;
  sin_theta = sqrt(1-cos_theta*cos_theta) ;
  phi = flat_random()*2*PI ;
  P.z = cos_theta ;
  P.x = sin_theta*cos(phi) ;
  P.y = sin_theta*sin(phi) ;
  return P ;
}


/* --------------------------------------------------------------------------------------------- */


float gauss_random(float mean, float sigma)   /* gaussian random numbers with given mean and sigma */
{
  float u , r , phi , x , y ;
  u = exp_random(1.) ;
  r = sqrt(u) ;
  phi = 2*PI*flat_random() ;
  x = r*cos(phi) ;   /* this is gaussian with mean=0 and sigma=1/sqrt(2) */
  y = r*sin(phi) ;   /* this is gaussian with mean=0 and sigma=1/sqrt(2) */
  /* x and y are gaussian with mean=0 and sigma=1/sqrt(2)
     ==> x+y is gaussian with mean=0 and sigma=1 */
  return (x+y) * sigma + mean ;
}


/* --------------------------------------------------------------------------------------------- */


point gauss_3d_random(point mean, float sigma)
{
  point p ;
  p.x = gauss_random(mean.x,sigma) ;
  p.y = gauss_random(mean.y,sigma) ;
  p.z = gauss_random(mean.z,sigma) ;
  return p ;
}


/* --------------------------------------------------------------------------------------------- */


float chi2_random(int Nsamples)
{
  float x , r ;
  int i ;
  if ( Nsamples <= 0 )  return 9999999. ;
  for ( i=0 ; i<Nsamples ; i++ )
    {
      r = gauss_random(0.,1.) ;
      x += r*r ;
    }
  x /= (float) Nsamples ;
  return x ;
}


/* --------------------------------------------------------------------------------------------- */


int binom_random(int N, float p)
{
  int n, k ;
  if ( N == 0 ) return 0 ;
  n = 0 ;
  for ( k=0 ; k<N ; k++ )
    {
      if ( flat_random()<p )  n++ ;
    }
  return n ;
}


/* --------------------------------------------------------------------------------------------- */