f7-5

290 Chapter7. RandomNumbersSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).FUNCTION irbit2(iseed) INTEGER irbit2,iseed,IB1,IB2,IB5,IB18,MASK PARAMETER (IB1=1,IB2=2,IB5=16,IB18=131072,MASK=IB1+IB2+IB5) Returns as an integer a random bit, based on the 18 low-significance bits in iseed(which is modified for the next call). if(iand(iseed,IB18).ne.0)then Change all masked bits, shift, and put 1 into bit 1. iseed=ior(ishft(ieor(iseed,MASK),1),IB1)irbit2=1 else Shift and put 0 into bit 1. iseed=iand(ishft(iseed,1),not(IB1)) irbit2=0 endifreturn END A word of caution is: Don’t use sequential bits from these routines as the bits ofalarge,supposedlyrandom,integer,orasthebits inthemantissaofasupposedly random floating-point number. They are not very random for that purpose; see Knuth[1]. Examples of acceptable uses of these random bits are: (i) multiplying a signal randomlyby ±1at a rapid“chiprate,” so as to spreadits spectrumuniformly (but recoverably) across some desired bandpass, or (ii) Monte Carlo exploration of a binary tree, where decisions as to whether to branch left or right are to bemade randomly. Now we do not want you to go through life thinking that there is something special about the primitive polynomial of degree 18 used in the above examples. (We chose 18 because 2 18is small enough for you to verify our claims directly by numerical experiment.) The accompanying table [2]lists one primitive polynomial for each degree up to 100. (In fact there exist many such for each degree. For example, see §7.7 for a complete table up to degree 10.) CITED REFERENCES AND FURTHER READING: Knuth,D.E.1981, SeminumericalAlgorithms ,2nded.,vol.2of TheArtofComputerProgramming (Reading, MA: Addison-Wesley), pp. 29ff. [1] Horowitz,P.,andHill,W.1989, TheArtofElectronics ,2nded.(Cambridge:CambridgeUniversity Press), §§9.32–9.37. Tausworthe, R.C. 1965, Mathematics of Computation , vol. 19, pp. 201–209. Watson, E.J. 1962, Mathematics of Computation , vol. 16, pp. 368–369. [2] 7.5 Random Sequences Based on Data Encryption InNumericalRecipes’ firstedition,wedescribedhowtousetheDataEncryptionStandard (DES)[1-3]for the generation of random numbers. Unfortunately, when implemented in software in a high-level language like FORTRAN, DES is very slow, so excruciatingly slow, in fact, that our previous implementation can be viewed as more mischievous than useful.Here we give a much faster and simpler algorithm which, though it may not be secure in thecryptographic sense, generates about equally good random numbers. DES, like its progenitor cryptographic system LUCIFER, is a so-called “block product cipher” [4]. Itactson64bitsofinputbyiterativelyapplying(16times,infact)akindofhighly 7.5RandomSequencesBasedonDataEncryption 291Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).32-bit XORright 32-bit word left 32-bit word right 32-bit word left 32-bit wordg32-bit XORright 32-bit word left 32-bit word g Figure 7.5.1. The Data Encryption Standard (DES) iterates a nonlinear function gon two 32-bit words, in the manner shown here (after Meyer and Matyas [4]). nonlinear bit-mixing function. Figure 7.5.1 shows the flow of information in DES during this mixing. The function g, which takes 32-bits into 32-bits, is called the “cipher function. ” Meyer and Matyas [4]discuss the importance of the cipher function being nonlinear, as well as other design criteria. DES constructs its cipher function gfrom an intricate set of bit permutations and table lookups acting on short sequences of consecutive bits. Apparently, this function was chosento be particularly strong cryptographically (or conceivably as some critics contend, to havean exquisitely subtle cryptographic flaw!). For our purposes, a different function gthat can be rapidly computed in a high-level computer language is preferable. Such a function mayweaken the algorithm cryptographically. Our purposes are not, however, cryptographic: Wewant tofind thefastest g,and smallestnumber of iterationsofthe mixing procedure inFigure 7.5.1, such that our output random sequence passes the standard tests that are customarilyapplied to random number generators. The resulting algorithm will not be DES, but rather akind of“pseudo-DES, ”better suited to the purpose at hand. Following the criterion, mentioned above, that gshould be nonlinear, we must give the integer multiply operation a prominent place in g. Because 64-bit registers are not generally accessible in high-level languages, we must con fine ourselves to multiplying 16-bit operands into a 32-bit result. So, the general idea of g, almost forced, is to calculate the three distinct 32-bit products of the high and low 16-bit input half-words, and then to combinethese, and perhaps additional fixed constants, by fast operations (e.g., add or exclusive-or) into a single 32-bit result. There are only a limited number of ways of effecting this general scheme, allowing systematic exploration of the alternatives. Experimentation, and tests of the randomness ofthe output, lead to the sequence of operations shown in Figure 7.5.2. The few new elementsin thefigure need explanation: The values C 1and C2arefixed constants, chosen randomly with the constraint that they have exactly 16 1-bits and 16 0-bits; combining these constants 292 Chapter7. RandomNumbersSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).lo2 hi2XOR C1 XOR C2NOT +hi  lo reverse half-words+ Figure 7.5.2. The nonlinear function gused by the routine psdes. via exclusive-or ensures that the overall ghas no bias towards 0 or 1 bits. The“reverse half-words ”operation in Figure 7.5.2 turns out to be essential; otherwise, the very lowest and very highest bits are not properly mixed by the three multiplications.The nonobvious choices in gare therefore: where along the vertical “pipeline”to do the reverse; in what order to combine the three products and C 2; and with which operation (add orexclusive-or)should eachcombining bedone? Wetestedthesechoices exhaustively beforesettling on the algorithm shown in the figure. Itremains todeterminethesmallestnumber ofiterations N itthatwecangetaway with. The minimum meaningful Nitis evidently two, since a single iteration simply moves one 32-bit word without altering it. One can use the constants C1and C2to help determine an appropriate Nit: When Nit=2and C1=C2=0(an intentionally very poor choice), the generator fails several tests of randomness by easily measurable, though not overwhelming,amounts. When N it=4, on the other hand, or with Nit=2but with the constants C1,C 2nonsparse, we have been unable to findanystatistical deviation from randomness in sequences of up to 109floating numbers riderived fromthis scheme. The combined strength ofNit=4and nonsparse C1,C 2should therefore give sequences that are random to tests even far beyond those that we have actually tried. These are our recommended conservativeparameter values, notwithstanding the fact that N it=2(which is, of course, twice as fast) has no nonrandomness discernible (by us). We turn now to implementation. The nonlinear function shown in Figure 7.5.2 is not implementableinstrictlyportable FORTRAN,foratleastthreereasons: (1)Theadditionoftwo 32-bit integers may over flow, and the multiplication of two 16-bit integers may not produce the correct 32-bit product because of sign-bit conventions. We intend that the over flow be ignored, and that the 16-bit integers be multiplied as if they are positive. It is possibleto force this behavior on most machines. (2) We assume 32-bit integers; however, there 7.5RandomSequencesBasedonDataEncryption 293Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).is no reason to believe that longerintegers would be in any way inferior (with suitable extensions of the constants C1,C 2). (3) Your compiler may require a different notation for hex constants (see below). Wehavebeenabletorunthefollowingroutine, psdes,successfullyonmachinesranging from PCs to VAXes and both “big-endian ”and“little-endian ”UNIX workstations. (Big- and little-endian refer to the order in which the bytes are stored in a word.) A strictly portableimplementation is possible in C. If all else fails, you can make a FORTRAN-callable version of the Croutine, found in Numerical Recipes in C . SUBROUTINE psdes(lword,irword) INTEGER irword,lword,NITERPARAMETER (NITER=4) “Pseudo-DES” hashing of the 64-bit word (lword,irword) . Both 32-bit arguments are returned hashed on all bits. NOTE: This routine assumes that arbitrary 32-bit integers can be added without overflow. To accomplish this, you may need to compile with a specialdirective (e.g., /check=nooverflow for VMS). In other languages, such as C, one can instead type the integers as “unsigned.” INTEGER i,ia,ib,iswap,itmph,itmpl,c1(4),c2(4)SAVE c1,c2DATA c1 /Z’BAA96887’,Z’1E17D32C’,Z’03BCDC3C’, Your compiler may use a differ- entnotationforhexconstants! * Z’0F33D1B2’/, c2 /Z’4B0F3B58’,Z’E874F0C3’, * Z’6955C5A6’, Z’55A7CA46’/ do 11i=1,NITER Perform niteriterations of DESlogic, using a simpler (non- cryptographic) nonlinear function instead of DES’s. iswap=irword ia=ieor(irword,c1(i)) The bit-rich constants c1and (below) c2guarantee lots of nonlinear mixing. itmpl=iand(ia,65535) itmph=iand(ishft(ia,-16),65535) ib=itmpl**2+not(itmph**2) ia=ior(ishft(ib,16),iand(ishft(ib,-16),65535))irword=ieor(lword,ieor(c2(i),ia)+itmpl*itmph)lword=iswap enddo 11 return END The routine ran4, listed below, uses psdesto generate uniform random deviates. We adopttheconventionthatanegativevalueoftheargument idumsetstheleft32-bitword,while a positive value isets the right 32-bit word, returns the ith random deviate, and increments idumtoi+1. This is no more than a convenient way of de fining many different sequences (negative values of idum), but still with random access to each sequence (positive values ofidum). For getting a floating-point number from the 32-bit integer, we like to do it by the masking trick described at the end of §7.1, above. The hex constants 3F800000 and 007FFFFF are the appropriate ones for computers using the IEEE representation for 32-bitfloating-point numbers (e.g., IBM PCs and most UNIX workstations). For DEC VAXes, the correct hex constants are, respectively, 00004080 and FFFF007F. Note that your compilermayrequireadifferentnotationforhexconstants,e.g., x’3f800000’ ,’3F800000’X ,ore ven 16#3F800000 . Forgreaterportability,youcaninsteadconstructa floatingnumberbymaking the (signed) 32-bit integer nonnegative (typically, you add exactly 2 31if it is negative) and then multiplying it by a floating constant (typically 2.−31). Aninteresting,andsometimesuseful,featureoftheroutine ran4,below,isthatitallows randomaccesstothe nthrandomvalueinasequence, withoutthenecessityof firstgenerating values 1··· n−1. Thispropertyissharedbyanyrandomnumbergeneratorbasedon hashing (the technique of mapping data keys, which may be highly clustered in value, approximatelyuniformly into a storage address space) [5,6]. One might have a simulation problem in which some certain rare situation becomes recognizable by its consequences only considerably afterithasoccurred. Onemaywishtorestartthesimulationbackatthatoccurrence,usingidenticalrandom values but, say, varying some other control parameters. The relevant question mightthen be something like “what random numbers were used in cycle number 337098901? ”It mightalreadybecycle number395100273 before thequestion comes up. Randomgeneratorsbased on recursion, rather than hashing, cannot easily answer such a question. 294 Chapter7. RandomNumbersSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Values for Verifying the Implementation of psdes idum before psdescall after psdescall (hex) ran4(idum) lword irword lword irword VAX PC –1 1 1 604D1DCE 509C0C23 0.275898 0.219120 99 1 99 D97F8571 A66CB41A 0.208204 0.849246 –99 99 1 7822309D 64300984 0.034307 0.375290 99 99 99 D7F376F0 59BA89EB 0.838676 0.457334 Successive calls to ran4with arguments −1, 99,−99, and 99 should produce exactly the lwordandirwordvalues shown. Maskingconversion toareturned floating randomvalue is allowed to be machine dependent; values for VAX and PC are shown. FUNCTION ran4(idum)INTEGER idumREAL ran4 C USES psdes Returns a uniform random deviate in the range 0.0 to 1.0, generated by pseudo-DES(DES -like)hashingofthe64-bitword (idums,idum) ,where idumswasset byapreviouscallwith negative idum. Alsoincrements idum. Routine canbe used togenerate arandom sequence by successive calls, leaving idumunaltered between calls; or it can randomly access the nth deviate in a sequence by calling with idum =n. Different sequences are initialized by calls with differing negative values of idum. INTEGER idums,irword,itemp,jflmsk,jflone,lword REAL ftemp EQUIVALENCE (itemp,ftemp)SAVE idums,jflone,jflmsk DATA idums /0/, jflone /Z’3F800000’/, jflmsk /Z’007FFFFF’/ Thehexadecimalconstants jfloneandjflmskareusedtoproduceafloatingnumberbetween 1. and 2. by bitwise masking. They are machine-dependent. See text. if(idum.lt.0)then Reset idumsand prepare to return the first devi- ate in its sequence. idums=-idum idum=1 endif irword=idum lword=idumscall psdes(lword,irword) “Pseudo-DES” encode the words. itemp=ior(jflone,iand(jflmsk,irword)) Mask to a floating number between 1 and 2. ran4=ftemp-1.0 Subtraction moves range to 0. to 1. idum=idum+1returnEND The accompanying table gives data for verifying that ran4andpsdeswork correctly on your machine. We do not advise the use of ran4unless you are able to reproduce the hex values shown. Typically, ran4is about 4 times slower than ran0(§7.1), or about 3 times slower than ran1. CITED REFERENCES AND FURTHER READING: Data Encryption Standard , 1977 January 15, Federal Information Processing Standards Publi- cation, number 46(Washington: U.S. Department of Commerce, NationalBureau of Stan-dards). [1] GuidelinesforImplementingandUsingtheNBSDataEncryptionStandard ,1981April1,Federal Information Processing Standards Publication, number 74 (Washington: U.S. Departmentof Commerce, National Bureau of Standards). [2] 7.6SimpleMonteCarloIntegration 295Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Validating the Correctness of Hardware Implementations of the NBS Data Encryption Stan- dard,1980,NBSSpecialPublication500 –20(Washington:U.S.DepartmentofCommerce, National Bureau of Standards). [3] Meyer,C.H.andMatyas,S.M.1982, Cryptography:ANewDimensioninComputerDataSecurity (New York: Wiley). [4] Knuth,D.E. 1973, SortingandSearching ,v ol.3of TheArtofComputerProgramming (Reading, MA: Addison-Wesley), Chapter 6. [5] Vitter,J.S.,andChen,W-C.1987, DesignandAnalysisofCoalescedHashing (NewYork:Oxford University Press). [6] 7.6 Simple Monte Carlo Integration Inspirationsfornumericalmethodscanspringfromunlikelysources. “Splines” first wereflexible strips of wood used by draftsmen. “Simulated annealing ”(we shall see in §10.9)is rooted in a thermodynamicanalogy. And who does not feel at least a faint echo of glamor in the name “Monte Carlo method ”? Supposethatwe pick Nrandompoints,uniformlydistributedina multidimen- sional volume V. Call them x1,...,x N. Then the basic theorem of Monte Carlo integrationestimates the integralofafunction foverthe multidimensionalvolume, /integraldisplay fd V ≈V/angbracketleftf/angbracketright± V/radicalBigg /angbracketleftf2/angbracketright−/angbracketleft f/angbracketright2 N(7.6.1 ) Heretheanglebracketsdenotetakingthearithmeticmeanoverthe Nsamplepoints, /angbracketleftf/angbracketright≡1 NN/summationdisplay i=1f(xi)/angbracketleftbig f2/angbracketrightbig ≡1 NN/summationdisplay i=1f2(xi)( 7.6.2 ) The“plus-or-minus ”term in (7.6.1) is a one standard deviation error estimate for the integral, not a rigorous bound; further, there is no guarantee that the error is distributedasaGaussian,sotheerrortermshouldbetakenonlyasaroughindication of probable error. Suppose that you want to integrate a function gover a region Wthat is not easy to sample randomly. For example, Wmight have a very complicated shape. No problem. Just find a region Vthatincludes Wand thatcaneasily be sampled (Figure 7.6.1), and then de finefto be equal to gfor points in Wand equal to zero for points outside of W(but still inside the sampled V). You want to try to make Venclose Was closely as possible, because the zero values of fwill increase the error estimate term of (7.6.1). And well they should: points chosen outside of W have no information content, so the effective value of N, the number of points, is reduced. The error estimate in (7.6.1) takes this into account. General purpose routines for Monte Carlo integration are quite complicated (see§7.8),butaworkedexamplewill showtheunderlyingsimplicityofthemethod. Supposethatwe want to findthe weightandthe positionof thecenterof mass ofan