Numerical Recipes
`in Fortran 77
`Second Edition
`
`Volume 1 of
`Fortran Numerical Recipes
`
`PGS Exhibit 2022
`WesternGeco v. PGS (IPR2015-00309, 310, 311)
`
`

`
`Numerical Recipes
`
`in Fortran 77
`The Art of Scientific Computing
`Second Edition
`
`Volume 1 of
`Fortran Numerical Recipes
`
`William H. Press
`Harvard-Smithsonian Center for Astrophysics
`
`Saul A. Teukolsky
`Department of Physics, Cornell University
`
`William T. Vetterling
`Polaroid Corporation
`
`Brian P. Flannery
`EXXON Research and Engineering Company
`
`~ C.Alv.IBRIDGE
`UNIVERSITY PRESS
`
`PGS Exhibit 2022
`WesternGeco v. PGS (IPR2015-00309, 310, 311)
`
`

`
`Published by the Press Syndicate of the University of Cambridge
`The Pitt Building, Trumpington Street, Cambridge CB2 1RP
`40 West 20th Street, New York, NY 10011-4211, USA
`10 Stamford Road, Oakleigh, Melbourne 3166, Australia
`
`Copyright © Cambridge University Press 1986, 1992
`except for §13.10, which is placed into the public domain,
`and except for all other computer programs and procedures, which are
`Copyright © Numerical Recipes Software 1986, 1992
`All Rights Reserved.
`
`Some sections of this book were originally published, in different form, in Computers
`in Physics magazine, Copyright© American Institute of Physics, 1988-1992.
`
`First Edition originally published 1986; Second Edition originally published 1992 as
`Numerical Recipes in FORTRAN: The Art of Scientific Computing
`Reprinted with corrections, 1993, 1994, 1995.
`Reprinted with corrections, 1996 as Numerical Recipes in Fortran 77: The Art of
`Scientific Computing (Vol. 1 of Fortran Numerical Recipes)
`
`This reprinting is corrected to software version 2.06
`
`Printed in the United States of America
`Typeset in 'IE;X
`
`Without an additional license to use the contained software, this book is intended as
`a text and reference book, for reading purposes only. A free license for limited use of the
`software by the individual owner of a copy of this book who personally types one or more
`routines into a single computer is granted under terms described on p. xxi. See the section
`"License Information" (pp. xx-xxiii) for information on obtaining more general licenses at
`low cost.
`Machine-readable media containing the software in this book, with included licenses
`for use on a single screen, are available from Cambridge University Press. See the
`order form at the back of the book, email to "orders@cup.org" (North America) or
`"trade@cup.cam.ac.uk" (rest of world), or write to Cambtidge University Press, 110
`Midland Avenue, Port Chester, NY 10573 (USA), for further information.
`The software may also be downloaded, with immediate purchase of a license
`also possible, from the Numerical Recipes Software Web Site (http: I /wwv. nr. com).
`Unlicensed transfer of Numerical Recipes programs to any other format, or to any computer
`except one that is specifically licensed, is strictly prohibited. Technical questions,
`corrections, and requests for iriformation should be addressed to Numerical Recipes
`Software, P.O. Box 243, Cambridge, MA 02238 (USA), email "info@nr.com", or fax
`617 863-1739.
`
`Library of Congress Cataloging in Publication Data
`Numerical recipes in Fortran 77 : the art of scientific computing I William H. Press
`. . . [et al.].
`- 2nd ed.
`Includes bibliographical references (p.
`ISBN 0-521-43064-X
`1. Numerical analysis-Computer programs. 2. Science-Mathematics-Computer programs.
`3. FORTRAN (Computer program language)
`I. Press, William H.
`QA297 .N866 1992
`519.4'0285'53-dc20
`
`) and index.
`
`92-8876
`
`A catalog record for this book is available from the British Library.
`
`ISBN
`ISBN
`ISBN
`ISBN
`ISBN
`ISBN
`
`0 52143064 X
`0 521 57439 0
`0 521 43721 0
`0 521 57440 4
`0 521 57608 3
`0 521 57607 5
`
`Volume 1 (this book)
`Volume 2
`Example book in FORTRAN
`FORTRAN diskette (IBM 3.5")
`CDROM (IBM PC/Macintosh)
`CDROM (UNIX)
`
`PGS Exhibit 2022
`WesternGeco v. PGS (IPR2015-00309, 310, 311)
`
`

`
`Chapter 7.
`
`Random Numbers
`
`7.0 Introduction
`
`It may seem perverse to use a computer, that most precise and deterministic of
`all machines conceived by the human mind, to produce "random" numbers. More
`than perverse, it may seem to be a conceptual impossibility. Any program, after all,
`will produce output that is entirely predictable, hence not truly "random."
`Nevertheless, practical computer "random number generators" are in common
`use. We will leave it to philosophers of the computer age to resolve the paradox in
`a deep way (see, e.g., Knuth [1] §3.5 for discussion and references). One sometimes
`hears computer-generated sequences termed pseudo-random, while the word random
`is reserved for the output of an intrinsically random physical process, like the elapsed
`time between clicks of a Geiger counter placed next to a sample of some radioactive
`element. We will not try to make such fine distinctions.
`A working, though_ imprecise, definition of randomness in the context of
`computer-generated sequences, is to say that the deterministic program that produces
`a random sequence should be different from, and -
`in all measurable respects -
`statistically uncorrelated with, the computer program that uses its output. In other
`words, any two different random number generators ought to produce statistically
`the same results when coupled to your particular applications program. If they don't,
`then at least one of them is not (from your point of view) a good generator.
`The above definition may seem circular, comparing, as it does, one generator to
`another. However, there exists a body of random number generators which mutually
`do satisfy the definition over a very, very broad class of applications programs.
`And it is also found empirically that statistically identical results are obtained from
`random numbers produced by physical processes. So, because such generators are
`known to exist, we can leave to the philosophers the problem of defining them.
`A pragmatic point of view, then, is that randomness is in the eye of the beholder
`(or programmer). What is random enough for one application may not be random
`enough for another. Still, one is not entirely adrift in a sea of incommensurable
`applications programs: There is a certain list of statistical tests, some sensible and
`some merely enshrined by history, which on the whole will do a very good job
`of ferreting out any correlations that are likely to be detected by an applications
`program (in this case, yours). Good random number generators ought to pass all of
`these tests; or at least the user had better be aware of any that they fail, so that he or
`she will be able to judge whether they are relevant to the case at hand.
`-
`
`266
`
`PGS Exhibit 2022
`WesternGeco v. PGS (IPR2015-00309, 310, 311)
`
`

`
`7. 1 Uniform Deviates
`
`267
`
`As for references on this subject, the one to tum to first is Knuth [1 ]. Then
`try [2]. Only a few of the standard books on numerical methods [3-4] treat topics
`relating to random numbers.
`
`CITED REFERENCES AND FURTHER READING:
`Knuth, D.E.1981, SBmlnumerlca/A/gorlthms, 2nded., vol. 2of TheArtofComputerProgrammlng
`(Reading, MA: Addison-Wesley), Chapter 3, especially §3.5. [1)
`Bratley, P., Fox, B.L., and Schrage, E.L. 1983, A Guide to Simulation (New York: Springer(cid:173)
`Verlag). [2]
`Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall),
`Chapter 11 . [3]
`Forsythe, G.E., Malcolm, M.A., and Moler, C.B. 19n, Computer Methods for Mathematical
`Computations (Englewood Cliffs, NJ: Prentice-Hall), Chapter 1 0. [4]
`
`7. 1 Uniform Deviates
`
`Uniform deviates are just random numbers that lie within a specified range
`(typically 0 to 1 ), with any one number in the range just as likely as any other. They
`are, in other words, what you probably think "random numbers" are. However,
`we want to distinguish uniform deviates from other sorts of random numbers, for
`example numbers drawn from a normal (Gaussian) distribution of specified mean
`and standard deviation. These other sorts of deviates are almost always generated by
`performing appropriate operations on one or more uniform deviates, as we will see
`in subsequent sections. So, a reliable source of random uniform deviates, the subject
`of this section, is an essential building block for any sort of stochastic modeling
`or Monte Carlo computer work.
`
`System-Supplied Random Number Generators
`
`Your computer very likely has lurking within it a library routine which is called
`a "random number generator." That routine typically has an unforgettable name like
`"ran," and a calling sequence like
`
`z-ran(ieeed)
`
`sets z to the next random number and updates iaeed
`
`You initialize iseed to a (usually) arbitrary value before the first call to ran.
`Each initializing value will typically return a different subsequent random sequence,
`or at least a different subsequence of some one enormously long sequence. The same
`initializing value of iseed will always return the same random sequence, however.
`Now our first, and perhaps most important, lesson in this chapter is: Be very,
`very suspicious of a system-supplied ran that resembles the one just described. If all
`scientific papers whose results are in doubt because of bad rans were to disappear
`from library shelves, there would be a gap on each shelf about as big as your
`fist. System-supplied rans are almost always linear congruential generators, which
`
`PGS Exhibit 2022
`WesternGeco v. PGS (IPR2015-00309, 310, 311)
`
`

`
`268
`
`Chapter 7.
`
`Random Numbers
`
`generate a sequence of integers It, /2, I a, ... , each between 0 and m - 1 (a large
`number) by the recurrence relation
`
`Ij+l = alj + c
`
`(mod m)
`
`(7.1.1)
`
`Herem is called the modulus, and a and care positive integers called the multiplier
`and the increment, respectively. The recurrence (7 .1.1) will eventually repeat itself,
`with a period that is obviously no greater than m. If m, a, and c are properly chosen,
`then the period will be of maximal length, i.e., of length m. In that case, all possible
`integers between 0 and m - 1 occur at some point, so any initial "seed" choice of I o
`is as good as any other: The sequence just takes off from that point. The real number
`between 0 and 1 which is returned is generally Ij+1 /m, so that it is strictly less than
`1, but occasionally (once in m calls) exactly equal to zero. iseed is set to Ij+l (or
`some encoding of it), so that it can be used on the next call to generate Ij+2, and so on.
`The linear congruential method has the advantage of being very fast, requiring
`only a few operations per call, hence its almost universal use. It has the disadvantage
`that it is not free of sequential correlation on successive calls. If k random numbers at
`a time are used to plot points in k dimensional space (with each coordinate between
`0 and 1 ), then the points will not tend to "fill up" the k-dimensional space, but
`rather will lie on (k- !)-dimensional "planes." There will be at most about m 1/k
`such planes. If the constants m, a, and c are not very carefully chosen, there will
`be many fewer than that. The number m is usually close to the machine's largest
`representable integer, e.g., ""'"' 2 32 • So, for example, the number of planes on which
`triples of points lie in three-dimensional space is usually no greater than about the
`cube root of 232 , about 1600. You might well be focusing attention on a physical
`process that occurs in a small fraction of the total volume, so that the discreteness
`of the planes can be very pronounced.
`Even worse, you might be using a ran whose choices of m, a, and c have
`been botched. One infamous such routine, RANDU, with a= 65539 and m = 2 31 ,
`was widespread on IBM mainframe computers for many years, and widely copied
`onto other systems [1 ]. One of us recalls producing a "random" plot with only 11
`planes, and being told by his computer center's programming consultant that he
`had misused the random number generator: "We guarantee that each number is
`random individually, but we don't guarantee that more than one of them is random."
`Figure that out.
`
`Correlation in k-space is not the only weakness oflinear congruential generators.
`Such generators often have their low-order (least significant) bits much less random
`than their high-order bits. If you want to generate a random integer between 1 and
`10, you should always do it using high-order bits, as in
`
`j•1+int(10.•ran(iseed))
`
`and never by anything resembling
`
`j~1+mod(int(1000000.•ran(iseed)),10)
`
`PGS Exhibit 2022
`WesternGeco v. PGS (IPR2015-00309, 310, 311)
`
`

`
`7. 1 Uniform Deviates
`
`269
`
`(which uses lower-order bits). Similarly you should never try to take apart a
`"ran" number into several supposedly random pieces. Instead use separate calls
`for every piece.
`
`Portable Random Number Generators
`
`Park and Miller [1 1 have surveyed a large number of random number generators
`that have been used over the last 30 years or more. Along with a good theoretical
`review, they present an anecdotal sampling of a number of inadequate generators
`that have come into widespread use. The historical record is nothing if not appalling.
`There is good evidence, both theoretical and empirical, that the simple mul(cid:173)
`tiplicative congruential algorithm
`
`Ij+l =ali
`
`(mod m)
`
`(7.1.2)
`
`can be as good as any of the more general linear congruential generators that have
`c =I= 0 (equation 7.1.1)- ifthe multiplier a and modulus mare chosen exquisitely
`carefully. Park and Miller propose a "Minimal Standard" generator based on the
`choices
`
`a= 75 = 16807
`
`m = 231
`
`- 1 = 2147483647
`
`(7 .1.3)
`
`First proposed by Lewis, Goodman, and Miller in 1969, this generator has in
`subsequent years passed all new theoretical tests, and (perhaps more importantly)
`has accumulated a large amount of successful use. Park and Miller do not claim that
`the generator is "perfect" (we will see below that it is not), but only that it is a good
`minimal standard against which other generators should be judged.
`It is not possible to implement equations (7 .1.2) and (7 .1.3) directly in a
`high-level language, since the product of a and m - 1 exceeds the maximum value
`for a 32-bit integer. Assembly language implementation using a 64-bit product
`register is straightforward, but not portable from machine to machine. A trick
`due to Schrage [2,3) for multiplying two 32-bit integers modulo a 32-bit constant,
`without using any intermediates larger than 32 bits (including a sign bit) is therefore
`extremely interesting: It allows the Minimal Standard generator to be implemented
`in essentially any programming language on essentially any machine.
`Schrage's algorithm is based on an approximate factorization of m,
`m = aq + r,
`
`i.e., q = [m/a], r = m mod a
`
`(7.1.4)
`
`with square brackets denoting integer part. If r is small, specifically r < q, and
`0 < z < m- 1, it can be shown that both a(z mod q) and r[z/q] lie in the range
`0, ... , m - 1, and that
`
`a(z mod q) - r(z/q]
`azmo m= a(z mod q) - r[z/q] + m
`d
`{
`
`if it is> 0,
`otherwise
`
`(7.1.5)
`
`The application of Schrage's algorithm to the constants (7.1.3) uses the values
`q = 127773 and r = 2836.
`Here is an implementation of the Minimal Standard generator:
`
`PGS Exhibit 2022
`WesternGeco v. PGS (IPR2015-00309, 310, 311)
`
`

`
`270
`
`Chapter 7. Random Numbers
`
`•
`
`FUNCTION ranO(idwa)
`INTEGER idum.IA.IH.IQ.IR.MASK
`REAL ranO • AM
`PARAMETER (IA•16807.IH•2147483647.AH•1./IH.
`IQ•127773.IR•2836.HASK•123469876)
`"Minimal" random number generator of Park and Miller. Returns a uniform random deviate
`between 0.0 and 1.0. Set or reset idum to any integer value (except the unlikely value MASK)
`to initialize the sequence; idum must not be altered between calls for successive deviates
`in a sequence.
`INTEGER k
`idum•ieor(idum.HASK)
`k•idum/IQ
`idum•IA•(idum-k•IQ)-IR•k
`if (idum.lt.O) idum•idum+IH
`ranO•AH•idum
`idum•ieor(idum.HASK)
`return
`END
`
`XORing with MASK allows use of zero and other simple
`bit patterns for idum.
`Compute idum-mod(IA•idum,IH) without overflows by
`Schrage's method.
`Convert idum to a floating result.
`Unmask before return.
`
`The period of ranO is 231 - 2 ~ 2.1 x 109 • A peculiarity of generators of
`the fonn (7 .1.2) is that the value 0 must never be allowed as the initial seed -
`it
`perpetuates itself- and it never occurs for any nonzero initial seed. Experience has
`shown that users always manage to call random number generators with the seed
`idum•O. That is why ranO performs its exclusive-or with an arbitrary constant both
`on entry and exit. If you are the first user in history to be proof against human error,
`you can remove the two lines with the ieor function.
`Park and Miller discuss two other multipliers a that can be used with the same
`m = 231 - 1. These are a = 48271 (with q = 44488 and r = 3399) and a = 69621
`(with q = 30845 and r -_23902). These can be substituted in the routine ranO
`if desired; they may be slightly superior to Lewis et al. 's longer-tested values. No
`values other than these should be used.
`The routine ranO is a Minimal Standard, satisfactory for the majority of appli(cid:173)
`cations, but we do not recommend it as the final word on random number generators.
`Our reason is precisely the simplicity of the Minimal Standard. It is not hard to think
`of situations where successive random numbers might be used in a way that acciden(cid:173)
`tally conflicts with the generation algorithm .. · For example, since successive numbers
`differ by a multiple of only 1.6 x 104 out of a modulus of more than 2 x 109 , very small
`random numbers will tend to be followed by smaller than average values. One time
`in 106 , for example, there will be a value < 10-6 returned (as there should be), but
`this will always be followed by a value less than about 0.0168. One can easily think
`of applications involving rare events where this property would lead to wrong results.
`There are other, more subtle, serial correlations present in ranO. For example,
`if successive points (Ji, /i+1) are binned into a two-dimensional plane for i =
`1, 2, ... , N, then the resulting distribution fails the x 2 test when N is greater than a
`few x 107 , much less than the period m - 2. Since low-order serial correlations have
`historically been such a bugaboo, and since there is a very simple way to remove
`them, we think that it is prudent to do so.
`The following routine, ran1, uses the Minimal Standard for its random value,
`but it shuffles the output to remove low-order serial correlations. A random deviate
`derived from the jth value in the sequence, I;, is output not on the jth call, but rather
`ori. a randomized later call, j + 32 on average. The shuffling algorithm is due to Bays
`and Durham as described in Knuth [4), and is illustrated in Figure 7 .1.1.
`
`PGS Exhibit 2022
`WesternGeco v. PGS (IPR2015-00309, 310, 311)
`
`

`
`7. 1 Uniform Deviates
`
`271
`
`•
`
`FUNCTION ranl (idum)
`INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV
`IlEAL ranl, AM, EPS, RNMX
`PARAMETER (IA•16807,IM•2147483647,AM•1./IM,IQ•127773,IR•2836,
`NTAB•32,NDIV•1+(IM-1)/NTAB,EPS•1.2e-7,RNMX•1.-EPS)
`"Minimal" random number generator of Park and Miller with Bays-Durham shuffle and
`added safeguards. Returns a uniform random deviate between 0.0 and 1.0 (exclusive of
`the endpoint values). Call with idum. a negative integer to initialize; thereafter, do not
`alter idum. between successive deviates in a sequence. RNMX should approximate the largest
`floating value that is less than 1.
`INTEGER j,k,iv(NTAB),iy
`SAVE iv,iy
`DATA iv /NTAB•O/, iy /0/
`if (idum .le. 0. or. iy. eq. 0) then Initialize.
`Be sure to prevent idum = 0.
`idum-maJt(-idum,l)
`do 11 j•NTAB+B,l,-1
`load the shuffle table {after 8 warm-ups).
`k•idum/IQ
`idum•IA•(idum-k•IQ)-IR•k
`if (idum.lt.O) idum•idum+IM
`if (j.le.NTAB) iv(j)•idum
`enddo 11
`iy•iv(l)
`end if
`k•idum/IQ
`idum•IA• (idum-k•IQ) -IR•k
`if (idum.lt.O) idum•idum+IM
`j•1+iy/NDIV
`iy-iv(j)
`iv(j)•idum
`ranl-min(AM•iy,RNMX)
`return
`END
`
`Start here when not initializing.
`Compute idum-mod(IA•idum, IM) without overflows. by
`Schrage's method.
`Will be in the range 1:NTAB.
`Output previously stored value and refill the shuffle ta-
`ble.
`Because users don't expect endpoint values.
`
`The routine rant passes those statistical tests that ranO is known to fail. In
`fact, we do not know of any statistical test that rant fails to pass, except when the
`number of calls starts to become on the order of the period m, say > 108 ~ m/20.
`For situations when even longer random sequences are needed, L'Ecuyer [6] has
`given a good way of combining two different sequences with different periods so
`as to obtain a new sequence whose period is the least common multiple of the two
`periods. The basic idea is simply to add the two sequences, modulo the modulus of
`either of them (call it m). A trick to avoid an intermediate value that overflows the
`integer wordsize is to subtract rather than add, and then add back the constant m - 1
`if the result is < 0, so as to wrap around into the desired interval 0, ... , m - 1.
`Notice that it is not necessary that this wrapped subtraction be able to reach all
`values 0, ... , m - 1 from every value of the first sequence. Consider the absurd
`extreme case where the value subtracted was only between 1 and 10: The resulting
`sequence would still be no less random than the first sequence by itself. As a
`practical matter it is only necessary that the second sequence have a range coveiing
`substantially all of the range of the first. L'Ecuyer recommends the use of the two
`generators m1 = 2147483563 (with a1 = 40014, Q1 = 53668, r1 = 12211) and
`m2 = 2147483399 (with a2 = 40692, Q2 = 52774, r2 = 3791). Both moduli
`are slightly less than 231 . The periods m1 - 1 = 2 x 3 x 7 x 631 x 81031 and
`m2 - 1 = 2 x 19 x 31 x 1019 x 1789 share only the factor 2, so the period of
`the combined generator is~ 2.3 x 1018 • For present computers, period exhaustion
`is a practical impossibility.
`
`PGS Exhibit 2022
`WesternGeco v. PGS (IPR2015-00309, 310, 311)
`
`

`
`272
`
`Chapter 7.
`
`Random Numbers
`
`OUTPUT
`
`CD
`
`RAN
`
`. . .
`
`Shuffling procedure used in ran1 to break up sequential correlations in the Minimal
`Figure 7 .1.1.
`Standard generator. Circled numbers indicate the sequence of events: On each call, the random number
`in iy is used to choose a random element in the array iv. That element becomes the output random
`number, and also is the next iy. Its spot in iv is refilled from the Minimal Standard routine.
`
`Combining the two generators breaks up serial correlations to a considerable
`extent. We nevertheless recommend the additional shuffle that is implemented in
`the following routine, ran2. We think that, within the limits of its floating-point
`precision, ran2 provides perfect random numbers; a practical definition of "perfect"
`is th~t we will pay $1 000 to the first reader who convinces us otherwise (by finding a
`statistical test that ran2 fails in a nontrivial way, excluding the ordinary limitations
`of a machine's floating-point representation).
`
`*
`*
`
`FUNCTION ran2(idum)
`INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
`REAL ran2,AM,EPS,RNMX
`PARAMETER (IM1=2147483563,IM2=2147483399,AM•1./IM1,IMM1~IM1-1,
`IA1•40014,IA2=40692,IQ1=53668,IQ2=62774,IR1•12211,
`IR2~3791,NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX•1.-EPS)
`Long period (> 2 x 1018 ) random number generator of L'Ecuyer with Bays-Durham shuffle
`and added safeguards. Returns a uniform random deviate between 0.0 and 1.0 (exclusive
`of the endpoint values). Call with idum a negative integer to initialize; thereafter, do not
`alter idum between successive deviates in a sequence. RNMX should approximate the largest
`floating value that is less than 1.
`INTEGER idum2,j,k,iv(NTAB),iy
`SAVE iv,iy,idum2
`DATA idum2/123466789/, iv/NTAB•O/, iy/0/
`Initialize.
`if (idum.1e.O) then
`Be sure to prevent idum = 0.
`idum=max(-idum, 1)
`idum2=idum
`do 11
`j=NTAB+S, 1,-1
`k=idum/IQ1
`
`Load the shuffle table (after 8 warm-ups).
`
`PGS Exhibit 2022
`WesternGeco v. PGS (IPR2015-00309, 310, 311)
`
`

`
`7. 1 Uniform Deviates
`
`273
`
`idum•IA1•(idum-k•IQ1)-k•IR1
`if (idum.lt.O) idum•idum+IH1
`if (j.le.NTAB) iv(j)•idum
`enddo 11
`iy•iv(1)
`end if
`k•idum/IQ1
`idum•IA1•(1dum-k•IQ1)-k•IR1
`if (idum.lt.O) idum•idum+IH1
`k•idum2/IQ2
`idum2•IA2•(idum2-k•IQ2)-k•IR2
`if (idum2.lt.O) idum2•idum2+IH2
`j•1+1y/NDIV
`iy•iv(j)-idum2
`iv(j)•idum
`1f(1y.lt.1)iy•1y+IHH1
`ran2-m1n(AH•iy,RNMX)
`return
`END
`
`Start here when not initializing.
`Compute 1dum-mod(IA1•idum, IH1) without over(cid:173)
`flows by Schrage's method .
`
`Compute idum2-mod(IA2•idum2, IH2) likewise.
`
`Will be in the range 1: NTAB.
`Here idum is shuffled, idum and idum2 are com(cid:173)
`bined to generate output.
`
`Because users don't expect endpoint values.
`
`L'Ecuyer [6] lists additional short generators that can be combined into longer
`ones, including generators that can be implemented in 16-bit integer arithmetic.
`Finally, we give you Knuth's suggestion [4] for a portable routine, which we
`have translated to the present conventions as ran3. This is not based on the linear
`congruential method at all, but rather on a subtractive method (see also [51). One
`might hope that its weaknesses, if any, are therefore of a highly different character
`from the weaknesses, if any, of ranl above. If you ever suspect trouble with one
`routine, it is a good idea to try the other in the same application. ran3 has one
`nice feature: if your machine is poor on integer arithmetic (i.e., is limited to 16-bit
`integers), substitution of the three "commented" lines for the ones directly preceding
`them will render the routine entirely floating-point.
`
`C
`
`C
`
`C
`
`Now initialize the rest of the table,
`in a slightly random order,
`with numbers that are not especially random.
`
`FUNCTION ran3(idum)
`Returns a uniform random deviate between 0.0 and 1.0. Set idum to any negative value
`to initialize or reinitialize the sequence.
`INTEGER idum
`INTEGER HBIG,HSEED,HZ
`REAL MBIG,HSEED,HZ
`REAL ran3,FAC
`PARAMETER (HBIG•1000000000,HSEED•161803398,HZ•O,FAC•1./HBIG)
`PARAMETER (HBIG•4000000.,HSEED•1618033.,HZ•O.,FAC•1./HBIG)
`According to Knuth, any large mbig, and any smaller (but still large) mseed can be sub(cid:173)
`stituted for the above values.
`INTEGER i,iff,ii,inext,inextp,k
`The value 55 is special and should not be modified; see
`INTEGER mj ,mk,ma(66)
`Knuth.
`REAL mj ,mk,ma(66)
`SAVE iff,inext,inextp,ma
`DATA iff /0/
`if(idum.lt .O.or. iff .eq.O)then
`iff•1
`mj•HSEED-iaba (idum)
`mj-mod(mj ,HBIG)
`ma(66)-mj
`mk•1
`do 11 1•1,64
`ii-mod(21*i ,66)
`ma(ii)-mk
`mk-mj-mk
`if(mk.lt.MZ)mk-mk+HBIG
`
`Initialization.
`
`Initialize ma(66) using the seed idum and the large num-
`ber maeed.
`
`PGS Exhibit 2022
`WesternGeco v. PGS (IPR2015-00309, 310, 311)
`
`

`
`274
`
`Chapter 7.
`
`Random Numbers
`
`We randomize them by "warming up the generator."
`
`Prepare indices for our first generated number.
`The constant 31 is special; see Knuth.
`
`mj-ma(ii)
`enddo 11
`do 13 k•1,4
`do 12 i•1,66
`ma(i)-ma(i)-ma(1+mod(i+30,66))
`if(ma(i).lt.MZ)ma(i)-ma(i)+MBIG
`enddo 12
`enddo 13
`inext•O
`inextp•31
`idum•1
`end it
`inext•inext+1
`it(inext.eq.66)inext•1
`inextp•inextp+1
`if(inextp.eq.66)inextp•1
`mj-ma(inext)-ma(inextp)
`if(mj.lt.MZ)mj-mj+MBIG
`ma(inext)-mj
`ran3-mj•FAC
`return
`END
`
`Here is where we start, except on initialization. Increment
`inext, wrapping around 56 to 1.
`Ditto for inextp.
`
`Now generate a new random number subtractively.
`Be sure that it is in range.
`Store it,
`and output the derived uniform deviate.
`
`Quick and Dirty Generators
`
`One sometimes would like a "quick and dirty" generator to embed in a program, perhaps
`taking only one or two lines of code, just to somewhat randomize things. One might wish to
`process data from an experiment not always in exactly the same order, for example, so that
`the first output is more "typical" than might otherwise be the case.
`For this kind of application, all we really need is a list of "good" choices for m, a, and
`c in equation (7.1.1). If we don't need a period longer than 104 to 106
`, say, we can keep the
`value of ( m - 1 )a + c small enough to avoid overflows that would otherwise mandate the
`extra complexity of Schrage's method (above). We can thus easily embed in our programs
`
`jran-mod(jran•ia+ic,im)
`ran•tloat(jran)/tloat(im)
`
`whenever we want a quick and dirty uniform deviate, or
`
`jran-mod(jran•ia+ic,im)
`j•jlo+((jhi-jlo+1)•jran)/im
`
`whenever we want an integer between jlo and jhi, inclusive. (In both cases jran was once
`initialized to any seed value between 0 and im-1.)
`Be sure to remember, however, that when im is small, the kth root of it, which is the
`number of planes in k-space, is even smaller! So a quick and dirty generator should never
`be used to select points in k-space with k > 1.
`With these caveats, some "good" choices for the constants are given in the accompanying
`table. These constants (i) give a period of maximal length im, and, more important, (ii) pass
`Knuth's "spectral test" for dimensions 2, 3, 4, S, and 6. The increment ic is a prime, close to
`the value ( ~ - 1v'3)im; actually almost any value of ic that is relatively prime to im will do
`just as well, but there is some "lore" favoring this choice (see [4], p. 84).
`
`PGS Exhibit 2022
`WesternGeco v. PGS (IPR2015-00309, 310, 311)
`
`

`
`7. 1 Uniform Deviates
`
`275
`
`Constants for Quick and Dirty Random Number Generators
`
`overftow at
`
`220
`
`221
`
`222
`
`223
`
`224
`
`226
`
`226
`
`im.
`6075
`
`ia
`106
`
`ic
`1283
`
`7875
`
`211
`
`1663
`
`7875
`
`421
`
`1663
`
`6075 1366
`6655
`936
`11979 430
`
`1283
`1399
`2531
`
`14406
`29282
`53125
`
`3041
`967
`419
`6173
`171 11213
`
`2731
`12960 1741
`2957
`14000 1541
`4621
`21870 1291
`31104
`625
`6571
`139968
`205 29573
`
`6173
`29282 1255
`421 17117
`81000
`281 28411
`134456
`
`overftow at
`
`227
`
`228
`
`229
`
`230
`
`231
`
`232
`
`ia
`im
`86436 1093
`121500 1021
`259200
`421
`
`117128 1277
`121500 2041
`312500
`741
`
`145800 3661
`175000 2661
`233280 1861
`244944 1597
`
`ic
`18257
`25673
`54773
`
`24749
`25673
`66037
`
`30809
`36979
`49297
`51749
`
`139968 3877
`29573
`214326 3613
`45289
`714025 1366 150889
`
`134456 8121
`259200 7141
`
`28411
`54773
`
`233280 9301
`49297
`714025 4096 150889
`
`An Even Quicker and Dirtier Generator
`
`Many FORTRAN compilers can be abused in such a way that they will multiply two 32-bit
`integers ignoring any resulting overflow. In such cases, on many machines, the value returned
`is predictably the low-order 32 bits of the true 64-bit product. (C compilers, incidentally,
`can do this without the requirement of abuse -
`it is guaranteed behavior for so-called
`unsigned long int integers. On VMS VAXes, the ne~ssary FORTRAN command is
`FORTRAN/CHECK•NOOVERFLOW.) If we now choose m = 232 , the "mod" in equation (7.1.1)
`is free, and we have simply
`
`1;+1 = al; + c
`(7.1.6)
`Knuth suggests a= 1664525 as a suitable multiplier for this value of m. H.W. Lewis
`has conducted extensive tests of this value of a with c = 1013904223, which is a prime close
`to ( v'5- 2)m. The resulting in-line generator (we will call it ranqdl) is simply
`
`idum•1664626•idum+1013904223
`
`This is about as good as any 32-bit linear congruential generator, entirely adequate for many
`uses. And, with only a single multiply and add, it is very fast.
`To check whether your compiler and machine have the desired overftow proper(cid:173)
`ties, see if you can generate the following sequence of 32-bit values (given here in
`hex): 00000000, 3C6EF35F, 47502932, D1CCF6E9, AAF95334, 6252E503, 9F2EC686,
`57FE6C2D, A3D95FA8, S1FDBEE7, 94FOAF1A, CBF633Bl.
`If you need floating-point values instead of 32-bit integers, and want to avoid a divide by
`floating-point 232 , a dirty trick is to mask in an exponent that makes the value lie between 1 and
`2, then subtract 1.0. The resulting in-line generator (call it ranqd2) will look something like
`
`PGS Exhibit 2022
`WesternGeco v. PGS (IPR2015-00309, 310, 311)
`
`

`
`276
`
`Chapter 7.
`
`Random Numbers
`
`INTEGER idum,itemp,jflone,jflmsk
`REAL ftemp
`EQUIVALENCE (itemp,ftemp)
`DATA jflone /Z'3F800000'/, jflmsk /Z'007FFFFF'/
`
`c
`
`idum•1664626•idum+1013904223
`itemp•ior(jflone,iand(jflmsk,idum))
`ran•ftemp-1. 0
`
`The hex constants 3F800000 and 007FFFFF are the appropriate ones for computers using
`the IEEE representation for 32-bit floating-point numbers (e.g., IBM PCs and most UNIX
`workstations). For DEC VAXes, the correct hex constants are, respectively, 00004080 and
`FFFF007F. Notice that the IEEE mask results in the floating-point number being constructed
`out of the 23 low-order bits of the integer, which is not ideal. Also notice that your compiler
`may require a different no