isisd/topology/random.c

   1 /*********************************************************************/
   2 /*                                                                   */
   3 /* current processor time in seconds                                 */
   4 /* difference between two calls is processor time spent by your code */
   5 /* needs: <sys/types.h>, <sys/times.h>                               */
   6 /* depends on compiler and OS                                        */
   7 /*                                                                   */
   8 /*********************************************************************/
   9
  10 #include <sys/types.h>
  11 #include <sys/times.h>
  12
  13 /*
  14  * Prototypes.
  15  */
  16 unsigned long timer(void);
  17 void init_rand(long);
  18 double rand01(void);
  19 double randg01(void);
  20 long nrand(long);
  21 void free_arc(void *);
  22
  23 unsigned long timer ()
  24    { struct tms hold;
  25
  26         times(&hold);
  27         return (unsigned long) ((float) (hold.tms_utime) / 60.0);
  28    }
  29
  30
  31 /*********************************************************************/
  32 /*                                                                   */
  33 /*            Family of random number generators                     */
  34 /*                                                                   */
  35 /*  Initialisation:                                                  */
  36 /*            void init_rand ( seed );                               */
  37 /*                 long seed     - any positive number               */
  38 /*                                 if seed<=0 init_rand takes time   */
  39 /*                                 from timer instead of seed        */
  40 /*                                                                   */
  41 /*  Whole number uniformly distributed on [0,n):                     */
  42 /*            long nrand (n);                                        */
  43 /*                 long n                                            */
  44 /*                                                                   */
  45 /*  Real number uniformly distributed on [0,1]                       */
  46 /*            double rand01();                                       */
  47 /*                                                                   */
  48 /*  Real number with Gauss(0,1) disitribution:                       */
  49 /*            double randg01();                                      */
  50 /*                                                                   */
  51 /*  Algorithm:                                                       */
  52 /*            x(n+1) = (x(n) * 5^13) mod 2^31                        */
  53 /*                                                                   */
  54 /*********************************************************************/
  55
  56 unsigned long internal_seed;
  57
  58 void init_rand ( init_seed )
  59
  60 long init_seed;
  61
  62 { internal_seed = ( init_seed > 0 )
  63                  ? (unsigned long) init_seed
  64                  : (unsigned long) timer();
  65
  66
  67   /* only odd numbers are acceptable */
  68   if ( internal_seed % 2 == 0 ) internal_seed --;
  69 }
  70
  71 /*********************************************************************/
  72 /*                                                                   */
  73 /* Internal function  irand  may depend on OS and compiler           */
  74 /*                                                                   */
  75 /* irand assumption:                                                 */
  76 /* unsigned long i,j;                                                */
  77 /*   if i*j > max(unsigned long)                                     */
  78 /*           1. No overflow interruption                             */
  79 /*           2. i*j = i*j mod max(unsigned long)                     */
  80 /*                                                                   */
  81 /* This assumption is true for a lot of computers.                   */
  82 /* If your computer fails:                                           */
  83 /*   rename: irand <---> xrand                                       */
  84 /*                                                                   */
  85 /*********************************************************************/
  86
  87 #define  A   1220703125
  88 #define  B   2147483647
  89 #define  BF  2147483647.
  90
  91 static long irand ()
  92
  93 { internal_seed = ( internal_seed * A ) & B;
  94   return (long) internal_seed ;
  95 }
  96
  97 #if 0 /* Not used. */
  98 /*********************************************************************/
  99 /*                                                                   */
 100 /*              computer independent variant of  irand               */
 101 /*                                                                   */
 102 /*********************************************************************/
 103
 104
 105 #define T15  32768
 106 #define T16  65536
 107 #define A1   37252
 108 #define A2   29589
 109
 110 static long xrand()
 111
 112 { unsigned long is1, is2;
 113
 114   is1 = internal_seed / T15;
 115   is2 = internal_seed % T15;
 116
 117   internal_seed = ( (((is2 * A1) + (is1 * A2))% T16 )* T15 + (is2 * A2) ) & B;
 118   return (long) ( internal_seed ) ;
 119 }
 120 #endif
 121
 122 /*********************************************************************/
 123
 124
 125 double rand01()
 126
 127 { return  (double) (irand() / BF) ;
 128 }
 129
 130 /*********************************************************************/
 131
 132 #define NK  12
 133
 134 double randg01()
 135
 136 { int i;
 137   double sum = 0;
 138
 139   for ( i = 0; i < NK; i++ ) sum += rand01();
 140   return sum - 6.;
 141
 142   /* if   NK != 12  then you must return (12/NK)*sum - (NK/2) */
 143 }
 144
 145 #undef NK
 146
 147
 148 /*********************************************************************/
 149
 150 long nrand ( n )
 151
 152 long n;
 153
 154 { return (long) ( rand01() * (double) n );
 155 }
 156
 157 /*********************************************************************/
 158
 159 #undef A
 160 #undef A1
 161 #undef A2
 162 #undef B
 163 #undef BF
 164 #undef T15
 165 #undef T16