LIBRARY OF THE UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN 510.84 The person charging this material is re- sponsible for its return to the library from which it was withdrawn on or before the Latest Date stamped below. Theft, mutilation, and underlining of books are reasons for disciplinary action and may result in dismissal from the University. UNIVERSITY OF ILLINOIS LIBRARY AT URBANA-CHAMPAIGN MAR 2 2 1976 OEC 1 4 1976 DEC 12^ L161 — O-1096 , , UIUCDCS-R-T3-56T fcAJ 5C1 0- y^t^z^ March 1973 CALCULATION OF GFSR PSEUDORANDOM NUMBER BINARY STARTING MATRIX by W. H. Payne %, w V'3 DEPARTMENT OF COMPUTER SCIENCE UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN URBANA, ILLINOIS ILLINOIS UNWERSlJ pA ,r,M UIUCDCS-R-73-567 CALCULATION OF GFSR PSEUDORANDOM NUMBER BINARY STARTING MATRIX by W. H. Payne* March 1973 DEPARTMENT OF COMPUTER SCIENCE UNIVERSITY OF ILLINOIS AT URBAN A-CH AMP AIGN URBANA, ILLINOIS 6l801 * Visiting Research Associate Professor of Computer Science on leave from the Computer Science Department and the Computing Center, Washington State University, Pullman, Washington 99163. Digitized by the Internet Archive in 2013 http://archive.org/details/calculationofgfs567payn 1. INTRODUCTION The generalized feedback shift register (GFSR) pseudorandom number generator produces the 'same' floating-point sequence of multi- dimensional distributed pseudorandom numbers on any computer. The 'same' sequence in the sense that high-order bits agree on all machines. Not only is the period of the pseudorandom sequence independent of the word size of the host machine, but the algorithm can be implemented with one exclusive-or (ffi) and two add (modulo p) instructions. This makes its execution extremely fast for mini- or microcomputers [l]. The pseudorandom number generator is based on a primitive trinomial modulo 2 of the form I + X + 1 and a resulting feedback shift register sequence. Statement of the GFSR algorithm is : GFSR Pseudorandom Number Algorithm Initialize integer J=0 Initialize table M(l) ,M(2) , . . . ,M(P) Step: 1. J -*■ J+l 2. If J > P, set J «- 1 3. K «- J+Q k. If K > P, set K «- 1 5. Exclusive-or, M(j) 9 M(K) 6. STORE, M(J) «- M(J) ® M(K) -2- The procedure of initializing the table, M(l) ,M(2) ,. . . ,M(P) of computer words (up to 98 bits) , is a nontrivial time-consuming task (about five minutes 360/75 CPU time). Specifically, a 98 x 98 binary matrix with linearly independent (modulo 2) columns is needed. Further, a judiciously selected delay of the basic feedback shift register sequence between columns to achieve excellent statistical performance of the generator must be sought. It is the purpose of this paper to present an initial M table, both in octal and hexadecimal, for the polynomial X + X + 1 Q ft (period 2 -l) computed from the initialization algorithm in [l]. For this specific polynomial, L-space distributed numbers to N-bit accuracy where LN < 98 will be produced by the GFSR algorithm. For example, it is possible to generate six-distributed random numbers to 15-bit accuracy since 6 x 15 = 90 < 98. 2. THE TABLE Alignment of both the octal and hexadecimal digits on the 98 bits was set so that the high-order octal digit attains a maximum value of three while the corresponding hexadecimal digit attains a maximum value of seven. This was done to insure that the sign digit was zero for those programmers who wish to implement the algorithm in compiler languages, -3- j r u HI l) M< 2) HI 3) Ml 41 Ml 5) Ml t>) Ml 7) Ml 8) Ml 9) MI10I Mill! Mt 12) HI13I Mil*) MI15I 11161 Mim HUB) MI19I Hlrtl MI21I MI22I M(23> M(2*l M125I MI26) HI27I MI28) M<29) MI3UI M131) MU2I MI33) M134) MI35I MI36I MI37i MI38J MI39) NHll MUD MI42I MI451 Ml**) MUSI KU6I HI47I HUB) MI49) MliOl M13U MI52) MI55I MI54) MISS) MIS61 M(SM MI58) MI59) MI60I MI61) MI62) NI63I MI64I MI6S) Mlbb) MlbT) MI68) MI69) MI70) MI71) M472I MI73) MI74) NI7SI MI76) MITM M(78) MI7V) Mtao) MI81) MI82) MI83J Mid,) Mlasi MIS6I MI87I MI88I MI8VI MI9QI MIV1I HI92) MI93I MIV4I MI94J Ml 96 I MI97) MI98I GFSR PSEUDORANDOM NUMBER GENERATOR 9B STARTING MATRIX 198 X 981 BITS). 27 UCTAL 04121 56 52 5*06*73 3326*6 12306 5362 10 106*1*30*72016615060*6 30533 43 0**/ 3*75 7250 73 7505 5115**5 7**5221172 3 213670*213710710*1*6*752 777133*30 3236*2216 7*11757*756122*60733*030 2053*36102*72 35233650652*7***716* 02**20350310170*76551* 7*72**65 70* 165*056*31266232221752 7*5 7000*362 2712216*012455*071120**1765*77335 2631*0160262 737*7*3 7171112*575317 115*210 75*700270*026152025362 77*0 300 733020170*51**52 765*0733 522363 11*2226176667 75526152*2 7671060101 25*5125766 7216*302**3612*****53*6 17*62666020127102561613310175*52* 116*40*2 750*2*302 751372*51*2*303* 1732 77276357*2633 7162 751735123 70* 0*12663*766*725 7050 7*050663315507 0*6 76216551*613505353*67*11521175 0556000*0*1*71*6602*65251**1*1172 11523*055557225612666 7*17*50750 72 5*3520716521360353265161663*32*01 1115 7211010304151131*442002112327 1**62710*06111664 0116115*0 5000*73 0*2656222127072*770221651*325*3** 0*3 5*25 7357510326*605*3141*5*7 76* 37705*2513113515*007721025440 3**0 167 7152722311326*602362271*163424 056404160710 7204 731 70 72230207 5 123 2123 174321 75602572 51626S 704032565 36240 7*15016305062532*03315674161 0* 16*60265635 70100*71 4 1*443 105*25 103333*1106*256461313432 531665374 016 700752 72741160542 700016 7340145 3726 54233600232650214 760662026600 55 745141 11002 354 76635 7065 70S 76636 35630316 73623544444741112615156 71 01172 7417500242241325505 71*0666** 2 7454 7006144406450600055 127134121 11505 7430343 151204055001 301672421 3*65 54 110 72*066 Oo*2U3 6*1 7*573 003* 136613600*115**2 7 75305*5012325406 1210655*15502 7637366*2 32613065112 32315055 3050160 766657 7506*21323 72 23 57*35123*16366500*2226010511061 33671577551116621501265*7163055*6 052416015057666 710406202705415 731 176456374546604532 14S472 724400205 11123552205513336*55*2 7*2 1707 7061 03 715752 732363243262 735612 7S42265 1237300026 17545 1000031546602 14700 13704125430160 714272346403422552 7 U251 76 706146 767TS4202300220001672 30355172 76462 70236 116OS5622210010 26742 7436104054 16 113737507 1363426 2371053564440234372433 767 72571060 1600 76365126230 7 77344673173062121 3 412 7 420 71*05006533*30052 5*2*0* 71 337 736007 11361*57026 7* 72 1025556*1 362353551345365022035 7525 7541 564 7 313 765*7652 5*6*5*012640643230**0 7 376727735370072*6*721311666 160515 115156*36 756**35565*465126*310651 2713720765256*2*633*4030335012550 220203122235176062 7 73320243466433 32560 70^4100523 7 100*2622*1361 1113 3307*1*522*3566 3603032 5*561163501 13150066 760316632*7226 7*6 735**025 0511*703*61076333 762*602*56206 753 10**6233655*605502 505155 736171656 303*7156 7232506242 571460336217310 24432312513 7 751013300763604412353 30067 050670615235556015454 706 7045 100564012574 71450536506 7424616632 23352703 721*033706 75216 7*65072*12 2600237*0335322312 755*7*111*1*755 2 5666 7051*51*63*46443532 75015 305 7 33347250412174240201523613562 7503 3065334 766263664550 '3 116760510 746 031675051250300616612602 30115 5 726 121537012054615161650256 756445540 153423043332256700171252201645106 0U36 7340345543 13106540 74032 146024 1561102030450245345404416 52651545 123616 7 7051145000 7603524505216403 2007016462 764206455542412635 7666 7 2 7332231506544152444432522*2**026 04624 7445412146452 4 7064 7066 52 5 406 11464366564144206504153135607 2431 1467 7 55 743050 725643 70 70*4*623*661 24061500027 737332 3 52140363374 3210 U 12036 11 12741 7051 54264 76 74231 6 755 351115607043575435057161142127456 16016 7051340507*3016313*17 7707073 1633 0*5 73 7716121517252261750366 3* 32370565 757*22*15**56 76206 7 73 3366 5131*7035*062 7100*130 755533211065 22 536*62156 76317022300602*3 73 7*63 HEXADECIMAL 14A3 7556003b60698A6357910 23*3189003B1A30998A0062*E 7301D510FA2D26C97C9522 7A6 45EE222F91C88669EAFF96E30 69E891UE13EF9EE294C3B7030 *2B8F10A7*EA6F51AA9E*9CEB 0A441D0C83C4FAC53CEA4D7BS 3AC17*656C9A*8FABCBC011E* 5CA*7*0S4B60E4A121FACF0BA 59980E082EFCFIF5C92A5F59E 26C**7B360B881635055E5FC0 60 76C20789*C957060EOO*9E6 262*B1FB6FE058DS17DC8C082 5652AF0BA3A30A*78A92495CC 3E63B60815C8571C5B20FB2A8 274A22F445185E9704A628C38 30AF07CEFBB37CE5E9EE94F88 10A09CFB*EAF1*782809B368E 137C8EB*CCS015D737S*0**FA 13E00410CE66C14D5552184F4 26A70SB6F4AE2B60E1F28F474 710459051783A06A7109C6A02 240E8904314B25996A0112 9AE 5265C8831276A09O082B02 76 1 1689 2*571D*FC2* 753 1AB1C8 1108AF77C21A030B29A6 59FE6 7F8B152C9 7*08O7EB85650E*C 3BF35 7*99206982 792E61CE28 17410E1C8EB4ECF102610F4A6 *S55E3*70C15EAVCB5E206AEA 79*1E1A0E62 8CAB50566EF0E2 13E9B20 73eC10273*CB59162A 21B6E123*5 7*C59 71AACE05F8 077030S078*E162E005BB80CA 70oB 13/80* J6AU9F0D905B00 77CA612*0*ECFB3BC6BC5FB3C 7730CEEF276492 78495 8D2F 72 04F5E1F4051285AB45E600B48 5E09COC64854A300202B970A2 268BE50E334A105A0160EEA22 755849 1041B0010 7A IF 2F6036 2F62F0109B22FEB0E5053568C 2B806C5685F5EF689AC 580494 699A206285B7DB3FEB01169F4 4EFBE94E1CF6A044960442462 6F7 3 7FB*93823*1SACE7316CC 1 5*381 A2F0B722OC82t2C5 7B2 3F*B9F966C2568CB3AEA*010A 2*A76A*2020B0208BC*78FC62 OF98EAEC3C046B2EEE280896A 29F60058FB2900066C0823380 2F885S8C1C598BA7340E256AE 0A9FB8C66FBFB104C048007T4 610A7AFA6SC2 789C20C 926010 38C5E3C44161C4BEF01CBCE2C *F915002*09C70*6FEF05E*60 380F9EA564C7F0C98B50BC8A2 74AF50E60A06A0t60556282 72 6FF780E*BC65E16F3A215B742 793AE02E57A8485BEA8EC5 74E 63F06 70459A580A006B03120E 7F75FBAF810403A2C90B1C29A 2698A30EE9107AC9A95A32 332 5CBEB7055014C0C8186082AD0 4820CA4903FOCBF60051COA36 6AE IF 3840A9F2 045928 5E2496 6C786 54A5BB3C186ACBB9CE82 2C0036F8335353A58C00D9026 1499C59B«F4B7F^9829 721B06 224C9B06CC200A8A60EF1E75C 61CE6EE9AA328AF3306F25090 52 34CAA3FF482081F3C242906 O0OE280C6 333B6E06CB380C46 20500 157CE651 56 A578A63B34 4005C3E8C00F1B04 7 79A8EA14 3804FC0006952B0B3C24C350A 3 76DC452999CM(98) = UABD. For an IBM 360 (32-bit fullword: 31-bit integer), M(l) = 1HA37556, M(2) = 23U3189D ,. . . , M(98) = 1+ABD322F. For a last example, the CDC 6000 series (word size, 60 bits: integer size, k8 bits; initialized to UT bits) requires M(l) = 00000512156525U06U73, M(2) = 0000106UlU30U720l66l,. . . , M(98) = 00022536U62136T631T. A multiple precision routine (of doubtful value) to 98 bits is easily implemented. 3. COMPILER IMPLEMENTATION Figure 1 presents an example of a non-ANSI standard FORTRAN coding of the GFSR algorithm. FUNCTICN RAND(M,P,Q, INTSIZ) C M(P) = TABLE CF P PREVIOUS RANDOM NUMBERS, C P, Q = POLVNCMIAL PARAMETERS: X**P+ X**Q-f 1 . C .NOT. OPERATOR IMPLEMENTED IN ARITHMETIC. C INTSIZ = INTEGER SIZE (BITS) OF HOST MACHINE: E.G., C IBM 360, 31; CDC 6000, 48; SRU 1100, 35; HP 2100, 15. LOGICAL AA,BB,LCCMPJ f LCOMPK INTEGER A, B,P,Q, INTSIZ, M(l ) EQUIVALENCE ( A A f A) , ( BB, B ) , ( MCOMPJ ,LC0MP J ) , ( MCOMPK , LCOMPK ) CATA J/0/ N*(2**(INTSIZ-l)-l)*2«-l J=J + 1 IF(J.GT.P) J=l K = J+Q IF(K.GT.P) K = K-P MCOMPjsN-M J) MCOMPK =N-M(K) A=M(K) B=M(J) BB = L CO MP J. AND. A A. OR. LCOMPK. AND. BB M( J)=B P.AND = FLOAT(M( J) |/FLOAT(N) RETURN END Figure 1. FORTRAN Implementation of GFSR Pseudorandom Number Generator -5- This particular coding has been certified for l) IBM 360 , 2) Sperry-Rand 1108. 3) Control Data Corporation .6h00 , and k) Hewlett-Packard 2116 computers. The first five GFSR pseudorandom numbers for each of these machines is given in Figure 2. Validity of coding any implementation can be checked against these values. 0.36963295936584470 0.40631365776062010 0.42877840995788570 0.47411382198333740 0.95315784215927120 (a) 0.36963297000000000 0.40631372000000000 0.42877845000000000 0.47411389000000000 0.95315778000000000 (b) 0.36963297409225149 0.40631371808778027 0.42877845193692465 0.47411388879095284 0.95315778681866803 (c) 0.36964017152786255 0.40632343292236328 0.42878508567810059 0.47410506010055542 0.95318460464477539 (d) Figure 2. The First Five Normalized Pseudorandom Numbers From GFSR: x 9 ° + x +1, delayed column=9o00, produced on the ft) IBM 360, b) SRU 1108, c) CDC 6k00 , and d) HP 2116 computers. -6- h. CONCLUSIONS The heretofore unachievable feat of essentially duplicating the same probabilistic simulation on any model computer using the fast, general* multidimensional GFSR pseudorandom number generator is now possible. REFERENCE [l] Lewis, T. G. and Payne, W. H. , "Generalized Feedback Shift Register Pseudorandom Number Algorithm" Journal of the ACM (to appear about August 1973) . -7- APPENDIX A Method ¥ Q Let X + X + 1 "be a primitive trinomial modulo 2 such that Q < P/2. Binary digits, a. , obey the recurrence a = a € a . The J? Q companion matrix for polynomial X + X + 1 is of the form thus , A = 1 ... 1 1 1 All zeros except a 1 in the Q-th column A x V. = l 10 10 a. i a i+l Vi i+1 = V2 = a i+2 a i+Q — a i+P-l a. <8a. i i+Q a i+P = V i+1 For this initialization application, let a_ - a. =. . .= a_ . = 1. Next form V = A x V n . Set V = U , and with V, form a N x 2 dimensioned matri: d 1 d such that V is the first column while V, is the second column: that is U d U 2 = V 1 x[l 0] + V d x[0 1] Form U = U n x[l 0] + A x U Q x 10 1 Continue this process. Matrix U„ will be a P x N matrix -with each of its columns of the form V_ , V,, V_ ...... ,V/. T ,\ ... Numerical values in the 0' d 2d (N-ljd table presented in this paper were to be both in octal and hexadecimal so a value of N , a multiple of 3 and k 3 was selected. Here d = 100 x P = 9800 and N = 2k. The first N bits of table M were the rows of matrix U, T . N Let U be the matrix of the second N bits. Clearly these are easily obtained from dN U N,2 = A XU N For statistical reasons, the generator was cycled another 5000 P times. Order of these operations does not affect the final result so II T was cycled 5000 P times before calculation of U . a 9 d -9- C GFSR PSEUDORANDOM NUMBER ALGORITHM INITIALIZATION TO 98 BITS. INTEGEP M<98,5) INTEGER ONE, DELAY, P,Q P*98 0=27 N = 24 DELAY=100*P 0NE=2**(N-1) DO I 1=1, P 1 M(I,1)=0NE DO 2 J = 1,N DO 3 K=l, DELAY 3 X=IXOP