random.c

#include <math.h>

#define IA 16807 
#define IM 2147483647 
#define AM (1.0/IM) 
#define IQ 127773 
#define IR 2836 
#define NTAB 32 
#define NDIV (1+(IM-1)/NTAB) 
#define EPS 1.2e-7 
#define RNMX (1.0-EPS) 

float ran1(long *idum)  
     /* Minimal  random number generator of Park and Miller with Bays-Durham shuffle and added safeguards.
    Returns a uniform random deviate between 0.0 and 1.0 (exclusive of the endpoint values). Call with idum 
    a negative integer to initialize; thereafter, do not alter idum between successive deviates in a sequence. 
    RNMX should approximate the largest floating value that is less than 1.
    Copyright: Numerical Recipes in C */               
{ 
  int j; 
  long k; 
  static long iy=0; 
  static long iv[NTAB]; 
  float temp; 
  
  if (*idum <= 0 || !iy) {       /* Initialize. */ 
    if (-(*idum) < 1) 
      *idum=1;			/* Be sure to prevent idum = 0. */ 
    else 
      *idum = -(*idum); 
    for (j=NTAB+7;j>=0;j--) {	/* Load the shuffle table (after 8 warm-ups). */
      k=(*idum)/IQ; 
      *idum=IA*(*idum-k*IQ)-IR*k; 
      if (*idum < 0) 
	*idum += IM; 
      if (j < NTAB) 
	iv[j] = *idum; 
    } 
    iy=iv[0]; 
  } 
  k=(*idum)/IQ;			/* Start here when not initializing. */
  *idum=IA*(*idum-k*IQ)-IR*k;	/* Compute idum=(IA*idum) % IM without overflows by Schrage's method. */
  if (*idum < 0) 
    *idum += IM; 
  j=iy/NDIV;			/* Will be in the range 0..NTAB-1. */
  iy=iv[j];			/* Output previously stored value and refill the shuffle table. */
  iv[j] = *idum; 
  if ((temp=AM*iy) > RNMX) 
    return RNMX;		/* Because users don t expect endpoint values. */
  else 
    return temp; 
}
	
void rand_perm(int n, long *idum, int *perm)
     /* Produces a random permutation of the elements 1,...,n using n swaps 
        perm is required to be an NR style ivector                          */
{
  int i, r, temp;

  for (i=1; i<=n; i++) 
        perm[i] = i;
  
  for (i=1; i<=n-1; i++) {
    /* swap i-th elt. with a random elt. from i to n */
      temp = perm[i];
      r = (int)floor((float)(n-i+1) * ran1(idum)) + i;
      perm[i] = perm[r];
      perm[r] = temp;
  }
}