random.c 5.5 KB

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