KISS4691, a potentially top-ranked RNG.

[Gib Bogle]

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.

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]

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..

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]

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?

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]

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?

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]

....
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.

 

[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]

| orz wrote:
|
| > 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.

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]

| orz wrote:
|
| > 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.

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.

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

 

[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).

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

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]

(snip of some suggestions, including a sign error in a shift)

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:

(snip)

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=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]

glen herrmannsfeldt wrote:

(snipped Fortran code)

That does it! Thanks a lot.

 

[geo]

I think there is a problem in the C code.

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 * 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 +
(Q<<13) + carry), that part is fine.
The new carry should be (Q >> 19), plus the carry in adding (Q +
(Q<<13) + carry) in 32 bit arithmetic.

The carry is found by adding all three items, and checking whether the
sum is less than Q .
That is incorrect when (Q<<13) + carry produces a carry in 32 bit
arithmetic.
And that happens when we start with carry = 0x800 and Q =
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.

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

 

[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]

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

 

[Gib Bogle]

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
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]

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.

 

[geo]

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

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

 

[glen herrmannsfeldt]

On Aug 18, 6:08 pm, "christian.bau" <[email protected]>
wrote:

(snip)

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).

That was my first thought, too, but then I wasn't so convinced.
It seems to me that will be wrong in the case that (x<<13) and c
are both zero, unless that can never happen.

This problem comes up in the S/370 AL (Add Logical) instruction
used for doing multiple precision addition. You have to check
for carry both when adding and when adding a previous carry.

(snip)

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 you could be (very) unlucky and it comes back such that
the period is shorter. I don't know the math well enough to
know that.

Or perhaps it could remain a mystery
and be further fodder for those who tend to
equate lack of understanding to "true" randomness.

There is a section in Knuth's TAOCP on randomly generated
random number generators.

-- glen