Author Message
George Marsaglia

Joined: 19 Nov 2009
Posts: 5
Location: United States Posted: Sun Nov 22, 2009 5:49 pm    Post subject: Help for a 64-bit RNG I have recently, Nov 3: RNGs: A Super KISS, sci.math, comp.lang.c, sci.crypt. developed a 32-bit KISS RNG with immense period, over 10^402575, great speed and excellent performance on tests of randomness. But it is not well-suited for extension to 64-bits or adaption to Fortran. As 64-bit RNGs are becoming more in demand, I have modified that recently proposed 32-bit one to produce a 64-bit RNG with similar period, this time a mere 10^397524, but with instructions that should serve in diverse programming languages. I list a C version below and ask for help from experts here on adapting it to a Fortran version. The basic idea is that we add each term from three sequences, Congruential, Xorshift and CMWC (Complimentary-Multiply-With-Carry). For the latter, we merely take the next Q element, except if they have all been used, the Q array is refilled, the index reset and the first Q element returned. The refilling process should work for signed or unsigned integers, as should the Congruential and Xorshift elements. The CMWC process requires that we extract the top and bottom 64 bits of a supposed 128-bit t=a*x+c, without actually forming t, using only 64-bit operations. The instructions, operations and assignments themselves all seem to have Fortran equivalents. The #include macros are merely sequences of assignments that can be transformed, and the Q() (nee Q[]) array must be accessible to all, which I would have done in the old days with common blocks. Here is the C listing: ----------------------------------------------------------- /* SUPRKISS64.c, period 5*2^1320480*(2^64-1) */ #include static unsigned long long Q,carry=36243678541LL, xcng=12367890123456LL,xs=521288629546311LL,indx=20632; #define CNG ( xcng=6906969069LL*xcng+123 ) #define XS ( xs^=xs<<13,xs^=xs>>17,xs^=xs<<43 ) #define SUPR ( indx<20632 ? Q[indx++] : refill() ) #define KISS SUPR+CNG+XS unsigned long long refill( ) {int i; unsigned long long z,h; for(i=0;i<20635;i++){ h=(carry&1); z=((Q[i]<<41)>>1)+((Q[i]<<39)>>1)+(carry>>1); carry=(Q[i]>>23)+(Q[i]>>25)+(z>>63); Q[i]=~((z<<1)+h); } indx=1; return (Q); } int main() {int i; unsigned long long x; for(i=0;i<20632;i++) Q[i]=CNG+XS; for(i=0;i<1000000000;i++) x=KISS; printf("Does x=-5061657174865032461\n x=%LLd.\n",x); } ----------------------------------------------------------- Given the initial seed values, a proper implementation of the mathematics behind the methods requires that the 10^9'th call return the 64-bit integer 13385086898844519155, which I have had printed as -5061657174865032461 for signed-integer implementations. Can I count on any of you to provide a Fortran version? George Marsaglia Here is the 32-bit version of the above: ---------------------------------------------------------- /* SUPRKISS32.c, period 5*2^1320481*(2^32-1) */ #include static unsigned long Q,indx=41265, carry=362,xcng=1236789,xs=521288629; #define CNG ( xcng=69609*xcng+123 ) /*Congruential*/ #define XS ( xs^=xs<<13,xs^=xs>>17,xs^=xs>>5 )/*Xorshift*/ #define SUPR ( indx<41265 ? Q[indx++] : refill() ) /*CMWC*/ #define KISS SUPR+CNG+XS int refill( ) {int i; unsigned long z,h; for(i=0;i<41265;i++){ h=(carry&1); z=((Q[i]<<9)>>1)+((Q[i]<<7)>>1)+(carry>>1); carry=(Q[i]>>23)+(Q[i]>>25)+(z>>31); Q[i]=~((z<<1)+h); } indx=1; return (Q); } int main( ) {unsigned long i,x; for(i=0;i<41265;i++) Q[i]=CNG+XS; for(i=0;i<1000000000;i++) x=KISS; printf("Does x=1493864468?\n x=%d.\n",x); } -----------------------------------------------------------    mecej4

Joined: 31 Oct 2006
Posts: 1084 Posted: Mon Nov 23, 2009 3:31 am    Post subject: Fortran-90 conversion of Marsaglia SuperKiss64 RNG Professor Marsaglia,

We appreciate your contributions to the art of RNGs.

This is a rather mechanical translation of your C code. It is likely to be slightly less efficient.

 Code: module suprkiss64_M integer,parameter :: I8=selected_int_kind(18)       ! 10^18 fits into 8 bytes integer(kind=I8) :: Q(20632),carry=36243678541_I8, &    xcng=12367890123456_I8,xs=521288629546311_I8,indx=20633_I8 contains function refill() result(s)    integer(kind=I8) :: s    integer :: i    integer(kind=I8) :: z,h    do i = 1,20635       h = iand(carry,1_I8)       z = ishft(ishft(Q(i),41),-1)+ &           ishft(ishft(Q(i),39),-1)+ &           ishft(carry,-1)       carry=ishft(Q(i),-23)+ishft(Q(i),-25)+ishft(z,-63)       Q(i)=not(ishft(z,1)+h)    end do    indx=2    s=Q(1)    return end function refill end module suprkiss64_M program tskiss64 use suprkiss64_M integer :: i integer(kind=I8) :: x,supr do i=1,20632    xcng=xcng*6906969069_I8+123    xs=ieor(xs,ishft(xs,13))    xs=ieor(xs,ishft(xs,-17))    xs=ieor(xs,ishft(xs,43))    Q(i)=xcng+xs end do do i=1,1000000000_I8    if(indx <= 20632)then       supr=Q(indx)       indx=indx+1    else       supr=refill()    endif    xcng=xcng*6906969069_I8+123    xs=ieor(xs,ishft(xs,13))    xs=ieor(xs,ishft(xs,-17))    xs=ieor(xs,ishft(xs,43))    x=xcng+xs+supr end do write(*,10)x 10 format(' Does x = -5061657174865032461 ?',/,6x,'x = ',I20) end program tskiss64

Last edited by mecej4 on Sun Dec 04, 2016 6:56 pm; edited 1 time in total   George Marsaglia

Joined: 19 Nov 2009
Posts: 5
Location: United States Posted: Mon Nov 23, 2009 2:49 pm    Post subject: Many thanks for that quick response. From different errors I got after trying to compile with both Silverfrost, (ftn95) and Nag (f95), I guessed at a possible solution and changed the three occurrences of (kind=8) to (kind=I8) (lines 7,9,28), and got successful compilations for both. The resulting NAG a.exe file ran and the 10^9th value was the same as that of the C program and consistent with the underlying mathematics. However, on this Windows Vista PC running in command mode, I was not able to link the .obj file produced by the Silverfrost compiler. Is there an easy way to ensure linkability? As for Fortran versions of this SuperKISS64, I would like to take advantage of statement functions---for example, for the Xorshift operations xs=ieor(xs,ishft(xs,13)) xs=ieor(xs,ishft(xs,-17)) xs=ieor(xs,ishft(xs,43)) I would use a statement function such as m(y,k)=ieor(ishft(y,k)) so that the code for the Xorshift operation would become xs=m(m(m(xs,13),-17),43) I liked statement functions in the old days when I used Fortran, but have read that they are currently frowned on. Would their use here be a bad choice? If I could link with the Silverfrost compiler I would have another comparison of speeds. The NAG a.exe file took around 33 seconds, so that RNGs were produced at around 30 million per second, whereas the gnu gcc compiler took 20 seconds, for a rate of 50 million 64-bit random integers per second. (The 32-bit RNG took 10 secs for a rate of 100 million/sec.) I plan to post a general description of this RNG at math.sci and sites that seem likely to make it a useful contribution, along with suggested code. I am already grateful for the help here, would welcome further comments that might lead to a better posting. George Marsaglia (And I was pleased with the output appearance of that newly downloaded Silverfrost FTN95,...if I could only figure out how to link the .obj file.)    mecej4

Joined: 31 Oct 2006
Posts: 1084 Posted: Mon Nov 23, 2009 7:54 pm    Post subject: Errors in first Fortran-9x version Mea culpa. I wrote up the code with "kind=8" for GFortran in Linux. I then did a manual conversion to kind=I8 for compilers whose kind number for 8 byte integers need not be 8. My conversion was incomplete. However, as you have already noted, the fix is easy -- replace all instances of integer(kind= to integer(kind=I . There is a problem with running the code using Salford Fortran. I have Windows XP rather than Vista. I can compile and link with no problems: s:\RNG> ftn95 supks64.f90 /optim /p6 s:\RNG> slink supks64.obj However, the program aborts right away with an integer overflow on the first executable statement xcng=xcng*6906969069_I8+123 Further examination showed that Salford Fortran-95 implements arithmetic on 8-byte integers using the FPU. Nor could I find the proper compiler switches to turn off the overflow checking. There are no problems running the code with other compilers. My best Fortran-90 timings are only 10 percent longer than my best C timings. As for using an ASF, I feel that they make the code more readable; as they are on their way out, I'd stay from them in writing new code.   mecej4

Joined: 31 Oct 2006
Posts: 1084 Posted: Mon Nov 23, 2009 8:50 pm    Post subject: Another minor bug in C and Fortran versions The array Q is dimensioned  in your original C version as well as my Fortran translation. The array is overrun in routine refill(), since the loop runs to 20635. The two limits need to be reconciled.Last edited by mecej4 on Sun Dec 04, 2016 7:08 pm; edited 1 time in total   George Marsaglia

Joined: 19 Nov 2009
Posts: 5
Location: United States Posted: Tue Nov 24, 2009 2:57 pm    Post subject: Yes, my typing slip, 20635 rather than 20632, needs correction. I have fixed it and got your Fortran version and the C version to give the same results. Learning of the -O4 option on NAG's f95, I get the same +%10% time cost versus C's gcc -O3 compiler. I am planning to post both C and Fortran versions on sites that seem suitable, but I am still puzzled on a compatible Fortran version of the 32-bit C listing. Here I need only ensure that integer arithmetic is the old standard mod 2^32, with positive residues for signed integers and least-absolute residues for signed integers. Can the default integer*4 still be used, and if not, which selected_int_kind(?) will---or is likely to---ensure compatibility?    mecej4

Joined: 31 Oct 2006
Posts: 1084 Posted: Tue Nov 24, 2009 8:07 pm    Post subject: Fortran version of 32-bit Super-KISS RNG Just as with the 64-bit RNG, the Fortran-90 version is straightforward. I have checked it out with both 32-bit and 64-bit Fortran compilers (both of which, however, have default 4-byte integers).

 Code: module suprkiss32_M integer,parameter :: I4=selected_int_kind(9)       ! 10^9 fits into 4 bytes integer(kind=I4),parameter :: QSIZ=41265_I4, &                               CMUL=69609_I4, &                               COFFS=123_I4 integer(kind=I4) :: Q(QSIZ),carry=362_I4, &    xcng=1236789_I4,xs=521288629_I4,       &    indx=QSIZ+1 contains function refill() result(s)    integer(kind=I4) :: s    integer :: i    integer(kind=I4) :: z,h    do i = 1,QSIZ       h = iand(carry,1_I4)       z = ishft(ishft(Q(i),9),-1)+ &           ishft(ishft(Q(i),7),-1)+ &           ishft(carry,-1)       carry=ishft(Q(i),-23)+ishft(Q(i),-25)+ishft(z,-31)       Q(i)=not(ishft(z,1)+h)    end do    indx=2    s=Q(1)    return end function refill end module suprkiss32_M program tskiss32 use suprkiss32_M integer :: i integer(kind=I4) :: x,supr do i=1,QSIZ    xcng=xcng*CMUL+COFFS    xs=ieor(xs,ishft(xs,13))    xs=ieor(xs,ishft(xs,-17))    xs=ieor(xs,ishft(xs,-5))    Q(i)=xcng+xs end do do i=1,1000000000_I4    if(indx <= QSIZ)then       supr=Q(indx)       indx=indx+1    else       supr=refill()    endif    xcng=xcng*CMUL+COFFS    xs=ieor(xs,ishft(xs,13))    xs=ieor(xs,ishft(xs,-17))    xs=ieor(xs,ishft(xs,-5))    x=xcng+xs+supr !   write(*,'(1x,i2,3I12,I14)')i,xcng,xs,supr,x !   if(i.gt.5)stop end do write(*,10)x 10 format(' Does x = 1493864468 ?',/,6x,'x =',I11) end program tskiss32   George Marsaglia

Joined: 19 Nov 2009
Posts: 5
Location: United States Posted: Wed Nov 25, 2009 3:14 pm    Post subject: I have fashioned a 64-bit version patterned after your model, with a few uses of semicolons to permit multiple expressions on a line, and with a KISS64 function in the module. The 32-bit version is similar, with a KISS32 function. Both are listed below. Your suggested use of QSIZ,CMUL,COFFS looked attractive as a way to make both versions look even more alike, but I encountered problems, so stuck with specific values for the 64 and 32 versions. I plan to make available both 32- and 64-bit versions for both C and Fortran to various newsgroups. May I acknowledge your contribution to the Fortran versions? Do you or other readers of this forum have any suggestions for improvements or otherwise? George Marsaglia The Fortran listings: [code:1:b0f6883a34] module suprkiss64_M ! period 5*2^1320480*(2^64-1) integer,parameter :: I8=selected_int_kind(18) integer(kind=I8) :: Q(20632),carry=36243678541_I8, & xcng=12367890123456_I8,xs=521288629546311_I8,indx=20633_I8 contains function KISS64() result(x) integer(kind=I8) :: x if(indx <= 20632)then; x=Q(indx); indx=indx+1 else; x=refill(); endif xcng=xcng*6906969069_I8+123 xs=ieor(xs,ishft(xs,13)) xs=ieor(xs,ishft(xs,-17)) xs=ieor(xs,ishft(xs,43)) x=x+xcng+xs return; end function KISS64 function refill() result(s) integer(kind=I8) :: i,s,z,h do i=1,20632 h=iand(carry,1_I8) z = ishft(ishft(Q(i),41),-1)+ & ishft(ishft(Q(i),39),-1)+ & ishft(carry,-1) carry=ishft(Q(i),-23)+ishft(Q(i),-25)+ishft(z,-63) Q(i)=not(ishft(z,1)+h) end do indx=2; s=Q(1) return; end function refill end module suprkiss64_M program testKISS64 use suprkiss64_M integer(kind=I8) :: i,x do i=1,20632 !fill Q with Congruential+Xorshift xcng=xcng*6906969069_I8+123 xs=ieor(xs,ishft(xs,13)) xs=ieor(xs,ishft(xs,-17)) xs=ieor(xs,ishft(xs,43)) Q(i)=xcng+xs end do do i=1,1000000000_I8; x=KISS64(); end do write(*,10) x 10 format(' Does x = 4013566000157423768 ?',/,6x,'x = ',I20) end program testKISS64 module suprkiss32_M ! period 5*2^1320481*(2^32-1) integer,parameter :: I4=selected_int_kind(9) integer(kind=I4) :: Q(41265),carry=362_I4, & xcng=1236789_I4,xs=521288629_I4,indx=41266_I4 contains function KISS32() result(x) integer(kind=I4):: x if(indx <= 41265)then;x=Q(indx); indx=indx+1 else; x=refill(); endif xcng=xcng*69069_I4+123 xs=ieor(xs,ishft(xs,13)) xs=ieor(xs,ishft(xs,-17)); xs=ieor(xs,ishft(xs,5)) x=x+xcng+xs return; end function KISS32 function refill() result(s) integer(kind=I4) :: i,s,z,h do i = 1,41265 h = iand(carry,1_I4) z = ishft(ishft(Q(i),9),-1)+ & ishft(ishft(Q(i),7),-1)+ & ishft(carry,-1) carry=ishft(Q(i),-23)+ishft(Q(i),-25)+ishft(z,-31) Q(i)=not(ishft(z,1)+h) end do indx=2; s=Q(1) return; end function refill end module suprkiss32_M program testKISS32 use suprkiss32_M integer(kind=I4) :: i,x do i=1,41265 !fill Q with Congruential+Xorshift xcng=xcng*69069_I4+123 xs=ieor(xs,ishft(xs,13)) xs=ieor(xs,ishft(xs,-17)[/c    mecej4

Joined: 31 Oct 2006
Posts: 1084 Posted: Sun Dec 04, 2016 8:36 pm    Post subject: George Marsaglia's post revisited, with corrected source cod Professor George Marsaglia (1924 - 2011), famous for his paper "Random numbers fall mainly in the planes" and other contributions, initiated this thread in 2009, when he was over 85 years old. See https://en.wikipedia.org/wiki/George_Marsaglia . Today, I was brought back to this thread after a search for some information on random number generators. I noticed that the Forum daemons had truncated his last post, in mid line-of-code even, and that there were still some typos in the programs. Here are corrected versions of his two programs. Both programs are written in the form of benchmark programs, so they need tweaking if you want to use the contained KISS RNGs in your work. He builds up a table using an RNG, and uses that table as a cache for generating the user's requested random numbers using the RNG along with lots of shuffling of bits and bytes. Note that these programs depend on integer overflow not being trapped, and are NOT standard-conforming. There are, however, useful under adult supervision. The 32-bit RNG, which uses only 4-byte integers, and can be used with FTN95-32 or FTN95-64, follows. [continued]Last edited by mecej4 on Mon Dec 05, 2016 4:56 pm; edited 3 times in total   mecej4

Joined: 31 Oct 2006
Posts: 1084 Posted: Sun Dec 04, 2016 8:37 pm    Post subject: [continuation] Here is the 32-bit RNG program of Marsaglia.
 Code: ! Posted by the late Professor George Marsaglia to the Silverfrost General forum ! in 2009, see http://forums.silverfrost.com/viewtopic.php?t=1480 . ! Works correctly with FTN95 32-bit and 64-bit. With Gfortran 4.9 and later, must ! use the -ftrapv flag to avoid an optimizer bug. !     Module suprkiss32_m ! period 5*2^1320481*(2^32-1)       Implicit None       Integer, Parameter :: i4 = selected_int_kind(9) ! 10^9 fits into 4 bytes       Integer (Kind=i4), Parameter :: QSIZ = 41265_i4, CMUL = 69609_i4, &         coffs = 123_i4       Integer (Kind=i4) :: q(41265), carry = 362_i4, xcng = 1236789_i4, &         xs = 521288629_i4, indx = QSIZ+1     Contains       Function kiss32() Result (x)         Implicit None         Integer (Kind=i4) :: x         If (indx<=qsiz) Then           x = q(indx)           indx = indx + 1         Else           x = refill()         End If         xcng = xcng*CMUL + COFFS         xs = ieor(xs, ishft(xs,13))         xs = ieor(xs, ishft(xs,-17))         xs = ieor(xs, ishft(xs,-5))         x = x + xcng + xs         Return       End Function kiss32       Function refill() Result (s)         Implicit None         Integer (Kind=i4) :: i, s, z, h         Do i = 1, qsiz           h = iand(carry, 1_i4)           z = ishft(ishft(q(i),9), -1) + ishft(ishft(q(i),7), -1) + &             ishft(carry, -1)           carry = ishft(q(i), -23) + ishft(q(i), -25) + ishft(z, -31)           q(i) = not(ishft(z,1)+h)         End Do         indx = 2         s = q(1)         Return       End Function refill     End Module suprkiss32_m     Program testkiss32       Use suprkiss32_m       Implicit None       Integer (Kind=i4) :: i, x, supr       Do i = 1, qsiz !fill Q with Congruential+Xorshift         xcng = xcng*CMUL + COFFS         xs = ieor(xs, ishft(xs,13))         xs = ieor(xs, ishft(xs,-17))         xs = ieor(xs, ishft(xs,-5))         q(i) = xcng + xs       End Do       Do i = 1, 1000000000_i4         x = kiss32()       End Do       Write (*, 100) x 100   Format (' Does x =  1493864468 ?', /, 6X, 'x =', I12)     End Program testkiss32

[continued]

Last edited by mecej4 on Mon Dec 05, 2016 12:53 pm; edited 2 times in total   mecej4

Joined: 31 Oct 2006
Posts: 1084 Posted: Sun Dec 04, 2016 8:42 pm    Post subject: [continuation]
Here is the 64 bit version of George Marsaglia's KISS RNG. It employs 8-byte integers, and works with FTN95-64. It does not work with FTN95-32 because INTEGER*8 arithmetic is implemented with X87 instructions and integer overflows cause FPU exceptions.
 Code: ! Posted by the late Professor George Marsaglia to the Silverfrost General forum ! in 2009, see http://forums.silverfrost.com/viewtopic.php?t=1480 . ! Does not work with FTN95 32-bit, since it uses the X87 to process 8-byte integers ! Works correctly with FTN95 64-bit. With Gfortran 4.9 and later, must ! use the -ftrapv flag to avoid an optimizer bug. !     Module suprkiss64_m ! period 5*2^1320480*(2^64-1)       Implicit None       Integer, Parameter :: i8 = selected_int_kind(18)       Integer, Parameter :: QSIZ = 20632       Integer (Kind=i8), Parameter :: CMUL = 6906969069_i8, COFFS = 123_i8       Integer (Kind=i8) :: q(QSIZ), carry = 36243678541_i8, &         xcng = 12367890123456_i8, xs = 521288629546311_i8, indx = QSIZ+1_i8     Contains       Function kiss64() Result (x)         Implicit None         Integer (Kind=i8) :: x         If (indx<=QSIZ) Then           x = q(indx)           indx = indx + 1         Else           x = refill()         End If         xcng = xcng*CMUL + COFFS         xs = ieor(xs, ishft(xs,13))         xs = ieor(xs, ishft(xs,-17))         xs = ieor(xs, ishft(xs,43))         x = x + xcng + xs         Return       End Function kiss64       Function refill() Result (s)         Implicit None         Integer (Kind=i8) :: i, s, z, h         Do i = 1, QSIZ           h = iand(carry, 1_i8)           z = ishft(ishft(q(i),41), -1) + ishft(ishft(q(i),39), -1) + &             ishft(carry, -1)           carry = ishft(q(i), -23) + ishft(q(i), -25) + ishft(z, -63)           q(i) = not(ishft(z,1)+h)         End Do         indx = 2         s = q(1)         Return       End Function refill     End Module suprkiss64_m     Program testkiss64       Use suprkiss64_m       Implicit None       Integer (Kind=i8) :: i, x       Do i = 1, qsiz !fill Q with Congruential+Xorshift         xcng = xcng*CMUL + COFFS         xs = ieor(xs, ishft(xs,13))         xs = ieor(xs, ishft(xs,-17))         xs = ieor(xs, ishft(xs,43))         q(i) = xcng + xs       End Do       Do i = 1, 1000000000_i8         x = kiss64()       End Do       Write (*, 100) x 100   Format (' Does x =  4013566000157423768 ?', /, 6X, 'x = ', I20)     End Program testkiss64   Display posts from previous: All Posts1 Day7 Days2 Weeks1 Month3 Months6 Months1 Year Oldest FirstNewest First
 All times are GMT + 1 Hour Page 1 of 1

 Jump to: Select a forum Admin----------------Announcements FTN95----------------GeneralKBaseSupportSuggestionsClearWin+Plato64-bit FTN77----------------Support
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum

Powered by phpBB © 2001, 2005 phpBB Group