
forums.silverfrost.com Welcome to the Silverfrost forums

View previous topic :: View next topic 
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 64bit RNG 


I have recently, Nov 3:
RNGs: A Super KISS, sci.math, comp.lang.c, sci.crypt.
developed a 32bit KISS RNG with immense period,
over 10^402575, great speed and excellent
performance on tests of randomness.
But it is not wellsuited for extension to
64bits or adaption to Fortran.
As 64bit RNGs are becoming more in demand, I
have modified that recently proposed 32bit one
to produce a 64bit 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
(ComplimentaryMultiplyWithCarry).
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 128bit t=a*x+c,
without actually forming t, using only 64bit 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^641) */
#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 64bit integer
13385086898844519155, which I have had printed as
5061657174865032461 for signedinteger implementations.
Can I count on any of you to provide a Fortran version?
George Marsaglia
Here is the 32bit version of the above:

/* SUPRKISS32.c, period 5*2^1320481*(2^321) */
#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);
}
 

Back to top 


mecej4
Joined: 31 Oct 2006 Posts: 828

Posted: Mon Nov 23, 2009 3:31 am Post subject: Fortran90 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 

Back to top 


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
functionsfor 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 64bit random
integers per second.
(The 32bit 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.) 

Back to top 


mecej4
Joined: 31 Oct 2006 Posts: 828

Posted: Mon Nov 23, 2009 7:54 pm Post subject: Errors in first Fortran9x 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 Fortran95 implements arithmetic on 8byte 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 Fortran90 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. 

Back to top 


mecej4
Joined: 31 Oct 2006 Posts: 828

Posted: Mon Nov 23, 2009 8:50 pm Post subject: Another minor bug in C and Fortran versions 


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.
Last edited by mecej4 on Sun Dec 04, 2016 7:08 pm; edited 1 time in total 

Back to top 


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 32bit C listing.
Here I need only ensure that integer arithmetic is the old
standard mod 2^32, with positive residues for signed integers
and leastabsolute residues for signed integers.
Can the default integer*4 still be used, and if not, which
selected_int_kind(?) willor is likely toensure compatibility? 

Back to top 


mecej4
Joined: 31 Oct 2006 Posts: 828

Posted: Tue Nov 24, 2009 8:07 pm Post subject: Fortran version of 32bit SuperKISS RNG 


Just as with the 64bit RNG, the Fortran90 version is straightforward. I have checked it out with both 32bit and 64bit Fortran compilers (both of which, however, have default 4byte 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



Back to top 


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 64bit 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 32bit 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 64bit 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^641)
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^321)
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 

Back to top 


mecej4
Joined: 31 Oct 2006 Posts: 828

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 lineofcode 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 standardconforming. There are, however, useful under adult supervision.
The 32bit RNG, which uses only 4byte integers, and can be used with FTN9532 or FTN9564, follows.
[continued]
Last edited by mecej4 on Mon Dec 05, 2016 4:56 pm; edited 3 times in total 

Back to top 


mecej4
Joined: 31 Oct 2006 Posts: 828

Posted: Sun Dec 04, 2016 8:37 pm Post subject: 


[continuation] Here is the 32bit 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 32bit and 64bit. With Gfortran 4.9 and later, must
! use the ftrapv flag to avoid an optimizer bug.
!
Module suprkiss32_m ! period 5*2^1320481*(2^321)
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 

Back to top 


mecej4
Joined: 31 Oct 2006 Posts: 828

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 8byte integers, and works with FTN9564. It does not work with FTN9532 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 32bit, since it uses the X87 to process 8byte integers
! Works correctly with FTN95 64bit. With Gfortran 4.9 and later, must
! use the ftrapv flag to avoid an optimizer bug.
!
Module suprkiss64_m ! period 5*2^1320480*(2^641)
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 


Back to top 




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
