Fortran and C: United with a KISS

485 views
Skip to first unread message

George Marsaglia

unread,
Jun 23, 2007, 3:10:23 PM6/23/07
to
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


Wade Ward

unread,
Jun 25, 2007, 3:50:53 PM6/25/07
to
On Jun 23, 3:10 pm, "George Marsaglia" <g...@stat.fsu.edu> wrote:

> /* 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
>

I enjoy posts like this that have code snippets done a couple
different ways, with at least one of them topical in the forum. With
my current capabilities in fortran, I wouldn't be able to pull a
multibyte number apart into bits. I can't wait to bring these home
and fed them to my compilers.
--
Wade Ward

George Marsaglia

unread,
Jun 27, 2007, 6:34:37 AM6/27/07
to
> 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

The above are the last four of a run
of 100000, not 10000. I apologize
for the goof.


Wade Ward

unread,
Jun 28, 2007, 6:31:45 PM6/28/07
to

"George Marsaglia" <g...@stat.fsu.edu> wrote in message
news:oIWdnXMK6MqooB_b...@comcast.com...
I'd like to understand either version of this program, but I'm not sure if
the variables are in the caller.
unsigned long KISS(void)
#include <stdio.h>
int main(void)
{
unsigned long g;
static unsigned long t,x=123456789,y=362436069,z=21288629,w=14921776,c=0;
g=KISS();
printf("%lu\n", g);
return 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); }
This draws about 200 errors from my C compiler. Unlike fortran, my C
compiler doesn't help me fix the errors. I was hoping to compare the two.
--
Wade Ward


mecej4

unread,
Jun 29, 2007, 11:21:09 AM6/29/07
to
Wade Ward wrote:
> "George Marsaglia" <g...@stat.fsu.edu> wrote in message
> news:oIWdnXMK6MqooB_b...@comcast.com...
>>> 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
>> The above are the last four of a run
>> of 100000, not 10000. I apologize
>> for the goof.
> I'd like to understand either version of this program, but I'm not sure if
> the variables are in the caller.

What variables? What caller?

> <<-- Removed incorrect C code -->>

> This draws about 200 errors from my C compiler. Unlike fortran, my C
> compiler doesn't help me fix the errors. I was hoping to compare the two.
> --
> Wade Ward
>
>

Try this:

unsigned long KISS(void);
static unsigned long t, x=123456789, y=362436069, z=21288629,
w=14921776, c=0;

#include <stdio.h>
int main(void){
unsigned long g;

g=KISS();
printf("%lu\n", g);
return 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);
}

-- mecej4

Wade Ward

unread,
Jun 29, 2007, 4:05:56 PM6/29/07
to

"mecej4" <mec...@operamail.com> wrote in message
news:46852369$0$4892$4c36...@roadrunner.com...

> Wade Ward wrote:
>> "George Marsaglia" <g...@stat.fsu.edu> wrote in message
>> news:oIWdnXMK6MqooB_b...@comcast.com...
>>>> 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
>>> The above are the last four of a run
>>> of 100000, not 10000. I apologize
>>> for the goof.
>> I'd like to understand either version of this program, but I'm not sure
>> if the variables are in the caller.
>
> What variables? What caller?
>
>> <<-- Removed incorrect C code -->>
unsigned long KISS(void);
static unsigned long t, x=123456789, y=362436069, z=21288629,
w=14921776, c=0;

#include <stdio.h>
int main(void){
unsigned long g, i;
for (i=1;i<=100000;++i)
{
g=KISS();
if(i > 99996){
printf("%lu %lu\n",i, g);
}
}
// system("pause");
return 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);
}

! end ssource begin output
99997 199275006
99998 86473693
99999 2209597521
100000 1298124039
! end output
Thanks for your reply. I had variables in main that needed to be at
different scope. Now I've got to figure out how to call it in fortran.
--
WW


e p chandler

unread,
Jun 29, 2007, 11:59:04 PM6/29/07
to

do 10 i=1,100000
k=KISS()
if (i>99996) print *,i,k
10 continue
end

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

Note: The 4th line of KISS contains a STATEMENT FUNCTION.

-- elliot


Wade Ward

unread,
Jun 30, 2007, 12:43:37 AM6/30/07
to

"e p chandler" <ep...@juno.com> wrote in message
news:1183175944.2...@k79g2000hse.googlegroups.com...
Oh goody gumdrops. With my anxiety that I have to perform miracles in
strange syntaxes now replaced by the anticipation of doing some
bit-twiddling in fortran, I sign off to look for some dutch girl's extra
time, and, of course, her KISS():-< .
--
WW


Reply all
Reply to author
Forward
0 new messages