home *** CD-ROM | disk | FTP | other *** search
/ Avalon - 3D Objects & Resources / Avalon.iso / objects / obj / terrain / grn.for < prev    next >
Encoding:
Text File  |  1995-01-01  |  5.1 KB  |  219 lines
  1.  
  2.  
  3.       FUNCTION GRN(KSEED)
  4.  
  5. C  GRN - Gaussian (Normal) Random Number Generator
  6. C
  7. C  Returns a normal random number with mean zero
  8. C  and standard deviation one.
  9. C
  10. C  Algorithm 488:  Collected Algorithms ACM
  11. C
  12. C  Reference:
  13. C
  14. C    Brent: A Gaussian pseudo-random number
  15. C    generator, Comm ACM, vol 17 #12, Dec 1974,
  16. C    704-706.
  17. C
  18. C  Implemented:
  19. C
  20. C    R Kamen, Code 3873, Sept 81
  21. C
  22. C  Modified:
  23. C
  24. C    K Lindblom, Code 3873, per
  25. C    L Lucas, Code 3873, 21 Jan 82
  26. C
  27. C  GRN is Brent's GRAND modified as follows:
  28. C
  29. C    1. The dummy argument N is replaced by
  30. C       the argument KSEED, which is used to
  31. C       hold the 'memory' of the uniform
  32. C       generator URN.
  33. C
  34. C    2. References to RAND(0) are replaced by
  35. C       references to URN(KSEED).
  36. C
  37. C    3. Dimension of D is reduced from 60 to 30.
  38. C
  39. C  Except on the first call to GRN...
  40. C  GRN returns a pseudo-random number having a
  41. C  Gaussian (i.e. normal) distribution with zero
  42. C  mean and unit variance.  Thus, the probability
  43. C  density is
  44. C
  45. C    f(x) = exp(- x**2 / 2.) / sqrt(2.*pi)
  46. C
  47. C  The first call initializes GRN and returns 0.
  48. C
  49. C  When GRN is called, KSEED must be the name of
  50. C  an integer variable in the calling program.
  51. C  Initially, KSEED must be an integer between
  52. C  1 and 2147483646.  KSEED should not be changed
  53. C  between calls to GRN.
  54. C
  55. C  GRN assumes that successive calls to URN yield
  56. C  independent, pseudo-random numbers distributed
  57. C  uniformly on (0,1) - possibly including 0 but
  58. C  NOT 1.
  59. C
  60. C  The method used was suggested by von Newmann
  61. C  and improved by Forsythe, Ahrens, Dieter and
  62. C  Brent.
  63. C
  64. C  On the average, there are 1.37746 calls to URN
  65. C  for each call to GRN.
  66. C
  67. C
  68. C
  69. C  Warning - dimensions and data statements below
  70. C            are machine-dependent.
  71. C
  72. C  Dimension of D must be at least the number of
  73. C  bits in the fraction of a floating-point number.
  74. C  Thus, on most machines the DATA statement below
  75. C  may be truncated.
  76. C
  77. C  If the integral of sqrt(2./pi)*exp(- x**2 / 2)
  78. C  from A(I) to infinity is 2**(-I), then
  79. C
  80. C       D(I) = A(I) - A(I-1)
  81. C
  82.       EXTERNAL URN
  83. c
  84. c     DIMENSION D(60)
  85.       DIMENSION  D(30)
  86.       DATA D(1)   /0.674489750/
  87.       DATA D(2)   /0.475859630/
  88.       DATA D(3)   /0.383771164/
  89.       DATA D(4)   /0.328611323/
  90.       DATA D(5)   /0.291142827/
  91.       DATA D(6)   /0.263684322/
  92.       DATA D(7)   /0.242508452/
  93.       DATA D(8)   /0.225567444/
  94.       DATA D(9)   /0.211634166/
  95.       DATA D(10)  /0.199924267/
  96.       DATA D(11)  /0.189910758/
  97.       DATA D(12)  /0.181225181/
  98.       DATA D(13)  /0.173601400/
  99.       DATA D(14)  /0.166841909/
  100.       DATA D(15)  /0.160796729/
  101.       DATA D(16)  /0.155349717/
  102.       DATA D(17)  /0.150409384/
  103.       DATA D(18)  /0.145902577/
  104.       DATA D(19)  /0.141770033/
  105.       DATA D(20)  /0.137963174/
  106.       DATA D(21)  /0.134441762/
  107.       DATA D(22)  /0.131172150/
  108.       DATA D(23)  /0.128125965/
  109.       DATA D(24)  /0.125279090/
  110.       DATA D(25)  /0.122610883/
  111.       DATA D(26)  /0.120103560/
  112.       DATA D(27)  /0.117741707/
  113.       DATA D(28)  /0.115511892/
  114.       DATA D(29)  /0.113402349/
  115.       DATA D(30)  /0.111402720/
  116. c     DATA D(31)  /0.109503852/
  117. c     DATA D(32)  /0.107697617/
  118. c     DATA D(33)  /0.105976772/
  119. c     DATA D(34)  /0.104334841/
  120. c     DATA D(35)  /0.102766012/
  121. c     DATA D(36)  /0.101265052/
  122. c     DATA D(37)  /0.099827234/
  123. c     DATA D(38)  /0.098448282/
  124. c     DATA D(39)  /0.097124309/
  125. c     DATA D(40)  /0.095851778/
  126. c     DATA D(41)  /0.094627461/
  127. c     DATA D(42)  /0.093448407/
  128. c     DATA D(43)  /0.092311909/
  129. c     DATA D(44)  /0.091215482/
  130. c     DATA D(45)  /0.090156838/
  131. c     DATA D(46)  /0.089133867/
  132. c     DATA D(47)  /0.088144619/
  133. c     DATA D(48)  /0.087187293/
  134. c     DATA D(49)  /0.086260215/
  135. c     DATA D(50)  /0.085361834/
  136. c     DATA D(51)  /0.084490706/
  137. c     DATA D(52)  /0.083645487/
  138. c     DATA D(53)  /0.082824924/
  139. c     DATA D(54)  /0.082027847/
  140. c     DATA D(55)  /0.081253162/
  141. c     DATA D(56)  /0.080499844/
  142. c     DATA D(57)  /0.079766932/
  143. c     DATA D(58)  /0.079053527/
  144. c     DATA D(59)  /0.078358781/
  145. c     DATA D(60)  /0.077681899/
  146. C
  147. C  End of machine-dependent statements.
  148. C
  149. C  U must be preserved between calls.
  150. C
  151.        DATA        U / 0.0 /
  152. C
  153. C  Initialize displacement A and counter I
  154. C
  155.        A = 0.0
  156.        I = 0
  157. C
  158. C  Increment counter and displacement if
  159. C  leading bit of U is one.
  160. C
  161.    10  U = U + U
  162.        IF ( U .LT. 1.0 )  GOTO 20
  163.        U = U - 1.0
  164.        I = I + 1
  165.        A = A - D(I)
  166.        GOTO 10
  167. C
  168. C  Form W uniform on 0 .LE. W .LT. D(I+1) from U
  169. C
  170.    20  W = D(I+1) * U
  171. C
  172. C  Form V = 0.5 * ((W - A)**2 - A**2)
  173. C
  174. C  Note:             0 .LE. V .LT. LOG(2)
  175. C
  176.        V = W * (0.5*W - A)
  177. C
  178. C  Generate new uniform U
  179. C
  180.    30  U = URN(KSEED)
  181. C
  182. C  Accept W as a random sample if V .LE. U
  183. C
  184.        IF ( V .LE. U )  GOTO 40
  185. C
  186. C  Generate random V
  187. C
  188.        V = URN(KSEED)
  189. C
  190. C  Loop if U .GT. V
  191. C
  192.        IF ( U .GT. V )  GOTO 30
  193. C
  194. C  Reject W and form a new uniform U from V and U
  195. C
  196.        U = (V - U) / (1.0 - U)
  197.        GOTO 20
  198. C
  199. C  Form new U (to be used on next call) from U and V
  200. C
  201.    40  U = (U - V) / (1.0 - V)
  202. C
  203. C  Use first bit of U for sign, return normal variate
  204. C
  205.        U = U + U
  206.        IF ( U .LT. 1.0 )  GOTO 50
  207.        U = U - 1.0
  208.        GRN = W - A
  209.        RETURN
  210. C
  211.    50  GRN = A - W
  212.        RETURN
  213.        END
  214.  
  215.  
  216.  
  217.  
  218.  
  219.