/*********************************************************************/
/*                                                                   */
/* current processor time in seconds                                 */
/* difference between two calls is processor time spent by your code */
/* needs: <sys/types.h>, <sys/times.h>                               */
/* depends on compiler and OS                                        */
/*                                                                   */
/*********************************************************************/

#include <sys/types.h>
#include <sys/times.h>

float timer()
   { struct tms hold;

        times(&hold);
        return  (float)(hold.tms_utime) / 60.0;
   }


/*********************************************************************/
/*                                                                   */
/*            Family of random number generators                     */
/*                                                                   */
/*  Initialisation:                                                  */
/*            void init_rand ( seed );                               */
/*                 long seed     - any positive number               */
/*                                 if seed<=0 init_rand takes time   */
/*                                 from timer instead of seed        */
/*                                                                   */
/*  Whole number uniformly distributed on [0,n):                     */
/*            long nrand (n);                                        */
/*                 long n                                            */
/*                                                                   */
/*  Real number uniformly distributed on [0,1]                       */
/*            double rand01();                                       */
/*                                                                   */
/*  Real number with Gauss(0,1) disitribution:                       */
/*            double randg01();                                      */
/*                                                                   */
/*  Algorithm:                                                       */
/*            x(n+1) = (x(n) * 5^13) mod 2^31                        */
/*                                                                   */
/*********************************************************************/

unsigned long internal_seed;

void init_rand ( init_seed )

long init_seed;

{ internal_seed = ( init_seed > 0 )
                 ? (unsigned long) init_seed
                 : (unsigned long) timer(); 
                

  /* only odd numbers are acceptable */
  if ( internal_seed % 2 == 0 ) internal_seed --;
}

/*********************************************************************/
/*                                                                   */
/* Internal function  irand  may depend on OS and compiler           */
/*                                                                   */
/* irand assumption:                                                 */
/* unsigned long i,j;                                                */
/*   if i*j > max(unsigned long)                                     */
/*           1. No overflow interruption                             */
/*           2. i*j = i*j mod max(unsigned long)                     */
/*                                                                   */
/* This assumption is true for a lot of computers.                   */
/* If your computer fails:                                           */
/*   rename: irand <---> xrand                                       */
/*                                                                   */
/*********************************************************************/
 
#define  A   1220703125
#define  B   2147483647
#define  BF  2147483647.

static long irand ()

{ internal_seed = ( internal_seed * A ) & B;
  return (long) internal_seed ;
}

/*********************************************************************/
/*                                                                   */
/*              computer independent variant of  irand               */
/*                                                                   */
/*********************************************************************/


#define T15  32768 
#define T16  65536
#define A1   37252
#define A2   29589

static long xrand()

{ unsigned long is1, is2;

  is1 = internal_seed / T15;
  is2 = internal_seed % T15;

  internal_seed = ( (((is2 * A1) + (is1 * A2))% T16 )* T15 + (is2 * A2) ) & B;
  return (long) ( internal_seed ) ;
}


/*********************************************************************/


double rand01()

{ return  (double) irand() / BF ;
}
  
/*********************************************************************/

#define NK  12

double randg01()

{ int i;
  double sum = 0;

  for ( i = 0; i < NK; i++ ) sum += rand01();
  return sum - 6.;

  /* if   NK != 12  then you must return (12/NK)*sum - (NK/2) */
}

#undef NK
     
  
/*********************************************************************/

long nrand ( n )

long n;

{ return (long) ( rand01() * (double) n );
}

/*********************************************************************/

#undef A
#undef A1
#undef A2
#undef B
#undef BF
#undef T15
#undef T16