home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 3 / PDCD_3.iso / utilities / utilsf / hugs / hs / minsrand < prev    next >
Text File  |  1995-02-14  |  4KB  |  126 lines
  1. -------------------------------------------------------------------------------
  2. -- The following random number generator is an implementation of the
  3. -- Minimum Standard generator recommended in
  4. --
  5. --    Random Number Generators: Good ones are hard to find
  6. --       Stephen K Park & Keith W Miller
  7. --       Communications of the ACM, Oct 88, Vol 31 No 10 1192 - 1201
  8. --
  9. -- Seeds must be in the range 1..2147483646, that is (1..(2**31)-2)
  10. -- Output will also be in that range. The generator is full period so that
  11. -- all 2147483646 values will be generated before the initial seed repeats.
  12. -- Dividing by 2147483647 (real) as in the Pascal code below will map it
  13. -- into the range (0..1) if required.
  14. --
  15. -- [This program assumes that you are working on a machine with (at least)
  16. -- 32 bit integers.  Folks using Hugs on a PC will have to stick with the
  17. -- less sophisticated random number generator in the file `randoms'.]
  18. -------------------------------------------------------------------------------
  19.  
  20. min_stand_test  :: Int -> Int
  21. min_stand_test n = if test > 0 then test else test + 2147483647
  22.            where test = 16807 * lo - 2836 * hi
  23.                  hi   = n `div` 127773
  24.                  lo   = n `rem` 127773
  25.  
  26. min_stand_randoms :: Int -> [Int]
  27. min_stand_randoms  = iterate min_stand_test
  28.  
  29. -- The article produced below also gives a test to check that the
  30. -- random number generator is working.  We can duplicate this test
  31. -- as follows:
  32. --
  33. --   ? strictIterate min_stand_test 1 !! 10000
  34. --   1043618065
  35. --   (149758 reductions, 240096 cells, 2 garbage collections)
  36. --
  37. -- Happily, this is the result that we expect to obtain.
  38. --
  39. -- The function strictIterate is defined below.  It is similar to the
  40. -- standard iterate function except that it forces the evaluation of
  41. -- each element in the list produced (except possibly the first).
  42. -- Had we instead tried to evaluate:
  43. --
  44. --   iterate min_stand_test 1 !! 10000
  45. --
  46. -- Hugs would have first constructed the expression graph:
  47. --
  48. --   min_stand_test (min_stand_test (... (min_stand_test 1) ...))
  49. --
  50. -- in which the min_stand_test function is applied 10000 times to 1
  51. -- and then attempted to evaluate this.  In either case, you'd need a
  52. -- large heap to represent the complete expression and a large stack so
  53. -- that you could handle 10000 levels of function calling.  Most standard
  54. -- configurations of Hugs aren't set up with sufficiently large defaults
  55. -- to make this possible, so the most likely outcome would be a runtime
  56. -- error of one kind or another!
  57.  
  58. strictIterate    :: (a -> a) -> a -> [a]
  59. strictIterate f x = x : strict (strictIterate f) (f x)
  60.  
  61. -------------------------------------------------------------------------------
  62. -- Some comments and code from:
  63. --
  64. -- Random Number Generators: Good ones are hard to find
  65. --    Stephen K Park & Keith W Miller
  66. --    Communications of the ACM, Oct 88, Vol 31 No 10 1192 - 1201
  67. -- 
  68. -- Minimum standard random number generator implementations
  69. -- 
  70. -- This version of Random will be correct if reals are represented
  71. -- with a 46-bit or larger mantissa (excluding the sign bit).
  72. -- For example, this version will be correct on all systems that support
  73. -- the IEEE 64-bit real arithmetic standard since the mantissa in that case
  74. -- is 53-bits.
  75. -- ... from page 1195 upper right quadrant
  76. -- 
  77. -- var seed : real;
  78. -- ...
  79. -- function Random : real;
  80. --     (* Real Version 1 *)
  81. -- const
  82. --    a = 16807.0;
  83. --    m = 2147483647.0;
  84. -- var
  85. --    temp : real;
  86. -- begin
  87. --    temp := a * seed;
  88. --    seed :=
  89. --       temp - m * Trunc(temp / m);
  90. --    Random := seed / m;
  91. -- end;
  92. -- 
  93. -- ... from page 1195 lower right quadrant, variant by L. Schrage, 1979, 1983
  94. --
  95. -- var seed : integer;
  96. -- ...
  97. -- function Random : real;
  98. --     (* Integer Version 2 *)
  99. -- const
  100. --    a = 16807;
  101. --    m = 2147483647;
  102. --    q = 127773;    (* m div a *)
  103. --    r = 2836;    (* m mod a *)
  104. -- var
  105. --    lo, hi, test : integer;
  106. -- begin
  107. --    hi := seed div q;
  108. --    lo := seed mod q;
  109. --    test := a * lo - r * hi;
  110. --    if test > 0 then
  111. --       seed := test
  112. --    else
  113. --       seed := test + m;
  114. -- 
  115. --    Random := seed / m;
  116. -- end;
  117. -- 
  118. -- From page 1195 lower left quadrant
  119. --
  120. -- seed := 1;
  121. -- for n := 1 to 10000 do
  122. --    u := Random;
  123. -- Writeln('The current value of seed is : ', seed);
  124. -- (* Expect 1043618065 *)
  125. -------------------------------------------------------------------------------
  126.