random.c

Go to the documentation of this file.

#include "random.h"

#include <sys/types.h>
#include <time.h>
#include <stdlib.h>
#include <math.h>

void randomReset()
{
  srand((unsigned int)time(NULL));
  rand();
  rand();
}

void randomSet(unsigned int seed)
{
  srand(seed);
  rand(); /*useful to avoid first rands to be almost equals*/
  rand();

}

double randomDouble()
{
  return((double)rand()/(double)RAND_MAX);
}

double randomNormal()
{
  static boolean oneCached=FALSE;
  double x1,x2,w;
  double y1;
  static double y2;

  if (oneCached)
    {
      oneCached=FALSE;
      return(y2);
    }
  else
    {
      do {
        x1=2.0*randomDouble()-1.0;
        x2=2.0*randomDouble()-1.0;
        w=x1*x1+x2*x2;
      } while (w>=1.0);
      
      w=sqrt((-2.0*log(w))/w);
      y1=x1*w;
      y2=x2*w;
 
      oneCached=TRUE;
      return(y1);
    }
}

double randomInInterval(Tinterval *t) 
{
  double a,b;

  a=LowerLimit(t);
  b=UpperLimit(t);

  return(a+randomDouble()*(b-a));
}

unsigned int randomMax(unsigned int m)
{
  if (m==0)
    return(0);
  else
    return((unsigned int)floor(randomDouble()*(m+1))); /* in [0,m+1) -> all values in 0,m have same chances */
}


unsigned int randomWithDistribution(unsigned int m,double s,double *d)
{
  double s1,r;
  unsigned int i;

  if (s==0.0)
    {
      for(i=0;i<m;i++)
        s+=d[i];
    }
 
  r=s*randomDouble();

  i=0;
  s1=d[0];
  while((r>s1)&&(i<(m-1)))
    {
      i++;
      s1+=d[i];
    }

  return(i);
}

void randomOnBall(double r,unsigned int k,double *p)
{
  unsigned int i;
  double n;

  n=0.0;
  for(i=0;i<k;i++)
    {
      p[i]=randomNormal();
      n+=(p[i]*p[i]);
    }
  n=r/sqrt(n);
  for(i=0;i<k;i++)
    p[i]*=n;
}

void randomInBall(double r,unsigned int k,double *p)
{
  unsigned int i;
  double n;
  
  /* random on surface ... */
  n=0.0;
  for(i=0;i<k;i++)
    {
      p[i]=randomNormal();
      n+=(p[i]*p[i]);
    }
  /* ... but scaling to select points inside the ball too. */
  n=pow(randomDouble(),1/((double)k))*(r/sqrt(n)); 
  for(i=0;i<k;i++)
    p[i]*=n;
}