Morgan - Numerical Methods (523161), страница 40
Текст из файла (страница 40)
Call congruent for a new number based upon the last seed.3. Use the lower byte of the MSW of that number as an index into that array.4. Get a new random number.5. Replace the previously selected number with this new number.6. Replace the seed with the previously selected number.7.
Scale the random number with rantop, a variable defining the maximumrandom number output by the routine.;; ******;irandom- generates random floats using the linear congruential methodirandomproc uses bx cx dx si dileasi, word ptr ran1moval, byte ptr initoral, aljnealready_initializedinvokerinit, xsubialready_initialized:invokecongruentandax, 0ffhshlax, 1shlax, 1addsi, axmovdi, siinvokecongruentmovbx, word ptr [si];check for initialization;default to 1;get a random number;every fourth byte, right?;multiply by four;point to number in array;so we can put one there too287NUMERICAL METHODSmovmovmovmovmovcx, word ptr [si][2]word ptr [di], axword ptr [di] [2], dxword ptr xsubi, bxword ptr xsubi[2], cxmovmulax, bxword ptr rantopmovax, dxretirandom;get number from array;replace it with another;seed for next random;scale output by rantop, the;maximum size of the;random number;if rantop were made 0ffffH,;the value could be used;directly as a fractionendpThe danger with a pseudo-random number generator is that it will look quiteacceptable on paper but may fail to produce good numbers.
Spectra1.c provides twoways to test irandom or any pseudo-random number generator. One is quite simple,allowing examination of the output in graphic format so that the numbers producedby irandom can be checked visually for any patterns or concentrations. Any serialcorrelations that might arise can be detected using this method, but it is no proof ofk-space, or multidimensional, randomness.The other test is the traditional Chi-square statistic. The output of this formulacan give a probabilistic indication as to whether your random number generator istruly random. The actual formula is stated:v=1(Ys-npS)2/npSkbut for the purposes of this algorithm is stated:v=288l/n(yS2/ps)-nA PSEUDO-RANDOM NUMBER GENERATORThis formula merely evaluates a sequence and produces a value indicating howmuch the sequence diverged from a probable or expected distribution.Say you generate 1,000 numbers, a, all of them less than 100, b.
You then dividea among 100 bins, c, based on the value of the number; in other words, a randomnumber of 55 would go into bin 55, and a number such as 32 would go in bin 32. Youwould probably expect 10 numbers in each bin. Of course, a random numbergenerator will seldom have an absolutely even distribution; it wouldn’t be randomif it did.In fact, this statistic is only an indicator and can vary from sampling to samplingon the same generator. Tables can be used to interpret the numbers output by thisformula. A good rule of thumb is that the statistic should be close, but not too close,2to the number of bins-probably within 2* ( b ) .This statistic can vary. While you could roll 10 sevens in a row, it simply won’thappen very often.
A statistic that varies widely from 2* (b), consistently producesthe same value, or is extremely close to b might be suspect. I tested Irandom —alongwith rand(), a “portable” pseudo-random number generator written in C and a thirdparty routine found little difference in this statistic. It always remained relativelyclose to b, only occasionally straying outside 2* (b).Both the visual and the Chi-square test are incorporated into a program calledspectral.c (no relation to Knuth’s spectral test; it is so called merely because of itsvisual aspect). The program is simple: Pairs of random numbers are scaled and usedas x and y coordinates for pixels on a graphics screen; 10,000 pixels are generatedthis way.
Serial correlations can show up as a sawtooth pattern or other concentrations in the display. Otherwise, the display should show a fairly uniform array ofwhite dots similar to a starry night.After painting the screen, program retrieves the seed to generate a sequence ofnumbers for the Chi-square statistic. The result is then displayed.spectral: Algorithm1. Prepare the screen, turn the cursor off, put the video in EGA graphicsmode, and retrieve a structure containing the current video configuration.289NUMERICAL METHODS2. a) Use the timer tick as a seed and generate 20,000 pseudo-random numbers,using pairs as x/y locations on the graphics screen to turn pixels on.b) Use the same seed to generate the sequence for Chi-square analysis;output the result to the screen.c) Print a message asking for a keystroke to continue, "q" to quit.3.
Return the screen to its previous state and exit to DOS.#include<conio.h>#include<stdio.h>#include<graph.h>#include<stdlib.h>#include<time.h>short modes[ ] = { _TEXTBW40,_TBXTC80,_MRESNOCOLOR,_HRESBW,_MRES16COLOR,_ERESNOCOLOR,_ERESCOLOR,_VRES16COLOR,_MRES256COLOR,};externexternexternextern_TBXTC40,_MRBS4COLOR,_TEXTBW80,_TEXTMONO,_HRESlGCOLOR,_HERCMONO,_VRES2COLOR,_ORESCOLORint irandom (void);void rinit (int);uraninit (long);double urand (void);/*this routine scales a random number to a maximum without using a modularoperation*/int get random (int max)unsigned long a, b;a=irandom();b = max*a;return(b/32768);290A PSEUDO-RANDOM NUMBER GENERATORvoid main ()short j, ch, x, Y, row, num = sizeof(modes) /size of (modes[0]);unsigned int i, e;long g, c:double rnum;double result;int n = 20000;int r = 100;insigned int f[l000];unsigned int a, b, d;int seed;float chi;float pi=22.0/7.0;struct videoconfig vc;_displaycursor(_GCURSOROFF);/*set up screen*/_setvideomode(_ERESNOCOLOR);/*EGA mode; change thisif you have somethingelse.
The table is above*/_getvideoconfig(&vc);/*get the video configuration*/do{do{seed=(unsigned)time(NULL);/*use the timer tick as theseed*/rinit(seed);_clearscreen(_GCLEARSCREEN);for(i=0;i<l0000;i++) {/*draw a starrynight on the screen*/x=getrandom(vc.numxpixels);y=getrandom(vc.numypixels);_setpixel(x,y);291NUMERICAL METHODSrinit(seed);/*calculate x-square basedupon same seed as display*/for (a = 0; a < r; a++) f[a] = 0;for (a = 0; a < n; a++) {f[getrandom(r)]++;for (a = 0, c = 0; a < r; a++)c += f[a] * f[a];chi= ((float)r * (float)c/(float)n)-(float)n;printf("\n(irandom) chi-square statistic for this set ofof random numbers is %f", chi);printf("\npress a key to continue...");)while((ch=getch()) != 'q');}while(ch != 'q');_displaycursor(_GCURSORON);_setvideomode(_DEFAULTMODE);292A PSEUDO-RANDOM NUMBER GENERATORlKnuth, D.
E. Seminumerical Algorithms. Reading, MA: Addison-Wesley Publishing Co., 1981, Pages 1-178.2Sedgewick, Robert. Algorithms in C. Reading, MA: Addison-Wesley Publishing Co., 1990, Page 517.293294Previous HomeAPPENDIX BTables and EquatesExtended Precision Values Used in Elementary Functionszeroone_halfone_over_pitwo_over_sihalf_pione_over_ln2ln2sqrt_halfexpepsepsymaxbig-xlittlexY0aY0bquarterqwordqwordqwordqwordqwordqwordqwordqwordqwordqwordqwordqwordqwordqwordqwordqword000000000000h3f0000000000h3ea2f9836e4eh3f22f9836e4eh3fc90fdaa221h3fb8aa3b295ch3f317217f7dlh3f3504f30000h338000000001h39fffff70000h45c90fdb0000h42a000000000h0c2a000000000h3ed5a9a80000h3f1714ba0000h3e8000000000hConstants for Cordic Functionscircularkhyperkqword 9b74eda7hqword 1351e8755hCommon Values Written for Quadword Fixed Point1/pp2÷pe1/e=====0.3183098869.8696044011.7724538512.7182818280.367879441=====517cc1b7h9de9e64dfh1c5bf891bh2b7e15163h5e2d58d9h295NextNUMERICAL METHODS======e2p/180÷2ln(p)÷3p7.3890560990.0174532931.4142135621.1447298861.7320508083.141592654======763992e35h477dla9h16a09e668h1250d048ehlbb67ae86h3243f6a89hNegative Powers of Two in Decimal2-12-22-32-42-52-62-72-82-92-102-l12-12============.5D.25D.125D.0625D.03125D.015625D.0078125D.00390625D.001953125D.0009765625D.00048828125D.000244140625DNegative Powers of Ten in 32 bit Hex Format10-l10-210 -310 -410 -510 -610 -7l0-810-9296=========.1999999aH.028f5c29H.00418037H.00068db9H.0000a7c5H.000010c6H.00000ladH.0000002aH.00000004HPrevious HomeAPPENDIX CFXMATH.ASM.dosseg.model small, c, os-dosinclude math.inc;.code; ******; add64 -Adds two fixed-point numbers.
The radix point lies between;word 2 and word 3.;the arguments are passed on the stack along with a pointer to;storage for the resultadd64 proc uses ax dx es di, addend0:qword, addendl:qword, result:wordmovdi,word ptr resultmovmovaddadcmovmovax, worddx, wordax, worddx, wordword ptrword ptrptr addend0[0]ptr addend0[2]ptr addend1[0]ptr addendl[2][di], ax[di][2], dx;ax = low word, addend0;dx = highword, addend0;add low word, addend1;add high word, addend1movmovadcadcmovmovretadd64 endpax, worddx, wordax, worddx, wordword ptrword ptrptr addend0[4]ptr addend0[6]ptr addend1[4]ptr addend1[6][di][4], ax[di][6], dx;ax = low word, addend0;dx = high word, addend0;add low word, addend1;add high word, addend1297NextNUMERICAL METHODS;; * sub64;arguments passed on the stack; pointer returned to resultsub64 proc uses dx es di,sub0:qword, sub1:qword, result:wordmovdi,word ptr resultmovmovsubsbbmovmovax, worddx, wordax, worddx, wordword ptrword ptrmovmovsbbsbbmovmovax, word ptr sub0 [4]dx, word ptr sub0 [6]ax, word ptr sub1 [4]dx, word ptr sub1 [6]word ptr [di][4],axword ptr [di][6],dxmovjncnotno-flag:ptr sub0 [0]ptr sub0 [2]ptr sub1 [0]ptr sub1 [2][di][0],ax[di] [2],dx;ax = low word, sub0;dx = high word, sub0;subtract low word, ;sub1;subtract high word, ;sub1;ax = low word, sub0;dx = high word, sub0;subtract low word,;sub1;subtract high word, sub1a,0no-flagax;result returned as dx:axretsub64 endp;;* sub128;arguments passed on the stack; pointer returned to resultsub128 proc uses ax dx es di,sub0:word, sub1:word, result:word298movmovdi,word ptr sub0si,word ptr sub1movmovsubax, word ptr [di] [0]dx, word ptr [di][2]ax, word ptr [si] [0];ax = low word, [di];dx = high word, [di];subtract low word, [si]FXMATH.ASMsbbmovmovdx, word ptr [si][2]word ptr [di],axword ptr [di][2],dx;subtract high word, [si]movmovsbbsbbmovmovax, word ptr [di][4]dx, word ptr [di][6]ax, word ptr [si][4]dx, word ptr [si] [6]word ptr [di][4],axword ptr [di][6],dx;ax = low word, [di];dx = high word, [di];subtract low word, [si];subtract high word, [si]movmovsbbsbbmovmovax, worddx, wordax, worddx, wordword ptrword ptrptr [di][8]ptr [di][10]ptr [si][8]ptr [si] [10][di][8],ax[di][10],dx;ax = low word, [di];dx = high word, [di];subtract low word, [si];subtract high word, [si]movmovsbbsbbmovmovax, word ptr [di][12]dx, word ptr [di][14]ax, word ptr [si][12]dx, word ptr [si][14]word ptr [di] [12],axword ptr [di][14],dx;ax = low word, [di];dx = high word, [di];subtract low word, [si];subtract high word, [si]movmovmovmovswrepretsub128 endpsi,didi,word ptr resultcx,8;result returned as dx:ax;*mullong - Multiplies two unsigned fixed point values.