Gmail Calendar Documents Reader Web more »
Recently Visited Groups | Help | Sign in
Google Groups Home
Message from discussion Fortran and C: United with a KISS
The group you are posting to is a Usenet group. Messages posted to this group will make your email address visible to anyone on the Internet.
Your reply message has not been sent.
Your post was successful
 
From:
To:
Cc:
Followup To:
Add Cc | Add Followup-to | Edit Subject
Subject:
Validation:
For verification purposes please type the characters you see in the picture below or the numbers you hear by clicking the accessibility icon. Listen and type the numbers you hear
 
George Marsaglia  
View profile  
 More options Jun 23 2007, 3:10 pm
Newsgroups: comp.lang.fortran
From: "George Marsaglia" <g...@stat.fsu.edu>
Date: Sat, 23 Jun 2007 15:10:23 -0400
Local: Sat, Jun 23 2007 3:10 pm
Subject: Fortran and C: United with a KISS
I have recently posed questions to this group:
Form t = a*x+c in 64 bits, get top and bottom 32 bits of t
Form k = i + j and test for overflow

Such problems arise in certain kinds of random number
generation---multiply-with-carry and add-with-carry.
http://tbf.coe.wayne.edu:16080/jmasm/vol2_no1.pdf

I was hoping to develop RNGs that would produce
exactly the same sequence of 32-bit random integers
in Fortran as in C, but using only simple code.

It seems that it is not yet feasible to carry out,
in most Fortrans, using 32-bit integers,
the t=a*x+c operation for  multiply-with-carry,
or the k=i+j overflow for add-with-carry.

However, add-with-carry is easy for 31-bit
integers.

I have therefore developed, and post here,
a version of my
   KISS (Keep It Simple, Stupid) RNG,
which uses only
add, shift, exclusive-or and 'and' operations
to produces exactly the same 32-bit integer
output, which C views as unsigned and
Fortran views as signed integers,
each made up of exactly the same pattern
of 32 0's and 1's.

Furthermore, this common output sequence
has period greater that 10**36 and seems
to pass all the  Diehard battery of tests
http://www.cs.hku.hk/~diehard/
and the tuftests in
http://www.jstatsoft.org/v07/i03/tuftests.pdf

/* C version */
static unsigned long t,x=123456789,y=362436069,z=21288629,w=14921776,c=0;
unsigned long KISS(){
x+=545925293;
y^=(y<<13); y^=(y>>17); y^=(y<<5);
t=z+w+c; z=w; c=(t>>31); w=t&2147483647;
return(x+y+w);      }

! Fortran version, (using several statements/line via semicolons.)
     function KISS()
      integer x,y,z,w,c,t
      data x,y,z,w,c/123456789,362436069,21288629,14921776,0/
      m(y,k)=ieor(y,ishft(y,k))
      x=x+545925293
      y=m(m(m(y,13),-17),5)
      t=z+w+c; z=w; c=ishft(t,-31); w=iand(t,2147483647)
      KISS=x+y+w
      return
     entry kisset(ix,iy,iz,iw,ic)
      x=ix; y=iy; z=iz; w=iw; c=ic
      kisset=1
      return
      end

The output of both KISS() functions will be
a sequence of 32-bit integers with period greater
than 2**121 or 10**36:
  period=576384491062058838*4294967296*4294967295
The value returned from KISS() can be any choice
of x op y op w, with op either '+' or 'exclusive-or'.
Each of the four possible choices produced good
results from the Diehard and tuftests, for a
variety of seed sets.

When b=2^31, b^2+b-1 is not prime, but factors into
          610092078393289*7559
so the full period from KISS() requires that
for the seeds x,y,z,w,c:
x can be any 32-bit integer,
y can be any 32-bit integer not 0,
z and w any 31-bit integers not multiples of 7559
c can be 0 or 1.

I was surprised that three such simple
components, individually likely to fail
various tests of randomness, could
collectively, via x op y op w, produce
sequences that do as well in stringent tests
of randomness as do any I have previously
encountered.
That, with their long periods, simple
operations and resulting speed,
(84 million per second on my PC),
makes them worth considering.

Furthermore, by applying that last component,
add-with-carry, to a suitable pair of elements
from an array  of 31-bit seeds z[1],z[2],...,z[n],
rather than just to the current z,w pair,
you can get periods on the order of 2**(31*n+64),
roughly the number of possible ways to initialize
x,y,c and the n 31-bit z's,
and have the same  Fortran-C compatibility,
simplicity and virtually the same speed.

You are welcome to use one or both of the above,
and perhaps try assessing the randomness yourself.
(If you generate 10000 integers with
either the C or Fortran version listed
above, the last four values should be
    Fortran        C
   199275006    199275006
    86473693     86473693
 -2085369775   2209597521
  1298124039   1298124039

George Marsaglia


    Reply to author    Forward  
You must Sign in before you can post messages.
To post a message you must first join this group.
Please update your nickname on the subscription settings page before posting.
You do not have the permission required to post.

Create a group - Google Groups - Google Home - Terms of Service - Privacy Policy
©2009 Google