\ r250d ranf r250d_init lcm_rand seed R250 Pseudo-Random number generator \ Forth Scientific Library Algorithm #23 \ algorithm from: \ Kirkpatrick, S., and E. Stoll, 1981; A Very Fast Shift-Register \ Sequence Random Number Generator, Journal of Computational Physics, \ V. 40. p. 517 \ \ see also: \ Maier, W.L., 1991; A Fast Pseudo Random Number Generator, \ Dr. Dobb's Journal, May, pp. 152 - 157 \ Uses the Linear Congruential Method, \ the "minimal standard generator" \ Park & Miller, 1988, Comm of the ACM, 31(10), pp. 1192-1201 \ for initialization \ For a review of BOTH of these generators, see: \ Carter, E.F, 1994; Generation and Application of Random Numbers, \ Forth Dimensions, Vol. XVI, Numbers 1,2 May/June, July/August \ This is an ANS Forth program requiring: \ 1. The Double-Number word set \ 2. The Floating-Point word set (for ranf) \ 3. The word '}' to dereference a one-dimensional array. \ 4. The words umd* umd/mod dxor dor and dand (see alternative note below) \ 5. Uses words 'Private:', 'Public:' and 'Reset_Search_Order' \ to control visibility of internal code \ 6. The Programming-Tools wordset is need to control generator compilation \ options (see note below). \ 7. The compilation of the test code is controlled by the VALUE TEST-CODE? \ and the conditional compilation words in the Programming-Tools wordset \ Note: \ Assumes that double integers are (at least) 32 bits. If regular integers \ are (at least) 32 bits then umd* and umd/mod will not be needed if \ lcm_rand is vectored to use lcm_rand32 (this will also be MUCH faster \ than using lcm_rand16 for such systems). The private word, vector_lcm attempts \ to do this automatically at load time. \ \ Alternatively, one can set the CONSTANT SINGLE-GENERATOR? to TRUE before loading, \ in this case the generator lcm_randB will be used. This will work with both 16 \ and 32 bit systems with a performance penalty of a little more than a factor of \ 2. This performance hit will have a small impact on the use of R250 (since \ lcm_rand is only used to initialize R250), but could be a disadvantage when \ using lcm_rand directly. When compiling this way, umd* and umd/mod are not \ needed. (all three lcm generator implementations give identical results). \ (c) Copyright 1994 Everett F. Carter. Permission is granted \ by the author to use this software for any application provided \ this copyright notice is preserved. \ CR .( R250.SEQ V1.5 16 December 1994 EFC ) Private: FALSE CONSTANT SINGLE-GENERATOR? DECIMAL 2147483647. 2CONSTANT max32 SINGLE-GENERATOR? [IF] \ multiply unsigned double by unsigned single, giving unsigned double result. : mu* ( ud1 u--ud2 ) TUCK * >R UM* R> + ; \ divide unsigned double by unsigned single, giving unsigned single result. : um/ ( ud u--u ) UM/MOD NIP ; \ divide unsigned double by unsigned single, giving unsigned double result. : mu/ ( ud u--ud ) >R 0 R@ UM/MOD R> SWAP >R um/ R> ; \ multiply unsigned single u1 by unsigned single u2, then divide by unsigned \ single u3, giving unsigned double result. : um*/ ( u1 u2 u3--ud ) >R UM* R> mu/ ; \ divide non-negative double +d1 by strictly positive double +d3, \ giving double quotient d3. \ \ The algorithm is described in "Long Divisors and Short Fractions \ by Prof. Nathaniel Grossman, in Forth Dimensions Volume VI No. 3. \ \ Grossman cites Abramowitz M and I A Stegun, Handbook of Mathematical \ Functions, National Bureau of Standards Applied Mathematics Series, 55. \ (Reprinted by Dover Publications) page 21 and Knuth, Seminumerical Algorithms \ as his references. : +d/ ( +d1 +d2--+d3 ) ?DUP IF DUP 1+ 0 1 ROT um/ DUP >R mu* >R OVER SWAP R@ um*/ D- 2R> M*/ NIP 0 ELSE mu/ THEN ; \ multiply unsigned double by unsigned double, giving double result. : ud* ( ud ud--ud ) DUP IF 2SWAP THEN DROP mu* ; \ divide non-negative double +d1 by strictly positive double +d3, \ giving double remainder d3 and double quotient d4. : +d/mod ( +d1 +d2--+d3 +d4 ) 4DUP +d/ 2DUP 2>R ud* D- 2R> ; [THEN] Public: 2VARIABLE seed 1234. seed 2! v: lcm_rand \ Linear Congruential Method, the "minimal standard generator" \ Park & Miller, 1988, Comm of the ACM, 31(10), pp. 1192-1201 SINGLE-GENERATOR? [IF] \ the following implementation works for both 16 and 32 bit systems : lcm_randB ( -- d ) seed 2@ 127773. +d/mod 2>R 16807 mu* 2R> 2836 mu* D- 2DUP D0< IF max32 D+ THEN 2DUP seed 2! ; USE( lcm_randB defines lcm_rand [ELSE] \ for a cell size of at least 16 bits and a double size of at least 32 bits : lcm_rand16 ( -- d ) seed 2@ 16807. umd* max32 umd/mod 2DROP 2DUP seed 2! ; \ the following implementation gives identical results and will \ run faster than lcm_rand16 for a system with a cell size of at least 32 bits : lcm_rand32 ( -- d ) seed 2@ 127773 UM/MOD >R 16807 UM* R> 2836 UM* D- DUP 0< IF max32 D+ THEN 2DUP seed 2! ; Private: \ +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ \ WARNING!!!!!! MAKE SURE THIS IS SET APPROPRIATELY FOR YOUR MACHINE/COMPILER !! : vector_lcm 1 CELLS 4 < IF \ for 16 BIT INTS USE( lcm_rand16 ELSE \ for 32 BIT (or larger) INTS USE( lcm_rand32 THEN defines lcm_rand ; vector_lcm \ +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ [THEN] Private: \ R250 code -- 31 bit version \ Kirkpatrick & Stoll, 1981; Jour. Computational Physics, \ 40, p. 517 2VARIABLE dmask 2VARIABLE dmsb VARIABLE front_index VARIABLE back_index 250 DOUBLE ARRAY r250d_buffer{ Public: : r250d_init ( d -- ) seed 2! 0 front_index ! 103 back_index ! \ initialize the buffer with random values 250 0 DO lcm_rand r250d_buffer{ I } 2! LOOP 250 0 [ HEX ] DO lcm_rand 020000000. D> IF r250d_buffer{ i } DUP >R 2@ 040000000. dor R> 2! THEN LOOP 040000000. dmsb 2! 07fffffff. dmask 2! [ DECIMAL ] 31 0 DO r250d_buffer{ I 7 * 3 + } DUP >R 2@ dmask 2@ dand dmsb 2@ dor R> 2! dmask DUP >R 2@ D2/ R> 2! dmsb DUP >R 2@ D2/ R> 2! LOOP ; : r250d ( -- d ) \ 32 bit positive (i.e. 31 bit) number r250d_buffer{ back_index @ } 2@ r250d_buffer{ front_index @ } DUP >R 2@ dxor 2DUP R> 2! \ save new value at front index location 1 front_index +! front_index @ 249 > IF 0 front_index ! THEN 1 back_index +! back_index @ 249 > IF 0 back_index ! THEN ; : ranf ( --, f: -- x) \ generate a random value from 0.0 to 1.0 r250d D>F max32 D>F F/ ; Reset_Search_Order TEST-CODE? [IF] \ test code ============================================== \ test the LCM generator : lcm_test ( -- ) 1. seed 2! \ set initial seed value 1. \ push a temporary value 10000 0 DO 2DROP lcm_rand LOOP ." final value: " D. ." should be 1043618065" CR ; \ test the R250 generator : r250_test ( -- ) 1. r250d_init 1. 10000 0 DO 2DROP r250d LOOP ." final value: " D. ." should be 267771767" CR ; [THEN]