Morgan - Numerical Methods (523161), страница 39
Текст из файла (страница 39)
The techniquesare given in books, but the art comes from experience with the arithmetic itself.Combining knowledge of how it behaves with science produces the best results.278THE ELEMENTARY FUNCTIONS1Horden, Ira. An FFT Algorithm For MCS-96 Products Including SupportingRoutines and Examples. Mt. Prospect, IL:Intel Corp., 1991, AP-275.2 Feynman, Richard P.
The Feynman Lectures On Physics. Reading, MA:Addison-Wesley Publishing Co., 1963, Volume I, Chapter 22.3Knuth, D. E. Fundamental Algorithms. Reading, MA: Addison-Wesley Publishing Co., 1973, Page 26, Exercise 28.4Jarvis, Pitts. Implementing CORDIC Algorithms. Dr. Dobb’s Journal, October1990, Pages 152-156.5Nielsen, Kaj L. Modern Trigonometry.
New York, NY: Barnes & Noble, 1966,Page 1696Acton, Forman S. Numerical Methods That Usually Work. Washington D.C.:Mathematical Association of America, 1990.Hamming, R. W. Numerical Methods for Scientists and Engineers. New York,NY: Dover Publications, 1973.7Sedgewick, Robert.
Algorithms in C. New York, NY: Addison-Welsley Publishing Co., 1990, Page 525.8Crenshaw, Jack W. Square Roots are Simple? Embedded Systems Programming, Nov. 1991, 4/11, Pages 30-52.9Cody, William J. and William Waite. Software Manual for the ElementaryFunctions. Englewood, NJ: Prentice-Hall, Inc., 1980.279280Previous Home NextAPPENDIX AA Pseudo-RandomNumber GeneratorTo test the floating-point routines in this book, I needed something that wouldgenerate an unpredictable and fairly uniform series of numbers. These routines arecomplex enough that a forgotten carry, incorrect two’s complement, or occasionaloverflow could easily hide from an ordinary “peek and poke” test.
Even with arandom number generator, it took many hours and tests with a number of data rangesto find some of the ugliest bugs.Of course, the standard C library has a random number generator, rand(), but thecode for it was unavailable and there were no guarantees as to how it worked. Somerandom number generators have such a high serial correlation (sequential dependence) that if a sequence of numbers was mapped to x/y locations on a monitor,patterns would appear. With others, users were warned that although each numbergenerated was guaranteed to be random individually, no sequence was guaranteedto be random.Generating random numbers isn’t as easy as it might sound.
Random numbersand arbitrary numbers are very different; if you asked a friend for a random number,you would really receive an arbitrarily chosen number. To be truly random, a numbermust have an equal chance of being chosen out of some known range and precision.Games of dice, cards, and the lottery all depend on a sequence of randomnumbers, and most use a means other than computers to generate them. People don’ttrust machines to generate random numbers because machines can become predictable and repetitive. But with the kind of simulations and testing needed to test thefloating-point routines in this book, drawing each number from a pot would take fartoo long.
Some other method had to be devised.281NUMERICAL METHODSOne of the first techniques for generating random numbers was originated byJohn Von Neumann and called the middle-square method. 1 It consisted of taking theseed, or previous random number, squaring it, and taking the middle digits.Unfortunately, this method had serious disadvantages that prevented it from beingwidely used. It didn’t take much for it to get into a rut; if a zero found its way intothese middle digits, for instance, there it would stay.A number of pseudo-random number generators are in general use, though notall of them are well tested and not all of them are good.
A good random numbergenerator is difficult to define exactly. The one quality that these generators mustpossess is randomness. An instance of this is given in the chi-square test, presentedlater. Given a uniformly distributed, pseudo-random sequence of a certain length, n,of numbers, all between 0 and some limit, l, divided among l bins, an equal numberof numbers in each bin would be highly suspicious.The most popular pseudo-random number generator in use, and the one chosenfor this book, is the multiplicative congruential method. This technique was firstproposed by D.
H. Lehmer1 in 1949. It is based on the formulaXn+l = (aXn + c) mod mEach new number is produced from a number, Xn, which is either the seed or theprevious number, through multiplication and modular division. It requires a multiplier, a, that must be equal to or greater than zero and less than the modulus, anadditive or increment, c, that must also be equal to or greater than zero and less thanthe modulus, and a modulus, m, that is greater than zero. Simply supplying numbersfor these variables won’t result in a good random number generator; the two “bad”generators described earlier were linear congruential generators.Here are a few guidelines, summarized from the materials of Donald Knuth:lThe seed, Xn, may be arbitrary and may, in fact, be the previously generatednumber in a pseudo-random sequence.
Irandom, the pseudo-random numbergenerator created for this book expects a double as the seed; in the demonstrationroutine spectral.c, the DOS timer tick is used.282A PSEUDO-RANDOM NUMBER GENERATORThe modulus, m, should be at least 230. Very often it is the word size (or a multiplethereof) of the computer, making division and extraction of the remainder trivial.The subroutine that actually produces the random number uses a modulus of 232.This means that after the seed is multiplied by a, the least significant doublewordis the new random number. The result would be the same if the product of a*Xwere divided by l00000000H.If you intend to run the random number generator on a binary computer, themultiplier, a, should be chosen; a mod 8=5.
If the target machine is decimal, thena mod 200 = 21. The multiplier and increment determine the period, or the lengthof the sequence before it starts again, and potency, or randomness, of the randomnumber generator. The multiplier, a, in irandom (presented later in this chapter)is 69069D, which is congruent to 5 mod 8.The multiplier should be between .0lm and .99m and should not involve a regularpattern.
The multiplier in irand is actually less than .0lm, but so was themultiplier in the original psuedo-random number generator proposed by Lehmer.In truth, it was chosen partly because of its size; the arithmetic was easier andfaster. In tests described later in this appendix, this multiplier performed as wellas those of two other generators.If you have a good multiplier, the value of c, the increment, is not important. Itmay be equal to one or even a. In irandom, c = 0.Beware of the least significant digits. They are not very random and should notbe used for decisions. Avoid methods of scaling random numbers that involvemodular operations, such as those found in the Microsoft C getrandom macro; themodular function will return the least random part of the number.
Instead, treatthe value returned by the random number generator as a fraction and use it to scalea user-determined maximum.The technique chosen for the random number generator here is a combination oflinear congruential and shuffling. In this sense, shuffling means that the randomnumbers are somehow moved around, or shuffled, before they’re generated. Thisbreaks up any serial correlation the sequence might have and provides a much longer,possibly infinite, period.283NUMERICAL METHODSThe pseudo-random number generator here comprises three routines. The firstis an initialization, rinit, in which an array of 256 doublewords is filled with numberscreated using the routine congruent and a seed value.The actual generation is done by irand.
First, this routine creates a new randomnumber based on the current seed, which is nothing more than the last numbergenerated. It then uses the lower byte of the upper word of this new random numberas an index into the array of 256 numbers created at initialization. A new randomnumber is created to replace the one selected and the routine exits, returning thenumber from the array. The initialization routine, rinit, must be called beforeirandom if the user wishes to select their own seeds; otherwise, the value 1 is chosen.Pseudocode for each of the routines is as followsrinit: Algorithm1.
Point to the double word array in RAM. This will be the initial list ofrandom numbers.2. Place the input seed in the seed variable. In these routines, the timertick is used as the seed.3. Call the routine congruent 256 times to fill the array.4. Exit.; ******;rinit - initializes random number generator based upon input seed.dataadwordIMAX equrantop wordran1 dwordxsubi dword2846906932767IMAX256 dup (0)lh;global iterative seed for;random number generator, change;this value to change defaultA PSEUDO-RANDOM NUMBER GENERATORinitbyteOh;global variable signaling;whether the generator has been;initialized or not.coderinit proc uses bx cx dx si di, seed:dwordleadi, word ptr ran1movmovmovmovax, wordword ptrax, wordword ptrmovcx, 256fill_array:invokemovmovaddlooprinit_exit:subnotmovretptr seed[2]xsubi[2], axptr seedxsubi, ax;put in seed variable;get seedcongruentword ptr [di], axword ptr [di][2], dxdi, 4fill_arrayax, axaxbyte ptr init, alrinit endpcongruent: Algorithm1.
Move the lower word of the seed, xsubi, into AX and multiply by the lowerword of the multiplier, a. This will produce a result in DX:AX, with theupperword of the product in DX. (This routine performs a multipleprecisionmultiply. This is a standard polynomial multiply; it is a bitsimpler and more direct because the multiplier is known.)2. Save the lower word of this product in BX and the upper word in CX.285NUMERICAL METHODS3. Place the upper word of the seed, xsubi, in AX and multiply by the lowerword of the multiplier, a.4. Add the lower word of the product of this last multiplication to the upperword of the product from the first multiplication, and propagate anycarries.5. Add to AX the lower word of xsubi, and to DX the upper word of xsubi.The multiplier used in this routine is 69069D, or 10dcdH. Themultiplications performed prior to this step all involved the lower word,0dcdH.
To multiply by 10000H, you need only shift the multiplicand 16places to the left and add it to the previous subproduct.6. Replace DX with BX, the LSW of the multiple-precision product. The MSWis discarded because it is purely overflow from any carries that havepropagated forward. Instead, the lesser words are used. They might beregarded as the fractional extension of any integer in the MSW.7. Write BX to the LSW of the seed and AX to the MSW.8.
Return.; ******;congruent -performs simple congruential algorithmcongruent proc uses bx cxmovax, wordmulword ptrmovbx, axmovcx, dxmovax, wordmulword ptraddax, cxadcdx, 0286ptr xsubia;a*seed (mod2ˆ32);lower word of result;upper wordptr xsubi[2]aaddax, word ptr xsubiadcmovmovdx, word ptr xsubi [2]dx, bxword ptr xsubi, bx;a multiplication by one is just;an add, right?A PSEUDO-RANDOM NUMBER GENERATORmovretcongruent endpword ptr xsubi [2], axirandom: Algorithm1. Point to the array of random numbers.2.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
