Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 72
Текст из файла (страница 72)
They are not very random for that purpose; seeKnuth [1]. Examples of acceptable uses of these random bits are: (i) multiplying asignal randomly by ±1 at a rapid “chip rate,” so as to spread its spectrum uniformly(but recoverably) across some desired bandpass, or (ii) Monte Carlo explorationof 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 somethingspecial about the primitive polynomial of degree 18 used in the above examples.(We chose 18 because 218 is small enough for you to verify our claims directly bynumerical experiment.) The accompanying table [2] lists one primitive polynomialfor each degree up to 100.
(In fact there exist many such for each degree. Forexample, see §7.7 for a complete table up to degree 10.)CITED REFERENCES AND FURTHER READING:Knuth, D.E. 1981, Seminumerical Algorithms, 2nd ed., vol. 2 of The Art of Computer Programming(Reading, MA: Addison-Wesley), pp. 29ff. [1]Horowitz, P., and Hill, W.
1989, The Art of Electronics, 2nd ed. (Cambridge: Cambridge UniversityPress), §§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 DataEncryptionIn Numerical Recipes’ first edition,we described how to use the Data Encryption Standard(DES) [1-3] for the generation of random numbers. Unfortunately, when implemented insoftware in a high-level language like C, DES is very slow, so excruciatingly slow, in fact, thatour previous implementation can be viewed as more mischievous than useful.
Here we givea much faster and simpler algorithm which, though it may not be secure in the cryptographicsense, generates about equally good random numbers.DES, like its progenitor cryptographic system LUCIFER, is a so-called “block productcipher” [4]. It acts on 64 bits of input by iteratively applying (16 times, in fact) a kind of highlynonlinear bit-mixing function. Figure 7.5.1 shows the flow of information in DES duringthis 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 wellas other design criteria.DES constructs its cipher function g from an intricate set of bit permutations and tablelookups 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 g that canbe rapidly computed in a high-level computer language is preferable. Such a function mayweaken the algorithm cryptographically.
Our purposes are not, however, cryptographic: Wewant to find the fastest g, and smallest number of iterations of the mixing procedure in Figure7.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 g should be nonlinear, we must givethe integer multiply operation a prominent place in g.
Because 64-bit registers are notgenerally accessible in high-level languages, we must confine ourselves to multiplying 16-bit7.5 Random Sequences Based on Data Encryptionleft 32-bit word301right 32-bit wordg32-bit XORleft 32-bit wordright 32-bit wordg32-bit XORleft 32-bit wordright 32-bit wordFigure 7.5.1. The Data Encryption Standard (DES) iterates a nonlinear function g on two 32-bit words,in the manner shown here (after Meyer and Matyas [4]).operands into a 32-bit result. So, the general idea of g, almost forced, is to calculate the threedistinct 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, allowingsystematic 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 the figure need explanation: The values C1 and C2 are fixed constants, chosen randomlywith the constraint that they have exactly 16 1-bits and 16 0-bits; combining these constantsvia exclusive-or ensures that the overall g has 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 g are therefore: where along the vertical “pipeline” to do thereverse; in what order to combine the three products and C2 ; and with which operation (addor exclusive-or) should each combining be done? We tested these choices exhaustively beforesettling on the algorithm shown in the figure.It remains to determine the smallest number of iterations Nit that we can get away with.The minimum meaningful Nit is evidently two, since a single iteration simply moves one32-bit word without altering it.
One can use the constants C1 and C2 to help determine anappropriate Nit : When Nit = 2 and C1 = C2 = 0 (an intentionally very poor choice), thegenerator fails several tests of randomness by easily measurable, though not overwhelming,amounts. When Nit = 4, on the other hand, or with Nit = 2 but with the constantsC1 , C2 nonsparse, we have been unable to find any statistical deviation from randomness insequences of up to 109 floating numbers ri derived from this scheme. The combined strengthof Nit = 4 and nonsparse C1 , C2 should therefore give sequences that are random to testseven far beyond those that we have actually tried.
These are our recommended conservative302Chapter 7.C1Random NumbersXORhi 2lo 2hi • loNOT+reversehalf-wordsC2XOR+Figure 7.5.2.The nonlinear function g used by the routine psdes.parameter values, notwithstanding the fact that Nit = 2 (which is, of course, twice as fast)has no nonrandomness discernible (by us).Implementation of these ideas is straightforward. The following routine is not quitestrictly portable, since it assumes that unsigned long integers are 32-bits, as is the caseon most machines. However, there is no reason to believe that longer integers would be inany way inferior (with suitable extensions of the constants C1 , C2 ). C does not provide aconvenient, portable way to divide a long integer into half words, so we must use a combinationof masking (& 0xffff) with left- and right-shifts by 16 bits (<<16 and >>16).
On somemachines the half-word extraction could be made faster by the use of C’s union construction,but this would generally not be portable between “big-endian” and “little-endian” machines.(Big- and little-endian refer to the order in which the bytes are stored in a word.)#define NITER 4void psdes(unsigned long *lword, unsigned long *irword)“Pseudo-DES” hashing of the 64-bit word (lword,irword). Both 32-bit arguments are returned hashed on all bits.{unsigned long i,ia,ib,iswap,itmph=0,itmpl=0;static unsigned long c1[NITER]={0xbaa96887L, 0x1e17d32cL, 0x03bcdc3cL, 0x0f33d1b2L};static unsigned long c2[NITER]={0x4b0f3b58L, 0xe874f0c3L, 0x6955c5a6L, 0x55a7ca46L};for (i=0;i<NITER;i++) {Perform niter iterations of DES logic, using a simpler (non-cryptographic) nonlinear function instead of DES’s.7.5 Random Sequences Based on Data Encryption303ia=(iswap=(*irword)) ^ c1[i];The bit-rich constants c1 and (below)itmpl = ia & 0xffff;c2 guarantee lots of nonlinear mixitmph = ia >> 16;ing.ib=itmpl*itmpl+ ~(itmph*itmph);*irword=(*lword) ^ (((ia = (ib >> 16) |((ib & 0xffff) << 16)) ^ c2[i])+itmpl*itmph);*lword=iswap;}}The routine ran4, listed below, uses psdes to generate uniform random deviates.
Weadopt the convention that a negative value of the argument idum sets the left 32-bit word, whilea positive value i sets the right 32-bit word, returns the ith random deviate, and incrementsidum to i + 1. This is no more than a convenient way of defining 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 themasking trick described at the end of §7.1, above. The hex constants 3F800000 and 007FFFFFare the appropriate ones for computers using the IEEE representation for 32-bit floating-pointnumbers (e.g., IBM PCs and most UNIX workstations). For DEC VAXes, the correct hexconstants are, respectively, 00004080 and FFFF007F.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
