KISS4691, a potentially top-ranked RNG.
in [Fortran]
|
Prev: solutions book
Next: real kind declaration
From: Uno on 30 Jul 2010 04:33 robin wrote: > "jacob navia" <jacob(a)spamsink.net> wrote in message news:i2fir2$op4$1(a)speranza.aioe.org... > > | This doesn't work with systems that have unsigned long as a 64 bit quantity. > | > | I obtain: > | > | MWC result=3740121002 ? > | 4169348530 > | KISS result=2224631993 ? > | 1421918629 > > For a 64-bit machine (using 64-bit integer arithmetic), > you'd need to truncate each result to 32 bits. That not > only applies to the multiplication, it also applies to addition, etc. > On a 32-bit machine, these extra bits are discarded, > but in 64-bit arithmetic, they are retained, > and unless they are similarly discarded, > you won't get the same results. > I suggest using IAND(k, 2*2147483647+1) > for the truncation. > > With such modifications in the program, > it should then produce the same results on both 32-bit and > 64-bit machines. > > P.S. the product 2*2... is best obtained using ISHFT. > > | Compiling with 32 bit machine yields: > | MWC result=3740121002 ? > | 3740121002 > | KISS result=2224631993 ? > | 2224631993 > > > First of all, I think we're talking to the actual George Marsaglia here. Thank you so much for posting. You may have displaced Terence as our senior member. Are you creating bigger numbers just to accomodate your age?:-) $ gcc -Wall -Wextra geo1.c -o out geo1.c: In function �MWC�: geo1.c:5: warning: type defaults to �int� in declaration of �c� geo1.c:5: warning: type defaults to �int� in declaration of �j� geo1.c:5: warning: unused variable �i� geo1.c: In function �main�: geo1.c:21: warning: format �%22u� expects type �unsigned int�, but argument 2 has type �long unsigned int� geo1.c:23: warning: format �%22u� expects type �unsigned int�, but argument 2 has type �long unsigned int� geo1.c:24: warning: control reaches end of non-void function $ ./out MWC result=3740121002 ? 3740121002 KISS result=2224631993 ? 2224631993 $ cat geo1.c static unsigned long xs=521288629,xcng=362436069,Q[4691]; unsigned long MWC(void) /*takes about 4.2 nanosecs or 238 million/ second*/ {unsigned long t,x,i; static c=0,j=4691; j=(j<4690)? j+1:0; x=Q[j]; t=(x<<13)+c+x; c=(t<x)+(x>>19); return (Q[j]=t); } #define CNG ( xcng=69069*xcng+123 ) #define XS ( xs^=(xs<<13), xs^=(xs>>17), xs^=(xs<<5) ) #define KISS ( MWC()+CNG+XS ) /*138 million/sec*/ #include <stdio.h> int main() {unsigned long i,x; for(i=0;i<4691;i++) Q[i]=CNG+XS; for(i=0;i<1000000000;i++) x=MWC(); printf(" MWC result=3740121002 ?\n%22u\n",x); for(i=0;i<1000000000;i++) x=KISS; printf("KISS result=2224631993 ?\n%22u\n",x); } // gcc -Wall -Wextra geo1.c -o out $ So, what is all this? In particular, is there something special about the value of 3.7 billion? -- Uno
From: Uno on 30 Jul 2010 04:58 glen herrmannsfeldt wrote: > In comp.lang.fortran Ron Shepard <ron-shepard(a)nospam.comcast.net> wrote: >> In article <i2ssst$nfi$1(a)speranza.aioe.org>, >> glen herrmannsfeldt <gah(a)ugcs.caltech.edu> wrote: > >>> In a large fraction of the cases, 2 billion different seeds >>> are enough, but one can still desire the appropriate randomness >>> from those different seeds. > >> My understanding is that pseudorandom number generators return a >> sequence of values with some period. The numbers in that sequence >> eventually repeat with some, hopefully long, cycle. The put=array >> argument gives the starting point within that sequence, but it doesn't >> affect the cycle length or the "randomness" of the values, those things >> are determined by the underlying algorithm, not by the initial seed. > >> The situation that needs to be avoided is to run the code once with one >> seed, and then run it again with another seed that results in an overlap >> of the sequences of values for the two runs. In some applications this >> is unimportant, but in other applications it would be bad for two >> supposedly independent runs to have a long sequence of identical values. >> It would be nice if there were some prescription to generate initial >> seeds that would avoid this situation. However, I don't really know how >> this could be specified in the fortran standard without specifying the >> algorithm itself (which is not done, apparently on purpose). Anyone >> have any ideas how this could be done? > > Here is the suggestion. Add to RANDOM_SEED a form that allows > a single integer variable to be supplied as a seed source. > > As far as I know, there is no requirement on period or otherwise > on the quality of the PRNG used, but the standard could suggest > that good seeds be selected from that integer. For example, > they could be chosen such that they have a fairly long sequence > before they overlap the sequence generated by another input > value. If the RNG only has 32 bits of state, or maybe less, then > it likely takes it directly. > > I don't know the math of the newer PRNGs enough to say, but I > think it is possible to do that. > > With the current standard there is no way to even suggest to > RANDOM_SEED that you want a good seed. > > OK, here is another suggestion that could be implemented without > any change in the standard. In the case where the supplied > seed array has all except the first element zero, then select > a good seed based on that first element. Ron asks an interesting question, and I think you shoot from the hip and may have missed on this one. First of all, when crossposted to clc and sci.math, there is no "standard" that does not require disambiquity. Given what it is, we might say "the C standard" when talking about C, or its lurking uncle, the fortran standard. I don't agree that you need to "add" to random_seed a form that allows a single integer variable supplied as a seed source. Any fortran program that implements this will be fortran-conforming by virtue of there being little specificity in either standard on this and no appetite to change it in either community. "disambiquity:" I kind of like that word. -- Uno
From: David Bernier on 30 Jul 2010 05:16 Ron Shepard wrote: [...] > The situation that needs to be avoided is to run the code once with one > seed, and then run it again with another seed that results in an overlap > of the sequences of values for the two runs. In some applications this > is unimportant, but in other applications it would be bad for two > supposedly independent runs to have a long sequence of identical values. I don't know how likely it would be for two seeds to produce overlap in a short time. But I like using seeds that I've almost surely never used before. I happen to use the Mersenne Twister algorithm, but getting new "random" seeds is always good for any algorithm. > It would be nice if there were some prescription to generate initial > seeds that would avoid this situation. However, I don't really know how > this could be specified in the fortran standard without specifying the > algorithm itself (which is not done, apparently on purpose). Anyone > have any ideas how this could be done? [...] I often want to do just that: use new seeds. I used to go to http://www.random.org/ and get random bits. The way I do things now is to have a file "data" and take a sha256 hash value, which is presented in hexadecimal. To avoid reusing the same seed, I then append the hash-value to the file "data" . This can be automated with a script which uses the "bash" shell (popular in Linux); for a C-shell, something similar should work. [The next step would be to be to access the script from a fortran or C program ... ] listing of rnd.sh file (3 lines long, made executable by: chmod u+x rnd.sh ) #! /bin/bash sha256sum wcisha >> wcisha sha256sum wcisha The rnd.sh file is in my home directory, as well as the "data" file 'wcisha'. In use: [ ~]$ ./rnd.sh [Enter] 52897abc80205035b14986dd8d4d7b1549b98dca8c959f944f43b7c535e4f059 wcisha [ ~]$ ./rnd.sh [Enter] 6822ad1bd1d1e453e66d2b65ea79f82a86a1ae53f90512b08f1feedcce6ae9f0 wcisha David Bernier
From: Uno on 30 Jul 2010 06:46 robin wrote: > "geo" <gmarsaglia(a)gmail.com> wrote in message news:a82cebe3-cdb9-48af-8080-bca935eeb9b1(a)l14g2000yql.googlegroups.com... > |I have been asked to recommend an RNG > | (Random Number Generator) that ranks > | at or near the top in all of the categories: > | performance on tests of randomness, > | length of period, simplicity and speed. > | The most important measure, of course, is > | performance on extensive tests of randomness, and for > | those that perform well, selection may well depend > | on those other measures. > > I have already posted a PL/I version using unsigned arithmetic. > > Here is another version, this time using signed arithmetic :-- > > (NOSIZE, NOFOFL): > RNG: PROCEDURE OPTIONS (MAIN, REORDER); > > declare (xs initial (521288629), xcng initial (362436069), > Q(0:4690) ) static fixed binary (31); > > MWC: procedure () returns (fixed binary (31)); > declare (t,x,i) fixed binary (31); > declare (c initial (0), j initial (4691) ) fixed binary (31) static; > declare (t1, t2, t3) fixed binary (31); > > if j < hbound(Q,1) then j = j + 1; else j = 0; > x = Q(j); > t = isll(x,13)+c+x; > t1 = iand(x, 3) - iand(t, 3); > t2 = isrl(x, 2) - isrl(t, 2); > if t2 = 0 then t2 = t1; > if t2 > 0 then t3 = 1; else t3 = 0; > c = t3 + isrl(x, 19); > Q(j)=t; > return (t); > end MWC; > > CNG: procedure returns (fixed binary (31)); > xcng=bin(69069)*xcng+bin(123); > return (xcng); > end CNG; > > XXS: procedure returns (fixed binary (31)); > xs = ieor (xs, isll(xs, 13) ); > xs = ieor (xs, isrl(xs, 17) ); > xs = ieor (xs, isll(xs, 5) ); > return (xs); > end XXS; > > KISS: procedure returns (fixed binary (31)); > return ( MWC()+CNG+XXS ); > end KISS; > > declare (i,x) fixed binary (31); > declare y fixed decimal (11); > > Q = CNG+XXS; /* Initialize. */ > do i = 1 to 1000000000; x=MWC(); end; > put skip edit (" Expected MWC result = 3740121002", 'computed =', x) > (a, skip, x(12), a, f(11)); > y = iand(x, 2147483647); > if x < 0 then y = y + 2147483648; > put skip edit (y) (x(11), f(22)); put skip; > do i = 1 to 1000000000; x=KISS; end; > put skip edit ("Expected KISS result = 2224631993", 'computed =', x) > (a, skip, x(12), a, f(11)); > y = iand(x, 2147483647); > if x < 0 then y = y + 2147483648; > put skip edit (y) (x(11), f(22)); > > end RNG; > > If you were to comment out the PL/I command line that compiled this, what would it be? -- Uno
From: robin on 30 Jul 2010 10:09
"glen herrmannsfeldt" <gah(a)ugcs.caltech.edu> wrote in message news:i2trdk$rip$1(a)speranza.aioe.org... | In comp.lang.fortran Gib Bogle <g.bogle(a)auckland.no.spam.ac.nz> wrote: | (snip) | | >> unsigned_result = 3740121002 | 1 | >> Error: Integer too big for its kind at (1). This check can be disabled | >> with the option -fno-range-check | | You can subtract 2**32, or find two values to IOR together. Subtracting 2**32 won't do anything because it is entirely put of range (i.e., requires 33 bits. The magic number is to subtract or add 2**31. |