Implementing the KISS4691 RNG

G

geo

For those interested in the KISS4691 RNG, I would
like to recommend a final form for the MWC
(Multiply-With-Carry) component, one that provides
a pattern for applications with signed integers,
takes care of the nuisance cases where (x<<13)+c
overflows, and, with different choice of multiplier
and prime but with the same mathematics, permits
going through a full cycle to verify the period.

First, a simple version for confirming the period:
Here we still have b=2^32, but multiplier a=5 and
prime p=ab-1, for which the order of b mod p is
(p-1)/2 = 5*2^31-1 = 10737418239

Note that, as with the full MWC() function below, we
deliberately seek multipliers 'a' of the form 2^i+1
to ensure the ability to carry out the MWC operations
in mod 2^32 arithmetic. This required a search for primes
of the form (2^i+1)*b^r-1, and I as able to find
the prime (2^13+1)*b^4691-1 for the KISS4691 version.

But the simpler version with p=(2^2+1)b-1 permits going
through a full cycle (using only 32-bit operations)
in less that a minute. Here we have a function

f([x,c])=[(5x+c) mod 2^32,trunc((5x+c)/2^32)]

with inverse

f^{-1}([x,c])=[trunc((c*2^32+x)/5),(c*2^32+x) mod 5]

Initial pairs [x,c]=[0,0] or [2^32-1,4] will make the
sequence f([x,c]),f(f([x,c])),f(f(f([x,c]))),... have
period 1. All others, with 0<=x<2^32, 0<=c<5, will have
period (p-1)/2=10737418239 and the (hexed) x's will form,
in reverse order, the hex expansion of j/p for some 0<j<p.

Running this C program:
-----------------------------------------------
#include <stdio.h>
int main()
{ unsigned long x=123456789,c=3,t;
unsigned long long k;
for(k=0;k<10737418239LL;k++)
{t=(x<<2)+c;
if(t<c) { c=(x>>30)+1; x=t+x;}
else {c=(x>>30); x=t+x; c+=(x<t);}
}
printf("%u,%u\n",x,c);
}
-----------------------------------------------
will go through a full cycle and, after some
35 seconds, reproduce the initial x,c pair:
123456789,3
Try it.

For the KISS4691 version, the MWC multiplier is a=2^13+1,
and I had a mental image of c taking 13 bits with the
rightmost 13 bits of (x<<13) all 0's, so that (x<<13)+c
could not overflow---indeed, that was the reason I sought
multipliers of the form a=2^i+1. But I overlooked the
rare case where c could occupy 14 bits: c=8193, and cause
an overflow when the rightmost 19 bits of x were all 1's.
Astute readers pointed that out. The above version, in
which we first test for overflow after t=(x<<2)+c, then,
if necessary, test with t=t+x, covers all cases and has the
advantage of permitting similar structures for programs
using signed integers---if we use a clever device advocated
by Glen Herrmannsfeldt when providing a Fortran version of
KISS4691 (comp.lang.fortran Aug 18, 2010).

With t=x+y for integers, the C version of the overflow for
signed integers, (t<x), may be replaced by (t'<x') for signed
integers, where the prime means: flip the sign bit. Thus,
in C, with m=(1<<31), the overflow is ((t^m)<(x^m)), while
for Fortran, with m=ishft(1,31), the logic is
ieor(t,m) .lt. ieor(x,m).


Thus I suggest using this listing for the MWC function
in KISS4691:
------------------------------------------------

unsigned long MWC(void)
{unsigned long t,x; static c=0,j=4691;
j=(j<4690)? j+1:0;
x=Q[j]; t=(x<<13)+c;
if(t<c) {c=(x>>19)+1; t=t+x;}
else {t=t+x; c=(x>>19)+(t<x);}
return (Q[j]=t);
}

------------------------------------------------

and recommend that it be the pattern for implementations in
PL/I, asm, Fortran or others. For signed integers, flip the
sign bits so that t<x becomes (t^m)<(x^m) or equivalents.
(C versions using signed integers should also avoid the
problem of right shifts, replacing (x>>19) by
((unsigned)x>>19) or ((unsigned)xs>>17) in the XS component.)


In case you do not have, or don't want to fetch, the original,
here is the entire listing, with the recommended form for MWC():

------------------------------------------------------------
static unsigned long xs=521288629,xcng=362436069,Q[4691];

unsigned long MWC(void)
{unsigned long t,x; static c=0,j=4691;
j=(j<4690)? j+1:0;
x=Q[j]; t=(x<<13)+c;
if(t<c) {c=(x>>19)+1; t=t+x;}
else {t=t+x; c=(x>>19)+(t<x);}
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 )

#include <stdio.h>
int main()
{unsigned long i,x;
for(i=0;i<4691;i++) Q=CNG+XS;
printf("Does MWC result=3740121002 ?\n");
for(i=0;i<1000000000;i++) x=MWC();
printf("%27u\n",x);
printf("Does KISS result=2224631993 ?\n");
for(i=0;i<1000000000;i++) x=KISS;
printf("%27u\n",x);
}
------------------------------------------------------------

Unfortunately, even at the rate of 238 million per second
for MWC(), it would take over 10^40407 years to verify the
period by running through a complete loop.

George Marsaglia
 
I

Ian Collins

On 09/ 6/10 09:04 AM, geo wrote:

Thus I suggest using this listing for the MWC function
in KISS4691:

static what?
j=(j<4690)? j+1:0;
x=Q[j]; t=(x<<13)+c;
if(t<c) {c=(x>>19)+1; t=t+x;}
else {t=t+x; c=(x>>19)+(t<x);}
return (Q[j]=t);
}

This code must be making assumptions about sizes, it gives different
results in 32 and 64 bit mode.
------------------------------------------------

and recommend that it be the pattern for implementations in
PL/I, asm, Fortran or others. For signed integers, flip the
sign bits so that t<x becomes (t^m)<(x^m) or equivalents.
(C versions using signed integers should also avoid the
problem of right shifts, replacing (x>>19) by
((unsigned)x>>19) or ((unsigned)xs>>17) in the XS component.)


In case you do not have, or don't want to fetch, the original,
here is the entire listing, with the recommended form for MWC():

------------------------------------------------------------
static unsigned long xs=521288629,xcng=362436069,Q[4691];

unsigned long MWC(void)
{unsigned long t,x; static c=0,j=4691;
j=(j<4690)? j+1:0;
x=Q[j]; t=(x<<13)+c;
if(t<c) {c=(x>>19)+1; t=t+x;}
else {t=t+x; c=(x>>19)+(t<x);}
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 )

#include<stdio.h>
int main()
{unsigned long i,x;
for(i=0;i<4691;i++) Q=CNG+XS;
printf("Does MWC result=3740121002 ?\n");
for(i=0;i<1000000000;i++) x=MWC();
printf("%27u\n",x);
printf("Does KISS result=2224631993 ?\n");
for(i=0;i<1000000000;i++) x=KISS;
printf("%27u\n",x);
}


Try compiling with warnings enabled!
 
D

David Bernier

geo wrote:
[...]
In case you do not have, or don't want to fetch, the original,
here is the entire listing, with the recommended form for MWC():

------------------------------------------------------------
static unsigned long xs=521288629,xcng=362436069,Q[4691];

unsigned long MWC(void)
{unsigned long t,x; static c=0,j=4691;
j=(j<4690)? j+1:0;
x=Q[j]; t=(x<<13)+c;
if(t<c) {c=(x>>19)+1; t=t+x;}
else {t=t+x; c=(x>>19)+(t<x);}
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 )

#include<stdio.h>
int main()
{unsigned long i,x;
for(i=0;i<4691;i++) Q=CNG+XS;
printf("Does MWC result=3740121002 ?\n");
for(i=0;i<1000000000;i++) x=MWC();
printf("%27u\n",x);
printf("Does KISS result=2224631993 ?\n");
for(i=0;i<1000000000;i++) x=KISS;
printf("%27u\n",x);
}
------------------------------------------------------------

Unfortunately, even at the rate of 238 million per second
for MWC(), it would take over 10^40407 years to verify the
period by running through a complete loop.

[...]

I thought my compiler was using 32-bit
unsigned longs, but it turns out it uses
64-bit unsigned longs ...

I realized by adding lines to print the size of
various data types that with the compiler and
system I had, "unsigned long"s were 8 bytes
and "unsigned int"s were 4 bytes.

So by using "unsigned int"s, it's now
working.

David Bernier
----------

$ time ./a.out
size of unsigned long is: 8
size of unsigned int is: 4
size of char is: 1

Does MWC result=3740121002 ?
3740121002
Does KISS result=2224631993 ?
2224631993

real 0m32.828s // 33 seconds

-------- C code for above: ----------

#include <stdio.h>

static unsigned int xs=521288629,xcng=362436069,Q[4691];

unsigned int MWC(void)
{unsigned int t,x; static unsigned int c=0,j=4691;
j=(j<4690)? j+1:0;
x=Q[j]; t=(x<<13)+c;
if(t<c) {c=(x>>19)+1; t=t+x;}
else {t=t+x; c=(x>>19)+(t<x);}
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 )


int main()
{unsigned int i,x;

printf("size of unsigned long is: %d\n", sizeof(unsigned long));
printf("size of unsigned int is: %d\n", sizeof(unsigned int));
printf("size of char is: %d\n\n", sizeof(char));

for(i=0;i<4691;i++) Q=CNG+XS;
printf("Does MWC result=3740121002 ?\n");
for(i=0;i<1000000000;i++) x=MWC();
printf("%27u\n",x);
printf("Does KISS result=2224631993 ?\n");
for(i=0;i<1000000000;i++) x=KISS;
printf("%27u\n",x);
}
 
R

robin

| On 09/ 6/10 09:04 AM, geo wrote:
|
| <snip>
|
| >
| > Thus I suggest using this listing for the MWC function
| > in KISS4691:
| > ------------------------------------------------
| >
| > unsigned long MWC(void)
| > {unsigned long t,x; static c=0,j=4691;
|
| static what?
|
| > j=(j<4690)? j+1:0;
| > x=Q[j]; t=(x<<13)+c;
| > if(t<c) {c=(x>>19)+1; t=t+x;}
| > else {t=t+x; c=(x>>19)+(t<x);}
| > return (Q[j]=t);
| > }
|
| This code must be making assumptions about sizes,

It's for 32 bit words.
As previously, for longer words, it's necessary to truncate.
 
R

robin

| On 09/ 6/10 09:22 PM, robin wrote:

| > It's for 32 bit words.
| > As previously, for longer words, it's necessary to truncate.
|
| Then it shouldn't make assumptions about the sizes of types.

It doesn't.
Find the phrase "in mod 2^32 arithmetic" in George's posting.
 
R

robin

| > It's for 32 bit words.
| > As previously, for longer words, it's necessary to truncate.
|
| Then it shouldn't make assumptions about the sizes of types. If fixed
| width types are required, use them.

You would need to truncate to 32 bits immediately after each right shift.
 
I

Ian Collins

|> It's for 32 bit words.
|> As previously, for longer words, it's necessary to truncate.
|
| Then it shouldn't make assumptions about the sizes of types. If fixed
| width types are required, use them.

You would need to truncate to 32 bits immediately after each right shift.

As I said, use fixed size types. Otherwise the code won't work in 64
bit builds.

#include <stdint.h>

uint32_t MWC(void)
{
uint32_t t,x;
static uint32_t c=0,j=4691;

j=(j<4690)? j+1:0;
x=Q[j];
t=(x<<13)+c;

if(t<c) {
c=(x>>19)+1;
t=t+x;
}
else {
t=t+x;
c=(x>>19)+(t<x);
}
return Q[j]=t;
}
 
G

geo

On 09/ 6/10 09:04 AM, geo wrote:

<snip>


Thus I suggest using this listing for the MWC function
in KISS4691:
------------------------------------------------
  unsigned long MWC(void)
  {unsigned long t,x; static c=0,j=4691;

static what?
    j=(j<4690)? j+1:0;
    x=Q[j]; t=(x<<13)+c;
      if(t<c) {c=(x>>19)+1; t=t+x;}
      else {t=t+x; c=(x>>19)+(t<x);}
    return (Q[j]=t);
  }

This code must be making assumptions about sizes, it gives different
results in 32 and 64 bit mode.




Try compiling with warnings enabled!

Evidently I had in mind: here is a list of unsigned longs,

unsigned long t,x; static c=0,j=4691;

among which I want c and j to be saved after each call.

Wrong thinking, but why did the compiler let me get away with it?

Sorry for that; I'll first try full warnings from now on.

George Marsaglia
 
M

Mark Bluemel

geo said:
Evidently I had in mind: here is a list of unsigned longs,

unsigned long t,x; static c=0,j=4691;

among which I want c and j to be saved after each call.

Wrong thinking, but why did the compiler let me get away with it?

implicit int...
 
R

robin

|
| For those interested in the KISS4691 RNG, I would
| like to recommend a final form for the MWC
| (Multiply-With-Carry) component, one that provides
| a pattern for applications with signed integers,
| takes care of the nuisance cases where (x<<13)+c
| overflows, and, with different choice of multiplier
| and prime but with the same mathematics, permits
| going through a full cycle to verify the period.
|
| First, a simple version for confirming the period:
| Here we still have b=2^32, but multiplier a=5 and
| prime p=ab-1, for which the order of b mod p is
| (p-1)/2 = 5*2^31-1 = 10737418239
|
| Note that, as with the full MWC() function below, we
| deliberately seek multipliers 'a' of the form 2^i+1
| to ensure the ability to carry out the MWC operations
| in mod 2^32 arithmetic. This required a search for primes
| of the form (2^i+1)*b^r-1, and I as able to find
| the prime (2^13+1)*b^4691-1 for the KISS4691 version.
|
| But the simpler version with p=(2^2+1)b-1 permits going
| through a full cycle (using only 32-bit operations)
| in less that a minute. Here we have a function
|
| f([x,c])=[(5x+c) mod 2^32,trunc((5x+c)/2^32)]
|
| with inverse
|
| f^{-1}([x,c])=[trunc((c*2^32+x)/5),(c*2^32+x) mod 5]
|
| Initial pairs [x,c]=[0,0] or [2^32-1,4] will make the
| sequence f([x,c]),f(f([x,c])),f(f(f([x,c]))),... have
| period 1. All others, with 0<=x<2^32, 0<=c<5, will have
| period (p-1)/2=10737418239 and the (hexed) x's will form,
| in reverse order, the hex expansion of j/p for some 0<j<p.
|
| Running this C program:
| -----------------------------------------------
| #include <stdio.h>
| int main()
| { unsigned long x=123456789,c=3,t;
| unsigned long long k;
| for(k=0;k<10737418239LL;k++)
| {t=(x<<2)+c;
| if(t<c) { c=(x>>30)+1; x=t+x;}
| else {c=(x>>30); x=t+x; c+=(x<t);}
| }
| printf("%u,%u\n",x,c);
| }
| -----------------------------------------------
| will go through a full cycle and, after some
| 35 seconds, reproduce the initial x,c pair:
| 123456789,3
| Try it.
|
| For the KISS4691 version, the MWC multiplier is a=2^13+1,
| and I had a mental image of c taking 13 bits with the
| rightmost 13 bits of (x<<13) all 0's, so that (x<<13)+c
| could not overflow---indeed, that was the reason I sought
| multipliers of the form a=2^i+1. But I overlooked the
| rare case where c could occupy 14 bits: c=8193, and cause
| an overflow when the rightmost 19 bits of x were all 1's.
| Astute readers pointed that out. The above version, in
| which we first test for overflow after t=(x<<2)+c, then,
| if necessary, test with t=t+x, covers all cases and has the
| advantage of permitting similar structures for programs
| using signed integers---if we use a clever device advocated
| by Glen Herrmannsfeldt when providing a Fortran version of
| KISS4691 (comp.lang.fortran Aug 18, 2010).
|
| With t=x+y for integers, the C version of the overflow for
| signed integers, (t<x), may be replaced by (t'<x') for signed
| integers, where the prime means: flip the sign bit. Thus,
| in C, with m=(1<<31), the overflow is ((t^m)<(x^m)), while
| for Fortran, with m=ishft(1,31), the logic is
| ieor(t,m) .lt. ieor(x,m).
|
|
| Thus I suggest using this listing for the MWC function
| in KISS4691:
| ------------------------------------------------
|
| unsigned long MWC(void)
| {unsigned long t,x; static c=0,j=4691;
| j=(j<4690)? j+1:0;
| x=Q[j]; t=(x<<13)+c;
| if(t<c) {c=(x>>19)+1; t=t+x;}
| else {t=t+x; c=(x>>19)+(t<x);}
| return (Q[j]=t);
| }
|
| ------------------------------------------------
|
| and recommend that it be the pattern for implementations in
| PL/I, asm, Fortran or others. For signed integers, flip the
| sign bits so that t<x becomes (t^m)<(x^m) or equivalents.
| (C versions using signed integers should also avoid the
| problem of right shifts, replacing (x>>19) by
| ((unsigned)x>>19) or ((unsigned)xs>>17) in the XS component.)
|
|
| In case you do not have, or don't want to fetch, the original,
| here is the entire listing, with the recommended form for MWC():
|
| ------------------------------------------------------------
| static unsigned long xs=521288629,xcng=362436069,Q[4691];
|
| unsigned long MWC(void)
| {unsigned long t,x; static c=0,j=4691;
| j=(j<4690)? j+1:0;
| x=Q[j]; t=(x<<13)+c;
| if(t<c) {c=(x>>19)+1; t=t+x;}
| else {t=t+x; c=(x>>19)+(t<x);}
| 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 )
|
| #include <stdio.h>
| int main()
| {unsigned long i,x;
| for(i=0;i<4691;i++) Q=CNG+XS;
| printf("Does MWC result=3740121002 ?\n");
| for(i=0;i<1000000000;i++) x=MWC();
| printf("%27u\n",x);
| printf("Does KISS result=2224631993 ?\n");
| for(i=0;i<1000000000;i++) x=KISS;
| printf("%27u\n",x);
| }
| ------------------------------------------------------------
|
| Unfortunately, even at the rate of 238 million per second
| for MWC(), it would take over 10^40407 years to verify the
| period by running through a complete loop.
|
| George Marsaglia


Here is the PL/I version, using 32-bit unsigned arithmetic :-

(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);
declare (t,x,i) fixed binary (32) unsigned;
declare (c initial (0), j initial (4691) ) fixed binary (32) unsigned static;

if j < hbound(Q,1) then j = j + 1; else j = 0;
x = Q(j); t = isll(x,13) + c;
if t<c then
do; c = isrl(x, 19) + 1; t = t + x; end;
else
do; t = t + x; c = isrl(x, 19) + (t<x); end;
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 );
end KISS;

declare (i,x) fixed binary (32) unsigned;

Q = CNG + XXS; /* Initialize: */
do i = 1 to 1000000000; x=MWC(); end;
put skip edit (" MWC result=3740121002",x) (a, skip, f(22));
do i = 1 to 1000000000; x=KISS; end;
put skip edit ("KISS result=2224631993",x) (a, skip, f(22));

end RNG;
 
K

Keith Thompson

geo said:
Evidently I had in mind: here is a list of unsigned longs,

unsigned long t,x; static c=0,j=4691;

among which I want c and j to be saved after each call.

The above is exactly equivalent to:

unsigned long t,x;
static c=0,j=4691;

There is no significance to the fact that they're on the same line.
Wrong thinking, but why did the compiler let me get away with it?

In C90 (but not in C99), a declaration with no explicit type defaults to
"int". This rule is widely considered to have been a bad idea (it goes
back to C's ancestor languages, B and/or BCPL), which is why it was
removed in C99. So in C90, those lines are equivalent to:

unsigned long t,x;
static int c=0,j=4691;

Or, in better style (IMHO):

unsigned long t;
unsigned long x;
static int c = 0;
static int j = 4691;

C lets you declare multiple objects in a single declaration, but it's a
rich source of errors. For one thing, you have to remember which things
are shared across everything in a declaration and which aren't:

int x, y = 0; /* y is initialized to 0, x is not */
static int x, y; /* both x and y are static */
int* x, y; /* x is a pointer, y is an int */
int *x, y; /* Again, x is a pointer, y is an int */

In the last, the spacing makes it clearer that the "*" applies to x,
but it's still much clearer to use two separate declarations:

int *x;
int y;

Due to C's declaration-follows-use syntax, the declaration "int *x;"
means "*x is of type int", which implies that x is of type int*.

[...]
 
G

glen herrmannsfeldt

(snip)
There is no significance to the fact that they're on the same line.
(snip)
C lets you declare multiple objects in a single declaration, but it's a
rich source of errors.

and Fortran lets you declare one object in multiple declarations:

INTEGER X
DIMENSION X(10)
POINTER X
TARGET X
INTENT(IN) X
PROTECTED X
SAVE X
VALUE X
VOLATILE X
ASYNCHRONOUS X

I am not so sure that you can use all those on a single variable,
but many of them.
For one thing, you have to remember which things
are shared across everything in a declaration and which aren't:
int x, y = 0; /* y is initialized to 0, x is not */

INTEGER X,Y
DATA X,Y/1,2/

initializes X to 1, and Y to 2.

Which reminds me of an IBM feature that never got into
standard Fortran. In IBM Fortran IV:

INTEGER X,Y/1/

initializes Y to 1. You can initialize more than one variable with
a list of values (in //) in DATA, but not in explicit type statements.

Many compilers might allow this for IBM compatability.

-- glen
 
R

Richard Maine

glen herrmannsfeldt said:
and Fortran lets you declare one object in multiple declarations:

which is also a source or errors. For example, it is common (no pun
intended, but...) advice to use include files for COMMON in f77 so as to
make sure that the COMMON declaration is the same in all scopes where it
is used. But the above-mentioned "feature" can defeat that and make the
commons disagree in spite of using include.

An include file might have something like

common /com/ x, y
real x, y(10)

with explicit type declrations carefully used. But there is no way that
it can ensure that one scoping unit doesn't have something like

dimension x(4)

buried away somewhere that is easy to overlook - heck maybe even in some
other include file. Some linkers might catch the resulting common size
mismatch, but there is no guarantee of that, and that's just an
after-the-fact catch. The "feature" is inherently vulnerable to errors.

Modules aren't subject to that same kind of error (well, ignoring the
fact that you can have a common block in a module) because the language
standard explicitly disallows USEing something from a module and adding
on attributes to it outside of the module.
 
R

robin

|
| For signed integers, flip the
| sign bits so that t<x becomes (t^m)<(x^m) or equivalents.
| (C versions using signed integers should also avoid the
| problem of right shifts, replacing (x>>19) by
| ((unsigned)x>>19) or ((unsigned)xs>>17) in the XS component.)

In the case of signed integers, both t<c AND t<x
must have their sign bits flipped, not just t<x.

Thus, the signed version in PL/I is:

(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 m fixed binary (31) static initial ((isll(1,31)));

if j < hbound(Q,1) then j = j + 1; else j = 0;
x = Q(j); t = isll(x,13) + c;
if ieor(t,m) < ieor(c,m) then
do; c = isrl(x, 19) + 1; t = t + x; end;
else
do; t = t + x; c = isrl(x, 19) + (ieor(t,m) < ieor(x,m)); end;
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 (10);

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;
 
U

Uno

geo said:
> Evidently I had in mind: here is a list of unsigned longs,
>
> unsigned long t,x; static c=0,j=4691;
>
> among which I want c and j to be saved after each call.
>
> Wrong thinking, but why did the compiler let me get away with it?
>
> Sorry for that; I'll first try full warnings from now on.

George,

The C mistakes in this source have been pointed out time and again and
seem to fall into categories:

1) You use a more traditional C where things like implicit int were
tolerated. Dann Corbitt published a version of your source with all
these shortcomings fixed, but you seem not to have noticed. Part of
posting is at least skimming *all* germane responses. If you can't read
these, then post less.

2) For your program to be useful to me, it needs a random seed, maybe
something like this for the *nix platform:

$ gcc -Wall -Wextra geo3.c -o out
$ ./out
198327871
$ ./out
2716291704
$ cat geo3.c
#include <stdio.h>

int main(void)
{
FILE *f = fopen("/dev/random", "rb");
if (f != NULL)
{
unsigned int data;
if (1 == fread(&data, sizeof data, 1, f))
{
printf("%u\n", data);
}
fclose(f);
}
return 0;
}

// gcc -Wall -Wextra geo3.c -o out
$

When you turn up warnings, these are, in my opinion, good suggestions.

3) My strongest criticism of this source is the lack of explanation for
hard-coded numbers. Why is j 4691 and not 4692? A comment for each of
these values goes a long way.

As it is, I don't know how to use your code effectively, because I can't
see what needed to start as something random.

Most people can hardly write. You, however, seem to be like Richard
Heathfield in that writing comes more easily than reading. When you
don't incorporate people's germane criticisms, people notice that you're
not learning simple, fixable things, and their attitude changes. I've
seen these same old implicit ints for your last several threads as OP.

Cheers,
 
D

Dann Corbit

For those interested in the KISS4691 RNG, I would
like to recommend a final form for the MWC
(Multiply-With-Carry) component, one that provides
a pattern for applications with signed integers,
takes care of the nuisance cases where (x<<13)+c
overflows, and, with different choice of multiplier
and prime but with the same mathematics, permits
going through a full cycle to verify the period.

First, a simple version for confirming the period:
Here we still have b=2^32, but multiplier a=5 and
prime p=ab-1, for which the order of b mod p is
(p-1)/2 = 5*2^31-1 = 10737418239

Note that, as with the full MWC() function below, we
deliberately seek multipliers 'a' of the form 2^i+1
to ensure the ability to carry out the MWC operations
in mod 2^32 arithmetic. This required a search for primes
of the form (2^i+1)*b^r-1, and I as able to find
the prime (2^13+1)*b^4691-1 for the KISS4691 version.

But the simpler version with p=(2^2+1)b-1 permits going
through a full cycle (using only 32-bit operations)
in less that a minute. Here we have a function

f([x,c])=[(5x+c) mod 2^32,trunc((5x+c)/2^32)]

with inverse

f^{-1}([x,c])=[trunc((c*2^32+x)/5),(c*2^32+x) mod 5]

Initial pairs [x,c]=[0,0] or [2^32-1,4] will make the
sequence f([x,c]),f(f([x,c])),f(f(f([x,c]))),... have
period 1. All others, with 0<=x<2^32, 0<=c<5, will have
period (p-1)/2=10737418239 and the (hexed) x's will form,
in reverse order, the hex expansion of j/p for some 0<j<p.

For many compilers {specifically those for which unsigned long is
exactly 32 bits}, this version will produce acceptable answers. For C99
compilers, we can use specifically sized types that will enable it to
run portably on any machine with a C99 compiler (no hen's teeth jokes,
please).

static unsigned long xs = 521288629;
static unsigned long xcng = 362436069;
static unsigned long Q[4691];

unsigned long MWC(void)
{
static unsigned long c = 0;
static unsigned long j = 4691;
unsigned long t;
unsigned long x;
j = (j < 4690) ? j + 1 : 0;
x = Q[j];
t = (x << 13) + c;
if (t < c) {
c = (x >> 19) + 1;
t = t + x;
} else {
t = t + x;
c = (x >> 19) + (t < x);
}
return (Q[j] = t);
}

void initMWC(void)
{
unsigned long i;
for (i = 0; i < sizeof Q / sizeof Q[0]; i++)
Q = (xcng = 69069 * xcng + 123) + (xs ^= (xs << 13),
xs ^= (xs >> 17), xs ^= (xs << 5));
}

#ifdef UNIT_TEST

#include <stdio.h>

int main()
{
unsigned long i;
unsigned long x;

initMWC();

printf("Does MWC result=3740121002 ?\n");
for (i = 0; i < 1000000000; i++)
x = MWC();
printf("%27u\n", x);

printf("Does KISS result=2224631993 ?\n");
for (i = 0; i < 1000000000; i++)
x = (MWC() + (xcng = 69069 * xcng + 123) + (xs ^= (xs << 13),
xs ^= (xs >> 17), xs ^= (xs << 5)));
printf("%27u\n", x);
return 0;
}
#endif
/*
c:\tmp>cl /W4 /Ox /DUNIT_TEST kiss7.c
Microsoft (R) C/C++ Optimizing Compiler Version 15.00.30729.01 for x64
Copyright (C) Microsoft Corporation. All rights reserved.

kiss7.c
c:\tmp\kiss7.c(46) : warning C4701: potentially uninitialized local
variable 'x' used
Microsoft (R) Incremental Linker Version 9.00.30729.01
Copyright (C) Microsoft Corporation. All rights reserved.

/out:kiss7.exe
kiss7.obj

c:\tmp>kiss7
Does MWC result=3740121002 ?
3740121002
Does KISS result=2224631993 ?
2224631993

P.S.
Is C4701 up above the dumbest warning you ever saw in your life?
*/
 
K

Keith Thompson

Dann Corbit said:
c:\tmp\kiss7.c(46) : warning C4701: potentially uninitialized local
variable 'x' used [...]
P.S.
Is C4701 up above the dumbest warning you ever saw in your life?
*/

Um, no. It merely means that the compiler wasn't able to prove that x
is always initialized before it's used. You can probably silence it by
providing an initializer on the declaration.

For example, if I wrote something like this:

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

int main(void) {
int x;
srand(time(NULL));
if (rand() % 2 == 0) {
puts("initializing x");
x = 42;
}
printf("x = %d\n", x);
return 0;
}

I'd certainly like my compiler to warn me about it.

(Interestingly, gcc doesn't warn with "-Wall -Wextra", but at -O1
and higher x is 42 whether the message is printed or not. Apparently
the compiler recognized the undefined behavior but didn't bother to
tell me about it.)

Why do you find the warning dumb?

Followups redirected to comp.lang.c.
 
D

Dann Corbit

Dann Corbit said:
c:\tmp\kiss7.c(46) : warning C4701: potentially uninitialized local
variable 'x' used [...]
P.S.
Is C4701 up above the dumbest warning you ever saw in your life?
*/

Um, no. It merely means that the compiler wasn't able to prove that x
is always initialized before it's used. You can probably silence it by
providing an initializer on the declaration.

For example, if I wrote something like this:

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

int main(void) {
int x;
srand(time(NULL));
if (rand() % 2 == 0) {
puts("initializing x");
x = 42;
}
printf("x = %d\n", x);
return 0;
}

I'd certainly like my compiler to warn me about it.

(Interestingly, gcc doesn't warn with "-Wall -Wextra", but at -O1
and higher x is 42 whether the message is printed or not. Apparently
the compiler recognized the undefined behavior but didn't bother to
tell me about it.)

Why do you find the warning dumb?

Because the loops absolutely, positively must execute at least once.
Followups redirected to comp.lang.c.

My attempt at a C++ version failed. I am not sure what is wrong with
the second method.


class MWC {

private:
unsigned long xs;
unsigned long xcng;
unsigned long Q[4691];
unsigned long t;
unsigned long x;
unsigned long c;
unsigned long j;

public:

// pseudo-random number method:
unsigned long kissRandom(void) {
j = (j < 4690) ? j + 1 : 0;
x = Q[j];
t = (x << 13) + c;
if (t < c) {
c = (x >> 19) + 1;
t = t + x;
} else {
t = t + x;
c = (x >> 19) + (t < x);
}
return (Q[j] = t);
}

// another flavor, with a bit of mixing:
unsigned long kissRandom2(void) {
return (kissRandom() + (xcng = 69069 * xcng + 123) +
(xs ^= (xs << 13), xs ^= (xs >> 17),
xs ^= (xs << 5)));
}

// Constructor:
MWC(void) {
unsigned long i;
c = 0;
j = 4691;
xs = 521288629;
long xcng = 362436069;
for (i = 0; i < sizeof Q / sizeof Q[0]; i++)
Q = (xcng = 69069 * xcng + 123) +
(xs ^= (xs << 13), xs ^= (xs >> 17),
xs ^= (xs << 5));
}
};

#ifdef UNIT_TEST

#include <iostream>

int main()
{
unsigned long i;
unsigned long x;
MWC prng;

std::cout << "Does MWC result=3740121002 ? " << std::endl;
for (i = 0; i < 1000000000; i++)
x = prng.kissRandom();
std::cout << x << std::endl;

std::cout << "Does KISS result=2224631993 ?" << std::endl;
for (i = 0; i < 1000000000; i++)
x = prng.kissRandom2();
std::cout << x << std::endl;
return 0;
}
#endif
/*
c:\tmp>cl /EHsc /DUNIT_TEST /W4 /Ox kiss7.cpp
Microsoft (R) C/C++ Optimizing Compiler Version 15.00.30729.01 for x64
Copyright (C) Microsoft Corporation. All rights reserved.

kiss7.cpp
c:\tmp\kiss7.cpp(60) : warning C4701: potentially uninitialized local
variable 'x' used
Microsoft (R) Incremental Linker Version 9.00.30729.01
Copyright (C) Microsoft Corporation. All rights reserved.

/out:kiss7.exe
kiss7.obj

c:\tmp>kiss7
Does MWC result=3740121002 ?
3740121002
Does KISS result=2224631993 ?
885060915
*/
 
G

glen herrmannsfeldt

In comp.lang.fortran Keith Thompson said:
Dann Corbit said:
c:\tmp\kiss7.c(46) : warning C4701: potentially uninitialized local
variable 'x' used [...]
P.S.
Is C4701 up above the dumbest warning you ever saw in your life?
*/
Um, no. It merely means that the compiler wasn't able to prove that x
is always initialized before it's used. You can probably silence it by
providing an initializer on the declaration.

(snip)

If you program Java for a while, you get used to it. Java compilers
are supposed to treat it as an error, not a warning, if the
compiler can't determine that a scalar variable is initialized.
(No such test for arrays.)

But, yes, I always think it is dumb when I know it is initialized,
but the compiler doesn't. With the Java exception model, though,
it isn't always so obvious. If a called method throws an exception,
it might be that the variable isn't assigned a value, even if it
is in an assignment statement. Or consider:

unsigned i,x;
for(i=-1;i<3;i++) x=3

If you try hard enough, you can always find a way to fool the
compiler initialization test logic.

The solution is to put an initializer on the declaration.
I sometimes add a comment about the dumb compiler as the
reason for the initializer.

-- glen
 

Ask a Question

Want to reply to this thread or ask your own question?

You'll need to choose a username for the site, which only take a couple of moments. After that, you can post your question and our members will help you out.

Ask a Question

Members online

No members online now.

Forum statistics

Threads
473,818
Messages
2,569,732
Members
45,691
Latest member
Dick331194

Latest Threads

Top