KISS4691, a potentially top-ranked RNG.
geo
(Random Number Generator) that ranks
at or near the top in all of the categories:
performance on tests of randomness,
length of period, simplicity and speed.
The most important measure, of course, is
performance on extensive tests of randomness, and for
those that perform well, selection may well depend
on those other measures.
The following KISS version, perhaps call it KISS4691,
seems to rank at the top in all of those categories.
It is my latest, and perhaps my last, as, at age 86,
I am slowing down.
Compiling and running the following commented
C listing should produce, within about four seconds,
10^9 calls to the principal component MWC(), then
10^9 calls to the KISS combination in another ~7 seconds.
Try it; you may like it.
George Marsaglia
/*
The KISS4691 RNG, a Keep-It-Simple-Stupid combination of
MWC (Multiply-With-Carry), Congruential and Xorshift sequences.
Expressed as a simple MWC() function and C macros
#define CNG ( xcng=69069*xcng+123 )
#define XS ( xs^=(xs<<13), xs^=(xs>>17), xs^=(xs<<5) )
#define KISS ( MWC()+CNG+XS )
but easily expressed in other languages, with a slight
complication for signed integers.
With its immense period, >10^45000, and speed: MWC()s at
around 238 million/sec or full KISSes at around 138 million,
there are few RNGs that do as well as this one
on tests of randomness and are comparable in even one
of the categories: period, speed, simplicity---not to
mention comparable in all of them.
The MWC4691 sequence x[n]=8193*x[n-4691]+carry mod b=2^32
is based on p=8193*b^4691-1, period ~ 2^150124 or 10^45192
For the MWC (Multiply-With-Carry) process, given a current
x and c, the new x is (8193*x+c) mod b,
the new c is the integer part of (8193*x+c)/b.
The routine MWC() produces, in reverse order, the base-b=2^32
expansion of some j/p with 0<j<p=8193*2^150112-1 with j
determined by the 64 bits of seeds xcng and xs, or more
generally, by 150112 random bits in the Q[] array.
*/
static unsigned long xs=521288629,xcng=362436069,Q[4691];
unsigned long MWC(void) /*takes about 4.2 nanosecs or 238 million/
second*/
{unsigned long t,x,i; static c=0,j=4691;
j=(j<4690)? j+1:0;
x=Q[j];
t=(x<<13)+c+x; c=(t<x)+(x>>19);
return (Q[j]=t);
}
#define CNG ( xcng=69069*xcng+123 )
#define XS ( xs^=(xs<<13), xs^=(xs>>17), xs^=(xs<<5) )
#define KISS ( MWC()+CNG+XS ) /*138 million/sec*/
#include <stdio.h>
int main()
{unsigned long i,x;
for(i=0;i<4691;i++) Q[i]=CNG+XS;
for(i=0;i<1000000000;i++) x=MWC();
printf(" MWC result=3740121002 ?\n%22u\n",x);
for(i=0;i<1000000000;i++) x=KISS;
printf("KISS result=2224631993 ?\n%22u\n",x);
}
/*
This RNG uses two seeds to fill the Q[4691] array by
means of CNG+XS mod 2^32. Users requiring more than two
seeds will need to randomly seed Q[] in main().
By itself, the MWC() routine passes all tests in the
Diehard Battery of Tests, but it is probably a good
idea to include it in the KISS combination.
Languages requiring signed integers should use equivalents
of the same operations, except that the C statement:
c=(t<x)+(x>>19);
can be replaced by that language's version of
if sign(x<<13+c)=sign(x) then c=sign(x)+(x>>19)
else c=1-sign(t)+(x>>19)
*/
Geoff
wrote:
The compiler also notices that i is defined but unused in MWC().
Purely cosmetic changes for readability and style:
/*
* The KISS4691 RNG, a Keep-It-Simple-Stupid combination of
* MWC (Multiply-With-Carry), Congruential and Xorshift sequences.
* Expressed as a simple MWC() function and C macros
* #define CNG ( xcng=69069*xcng+123 )
* #define XS ( xs^=(xs<<13), xs^=(xs>>17), xs^=(xs<<5) )
* #define KISS ( MWC()+CNG+XS )
* but easily expressed in other languages, with a slight
* complication for signed integers.
*
* With its immense period, >10^45000, and speed: MWC()s at
* around 238 million/sec or full KISSes at around 138 million,
* there are few RNGs that do as well as this one
* on tests of randomness and are comparable in even one
* of the categories: period, speed, simplicity---not to
* mention comparable in all of them.
*
* The MWC4691 sequence x[n]=8193*x[n-4691]+carry mod b=2^32
* is based on p=8193*b^4691-1, period ~ 2^150124 or 10^45192
* For the MWC (Multiply-With-Carry) process, given a current
* x and c, the new x is (8193*x+c) mod b,
* the new c is the integer part of (8193*x+c)/b.
*
* The routine MWC() produces, in reverse order, the base-b=2^32
* expansion of some j/p with 0<j<p=8193*2^150112-1 with j
* determined by the 64 bits of seeds xcng and xs, or more
* generally, by 150112 random bits in the Q[] array.
* George Marsaglia
*/
static unsigned long xs = 521288629;
static unsigned long xcng = 362436069;
static unsigned long Q[4691];
/* takes about 4.2 nanosecs or 238 million/second */
unsigned long MWC(void)
{
unsigned long t, x, i;
static c = 0, j = 4691;
j = (j < 4690) ? j + 1 : 0;
x = Q[j];
t = (x<<13) + c + x;
c = (t < x) + (x>>19);
return (Q[j] = t);
}
#define CNG ( xcng = 69069 * xcng + 123 )
#define XS ( xs^=(xs<<13), xs^=(xs>>17), xs^=(xs<<5) )
#define KISS ( MWC() + CNG + XS ) /*138 million/sec*/
#include <stdio.h>
int main()
{
unsigned long i, x;
for(i = 0; i < 4691; i++)
Q[i] = CNG + XS;
for(i = 0; i < 1000000000; i++)
x = MWC();
printf(" MWC result=3740121002 ?\n%22u\n",x);
for(i = 0; i < 1000000000; i++)
x = KISS;
printf("KISS result=2224631993 ?\n%22u\n",x);
return 0;
}
/*
* This RNG uses two seeds to fill the Q[4691] array by
* means of CNG+XS mod 2^32. Users requiring more than two
* seeds will need to randomly seed Q[] in main().
* By itself, the MWC() routine passes all tests in the
* Diehard Battery of Tests, but it is probably a good
* idea to include it in the KISS combination.
*
* Languages requiring signed integers should use equivalents
* of the same operations, except that the C statement:
* c=(t<x)+(x>>19);
* can be replaced by that language's version of
* if sign(x<<13+c)=sign(x) then c=sign(x)+(x>>19)
* else c=1-sign(t)+(x>>19)
*/
jacob navia
This doesn't work with systems that have unsigned long as a 64 bit quantity.
I obtain:
MWC result=3740121002 ?
4169348530
KISS result=2224631993 ?
1421918629
Compiling with 32 bit machine yields:
MWC result=3740121002 ?
3740121002
KISS result=2224631993 ?
2224631993
Gene
Here's an implementation in Ada, verified to produce the same answers
as the C code.
The Ada version of gcc compiles this to nearly the same code as
produced
for the C. Function call overhead adds about 40% to the run time.
What's a good way to take an arbitrary, single 32-bit seed value to
a complete state initialization? Accept that lots of people will
pick very small numbers.
-- kiss_random_numbers.ads
-- Specification for George Marsaglia's KISS random number generator.
package KISS_Random_Numbers is
type Result_Type is mod 2 ** 32;
type State_Type is private;
function Next_Random_Value(State : access State_Type) return
Result_Type;
private
type State_Index_Type is mod 4691;
type State_Vector_Type is array (State_Index_Type) of Result_Type;
Initial_Q : State_Vector_Type; -- set in package body
type Substate_Type is
record
Xs : Result_Type := 521288629;
Xcng : Result_Type := 362436069;
end record;
type State_Type is
record
Sub : aliased Substate_Type;
C : Result_Type := 0;
Q : State_Vector_Type := Initial_Q;
J : State_Index_Type := State_Index_Type'Last;
end record;
end KISS_Random_Numbers;
-- kiss_random_numbers.adb
-- Implementation of George Marsaglia's KISS random number generator.
package body KISS_Random_Numbers is
function Mwc (State : access State_Type) return Result_Type is
T, X : Result_Type;
begin
State.J := State.J + 1;
X := State.Q(State.J);
T := X * 2**13 + State.C + X;
if T < X then
State.C := X / 2**19 + 1;
else
State.C := X / 2**19;
end if;
State.Q(State.J) := T;
return T;
end Mwc;
pragma Inline(Mwc);
function Cng (State : access Substate_Type) return Result_Type is
begin
State.Xcng := 69069 * State.Xcng + 123;
return State.Xcng;
end Cng;
pragma Inline(Cng);
function Xs (State : access Substate_Type) return Result_Type is
Xs : Result_Type;
begin
Xs := State.Xs;
Xs := Xs xor (Xs * 2**13);
Xs := Xs xor (Xs / 2**17);
Xs := Xs xor (Xs * 2**5);
State.Xs := Xs;
return Xs;
end Xs;
pragma Inline(Xs);
function Next_Random_Value(State : access State_Type) return
Result_Type is
begin
return Mwc(State) + Cng(State.Sub'Access) +
Xs(State.Sub'Access);
end Next_Random_Value;
begin
declare
S : aliased Substate_Type;
begin
for I in Initial_Q'Range loop
Initial_Q(I) := Cng(S'access) + Xs(S'Access);
end loop;
end;
end KISS_Random_Numbers;
-- test_kiss.adb
-- Simple test of George Marsaglia's KISS random number generator.
with Ada.Text_IO; use Ada.Text_IO;
with KISS_Random_Numbers;
use KISS_Random_Numbers;
procedure Test_Kiss is
X : Result_Type;
S : aliased State_Type;
begin
for I in 1 .. 1000000000 loop
X := Next_Random_Value(S'Access);
end loop;
Put(Result_Type'Image(X));
New_Line;
end;
Gib Bogle
>
> Languages requiring signed integers should use equivalents
> of the same operations, except that the C statement:
> c=(t<x)+(x>>19);
> can be replaced by that language's version of
> if sign(x<<13+c)=sign(x) then c=sign(x)+(x>>19)
> else c=1-sign(t)+(x>>19)
> */
Hi George,
I translated this into Fortran, and found that I get different results than with
C. I've tracked the difference into MWC(). The following Fortran code, with my
laborious comparison of two signed integers treating them as unsigned, gives
correct results. If I comment out the line
c = tLTx + SHIFTR(x,19)
and uncomment the following lines that implement your suggestion above to
compute c, I get different results.
integer function MWC()
integer :: t, x, i
integer, save :: c = 0, j = 4691
integer :: tLTx
if (j < 4690) then
j = j + 1
else
j = 0
endif
x = Q(j)
t = SHIFTL(x,13) + c + x
if ((t >= 0 .and. x >= 0) .or. (t < 0 .and. x < 0)) then
if (t < x) then
tLTx = 1
else
tLTx = 0
endif
elseif (x < 0) then
tLTx = 1
elseif (t < 0) then
tLTx = 0
endif
c = tLTx + SHIFTR(x,19)
!if (sign(1,SHIFTL(x,13)+c) == sign(1,x)) then
! c = sign(1,x) + SHIFTR(x,19)
!else
! c = 1 - sign(1,t) + SHIFTR(x,19)
!endif
Q(j) = t
MWC = t
end function
Is it the case that although your suggested workaround gives different results
from the C code in this case, it is still equivalent as a RNG?
Cheers
Gib
geo
Thanks very much for the Fortran version.
I made a mistake in the comment on versions
for signed integers. This:
Languages requiring signed integers should use equivalents
of the same operations, except that the C statement:
c=(t<x)+(x>>19);
can be replaced by that language's version of
if sign(x<<13+c)=sign(x) then c=sign(x)+(x>>19)
else c=1-sign(t)+(x>>19)
should have been:
Languages requiring signed integers should use equivalents
of the same operations, except that the C statement:
c=(t<x)+(x>>19);
can be replaced by that language's version of
if sign(x<<13+c)=sign(x) then c=sign(x)+(x>>19)
else c=1-sign(x<<13+c)+(x>>19)
Sorry for that error.
I still like inline functions in Fortan,
so would tend to define
isign(x)=ishft(x,-31)
and
m=ishft(x,13)+c
if(isign(m).eq.isign(x)) then c=isign(x)+ishft(x,-19)
else c=1-isign(m)+ishft(x,-19)
and
Q[j]=m+x
If calculating the carry c of the MWC operation fails
to fix that extra increment properly, then rather than
a systematic expansion, in reverse order, 32 bits at a time,
of the binary representation of j/(1893*2^150112-1) for some
j determined by the random seeds, we will be jumping around
in that expansion, and we can't be sure that the period will
still be the order of b=2^32 for the prime p=1893*b^4196-1.
gm
Dick Hendrickson
[snip]
>
> I still like inline functions in Fortan,
> so would tend to define
> isign(x)=ishft(x,-31)
function in Fortran and compilers should generate near perfect code
inline for things like sign(1,x). The generic sign function was added
to Fortran 90. In FORTRAN 77, you'd need to use the specific ISIGN(1,X)
function.
Dick Hendrickson
[snip]
>
> gm
Gib Bogle
> Thanks very much for the Fortran version.
> I made a mistake in the comment on versions
> for signed integers. This:
>
> Languages requiring signed integers should use equivalents
> of the same operations, except that the C statement:
> c=(t<x)+(x>>19);
> can be replaced by that language's version of
> if sign(x<<13+c)=sign(x) then c=sign(x)+(x>>19)
> else c=1-sign(t)+(x>>19)
>
> should have been:
>
>
> Languages requiring signed integers should use equivalents
> of the same operations, except that the C statement:
> c=(t<x)+(x>>19);
> can be replaced by that language's version of
> if sign(x<<13+c)=sign(x) then c=sign(x)+(x>>19)
> else c=1-sign(x<<13+c)+(x>>19)
>
> Sorry for that error.
That produces different c values from those generated by the method based on the
value of (t<x), therefore it deviates from the C code. This is what I used:
m = shiftl(x,13) + c
if (sign(1,m) == sign(1,x)) then
c = sign(1,x) + shiftr(x,19)
else
c = 1 - sign(1,m) + shiftr(x,19)
endif
robin
|I have been asked to recommend an RNG
| (Random Number Generator) that ranks
| at or near the top in all of the categories:
| performance on tests of randomness,
| length of period, simplicity and speed.
| The most important measure, of course, is
| performance on extensive tests of randomness, and for
| those that perform well, selection may well depend
| on those other measures.
|
| The following KISS version, perhaps call it KISS4691,
| seems to rank at the top in all of those categories.
| It is my latest, and perhaps my last, as, at age 86,
| I am slowing down.
|
| Compiling and running the following commented
| C listing should produce, within about four seconds,
| 10^9 calls to the principal component MWC(), then
| 10^9 calls to the KISS combination in another ~7 seconds.
|
| Try it; you may like it.
|
| George Marsaglia
Here's the PL/I version:
(NOSIZE, NOFOFL):
RNG: PROCEDURE OPTIONS (MAIN, REORDER);
declare (xs initial (521288629), xcng initial (362436069),
Q(0:4690) ) static fixed binary (32) unsigned;
MWC: procedure () returns (fixed binary (32) unsigned);
/*takes about 4.2 nanosecs or 238 million/second*/
declare (t,x,i) fixed binary (32) unsigned;
declare (c initial (0), j initial (4691) ) fixed binary (32) unsigned static;
if j < 4690 then j = j + 1; else j = 0;
x = Q(j);
t = isll(x,13)+c+x; c = (t<x)+isrl(x,19);
Q(j)=t;
return (t);
end MWC;
CNG: procedure returns (fixed binary (32) unsigned);
xcng=bin(69069)*xcng+bin(123);
return (xcng);
end CNG;
XXS: procedure returns (fixed binary (32) unsigned);
xs = ieor (xs, isll(xs, 13) );
xs = ieor (xs, isrl(xs, 17) );
xs = ieor (xs, isll(xs, 5) );
return (xs);
end XXS;
KISS: procedure returns (fixed binary (32) unsigned);
return ( MWC()+CNG+XXS ); /*138 million/sec*/
end KISS;
declare (i,x) fixed binary (32) unsigned;
/* Initialize: */
do i = 0 to 4691-1; Q(i) = CNG+XXS; end;
put skip list (q(0), q(4690));
put skip list ('initialized'); put skip;
do i = 0 to 1000000000-1; x=MWC(); end;
put skip edit (" MWC result=3740121002 ",x) (a, f(23));
do i = 0 to 1000000000-1; x=KISS; end;
put skip edit ("KISS result=2224631993 ",x) (a, f(23));
end RNG;
robin
| geo wrote:
| > Thanks very much for the Fortran version.
| > I made a mistake in the comment on versions
| > for signed integers. This:
| >
| > Languages requiring signed integers should use equivalents
| > of the same operations, except that the C statement:
| > c=(t<x)+(x>>19);
| > can be replaced by that language's version of
| > if sign(x<<13+c)=sign(x) then c=sign(x)+(x>>19)
| >
| > Sorry for that error.
|
| That produces different c values from those generated by the method based on the
| value of (t<x), therefore it deviates from the C code. This is what I used:
|
| m = shiftl(x,13) + c
| if (sign(1,m) == sign(1,x)) then
| c = sign(1,x) + shiftr(x,19)
| else
| c = 1 - sign(1,m) + shiftr(x,19)
| endif
Maybe I missed something,
but isn't this exactly equivalent to what George wrote?
Just substitute x<<13+c for m in your last two assignments ...
Geoff
wrote:
I am no FORTRAN hacker but I think there's a difference between
sign(x) and sign(1,x).
Harald Anlauf
> This RNG uses two seeds to fill the Q[4691] array by
> means of CNG+XS mod 2^32. Users requiring more than two
> seeds will need to randomly seed Q[] in main().
> By itself, the MWC() routine passes all tests in the
> Diehard Battery of Tests, but it is probably a good
> idea to include it in the KISS combination.
Does this mean that using different seeds will lead to
streams that are always statistically independent
(as long as one does not exhaust the RNG's period)?
Or are there restrictions on the possible combinations
of seeds?
I am currently using D.E.Knuth's generator from TAOCP,
which IIRC allows for 2^30-2 independent streams, and
which are asserted to be independent, but being able to
extend the "limit" would be nice.
Harald
Gib Bogle
I hope so. Maybe I didn't express myself clearly enough. I'll try again.
Using my implementation of George's corrected code, I get results from the
Fortran code that differ from those generated by his C code.
Gib Bogle
> I am no FORTRAN hacker but I think there's a difference between
> sign(x) and sign(1,x).
sign(1,x) returns the sign of x.
nm...@cam.ac.uk
Harald Anlauf <anlau...@arcor.de> wrote:
>On Jul 24, 3:02=A0pm, geo <gmarsag...@gmail.com> wrote:
>> This RNG uses two seeds to fill the Q[4691] array by
>> means of CNG+XS mod 2^32. Users requiring more than two
>> seeds will need to randomly seed Q[] in main().
>> By itself, the MWC() routine passes all tests in the
>> Diehard Battery of Tests, but it is probably a good
>> idea to include it in the KISS combination.
>
>Does this mean that using different seeds will lead to
>streams that are always statistically independent
>(as long as one does not exhaust the RNG's period)?
>Or are there restrictions on the possible combinations
>of seeds?
No. Without checking it more carefully, I can't say definitely,
but it looks as if you would be very unlikely to notice a PRACTICAL
association with two different seeds, provided that you throw
away the first 10,000 or so numbers - and even that qualification
may be unnecessary. But, unless I have some missed some subtlety,
the sequences cannot be guaranteed to be pseudo-independent.
The only two methods I know of of guaranteeing pseudo-independence
are using coprime sequences and by choosing them using the spectral
test or equivalent. Even then, there are some qualifications on
what is meant by pseudo-independence. However, in practice, it's
rare to have trouble with high-quality generators.
>I am currently using D.E.Knuth's generator from TAOCP,
>which IIRC allows for 2^30-2 independent streams, and
>which are asserted to be independent, but being able to
>extend the "limit" would be nice.
Eh? Which version? There was assuredly nothing with those
properties in edition 2. The first papers on the topic were later.
You need to be very careful to distinguish separate (i.e. disjoint)
sequences from pseudo-independent ones, and FAR too many papers
written by people who ought to know better confuse the two. Doing
that is a common cause of seriously wrong answers in some types of
calculation.
Regards,
Nick Maclaren.
Harald Anlauf
> In article <fa9dd141-823a-4179-a80b-18d214a3e...@t2g2000yqe.googlegroups.com>,
> Harald Anlauf <anlauf.2...@arcor.de> wrote:
>
> >Does this mean that using different seeds will lead to
> >streams that are always statistically independent
> >(as long as one does not exhaust the RNG's period)?
> >Or are there restrictions on the possible combinations
> >of seeds?
>
> No. Without checking it more carefully, I can't say definitely,
> but it looks as if you would be very unlikely to notice a PRACTICAL
> association with two different seeds, provided that you throw
> away the first 10,000 or so numbers - and even that qualification
> may be unnecessary. But, unless I have some missed some subtlety,
> the sequences cannot be guaranteed to be pseudo-independent.
>
> The only two methods I know of of guaranteeing pseudo-independence
> are using coprime sequences and by choosing them using the spectral
> test or equivalent. Even then, there are some qualifications on
> what is meant by pseudo-independence. However, in practice, it's
> rare to have trouble with high-quality generators.
I was shown funny experiences with simple generators... =:-o
> >I am currently using D.E.Knuth's generator from TAOCP,
> >which IIRC allows for 2^30-2 independent streams, and
> >which are asserted to be independent, but being able to
> >extend the "limit" would be nice.
>
> Eh? Which version? There was assuredly nothing with those
> properties in edition 2. The first papers on the topic were later.
http://www-cs-faculty.stanford.edu/~uno/news02.html#rng
> You need to be very careful to distinguish separate (i.e. disjoint)
> sequences from pseudo-independent ones, and FAR too many papers
> written by people who ought to know better confuse the two. Doing
> that is a common cause of seriously wrong answers in some types of
> calculation.
For an application which needs to distribute a large problem over many
processors it is very useful to have many (pseudo-)independent
streams,
at least they should be practically indistinguishable from independent
ones. An accidental correlation is likely worse than a short period.
Harald
glen herrmannsfeldt
> In article <fa9dd141-823a-4179...@t2g2000yqe.googlegroups.com>,
> Harald Anlauf <anlau...@arcor.de> wrote:
>>Does this mean that using different seeds will lead to
>>streams that are always statistically independent
>>(as long as one does not exhaust the RNG's period)?
>>Or are there restrictions on the possible combinations
>>of seeds?
> No. Without checking it more carefully, I can't say definitely,
> but it looks as if you would be very unlikely to notice a PRACTICAL
> association with two different seeds, provided that you throw
> away the first 10,000 or so numbers - and even that qualification
> may be unnecessary. But, unless I have some missed some subtlety,
> the sequences cannot be guaranteed to be pseudo-independent.
My biggest complaint about the current standard RANDOM_SEED
is that it doens't provide a way to get a reliably good
seed from a default (likely 32 bit) integer.
There are many generators with extrememly long periods,
and correspondingly long state. As the designers of the RNG
are the ones likely to know how to choose a good seed, it
would seem they would be the best ones to supply a good
seed generator.
> The only two methods I know of of guaranteeing pseudo-independence
> are using coprime sequences and by choosing them using the spectral
> test or equivalent. Even then, there are some qualifications on
> what is meant by pseudo-independence. However, in practice, it's
> rare to have trouble with high-quality generators.
Now, one can supply an array of the appropriate length to
RANDOM_SEED(PUT=...), but how to generate such an array
from a smaller seed? There is no way to know.
(snip)
> You need to be very careful to distinguish separate (i.e. disjoint)
> sequences from pseudo-independent ones, and FAR too many papers
> written by people who ought to know better confuse the two. Doing
> that is a common cause of seriously wrong answers in some types of
> calculation.
-- glen
nm...@cam.ac.uk
Harald Anlauf <anlau...@arcor.de> wrote:
>>
>> The only two methods I know of of guaranteeing pseudo-independence
>> are using coprime sequences and by choosing them using the spectral
>> what is meant by pseudo-independence. =A0However, in practice, it's
>> rare to have trouble with high-quality generators.
>
Yup. Very common.
>> >I am currently using D.E.Knuth's generator from TAOCP,
>> >which IIRC allows for 2^30-2 independent streams, and
>> >which are asserted to be independent, but being able to
>> >extend the "limit" would be nice.
>>
>> Eh? =A0Which version? =A0There was assuredly nothing with those
>> properties in edition 2. =A0The first papers on the topic were later.
>
>http://www-cs-faculty.stanford.edu/~uno/news02.html#rng
Ah. Thanks. It's new, then. I haven't been active in this area
ina couple of decades, so haven't been keeping track.
>> You need to be very careful to distinguish separate (i.e. disjoint)
>> sequences from pseudo-independent ones, and FAR too many papers
>> written by people who ought to know better confuse the two. =A0Doing
>> that is a common cause of seriously wrong answers in some types of
>> calculation.
>
>For an application which needs to distribute a large problem over many
>processors it is very useful to have many (pseudo-)independent
>streams,
>at least they should be practically indistinguishable from independent
>ones. An accidental correlation is likely worse than a short period.
Right. That's precisely why I got involved.
Regards,
Nick Maclaren.
robin
George gave general advice on how to do it.
That advice wasn't specific to Fortran.
It's necessary to use the appropriate Fortran construct --
and sign(1,x) is the only way to do that.
Geoff
wrote:
Thanks for that clarification. I have not done any Fortran code since
1975 and it looked a whole lot different than it does today.
ss3329
’s shoe is the most famous product [url=http://www.b2bsharing.com/footwear-nike-shox-shoes-c-3502_3734.html], etc. [b]Nike
Shox running shoes[/b][/url] produced material is comfortable and airy. [url=http://www.b2bsharing.com/nike-shox-shoes-
mens-shox-r3-shoes-c-3502_3734_3739.html][b]Shox r3[/b][/url] has managed to invite the awareness of Disney, For demand. As
a part of the [url=http://www.b2bsharing.com/footwear-nike-shox-shoes-c-3502_3734.html][b]Nike Shox[/b][/url] family, this
stylish shoe is sure to be a hit amongst athletes and others looking for the latest in style and fashion. There are more and
more people to buy [url=http://www.b2bsharing.com/footwear-nike-shox-shoes-c-3502_3734.html][b]wholesale nike shox [/b][/url]
Nike shox r3 comes equipped with extra padding built into the heel for athletic purposes as well the [b]
[url=http://www.b2bsharing.com/nike-shox-shoes-womens-shox-turb-shoes-c-3502_3734_3750.html]Women's Shox Turb Shoes[/url]
[/b]is designed streamlined for style, comfort, as well as for practical reasons as arunning shoe or trainer.
nm...@cam.ac.uk
Geoff <ge...@invalid.invalid> wrote:
>>|
>>| I am no FORTRAN hacker but I think there's a difference between
>>| sign(x) and sign(1,x).
>>
>>George gave general advice on how to do it.
>>That advice wasn't specific to Fortran.
>>It's necessary to use the appropriate Fortran construct --
>>and sign(1,x) is the only way to do that.
>>
>Thanks for that clarification. I have not done any Fortran code since
>1975 and it looked a whole lot different than it does today.
Actually, that aspect hasn't changed since :-) SIGN(X) always was
an implementation-dependent extension (now called processor-dependent).
Regards,
Nick Maclaren.
robin
| (Random Number Generator) that ranks
| at or near the top in all of the categories:
| performance on tests of randomness,
| length of period, simplicity and speed.
| The most important measure, of course, is
| performance on extensive tests of randomness, and for
| those that perform well, selection may well depend
| on those other measures.
I have already posted a PL/I version using unsigned arithmetic.
Here is another version, this time using signed arithmetic :--
(NOSIZE, NOFOFL):
RNG: PROCEDURE OPTIONS (MAIN, REORDER);
declare (xs initial (521288629), xcng initial (362436069),
Q(0:4690) ) static fixed binary (31);
MWC: procedure () returns (fixed binary (31));
declare (t,x,i) fixed binary (31);
declare (c initial (0), j initial (4691) ) fixed binary (31) static;
declare (t1, t2, t3) fixed binary (31);
if j < hbound(Q,1) then j = j + 1; else j = 0;
x = Q(j);
t = isll(x,13)+c+x;
t1 = iand(x, 3) - iand(t, 3);
t2 = isrl(x, 2) - isrl(t, 2);
if t2 = 0 then t2 = t1;
if t2 > 0 then t3 = 1; else t3 = 0;
c = t3 + isrl(x, 19);
Q(j)=t;
return (t);
end MWC;
CNG: procedure returns (fixed binary (31));
xcng=bin(69069)*xcng+bin(123);
return (xcng);
end CNG;
XXS: procedure returns (fixed binary (31));
xs = ieor (xs, isll(xs, 13) );
xs = ieor (xs, isrl(xs, 17) );
xs = ieor (xs, isll(xs, 5) );
return (xs);
end XXS;
KISS: procedure returns (fixed binary (31));
return ( MWC()+CNG+XXS );
end KISS;
declare (i,x) fixed binary (31);
declare y fixed decimal (11);
Q = CNG+XXS; /* Initialize. */
do i = 1 to 1000000000; x=MWC(); end;
put skip edit (" Expected MWC result = 3740121002", 'computed =', x)
(a, skip, x(12), a, f(11));
y = iand(x, 2147483647);
if x < 0 then y = y + 2147483648;
put skip edit (y) (x(11), f(22)); put skip;
do i = 1 to 1000000000; x=KISS; end;
put skip edit ("Expected KISS result = 2224631993", 'computed =', x)
(a, skip, x(12), a, f(11));
y = iand(x, 2147483647);
if x < 0 then y = y + 2147483648;
put skip edit (y) (x(11), f(22));
end RNG;
robin
| This doesn't work with systems that have unsigned long as a 64 bit quantity.
|
| I obtain:
|
| MWC result=3740121002 ?
| 4169348530
| KISS result=2224631993 ?
| 1421918629
For a 64-bit machine (using 64-bit integer arithmetic),
you'd need to truncate each result to 32 bits. That not
only applies to the multiplication, it also applies to addition, etc.
On a 32-bit machine, these extra bits are discarded,
but in 64-bit arithmetic, they are retained,
and unless they are similarly discarded,
you won't get the same results.
I suggest using IAND(k, 2*2147483647+1)
for the truncation.
With such modifications in the program,
it should then produce the same results on both 32-bit and
64-bit machines.
P.S. the product 2*2... is best obtained using ISHFT.
Uno
So for the put= values in fortran, you need a vector of pseudorandom
integers, which is as good as it gets without truly random devices,
making--one hopes-a period that is large with respect to the interval
you're interested in.
It doesn't seem like a problem with epistemology as much a mathematical
ceiling on how much randomness you can create by a handful of values.
--
Uno
glen herrmannsfeldt
(snip, I wrote)
>> My biggest complaint about the current standard RANDOM_SEED
>> is that it doens't provide a way to get a reliably good
>> seed from a default (likely 32 bit) integer.
>> There are many generators with extrememly long periods,
>> and correspondingly long state. As the designers of the RNG
>> are the ones likely to know how to choose a good seed, it
>> would seem they would be the best ones to supply a good
>> seed generator.
> So for the put= values in fortran, you need a vector of pseudorandom
> integers, which is as good as it gets without truly random devices,
> making--one hopes-a period that is large with respect to the interval
> you're interested in.
In a large fraction of the cases, 2 billion different seeds
are enough, but one can still desire the appropriate randomness
from those different seeds.
> It doesn't seem like a problem with epistemology as much
> a mathematical ceiling on how much randomness you can create
> by a handful of values.
Given a default integer, one might fill an array with that
integer and use that as a seed. That might be good for some,
not so good for others. Even more, for standard Fortran such
should be done without knowing the range of values of an integer
variable.
R has two seeding functions, one that takes a full length state
array, and the other takes a single integer. That makes sense
to me.
-- glen
Uno
Gib, can you post your entire fortran version?
--
Uno
Gib Bogle
> Gib, can you post your entire fortran version?
OK, here's the Fortran. I've removed almost all the comments to save space.
If you uncomment the line in MWC():
c = c_this_works
you'll see how to reproduce the C code results.
Gib
module kiss_mod
implicit none
integer :: xs = 521288629
integer :: xcng = 362436069
integer :: Q(0:4690)
contains
integer function MWC()
integer :: m, x, i
integer, save :: c = 0, j = 4691
integer :: t, tLTx, c_this_works
if (j < 4690) then
j = j + 1
else
j = 0
endif
x = Q(j)
m = SHIFTL(x,13) + c
Q(j) = m + x
t = Q(j)
! This is the laborious determination of c
if ((t >= 0 .and. x >= 0) .or. (t < 0 .and. x < 0)) then
if (t < x) then
tLTx = 1
else
tLTx = 0
endif
elseif (x < 0) then
tLTx = 1
elseif (t < 0) then
tLTx = 0
endif
c_this_works = tLTx + SHIFTR(x,19)
! This c gives results inconsistent with the C code
if (sign(1,m) == sign(1,x)) then
c = sign(1,x) + SHIFTR(x,19)
else
c = 1 - sign(1,m) + SHIFTR(x,19)
endif
! c = c_this_works
MWC = Q(j)
end function
subroutine test_RNG4691
integer :: i, x, unsigned_result
integer :: n = 1000000000
do i = 0, 4690
xcng = 69069 * xcng + 123
xs = XOR(xs,SHIFTL(xs,13))
xs = XOR(xs,SHIFTR(xs,17))
xs = XOR(xs,SHIFTL(xs,5))
Q(i) = xcng + xs
enddo
do i = 1, n
x = MWC()
enddo
unsigned_result = 3740121002
write(*,*) "MWC result=3740121002 ?: ",unsigned_result,x
do i = 1,n
xcng = 69069 * xcng + 123
xs = XOR(xs,SHIFTL(xs,13))
xs = XOR(xs,SHIFTR(xs,17))
xs = XOR(xs,SHIFTL(xs,5))
x = MWC() + xcng + xs
enddo
unsigned_result = 2224631993
write(*,*) "KISS result=2224631993 ?: ",unsigned_result,x
end subroutine
end module
program main
use kiss_mod
call test_RNG4691
end program
Uno
> OK, here's the Fortran. I've removed almost all the comments to save
> space.
> If you uncomment the line in MWC():
> c = c_this_works
> you'll see how to reproduce the C code results.
Alright, thx, gib, looks like there's a couple rough corners:
$ gfortran geo1.f90 -o out
geo1.f90:63.32:
unsigned_result = 3740121002
1
Error: Integer too big for its kind at (1). This check can be disabled
with the option -fno-range-check
So, do I ask for something bigger than the default integer, use the
iso_c_binding, or disable the warning?
I think of 3.5 billion as a value that busts types in C. I would have
thought that was going on except for the second error:
geo1.f90:72.32:
unsigned_result = 2224631993
1
Error: Integer too big for its kind at (1). This check can be disabled
with the option -fno-range-check
Then there's this. It looks like the bitshifting capabilities that
fortran and C have, but neither MR&C nor my slightly-inappropriate C
reference have a SHIFTL. Do you have a definition for it?
geo1.f90:55.16:
xs = XOR(xs,SHIFTL(xs,13))
1
Error: Function 'shiftl' at (1) has no IMPLICIT type
Beautiful night with the cicadas humming. I don't know that I've ever
seen one which makes them, aesthetically, my favorite bug. Cheers,
--
Uno
Ron Shepard
glen herrmannsfeldt <g...@ugcs.caltech.edu> wrote:
> In a large fraction of the cases, 2 billion different seeds
> are enough, but one can still desire the appropriate randomness
> from those different seeds.
My understanding is that pseudorandom number generators return a
sequence of values with some period. The numbers in that sequence
eventually repeat with some, hopefully long, cycle. The put=array
argument gives the starting point within that sequence, but it doesn't
affect the cycle length or the "randomness" of the values, those things
are determined by the underlying algorithm, not by the initial seed.
The situation that needs to be avoided is to run the code once with one
seed, and then run it again with another seed that results in an overlap
of the sequences of values for the two runs. In some applications this
is unimportant, but in other applications it would be bad for two
supposedly independent runs to have a long sequence of identical values.
It would be nice if there were some prescription to generate initial
seeds that would avoid this situation. However, I don't really know how
this could be specified in the fortran standard without specifying the
algorithm itself (which is not done, apparently on purpose). Anyone
have any ideas how this could be done?
$.02 -Ron Shepard
glen herrmannsfeldt
Here is the suggestion. Add to RANDOM_SEED a form that allows
a single integer variable to be supplied as a seed source.
As far as I know, there is no requirement on period or otherwise
on the quality of the PRNG used, but the standard could suggest
that good seeds be selected from that integer. For example,
they could be chosen such that they have a fairly long sequence
before they overlap the sequence generated by another input
value. If the RNG only has 32 bits of state, or maybe less, then
it likely takes it directly.
I don't know the math of the newer PRNGs enough to say, but I
think it is possible to do that.
With the current standard there is no way to even suggest to
RANDOM_SEED that you want a good seed.
OK, here is another suggestion that could be implemented without
any change in the standard. In the case where the supplied
seed array has all except the first element zero, then select
a good seed based on that first element.
-- glen
orz
wrote:
> In article <i2ssst$nf...@speranza.aioe.org>,
On Jul 29, 8:13 pm, Ron Shepard <ron-shep...@NOSPAM.comcast.net>
wrote:
> In article <i2ssst$nf...@speranza.aioe.org>,
If an RNG has a large number of possible states (say, 2^192) then it's
unlikely for any two runs to EVER reach identical states by any means
other than starting from the same seed. This RNG has a period vastly
in excess of "large", so for practical purposes it just won't
happen.
However, because this is a combination of multiple component RNGs
there is the related issue of overlapping within the period of one or
more component RNGs. MWC has essentially zero chance of being in an
identical state for any reason other than identical seeding, but both
XS and CNG will frequently have identical states in two different
runs. There does not seem to be any detectable correlation however
unless (either MWC or (both XS and CNG)) have identical states between
two different runs. Which is rather uncommon, though it does happen
often enough to see it in the real world on occasion.
Gib Bogle
Uno, the only reason for putting those lines in with unsigned_integer is to show
that the result that you see when x is written is actually the same as the
C-generated number. The Intel Fortran compiler is happy with the code, not even
issuing a warning. I suggest using -fno-range-check, or alternatively subtract
2^31 from the big numbers (I think). Or you could just take my word for it.
I didn't realize that SHIFTL and SHIFTR were not standard Fortran. You can use
ISHFT:
result = ISHFT (i,shift)
i (Input) Must be of type integer.
shift (Input) Must be of type integer. The absolute value for shift must be less
than or equal to BIT_SIZE( i).
Results
The result type is the same as i. The result has the value obtained by shifting
the bits of i by shift positions. If shift is positive, the shift is to the
left; if shift is negative, the shift is to the right. If shift is zero, no
shift is performed.
Bits shifted out from the left or from the right, as appropriate, are lost.
Zeros are shifted in from the opposite end.
ISHFT with a positive shift can also be specified as LSHIFT (or LSHFT). ISHFT
with a negative shift can also be specified as RSHIFT (or RSHFT)with | shift |.
I'm very surprised to hear you say you've never seen a cicada. In the summer
they are impossible to miss here. Sparrows catch them on the wing. But what
you are hearing at night are not cicadas, I think. As far I know they are
active only in daylight hours. The night shift (here) is taken by crickets,
whose song I find very musical. You don't get to see them much. They tend to
hide in cracks in the ground, which open up in the dry season. Sometimes
crickets come inside, black ones here.
glen herrmannsfeldt
(snip)
>> unsigned_result = 3740121002
1
>> Error: Integer too big for its kind at (1). This check can be disabled
>> with the option -fno-range-check
You can subtract 2**32, or find two values to IOR together.
(And remember that this is all twos-complement specific.)
-- glen
glen herrmannsfeldt
(snip)
> If an RNG has a large number of possible states (say, 2^192) then it's
> unlikely for any two runs to EVER reach identical states by any means
> other than starting from the same seed. This RNG has a period vastly
> in excess of "large", so for practical purposes it just won't
> happen.
That is true, but there are some assumptions. Since the
standard doesn't say much at all about the underlying algorithm,
it isn't so reliable a statement as it seems.
For one, many PRNGs have poor performance with some starting seeds,
or poor performance for some cycles, especially as viewed from
the low bits.
> However, because this is a combination of multiple component RNGs
> there is the related issue of overlapping within the period of one or
> more component RNGs. MWC has essentially zero chance of being in an
> identical state for any reason other than identical seeding, but both
> XS and CNG will frequently have identical states in two different
> runs. There does not seem to be any detectable correlation however
> unless (either MWC or (both XS and CNG)) have identical states between
> two different runs. Which is rather uncommon, though it does happen
> often enough to see it in the real world on occasion.
My previous comments were regarding the Fortran standard
RANDOM_SEED routine. I haven't thought about the seeding of
the RNG recently mentioned here.
-- glen
Gib Bogle
Ah, of course.
nm...@cam.ac.uk
Ron Shepard <ron-s...@NOSPAM.comcast.net> wrote:
>In article <i2ssst$nfi$1...@speranza.aioe.org>,
> glen herrmannsfeldt <g...@ugcs.caltech.edu> wrote:
>
>> In a large fraction of the cases, 2 billion different seeds
>> are enough, but one can still desire the appropriate randomness
>> from those different seeds.
>
>My understanding is that pseudorandom number generators return a
>sequence of values with some period. The numbers in that sequence
>eventually repeat with some, hopefully long, cycle. The put=array
>argument gives the starting point within that sequence, but it doesn't
>affect the cycle length or the "randomness" of the values, those things
>are determined by the underlying algorithm, not by the initial seed.
>
>The situation that needs to be avoided is to run the code once with one
>seed, and then run it again with another seed that results in an overlap
>of the sequences of values for the two runs. In some applications this
>is unimportant, but in other applications it would be bad for two
>supposedly independent runs to have a long sequence of identical values.
No, no, a thousand times, NO! That is NOT enough, though FAR too many
Web pages, published papers and books claim that it is. Disjointness
isn't even a poor relation of randomness.
This is easiest to see with a simple multiplicative congruential
generator, but many classes of generator have similar properties.
A sequence can be written as A.B^k, for k = 1...n, where B is fixed
by the properties of the generator but A also depends on the seed.
Two sequences X1 and X2 will always have the property that
A2.X1+A1.X2 = 0 modulo M.
So A1 and A2 need to be large and different, but you also need to
consider offset sequences - e.g. with A2 replaced by A2.B^p, for
all values of p < n. The spectral test theory tells you that you
will always have some values of p for which Q2.X1+Q1.X2 = 0 modulo M
and both Q1 and Q2 are very small. That's Bad News for many analyses.
This problem was a significant driving force behind the development
of the highly non-linear algorithms and composite algorithms. Both
are resistant to it, and the combination is very resistant. KISS4691
is one such.
Parallel random number generation is not easy, and 99% of the stuff
published on it is somewhere between grossly misleading and complete
nonsense.
Regards,
Nick Maclaren.
Uno
First of all, I think we're talking to the actual George Marsaglia here.
Thank you so much for posting. You may have displaced Terence as our
senior member.
Are you creating bigger numbers just to accomodate your age?:-)
$ gcc -Wall -Wextra geo1.c -o out
geo1.c: In function ‘MWC’:
geo1.c:5: warning: type defaults to ‘int’ in declaration of ‘c’
geo1.c:5: warning: type defaults to ‘int’ in declaration of ‘j’
geo1.c:5: warning: unused variable ‘i’
geo1.c: In function ‘main’:
geo1.c:21: warning: format ‘%22u’ expects type ‘unsigned int’, but
argument 2 has type ‘long unsigned int’
geo1.c:23: warning: format ‘%22u’ expects type ‘unsigned int’, but
argument 2 has type ‘long unsigned int’
geo1.c:24: warning: control reaches end of non-void function
$ ./out
MWC result=3740121002 ?
3740121002
KISS result=2224631993 ?
2224631993
$ cat geo1.c
static unsigned long xs=521288629,xcng=362436069,Q[4691];
unsigned long MWC(void) /*takes about 4.2 nanosecs or 238 million/
second*/
{unsigned long t,x,i; static c=0,j=4691;
j=(j<4690)? j+1:0;
x=Q[j];
t=(x<<13)+c+x; c=(t<x)+(x>>19);
return (Q[j]=t);
}
#define CNG ( xcng=69069*xcng+123 )
#define XS ( xs^=(xs<<13), xs^=(xs>>17), xs^=(xs<<5) )
#define KISS ( MWC()+CNG+XS ) /*138 million/sec*/
#include <stdio.h>
int main()
{unsigned long i,x;
for(i=0;i<4691;i++) Q[i]=CNG+XS;
for(i=0;i<1000000000;i++) x=MWC();
printf(" MWC result=3740121002 ?\n%22u\n",x);
for(i=0;i<1000000000;i++) x=KISS;
printf("KISS result=2224631993 ?\n%22u\n",x);
}
// gcc -Wall -Wextra geo1.c -o out
$
So, what is all this? In particular, is there something special about
the value of 3.7 billion?
--
Uno
Uno
Ron asks an interesting question, and I think you shoot from the hip and
may have missed on this one.
First of all, when crossposted to clc and sci.math, there is no
"standard" that does not require disambiquity. Given what it is, we
might say "the C standard" when talking about C, or its lurking uncle,
the fortran standard.
I don't agree that you need to "add" to random_seed a form that allows a
single integer variable supplied as a seed source. Any fortran program
that implements this will be fortran-conforming by virtue of there being
little specificity in either standard on this and no appetite to change
it in either community.
"disambiquity:" I kind of like that word.
--
Uno
David Bernier
[...]
> The situation that needs to be avoided is to run the code once with one
> seed, and then run it again with another seed that results in an overlap
> of the sequences of values for the two runs. In some applications this
> is unimportant, but in other applications it would be bad for two
> supposedly independent runs to have a long sequence of identical values.
I don't know how likely it would be for two seeds to produce
overlap in a short time. But I like using seeds that I've
almost surely never used before. I happen to use
the Mersenne Twister algorithm, but getting new
"random" seeds is always good for any algorithm.
> It would be nice if there were some prescription to generate initial
> seeds that would avoid this situation. However, I don't really know how
> this could be specified in the fortran standard without specifying the
> algorithm itself (which is not done, apparently on purpose). Anyone
> have any ideas how this could be done?
[...]
I often want to do just that: use new seeds. I used to go to
http://www.random.org/ and get random bits.
The way I do things now is to have a file "data" and take a sha256 hash value,
which is presented in hexadecimal. To avoid reusing the same seed,
I then append the hash-value to the file "data" .
This can be automated with a script which uses the "bash" shell (popular
in Linux); for a C-shell, something similar should work.
[The next step would be to be to access the script from
a fortran or C program ... ]
listing of rnd.sh file (3 lines long, made executable by: chmod u+x rnd.sh )
#! /bin/bash
sha256sum wcisha >> wcisha
sha256sum wcisha
The rnd.sh file is in my home directory, as well as the "data" file 'wcisha'.
In use:
[ ~]$ ./rnd.sh [Enter]
52897abc80205035b14986dd8d4d7b1549b98dca8c959f944f43b7c535e4f059 wcisha
[ ~]$ ./rnd.sh [Enter]
6822ad1bd1d1e453e66d2b65ea79f82a86a1ae53f90512b08f1feedcce6ae9f0 wcisha
David Bernier
Uno
If you were to comment out the PL/I command line that compiled this,
what would it be?
--
Uno
robin
| In comp.lang.fortran Gib Bogle <g.b...@auckland.no.spam.ac.nz> wrote:
| (snip)
|
| >> unsigned_result = 3740121002
| 1
| >> Error: Integer too big for its kind at (1). This check can be disabled
| >> with the option -fno-range-check
|
| You can subtract 2**32, or find two values to IOR together.
Subtracting 2**32 won't do anything because it is entirely put of range
(i.e., requires 33 bits.
The magic number is to subtract or add 2**31.
robin
| I think of 3.5 billion as a value that busts types in C. I would have
| thought that was going on except for the second error:
It doesn't "bust" types in C, because it is a valid unsigned number
for a 32-bit machine.
| geo1.f90:72.32:
|
| unsigned_result = 2224631993
| 1
| Error: Integer too big for its kind at (1). This check can be disabled
| with the option -fno-range-check
Same as above. It's valid unsigned number.
Ron Shepard
wrote:
> In article <ron-shepard-969C...@news60.forteinc.com>,
> Ron Shepard <ron-s...@NOSPAM.comcast.net> wrote:
> >In article <i2ssst$nfi$1...@speranza.aioe.org>,
> > glen herrmannsfeldt <g...@ugcs.caltech.edu> wrote:
> >
> >> In a large fraction of the cases, 2 billion different seeds
> >> are enough, but one can still desire the appropriate randomness
> >> from those different seeds.
> >
> >My understanding is that pseudorandom number generators return a
> >sequence of values with some period. The numbers in that sequence
> >eventually repeat with some, hopefully long, cycle. The put=array
> >argument gives the starting point within that sequence, but it doesn't
> >affect the cycle length or the "randomness" of the values, those things
> >are determined by the underlying algorithm, not by the initial seed.
> >
> >The situation that needs to be avoided is to run the code once with one
> >seed, and then run it again with another seed that results in an overlap
> >of the sequences of values for the two runs. In some applications this
> >is unimportant, but in other applications it would be bad for two
> >supposedly independent runs to have a long sequence of identical values.
>
> No, no, a thousand times, NO! That is NOT enough, though FAR too many
> Web pages, published papers and books claim that it is. Disjointness
> isn't even a poor relation of randomness.
That may be correct, but without specifying the PRNG algorithm, that
may be all that could be expected for a "good" choice of initial
seeds.
I guess my earlier question wasn't clear enough, so I will try to
rephrase it. What I'm looking for is whether there is a practical
way to define a subroutine like
call RANDOM_NEW_SEED( old_seed, n, new_seed )
new_seed(:) and old_seed(:) are integer arrays of the appropriate
length for the (unspecified) PRNG algorithm. old_seed(:) was some
old seed that was used elsewhere, perhaps in some previous execution
of the same code, or perhaps in some other initialization of the
PRNG on another node of a parallel machine. The integer n is the
number of values that will be generated, or some estimate thereof.
new_seed(:) is the new seed that is guaranteed to result in a
sequence of values that do not overlap with the previous run
(assuming the cycle length is long enough, of course).
I guess one might be able to impose other conditions on the values
that will be generated, such as those you mention, but again,
without actually specifying the PRNG algorithm, it isn't clear to me
how that could be done generally.
> Parallel random number generation is not easy, and 99% of the stuff
> published on it is somewhere between grossly misleading and complete
> nonsense.
I think in the parallel case, one would want to be able to generate
a seed to produce values that are guaranteed not to overlap with any
other node. Maybe something like
call RANDOM_NEW_SEED( old_seed, n, new_seed, my_node )
would be a sufficient interface. new_seed(:) would depend on
my_node in such a way that the generated sequence would not overlap
with that produced by any other possible value of my_node (again,
assuming the cycle length is long enough to satisfy that request).
$.02 -Ron Shepard
Gib Bogle
difference between the C-generated and Fortran-generated sequences is immaterial
(from the randomness point of view), or how to make the Fortran work in exactly
the same way as the C code. Anyone else is free to contribute in his place ...
orz
> >> In article <fa9dd141-823a-4179-a80b-18d214a3e...@t2g2000yqe.googlegroups.com>,
I think you have part of it backwards:
elseif (x < 0) then
tLTx = 1
elseif (t < 0) then
tLTx = 0
endif
should be:
elseif (x < 0) then
tLTx = 0
elseif (t < 0) then
tLTx = 1
endif
...if you want to produce an sequence of numbers identical to his C
code. Though the way you have it written produces significantly
higher quality output according to some statistical tests.
Another issue to worry about is the rightshift must be an unsigned
rightshift, but a quick googling suggests that fortran ISHFT is
unsigned anyway.
orz
intended to be a reply to Gib Bogle's most recent post.
Gib Bogle
> I think you have part of it backwards:
> elseif (x < 0) then
> tLTx = 1
> elseif (t < 0) then
> tLTx = 0
> endif
>
> should be:
> elseif (x < 0) then
> tLTx = 0
> elseif (t < 0) then
> tLTx = 1
> endif
>
> ...if you want to produce an sequence of numbers identical to his C
> code. Though the way you have it written produces significantly
> higher quality output according to some statistical tests.
>
> Another issue to worry about is the rightshift must be an unsigned
> rightshift, but a quick googling suggests that fortran ISHFT is
> unsigned anyway.
I think you are wrong. The code as I have it is correct, as far as I'm
concerned, and it produces identical results to the C code, when the signed
values are treated as unsigned (or vice versa). Remember, a signed negative
(e.g. x < 0) is representing an unsigned value greater than 2^31 - 1. In the
conditional code you refer to, signed x < 0 and signed t > 0 => unsigned x >=
2^31 and unsigned t < 2^31. Q.E.D.
As Glen pointed out, this all assumes a two's complement number representation,
such as is used by my Intel CPU/Intel compiler combination.
orz
Whoops, you are correct, I somehow got the two mixed up even when
testing them in code.
But in that case I don't know what you meant in your previous post
when you asked about differences between your output and the original
output. If the information produced is identical, then the
information produced is identical.
Gib Bogle
> Whoops, you are correct, I somehow got the two mixed up even when
> testing them in code.
>
> But in that case I don't know what you meant in your previous post
> when you asked about differences between your output and the original
> output. If the information produced is identical, then the
> information produced is identical.
Please look at the code. You'll see that tLTx is used only to compute
c_this_works. The code as I posted it contains the method of computing c
suggested by George. This generates results different from the C code. If you
uncomment the line c = c_this_works you get identical results to the C code.
I'm sure I already said this.
orz
Yes. Sorry. I was reading backwards from your last post and ended up
missing the point. And getting confused on the sign.
Anyway, the issue is that Georges code uses a different definition of
sign than your implementation of it - his code is actually correct if
sign(x) is 1 if x is positive and 0 if x is negative. Since your sign
function returns -1 on negative, using it produces the wrong
results.
side note: The incorrect results produced that way at a appear to have
vaguely similar statistical properties as the original C codes output,
passing and failing the same tests that the original C code does in my
brief tests.
Gib Bogle
> Yes. Sorry. I was reading backwards from your last post and ended up
> missing the point. And getting confused on the sign.
>
> Anyway, the issue is that Georges code uses a different definition of
> sign than your implementation of it - his code is actually correct if
> sign(x) is 1 if x is positive and 0 if x is negative. Since your sign
> function returns -1 on negative, using it produces the wrong
> results.
>
> side note: The incorrect results produced that way at a appear to have
> vaguely similar statistical properties as the original C codes output,
> passing and failing the same tests that the original C code does in my
> brief tests.
Interesting, who would have guessed that there is a language in which sign(-1) = 0.
Ilmari Karonen
On 2010-07-30, Ron Shepard <ron-s...@NOSPAM.comcast.net> wrote:
> In article <i2tvn7$ldf$1...@smaug.linux.pwf.cam.ac.uk>, nm...@cam.ac.uk
> wrote:
>> In article <ron-shepard-969C...@news60.forteinc.com>,
>> Ron Shepard <ron-s...@NOSPAM.comcast.net> wrote:
>> >
>> >The situation that needs to be avoided is to run the code once with one
>> >seed, and then run it again with another seed that results in an overlap
>> >of the sequences of values for the two runs. In some applications this
>> >is unimportant, but in other applications it would be bad for two
>> >supposedly independent runs to have a long sequence of identical values.
>>
>> No, no, a thousand times, NO! That is NOT enough, though FAR too many
>> Web pages, published papers and books claim that it is. Disjointness
>> isn't even a poor relation of randomness.
>
>
>> Parallel random number generation is not easy, and 99% of the stuff
>> published on it is somewhere between grossly misleading and complete
>> nonsense.
>
> I think in the parallel case, one would want to be able to generate
> a seed to produce values that are guaranteed not to overlap with any
> other node. Maybe something like
>
> call RANDOM_NEW_SEED( old_seed, n, new_seed, my_node )
>
> would be a sufficient interface. new_seed(:) would depend on
> my_node in such a way that the generated sequence would not overlap
> with that produced by any other possible value of my_node (again,
> assuming the cycle length is long enough to satisfy that request).
Let me try to demonstrate, using a simplified (and admittedly somewhat
artificial) counterexample, why lack of overlap is not a sufficient
condition for independence:
Assume we have a "random oracle" R: Z -> [0,1) which takes as its
input an integer and returns a uniformly distributed random real
number between 0 and 1, such that the same input always produces the
same output, but that the outputs for different inputs are completely
independent.
Given such an oracle, we can construct a perfect PRNG P: Z -> [0,1)^N
which takes as its input an integer k, and returns the sequence <R(k),
R(k+1), R(k+2), ...>. Obviously, the sequence generated by P would be
indistinguishable from random (since it is indeed, by definition,
random) and non-repeating.
Now, what if we wanted several independent streams of random numbers?
Obviously, we can't just pass different seed values to P, since we
know that the streams would eventually overlap. We could solve the
problem by modifying P, e.g. to make it return the sequence
<R(f(k,0)), R(f(k,1)), R(f(k,2)), ...> instead, where f is an
injective function from Z^2 to Z. But what if we wanted to do it
_without_ modifying P?
One _bad_ solution would be to define a new generator Q: Z -> [0,1)^N
as Q(k)_i = frac(P(0)_i + ak), where a is some irrational number and
frac: R -> [0,1) returns the part of a real number "after the decimal
point" (i.e. frac(x) = x - floor(x)).
Clearly, the sequences returned by Q for different seed values would
never overlap, and, individually, they would each still be perfectly
random. Yet, just as clearly, all those sequences would be related by
a very trivial linear relation that would be blatantly obvious if you,
say, plotted two of them against each other.
The _right_ solution, as I suggested above, would've been to redesign
the underlying generator so that the different streams would be not
just non-overlapping but actually statistically independent. The same
general advice holds for practical PRNGs too, not just for my
idealized example.
You can't just take any arbitrary PRNG, designed for generating a
_single_ stream of random-seeming numbers, and expect the output to
still look random if you get to compare several distinct output
streams generated from related seeds. It might turn out that it does
look so, if the designer of the PRNG was careful or just lucky. But
designing a PRNG to satisfy such a requirement is a strictly harder
problem than simply designing it to generate a single random-looking
stream.
--
Ilmari Karonen
To reply by e-mail, please replace ".invalid" with ".net" in address.
Uno
Now that's an interesting read. Ron is a person who makes "mistakes" in
contemporary scientific programming, but he always moves the scrum
forward in a useful way.
I notice that follow-ups are set to sci.math, which I've never really
read much. When I have needed mathematical help on usenet, I've turned
to the germans in de.sci.mathematik, who have a very well-stylized notation.
I'm used to seeing the notation of functions from the point of abstract
algebra (in particular the Gallian text, which follows me everywhere).
But there's been a lot of water under the bridge (and one concussion)
since I had the leisure of reading college textbooks.
So my question for you is going to be more of a request, namely, that
you say a few more words about the mappings, what sets are involved, and
what things like Z^2 mean exactly.
To reveal my naivete I would think a perfect prng would look like:
P: Z^N -> [0,1)^N
By the way, nice job to get the interval [0,1) correct. That's a
fortranism.
Cheers,
--
Uno
sturlamolden
> The _right_ solution, as I suggested above, would've been to redesign
> the underlying generator so that the different streams would be not
> just non-overlapping but actually statistically independent.
Yes. For the Mersenne Twister it has been suggested to use different
characterisic polynomials that are "relatively prime" to each other.
Another suggestion is to derive a fast "jump ahead" algorithm allowing
parallel work share with the same polynomial. Just picking different
seeds and praying they will not overlap/interact is risky.
http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/DC/dgene.pdf
http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/ARTICLES/jump-seta-lfsr.pdf
The other option is to observe that random numbers can usually be
produced much faster than they are consumed. That is because the
simple integer manipulation in the PRNG is much less expensive than
the complicated floating point maths in numerical algorithms. If we
can only consume a fraction of the 138 million random numbers per
second from KISS4691, we don't need a parallel PRNG. We can let one
processor fill N (spinlock synchronized) queues (fifos) with
pseudorandom integers from one KISS4691 sequence. The queues will
serve as sources of random bits for a group of N relatively/mutually
independent PRNGs. Amortizing the synchronization overhead comes down
to buffering inside the prng frontend harvesting the queues. E.g. with
an internal buffer of 256 integers, one only has to aquire the
spinlock and harvest a queue on every 256 call to the prng. The best
buffer size is dependent on the hardware (e.g. hierarchical memory).
Here we could note that "Keep it simple stupid" is fast compared to
Mersenne Twister, which makes it more suitable for this kind of
parallel algorithms.
Ron Shepard
<e5f36412-ddf2-4ed5...@u4g2000prn.googlegroups.com>,
orz <cdh...@gmail.com> wrote:
> If an RNG has a large number of possible states (say, 2^192) then it's
> unlikely for any two runs to EVER reach identical states by any means
> other than starting from the same seed. This RNG has a period vastly
> in excess of "large", so for practical purposes it just won't
> happen.
If you think about it for a second, you can see that this is not true
for a PRNG with the standard fortran interface. This interface allows
the seed values to be both input and output. The purpose (presumably,
although this is not required by the standard) for this is to allow the
seed to be extracted at the end of one run, and then used to initiate
the seed for the next run. This eliminates the possibility of the two
runs having overlapping sequences (provided the PRNG cycle is long
enough, of course).
But, with the ability to read the internal state of the PNRG, this would
allow the user to generate the exact same sequence twice (useful for
debugging, for example). Or the user could generate two sequences that
have arbitrary numbers of overlapping values (the sequences offset by 1,
offset by 2, offset by 100, or whatever). I can't think of any
practical reason for doing this. However, since it can be done with
100% certainty, it is not "unlikely" for it to occur, assuming the user
wants it to happen.
In a previous post, I suggested a new (intrinsic fortran) subroutine of
the form
call RANDOM_NEW_SEED( old_seed, n, new_seed )
One way to implement this subroutine would be for new_seed(:) to be the
internal state corresponding to the (n+1)-th element of the sequence
that results (or would result) from initialization with old_seed(:). Of
course, for this to be practical, especially for large values of n, this
would require the PRNG algorithm to allow the internal state for
arbitrary elements within a sequence to be determined in this way. That
would place practical restrictions on which PRNG algorithms can be used
with this interface, especially if additional features such as
independence of the separate sequences are also required.
$.02 -Ron Shepard
sturlamolden
I did a speed comparison on my laptop of KISS4691 against Mersenne
Twister 19937. KISS4691 produced about 110 million random numbers per
second, whereas MT19937 produced 118 million. If the compiler was
allowed to inline KISS4691, the performance increased to 148 millon
per second. The function call overhead is imporant, and without taking
it into consideration, MT19937 will actually be faster than KISS4691.
Also note that a SIMD oriented version of MT19937 is twice as fast as
the one used in my tests. There is also a version of Mersenne Twister
suitable for parallel processing and GPUs. So a parallel or SIMD
version of Mersenne Twister is likely to be a faster PRNG than
KISS4691.
Then for the question of numerical quality: Is KISS4691 better than
MT19937? I don't feel qualified to answer this. Marsaglia is the
world's foremost authority on random number generators, so I trust his
advice, but MT19937 is claimed to give impeckable quality except for
crypto.
For those who need details:
The speed test was basically running the PRNGs for five seconds or
more, querying Windows' performance counter on every ten million call,
and finally subtracting an estimate of the timing overhead. The C
compiler was Microsoft C/C++ version 15.00.30729.01 (the only
optimization flag I used was /O2, i.e. maximize for speed). The
processor is AMD Phenom N930 @ 2.0 Gz.
Sturla Molden
glen herrmannsfeldt
> In a previous post, I suggested a new (intrinsic fortran) subroutine of
> the form
> call RANDOM_NEW_SEED( old_seed, n, new_seed )
> One way to implement this subroutine would be for new_seed(:) to be the
> internal state corresponding to the (n+1)-th element of the sequence
> that results (or would result) from initialization with old_seed(:). Of
> course, for this to be practical, especially for large values of n, this
> would require the PRNG algorithm to allow the internal state for
> arbitrary elements within a sequence to be determined in this way. That
> would place practical restrictions on which PRNG algorithms can be used
> with this interface, especially if additional features such as
> independence of the separate sequences are also required.
The standard mostly does not specify what is, or is not, practical.
DO I=1,INT(1000000,SELECTED_INT_KIND(20))**3
ENDDO
For many generators, one should not specify too large a value
for the above RANDOM_NEW_SEED routine.
-- glen
Noob
> [...] MT19937 is claimed to give impeckable quality [...]
It cannot be eaten in small bites? ;-)
Uno
I'd be curious to see the source you used for this purpose. This is
very-slightly adapted from Dann Corbitt in ch 13 of unleashed.
$ gcc -Wall -Wextra cokus2.c -o out
$ ./out
3510405877 4290933890 2191955339 564929546 152112058 4262624192
2687398418
268830360 1763988213 578848526 4212814465 3596577449 4146913070
950422373
1908844540 1452005258 3029421110 142578355 1583761762 1816660702
2530498888
1339965000 3874409922 3044234909 1962617717 2324289180 310281170
981016607
908202274 3371937721 2244849493 675678546 3196822098 1040470160
3059612017
3055400130 2826830282 2884538137 3090587696 2262235068 3506294894
2080537739
1636797501 4292933080 2037904983 2465694618 1249751105 30084166
112252926
1333718913 880414402 334691897 3337628481 17084333 1070118630
...
3009858040 3815089086 2493949982 3668001592 1185949870 2768980234
3004703555
1411869256 2625868727 3108166073 3689645521 4191339889 1933496174
1218198213
3716194408 1148391246 1345939134 3517135224 3320201329 4292973312
3428972922
1172742736 275920387 617064233 3754308093 842677508 529120787
1121641339
$ cat cokus2.c
#include <stdio.h>
#include <stdlib.h>
#include "mtrand.h"
/*
uint32 must be an unsigned integer type capable of holding at least 32
bits; exactly 32 should be fastest, but 64 is better on an Alpha with
GCC at -O3 optimization so try your options and see what's best for you
*/
typedef unsigned long uint32;
/* length of state vector */
#define N (624)
/* a period parameter */
#define M (397)
/* a magic constant */
#define K (0x9908B0DFU)
/* mask all but highest bit of u */
#define hiBit(u) ((u) & 0x80000000U)
/* mask all but lowest bit of u */
#define loBit(u) ((u) & 0x00000001U)
/* mask the highest bit of u */
#define loBits(u) ((u) & 0x7FFFFFFFU)
/* move hi bit of u to hi bit of v */
#define mixBits(u, v) (hiBit(u)|loBits(v))
/* state vector + 1 extra to not violate ANSI C */
static uint32 state[N + 1];
/* next random value is computed from here */
static uint32 *next;
/* can *next++ this many times before reloading */
static int left = -1;
/*
**
** We initialize state[0..(N-1)] via the generator
**
** x_new = (69069 * x_old) mod 2^32
**
** from Line 15 of Table 1, p. 106, Sec. 3.3.4 of Knuth's
** _The Art of Computer Programming_, Volume 2, 3rd ed.
**
** Notes (SJC): I do not know what the initial state requirements
** of the Mersenne Twister are, but it seems this seeding generator
** could be better. It achieves the maximum period for its modulus
** (2^30) iff x_initial is odd (p. 20-21, Sec. 3.2.1.2, Knuth); if
** x_initial can be even, you have sequences like 0, 0, 0, ...;
** 2^31, 2^31, 2^31, ...; 2^30, 2^30, 2^30, ...; 2^29, 2^29 + 2^31,
** 2^29, 2^29 + 2^31, ..., etc. so I force seed to be odd below.
**
** Even if x_initial is odd, if x_initial is 1 mod 4 then
**
** the lowest bit of x is always 1,
** the next-to-lowest bit of x is always 0,
** the 2nd-from-lowest bit of x alternates ... 0 1 0 1 0 1 0 1 ... ,
** the 3rd-from-lowest bit of x 4-cycles ... 0 1 1 0 0 1 1 0 ... ,
** the 4th-from-lowest bit of x has the 8-cycle ... 0 0 0 1 1 1 1 0 ... ,
** ...
**
** and if x_initial is 3 mod 4 then
**
** the lowest bit of x is always 1,
** the next-to-lowest bit of x is always 1,
** the 2nd-from-lowest bit of x alternates ... 0 1 0 1 0 1 0 1 ... ,
** the 3rd-from-lowest bit of x 4-cycles ... 0 0 1 1 0 0 1 1 ... ,
** the 4th-from-lowest bit of x has the 8-cycle ... 0 0 1 1 1 1 0 0 ... ,
** ...
**
** The generator's potency (min. s>=0 with (69069-1)^s = 0 mod 2^32) is
** 16, which seems to be alright by p. 25, Sec. 3.2.1.3 of Knuth. It
** also does well in the dimension 2..5 spectral tests, but it could be
** better in dimension 6 (Line 15, Table 1, p. 106, Sec. 3.3.4, Knuth).
**
** Note that the random number user does not see the values generated
** here directly since reloadMT() will always munge them first, so maybe
** none of all of this matters. In fact, the seed values made here could
** even be extra-special desirable if the Mersenne Twister theory says
** so-- that's why the only change I made is to restrict to odd seeds.
*/
void mtsrand(uint32 seed)
{
register uint32 x = (seed | 1U) & 0xFFFFFFFFU,
*s = state;
register int j;
for (left = 0, *s++ = x, j = N; --j;
*s++ = (x *= 69069U) & 0xFFFFFFFFU);
}
uint32
reloadMT(void)
{
register uint32 *p0 = state,
*p2 = state + 2,
*pM = state + M,
s0,
s1;
register int j;
if (left < -1)
mtsrand(4357U);
left = N - 1, next = state + 1;
for (s0 = state[0], s1 = state[1], j = N - M + 1; --j; s0 = s1, s1
= *p2++)
*p0++ = *pM++ ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0U);
for (pM = state, j = M; --j; s0 = s1, s1 = *p2++)
*p0++ = *pM++ ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0U);
s1 = state[0], *p0 = *pM ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K
: 0U);
s1 ^= (s1 >> 11);
s1 ^= (s1 << 7) & 0x9D2C5680U;
s1 ^= (s1 << 15) & 0xEFC60000U;
return (s1 ^ (s1 >> 18));
}
uint32 mtrand(void)
{
uint32 y;
if (--left < 0)
return (reloadMT());
y = *next++;
y ^= (y >> 11);
y ^= (y << 7) & 0x9D2C5680U;
y ^= (y << 15) & 0xEFC60000U;
return (y ^ (y >> 18));
}
#define UNIT_TEST
#ifdef UNIT_TEST
int main(void)
{
int j;
/* you can seed with any uint32, but the best are odds in 0..(2^32
- 1) */
mtsrand(4357U);
/* print the first 2,002 random numbers seven to a line as an
example */
for (j = 0; j < 2002; j++)
printf(" %10lu%s", (unsigned long) mtrand(), (j % 7) == 6 ?
"\n" : "");
return (EXIT_SUCCESS);
}
#endif
// gcc -Wall -Wextra cokus2.c -o out
$ cat mtrand.h
/*
** The proper usage and copyright information for
** this software is covered in DSCRLic.TXT
** This code is Copyright 1999 by Dann Corbit
*/
/*
** Header files for Mersenne Twister pseudo-random number generator.
*/
extern void mtsrand(unsigned long seed);
extern unsigned long reloadMT(void);
extern unsigned long mtrand(void);
$
--
Uno
sturlamolden
> I'd be curious to see the source you used for this purpose. This is
> very-slightly adapted from Dann Corbitt in ch 13 of unleashed.
Dann Corbit
6f0708...@q35g2000yqn.googlegroups.com>, sturla...@yahoo.no
says...
> http://folk.uio.no/sturlamo/prngtest.zip
>
I got this:
Mersenne Twister 19937 speed, single (Hz): 1.63278e+008
Mersenne Twister 19937 speed, array (Hz): 1.3697e+008
KISS 4691 speed, single (Hz): 1.86338e+008
KISS 4691 speed, array (Hz): 1.87675e+008
robin
James Waldby
> "Uno" <merrilljensen> wrote:
[snip code]
>> If you were to comment out the PL/I command line that compiled this,
>> what would it be?
>
> ???
Does that mean you don't understand Uno's question,
or don't know the answer?
In case you don't understand the question, it appears
to be: "What command is used to compile the code?"
--
jiw
Dann Corbit
It will depend on the operating system.
Probably JCL along the lines of:
// EXEC PL1LFCLG,REGION.PL1L=256K
Peter Flass
or "plic -C" <filename>
sturlamolden
> Mersenne Twister 19937 speed, single (Hz): 1.63278e+008
> Mersenne Twister 19937 speed, array (Hz): 1.3697e+008
> KISS 4691 speed, single (Hz): 1.86338e+008
> KISS 4691 speed, array (Hz): 1.87675e+008
Those numbers are in 100 million samples per second (10^8), so you
have the same order of magnitude I reported. The array version of KISS
is inlined (by macro expansion).
These numbers are likely dependent on compiler and hardware.
In you case, KISS4691 is always faster than MT19937. That is what I
expected to see on my laptop as well, but did not. The speed
difference is not very substantial, though, less than a factor of 2.
I am more concerned about numerical quality. Which one should we use
based on that?
Sturla
Uno
I'll restate the question, and I'm sure you'll get my drift. When I
compile off a command line, I keep the command lines I used as the final
comments in that file. So there might, in fortran, exist
implicit real
pi = 4.0 * atan(1.0)
print *, pi
endprogram
!here it comes, the goocher:
! gfortran pi1.f90 -o out
1) What did you name this pli thing?
2) What command compiled it?
3) How does one comment in pli?
4) How does one caquire a pli facilty on ubuntu?
Thanks for your comment,
and holy balls did I get healthy today,
--
Uno
orz
MT19937 fails only a few relatively obscure empirical tests. There's
a widespread belief that it's a good RNG and its few failed tests have
little real world consequence. Reduced strength versions of MT19937
have a strong tendency to fail additional tests.
KISS4691 fails no single-seed empirical tests, but it does fail some
empirical tests for correlation between different seeds. Reduced
strength versions tend to do better than reduced strength versions of
MT19937, but not by much.
If you want to go all out for quality over speed, the standard method
is to simply encrypt the output of some RNG with AES or some
comparable algorithm. That's much much slower than MT19937 or
KISS4691, but pretty much perfect quality. Or you could xor the
output of KISS4691 and MT19937.
If you're willing to go with unknown-shortest-cycle-length RNGs then
there are other options, many of which do perfectly on all empirical
tests, and some of those also offer other advantages (such as being
substantially faster than either KISS4691 or MT19937, or offering some
degree of cryptographic security, or passing all empirical tests even
in drastically reduced strength versions). But the famous people in
RNG theory generally seem to think that unknown-shortest-cycle-length
RNGs should not be trusted.
orz
I have to correct myself for swapping 0 and 1 *again*. And I'm not
even dyslexic, so far as I know.
His code assumed sign returned 1 on negative, and 0 otherwise, as in a
simple unsigned 31 bit rightshift. The exact opposite of what I
said.
robin
| It will depend on the operating system.
| Probably JCL along the lines of:
| // EXEC PL1LFCLG,REGION.PL1L=256K
or
PLI <filename>
on the PC for various compilers.
robin
| On Tue, 03 Aug 2010 20:41:15 +1000, robin wrote:
| > "Uno" <merrilljensen> wrote:
| [snip code]
| >> If you were to comment out the PL/I command line that compiled this,
| >> what would it be?
| >
| > ???
|
| Does that mean you don't understand Uno's question,
| or don't know the answer?
It means that the question makes no sense.
Uno
I'll restate the question, and I'm sure you'll get my drift. When I
compile off a command line, I keep the command lines I used as the final
comments in that file. So there might, in fortran, exist
implicit real
pi = 4.0 * atan(1.0)
print *, pi
endprogram
!here it comes, the goocher:
! gfortran pi1.f90 -o out
1) What did you name this pli thing?
2) What command compiled it?
3) How does one comment in pli?
4) How does one acquire a pli facilty on ubuntu?
--
Uno
Uno
Zero: the other one.
Zero: One-Lite.
Telling left from right is sometimes the hardest thing.
--
Uno
robin
| I'll restate the question, and I'm sure you'll get my drift. When I
| compile off a command line, I keep the command lines I used as the final
| comments in that file. So there might, in fortran, exist
|
| implicit real
| pi = 4.0 * atan(1.0)
| print *, pi
| endprogram
|
| !here it comes, the goocher:
|
| ! gfortran pi1.f90 -o out
|
| 1) What did you name this pli thing?
RNG-2010.PLI
| 2) What command compiled it?
PL/I RNG-2010
| 3) How does one comment in pli?
/* Stuff */
mecej4
Those kinds of basic questions are mostly covered by the PL/I FAQ, which is
posted quite regularly in this newsgroup.
As far as I am aware, there is no production quality native PL/I compiler
available for Linux. There is VisualAge PL/I for Windows, which IBM makes
available through its Scholars Program to those who qualify or which may be
purchased (at significant cost) as part of the Rational Developer for
System Z product.
-- mecej4
Peter Flass
>
> As far as I am aware, there is no production quality native PL/I compiler
> available for Linux.
I wish everyone would stop repeating this.
Not that I want to plug a competing product, but Micro Focus Open PL/I
runs on Linux. (http://www.microfocus.com/products/studio/open-pli.aspx)
What's your definition of "production quality?" I haven't tried this,
but it at least sounds good.
Of course Micro Focus invites this anonymity, since their marketing is
probably second only to IBM's in its badness: no pricing, no demo, no
advertising or promotion, etc.
If you want something to try *now*, look at Iron Spring
(http://www.iron-spring.com/) - possibly not yet "production quality"
but getting there. If it doesn't do what you want, just ask.
Selling commercial software for Linux presents some challenges.
Gib Bogle
I tried your suggestion, replacing the sign() function with one based on your
rule. It doesn't work, i.e. it gives results different from the C code.
orz
If I made yet another mistake on that then I'm going to throw
something.
It produced identical results for me when I tried it in C using this
code:
typedef unsigned int Uint32;
Uint32 sign(Uint32 value) {return value>>31;}
Uint32 index, carry, data[4691];
Uint32 mwc4691() {
index = (index < 4690) ? index + 1 : 0;
Uint32 x = data[index];
Uint32 t = (x << 13) + carry + x;
if (sign((x<<13)+carry) == sign(x))
carry = sign(x) + (x>>19);
else carry = 1 - sign(t) + (x>>19);
data[index] = t;
return t;
}
I think the only difference between that and the algorithm he
expressed in pseudocode is that I added parentheses on the first left-
shift (he was apparently assuming higher precedence on the leftshift
operator than C uses). Maybe that's why you're getting different
results?
Gib Bogle
> If I made yet another mistake on that then I'm going to throw
> something.
>
> It produced identical results for me when I tried it in C using this
> code:
>
> typedef unsigned int Uint32;
> Uint32 sign(Uint32 value) {return value>>31;}
> Uint32 index, carry, data[4691];
> Uint32 mwc4691() {
> index = (index < 4690) ? index + 1 : 0;
> Uint32 x = data[index];
> Uint32 t = (x << 13) + carry + x;
> if (sign((x<<13)+carry) == sign(x))
> carry = sign(x) + (x>>19);
> else carry = 1 - sign(t) + (x>>19);
> data[index] = t;
> return t;
> }
>
> I think the only difference between that and the algorithm he
> expressed in pseudocode is that I added parentheses on the first left-
> shift (he was apparently assuming higher precedence on the leftshift
> operator than C uses). Maybe that's why you're getting different
> results?
You do realize that there's no unsigned integer type in Fortran? I don't think
there's much point in trying to predict what the Fortran code does without using
a Fortran compiler. Don't you have one?
Dann Corbit
g.b...@auckland.no.spam.ac.nz says...
Maybe he is using Sun's Fortran 95 compiler.
From:
http://docs.sun.com/source/819-0492/4_f95.html
"4.5 Unsigned Integers
The Fortran 95 compiler accepts a new data type, UNSIGNED, as an
extension to the language. Four KIND parameter values are accepted with
UNSIGNED: 1, 2, 4, and 8, corresponding to 1-, 2-, 4-, and 8-byte
unsigned integers, respectively.
The form of an unsigned integer constant is a digit-string followed by
the upper or lower case letter U, optionally followed by an underscore
and kind parameter. The following examples show the maximum values for
unsigned integer constants:
255u_1
65535u_2
4294967295U_4
18446744073709551615U_8
Expressed without a kind parameter (12345U), the default is the same as
for default integer. This is U_4 but can be changed by the -xtypemap
option, which will change the kind type for default unsigned integers.
Declare an unsigned integer variable or array with the UNSIGNED type
specifier:
UNSIGNED U
UNSIGNED(KIND=2) :: A
UNSIGNED*8 :: B"
orz
Let me add another qualifier. The psuedocode he posted is correct
given the following assumptions:
1. The integers are 32 bit integers.
2. The "sign" function returns the value right-shifted by 31 bits (0
for non-negative, 1 for negative).
3. Shift operators are treated as higher precedence or parentheses
are added around that shift.
4. The right-shifts are treated as unsigned right-shifts.
Aside from those qualifiers, the results are independent of whether
the numbers are considered to be signed or unsigned. I probably
should have mentioned #4, but since I believe Fortran treats all right-
shifts as unsigned, it should not matter in Fortran.
Some C code which produces identical results and illustrates that:
typedef unsigned int Uint32;
typedef signed int Sint32;
Sint32 index, carry, data[4691];
Sint32 rshift(Sint32 value, int bits) {return Uint32(value)>>bits;}
Sint32 sign(Sint32 value) {return rshift(value,31);}
Sint32 mwc4691() {
index = (index < 4690) ? index + 1 : 0;
Sint32 x = data[index];
Sint32 t = (x << 13) + carry + x;
if (sign((x<<13)+carry) == sign(x)) carry = sign(x) + rshift(x,19);
else carry = 1 - sign(t) + rshift(x,19);
data[index] = t;
return t;
}
Perhaps you could post the Fortran code that is producing incorrect
results?
Richard Maine
> In article <i3ssd1$be8$1...@speranza.aioe.org>,
> g.b...@auckland.no.spam.ac.nz says...
> > You do realize that there's no unsigned integer type in Fortran?
...
> Maybe he is using Sun's Fortran 95 compiler.
Though if my admittedly fallible memory serves me, the unsigned integer
in Sun's f95 has restrictions that are likely to surprise many casual
users. Being more specific might stretch my memory a bit farther than I
trust, but I think I recall a lack of mixed-mode operations, for
example. So one would have to tread pretty carefully in order to
correctly transcribe pseudocode into working code. For example, I'm not
at all sure that one can do such things as adding 1 to an unsigned
integer; I think you might have to add 1U instead.
--
Richard Maine | Good judgment comes from experience;
email: last name at domain . net | experience comes from bad judgment.
domain: summertriangle | -- Mark Twain
glen herrmannsfeldt
(snip)
> 3. Shift operators are treated as higher precedence or parentheses
> are added around that shift.
The C precedence rules for shifts are widely believed to
be wrong, including by the developers of C, but it is considered
too late to change. Best to always use parentheses around them.
(Many Fortran rules haven't been changed for the same reason.)
-- glen
robin
It's perfectly possible to predict how to program the operation
(unsigned arithmetic) in Fortran provided that we make the
assumption that the hardware provides twos complement arithmetic.
Gib Bogle
It's possible to predict, but more reliable to test. I can use Fortran to
reproduce the results generated by the C code, but not using the pseudo-code
that George provided. Unfortunately he seems to have lost interest in this
thread. Perhaps you'd like to suggest the Fortran code to implement the method
that George provided in C (assuming twos complement).
glen herrmannsfeldt
> It's possible to predict, but more reliable to test.
> I can use Fortran to reproduce the results generated by the
> C code, but not using the pseudo-code that George provided.
> Unfortunately he seems to have lost interest in this thread.
> Perhaps you'd like to suggest the Fortran code to implement
> the method that George provided in C (assuming twos complement).
You also have to assume 32 bit, or IAND() off the extra bits
(of more than 32).
There is no guarantee that twos complement overflow gives
the low bits of the full result, but it usually does.
(Fortran allows for almost anything to happen.)
> static unsigned long xs=521288629,xcng=362436069,Q[4691];
> unsigned long MWC(void) /*takes about 4.2 nanosecs or 238 million/
> second*/
> {unsigned long t,x,i; static c=0,j=4691;
c, j, should have the SAVE attribute.
> j=(j<4690)? j+1:0;
> x=Q[j];
> t=(x<<13)+c+x; c=(t<x)+(x>>19);
Fortran's ISHFT guarantees a logical shift (shift in zeros).
To do the equivalent of an unsigned less than in twos-complement
invert the sign bit before comparing.
j=j+1
if(j>4691) j=0
x=Q(j)
t=ishft(x,13)+c+x
c=ishft(x,-19)
if(ieor(t,ishft(1,-31)).lt.ieor(x,ishft(1,-31))) c=c+1
(Hopefully ishft(1,-31) is a compile-time constant. Otherwise,
substitute the appropriate constant.)
> return (Q[j]=t);
> }
Q(j)=t
MWC=t
end
> #define CNG ( xcng=69069*xcng+123 )
This can be done with shift and add if multiply doesn't give
the appropriate result.
> #define XS ( xs^=(xs<<13), xs^=(xs>>17), xs^=(xs<<5) )
> #define KISS ( MWC()+CNG+XS ) /*138 million/sec*/
-- glen
christian.bau
MWC () is supposed to implement multiplication of a huge number by
8193, with the huge number split into 32 bit items.
We could implement this by calculating (Q[i] * 8193 + carry) with 64
bit arithmetic, storing 32 bits, moving the upper 32 bits into carry.
In 32 bit arithmetic, we calculate the 32 bits to be stored as (Q[i] +
(Q[i]<<13) + carry), that part is fine.
The new carry should be (Q[i] >> 19), plus the carry in adding (Q[i] +
(Q[i]<<13) + carry) in 32 bit arithmetic.
The carry is found by adding all three items, and checking whether the
sum is less than Q [i].
That is incorrect when (Q[i]<<13) + carry produces a carry in 32 bit
arithmetic.
And that happens when we start with carry = 0x800 and Q [i] =
0xffffffff.
I have no idea how this would influence the randomness, but it changes
the random number generator from something that can be examined
mathematically to something that is quite random.
Gib Bogle
I guess you didn't test this. I did, and it gives results that are different
from those given by the C code that George posted. The workaround that I
posted, in which the comparison of t and x is done laboriously:
if ((t >= 0 .and. x >= 0) .or. (t < 0 .and. x < 0)) then
if (t < x) then
tLTx = 1
else
tLTx = 0
endif
elseif (x < 0) then
tLTx = 1
elseif (t < 0) then
tLTx = 0
endif
c = tLTx + SHIFTR(x,19)
gives results identical to the C code. Can you supply Fortran code that
reproduces George's numbers?
glen herrmannsfeldt
> glen herrmannsfeldt wrote:
>> In comp.lang.fortran Gib Bogle <g.b...@auckland.no.spam.ac.nz> wrote:
>>> It's possible to predict, but more reliable to test.
>>> I can use Fortran to reproduce the results generated by the
>>> C code, but not using the pseudo-code that George provided.
>>> Unfortunately he seems to have lost interest in this thread.
>>> Perhaps you'd like to suggest the Fortran code to implement
>>> the method that George provided in C (assuming twos complement).
> I guess you didn't test this. I did, and it gives results that are different
> from those given by the C code that George posted. The workaround that I
> posted, in which the comparison of t and x is done laboriously:
> gives results identical to the C code. Can you supply Fortran
> code that reproduces George's numbers?
Trying it with gcc, I don't get the same answers as George.
The following Fortran, modulo 2**32, gives the right answers:
The C code has been left in as comments. The before and after
comments removed, but one could add them back in from the original
post.
As mentioned previously, this assumes 32 bit, twos complement, and
that on overflow the proper low order bits are returned.
I haven't done any timing, or time comparisons to C.
integer function MWC()
! when did zero argument functions get added to Fortran?
integer q
common q(0:4690)
integer t,x,c,j
save c
data j/4691/
! unsigned long t,x,i; static c=0,j=4691;
j=j+1
if(j>4690) j=0
! j=(j<4690)? j+1:0;
x=q(j)
! x=Q[j];
t=ishft(x,13)+c+x
c=ishft(x,-19)
if(ieor(t,ishft(1,31)).lt.ieor(x,ishft(1,31))) c=c+1
! t=(x<<13)+c+x; c=(t<x)+(x>>19);
q(j)=t
MWC=t
return
end
! return (Q[j]=t);
! int main()
program main
integer xs,xcng
data xs,xcng/521288629,362436069/
! static unsigned long xs=521288629,xcng=362436069,Q[4691];
integer i,x
integer q
common q(0:4690)
! {unsigned long i,x;
do i=0,4690
xcng=69069*xcng+123
! #define CNG ( xcng=69069*xcng+123 )
xs=ieor(xs,ishft(xs, 13))
xs=ieor(xs,ishft(xs,-17))
xs=ieor(xs,ishft(xs, 5))
! #define XS ( xs^=(xs<<13), xs^=(xs>>17), xs^=(xs<<5) )
q(i)=xcng+xs
enddo
! for(i=0;i<4691;i++) Q[i]=CNG+XS;
do i=1,1000000000
x=MWC()
enddo
! for(i=0;i<1000000000;i++) x=MWC();
print *," MWC result=3740121002 (or -554846295) ?",x
! printf(" MWC result=3740121002 ?\n%22u\n",x);
! for(i=0;i<1000000000;i++) x=KISS;
do i=1,1000000000
xcng=69069*xcng+123
! #define CNG ( xcng=69069*xcng+123 )
xs=ieor(xs,ishft(xs, 13))
xs=ieor(xs,ishft(xs,-17))
xs=ieor(xs,ishft(xs, 5))
! #define XS ( xs^=(xs<<13), xs^=(xs>>17), xs^=(xs<<5) )
x=MWC()+xcng+xs
enddo
! #define KISS ( MWC()+CNG+XS ) /*138 million/sec*/
print *,"KISS result=2224631993 (or -2070335303) ?",x;
! printf("KISS result=2224631993 ?\n%22u\n",x);
! }
end
-- glen
Gib Bogle
(snipped Fortran code)
That does it! Thanks a lot.
geo
wrote:
I commend you on your analysis rather than concern over
programming style, but in the MWC method, with modulus b,(=2^32)
multiplier a<b, prime p=ab^n-1, (n=4196)
and 0 <= x < b and 0 <= c < a,
new x = a*x+c mod b, new c = integerpartof((a*x+c)/b)
and the new c can never reach, or exceed, a.
Formally, we may consider the set Z of a*b^4169 "vectors" z:
z=[c;x_1,x_2,...,x_{4169}], 0 <= c < a and 0 <= x_i < 2^32
and a function f(z): with t=a*x+c,
f([c;x_1,x_2,...,x_{4169}])=
[trunc(t/b);x_2,x_3,...,x_{4169},t mod b],
If k is the order of b mod p, the directed graph:z -> f(z)
has (p-1)/k loops of size k and two "loops" of size 1:
f(z)=z for z=[0;0,0,...,0] and z=[b-1,a-1,a-1,...,a-1].
The concatenated sequence made by taking the rightmost x
from each z forms, in reverse order,
the base-b expansion of j/p for some 0 <= j <= p.
When z=[0;0,0,...,0]
we get 0/p=.0000000...
When z=[a-1;d,d,d,...,d] with d=b-1,
we get p/p=.dddd...
In our case, we have four possible loops,
two of size (p-1)/2, two of size 1,
for a total of 2(p-1)/2 + 2 = p+1 = ab^4169,
related by the mappings
c <-> a-1-c, x <-> b-1-x, j/p <-> (p-j)/p.
The function f has an inverse: with s=c*b+x_{4169},
f^{-1}(z)=[s mod a;trunc(s/a),x_1,x_2,...,x_{4168}],
so that with g(z)=f^{-1}(z),
the trailing element in each of the vectors
g(z),g(g(z)),g(g(g(z))),...
forms, in normal order, the base-b expansion of some j/p.
But arithmetic in g(z) is not well-suited for computer
implementation, while that of f(z) is eminently so, the
reason I created it---the method, not the arithmetic.
By finding the prime p=ab^n-1 with
b=2^32, a=8193=2^13+1 and n=4169,
I as able to ensure that the necessary arithmetic
could be simply---and quickly--carried out with
the base b=2^32 arithmetic available to most of us.
In my original post I neglected to point out that the
initial choice of c should satisfy 0 <= c < a=8193
and that, choosing
c=0 and all the Q's zero, or
c=a-1 and all the Q's b-1
would makes the MWC generator have period 1
and the resulting KISS period a mere 2^64-2^32
rather than one exceeding 10^40000.
Please forgive me.
George Marsaglia
christian.bau
the problem is not that the carry could match or exceed the multiplier
8193.
The problem is that the carry can be as large as 8192.
In your code you write
t=(x<<13)+c+x; c=(t<x)+(x>>19);
To make it a bit more obvious what happens, I can change this to the
100% equivalent
t = (x << 13) + c;
c = x >> 19;
t += x; if (t < x) ++c;
The last line is the usual trick to add numbers and check for carry: t
will be less than x exactly if the addition t += x wrapped around and
produced a carry.
The problem is that (x << 13) + c can also produce a carry, and that
goes unnoticed. If c = 8192, and the lowest 19 bits in x are all set,
then (x << 13) has the highest 19 bits set and 13 bits zero. Adding c
"overflows" and gives a result of 0, and adding x leaves x. The
comparison (t < x) is false and produces zero, even though (x<<13)+c+x
should have produced a carry.
This is very, very rare. It doesn't actually happen in your original
source code with a total of 2 billion random numbers. If you produce
four billion random numbers, then it happens at i = 3596309491 and
again at i = 3931531774.
I am more interested in 64 bit performance, so I just made c, t, and x
64 bit numbers and wrote
uint64_t x = Q [i];
uint64_t t = (x << 13) + x + c;
c = (t >> 32);
That might help in Fortran as well, assuming 64 bit integers are
available, signed or unsigned doesn't matter.
Gib Bogle
> I haven't done any timing, or time comparisons to C.
On my system, using the Intel compilers (Windows), the C version executes in 9
sec, the Fortran in 28 sec.
baf
Odd. With gfortran 4.6, I get
MWC result=3740121002 (or -554846295) ? -554846294
KISS result=2224631993 (or -2070335303) ? -2070335303
Notice that the MWC result is off by 1
Gib Bogle
Presumably that's a typo. The signed value of 3740121002 is -554846294
according to me. Do this:
integer :: u = 3740121002
write(*,*) u
I'd be interested to hear how the C and Fortran times compare for you.
baf
> baf wrote:
>> On 8/17/2010 11:35 PM, Gib Bogle wrote:
>>> glen herrmannsfeldt wrote:
>>>
>>> (snipped Fortran code)
>>>
>>> That does it! Thanks a lot.
>>
>> Odd. With gfortran 4.6, I get
>>
>> MWC result=3740121002 (or -554846295) ? -554846294
>> KISS result=2224631993 (or -2070335303) ? -2070335303
>>
>> Notice that the MWC result is off by 1
>
> Presumably that's a typo. The signed value of 3740121002 is -554846294
Yes, that seems to be correct. Just a typo in the code.
> according to me. Do this:
>
> integer :: u = 3740121002
> write(*,*) u
>
> I'd be interested to hear how the C and Fortran times compare for you.
>
gcc -O2 27.7 sec
gfortran -O2 31.9
I am sure both of these times could be improved with some better choices
in compile options.
io_x
"christian.bau" <christ...@cbau.wanadoo.co.uk> ha scritto nel messaggio
news:9309cd00-acf9-4908...@x21g2000yqa.googlegroups.com...
> Hi George,
>
> the problem is not that the carry could match or exceed the multiplier
> 8193.
> The problem is that the carry can be as large as 8192.
>
> In your code you write
>
> t=(x<<13)+c+x; c=(t<x)+(x>>19);
>
> To make it a bit more obvious what happens, I can change this to the
> 100% equivalent
>
> t = (x << 13) + c;
> c = x >> 19;
> t += x; if (t < x) ++c;
>
> The last line is the usual trick to add numbers and check for carry: t
> will be less than x exactly if the addition t += x wrapped around and
> produced a carry.
>
> The problem is that (x << 13) + c can also produce a carry, and that
> goes unnoticed.
but if Compile it using RosAsm,
the result adding that carry is the same
---------------------------------------------
[mesA: D$ 0
mesB: D$ 0
StatC: D$ 0
StatJ: D$ 4691
xs: D$ 521288629
xcng: D$ 362436069
Q: D$ 0 #2346
secondQ: D$ 0 #2346 ; 2346*2=4692
; 0..4690
IMWCIIrisultatoI3740121002I: B$ "MWC risultato 3740121002?", 0
IKISSIrisultatoI2224631993I: B$ "KISS risultato 2224631993?", 0
msg: B$ "On each window press ok and wait for the next window", 0 , 0
]
[ArrForAll: D$ ? #2048
]
mwc:
push esi |push edi
mov eax, D$StatJ |cmp eax, 4690 |jae @1 |inc eax |jmp @2
@1: xor eax, eax
@2: mov D$StatJ , eax
mov edx, D$Q+eax*4 |mov ecx, D$StatC |mov edi, edx |xor esi, esi
shl edx, 13 |add edx, ecx |adc esi, 0
add edx, edi |adc esi, 0
shr edi, 19 |lea ecx, D$edi+esi
mov D$StatC , ecx
mov D$Q+eax*4 , edx
mov eax, edx
pop edi |pop esi
ret
CNGmXS:
mov eax, D$xcng |mov edx, 69069 |mul edx |add eax, 123 |mov D$xcng ,
eax
mov eax, D$xs |mov ecx, D$xs |shl eax, 13 |xor ecx, eax
mov eax, ecx |shr eax, 17 |xor ecx, eax
mov eax, ecx |shl eax, 5 |xor ecx, eax
mov D$xs , ecx |mov eax, D$xcng |add eax, ecx
ret
KISS:
push esi |push edi
mov eax, D$StatJ |cmp eax, 4690 |jae @1 |inc eax |jmp @2
@1: xor eax, eax
@2: mov D$StatJ , eax
mov edx, D$Q+eax*4 |mov ecx, D$StatC |mov edi, edx |xor esi, esi
shl edx, 13 |add edx, ecx |adc esi, 0
add edx, edi |adc esi, 0
shr edi, 19 |lea ecx, D$edi+esi
mov D$StatC , ecx
mov D$Q+eax*4 , edx
pop edi |pop esi
push edx
mov eax, D$xcng |mov edx, 69069 |mul edx |add eax, 123 |mov D$xcng ,
eax
mov eax, D$xs |mov ecx, D$xs |shl eax, 13 |xor ecx, eax
mov eax, ecx |shr eax, 17 |xor ecx, eax
mov eax, ecx |shl eax, 5 |xor ecx, eax
mov D$xs , ecx |mov eax, ecx |mov eax, D$xcng |add eax, ecx
pop edx
add eax, edx
ret
Main:
pushad
;iint3
push &MB_OK |push msg |push msg |push 0 |call 'user32.MessageBoxA'
mov esi, 0
@1: call CNGmXS |mov D$Q+esi*4 , eax |inc esi |cmp esi, 4691 |jb @1
;iint3
mov esi, 0
@2: call mwc |inc esi |cmp esi, 1000000000 |jb @2
mov D$mesA , eax
mov edx, esp |mov D$esp , "%u" |push eax |push edx |push ArrForAll
call 'user32.wsprintfA' |add esp, 12
push &MB_OK |push IMWCIIrisultatoI3740121002I |push ArrForAll |push 0
call 'user32.MessageBoxA'
;iint3
mov esi, 0
@3: call KISS |inc esi |cmp esi, 1000000000 |jb @3
mov D$mesB , eax
mov edx, esp |mov D$esp , "%u" |push eax |push edx |push ArrForAll
call 'user32.wsprintfA' |add esp, 12
push &MB_OK |push IKISSIrisultatoI2224631993I |push ArrForAll |push 0
call 'user32.MessageBoxA'
;iint3
popad
mov eax, 0
ret
------------------End RosAsm
-------bigin the language i use
PREPARSE MACRO2D
section DATA
%define MsgBoxA 'user32.MessageBoxA'
%define wsprintfA 'user32.wsprintfA'
mesA dd 0
mesB dd 0
StatC dd 0
StatJ dd 4691
xs dd 521288629
xcng dd 362436069
Q: times 2346 dd 0
secondQ times 2346 dd 0 ; 2346*2=4692
; 0..4690
"MWC risultato 3740121002?" db "MWC risultato 3740121002?" , 0
"KISS risultato 2224631993?" db "KISS risultato 2224631993?", 0
msg db "On each window press ok and wait for the next window", 0, 0
section BSS
ArrForAll: resd 2048
section TEXT
mwc:
<i,j
a=*StatJ|a<4690!#.1|++a|#.2
.1: a^=a
.2: *StatJ=a
r=*Q+a*4|c=*StatC|j=r|i^=i
r<<=13 |r+=c |i++=0
r+=j |i++=0
j>>=19 |c=&*j+i
*StatC=c
*Q+a*4=r
a=r
>i,j
ret
CNGmXS:
a=*xcng|r=69069|mul r|a+=123|*xcng=a
a=*xs|c=*xs|a<<=13|c^=a
a=c |a>>=17|c^=a
a=c |a<<=5 |c^=a
*xs=c|a=*xcng|a+=c
ret
KISS:
<i,j
a=*StatJ|a<4690!#.1|++a|#.2
.1: a^=a
.2: *StatJ=a
r=*Q+a*4|c=*StatC|j=r|i^=i
r<<=13 |r+=c |i++=0
r+=j |i++=0
j>>=19 |c=&*j+i
*StatC=c
*Q+a*4=r
>i,j
<r
a=*xcng|r=69069|mul r|a+=123|*xcng=a
a=*xs|c=*xs|a<<=13|c^=a
a=c |a>>=17|c^=a
a=c |a<<=5 |c^=a
*xs=c|a=c |a=*xcng|a+=c
>r
a+=r
ret
Main:
pushad
;iint3
MsgBoxA(0, msg, msg, &MB_OK )
i=0
.1: CNGmXS()|*Q+i*4=a|++i|i<4691#.1
;iint3
i=0
.2: mwc() |++i|i< 1000000000#.2
*mesA=a
r=s|D*s="%u"|wsprintfA<(ArrForAll, r, a)
MsgBoxA(0, ArrForAll, "MWC risultato 3740121002?" , &MB_OK )
;iint3
i=0
.3: KISS()|++i|i< 1000000000#.3
*mesB=a
r=s|D*s="%u"|wsprintfA<(ArrForAll, r, a)
MsgBoxA(0, ArrForAll, "KISS risultato 2224631993?", &MB_OK )
;iint3
popad
a=0
ret
-------end the language i use
christian.bau
> <code snipped>
I have no idea what that code is doing, but the situation where a
carry is missed is very rare, and the first time the result is
different is after more than 3 billion iterations. Not visible in the
1 billion you did.
geo
wrote:
Thanks for pointing that out.
The correction (t<x) should be (t<=x),
to take care of that rare, but inevitable, event
when the last 19 bits of x are 1's
and the carry c is a-1=2^13, making (x<<13)+c = 0,
and thus t=(x<<13)+c+x = x,
(mod 2^32, of course, the assumed underlying
integer arithmetic of the processor).
The entire MWC() routine should have that
(t<x) to (t<=x) correction:
unsigned long MWC(void)
{ unsigned long t,x,i; static c=0,j=4691;
j=(j<4690)? j+1:0;
x=Q[j];
t=(x<<13)+c+x; c=(t<=x)+(x>>19);
return (Q[j]=t);
}
and for interested readers who may have pasted
the original, please change that (t<x) to (t<=x).
The test values from 10^9 calls will still be as indicated,
but strict adherence to the underlying mathematics requires
the change from (t<x) to (t<=x) to ensure the full period.
Keeping that (t<x) may be viewed as, rather than making a
circular transition through the 8193*2^133407-1 base-b digits
of the expansion of some j/p, we jump---after an average of
b=2^32 steps---to another point in the immense circle.
It could be that the period is increased rather than
decreased, but it remains a curiosity. Perhaps light could
be shed with primes p=ab^n-1 for which the actual distribution
of such jumps, averaging every b steps, could be examined.
Or perhaps it could remain a mystery
and be further fodder for those who tend to
equate lack of understanding to "true" randomness.
George Marsaglia