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 52^1320480(2^64-1) */ #include <stdio.h> static unsigned long long Q[20632],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[0]); }

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 52^1320481(2^32-1) */ #include <stdio.h> static unsigned long Q[41265],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[0]); }

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); }

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.

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,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 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 = 4013566000157423768 ?',/,6x,'x = ',I20)
end program tskiss64

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.)

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=8) to integer(kind=I8).

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.

The array Q is dimensioned [20632] 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.

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?

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).

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

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:

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)

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]

[continuation] Here is the 32-bit RNG program of Marsaglia.

! Posted by the late Professor George Marsaglia to the Silverfrost General forum
! in 2009, see https://forums.silverfrost.com/Forum/Topic/1217 .
! 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]

[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.

! Posted by the late Professor George Marsaglia to the Silverfrost General forum
! in 2009, see https://forums.silverfrost.com/Forum/Topic/1217 .
! 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