* Re: RNGs with periods exceeding 10^(40million).
[not found] <603ebe15-a32f-4fbb-ba44-6c73f7919a33@t35g2000yqj.googlegroups.com>
@ 2011-01-17 16:10 ` robin
2011-01-18 10:09 ` robin
1 sibling, 0 replies; 2+ messages in thread
From: robin @ 2011-01-17 16:10 UTC (permalink / raw)
/*** 32-bit version in PL/I ***/
(nosize):
rng: procedure options (main);
declare Q(0:4194303) unsigned fixed binary (32),
carry unsigned fixed binary (32) static initial (0);
b32MWC: procedure () returns (unsigned fixed binary (32));
declare (t,x) unsigned fixed binary (32);
declare j fixed binary (31) static initial (4194303);
j=iand(j+1, 4194303);
x=Q(j); t=isll(x, 28)+carry;
carry=isrl(x, 4)-(t<x);
Q(j)=t-x;
return (Q(j));
end b32MWC;
KISS: procedure() returns (unsigned fixed binary (32));
declare t unsigned fixed binary (32);
t = b32MWC();
cng = bin(69069)*cng+(13579);
xs = ieor(xs, isll(xs, 13));
xs = ieor(xs, isrl(xs, 17));
xs = ieor(xs, isll(xs, 5));
return (t + cng + xs);
end KISS;
declare (i,x,cng static initial (123456789),
xs static initial (362436069)) unsigned fixed binary (32);
declare (start_time, finish_time) float (18);
start_time = secs();
/* First seed Q with CNG+XS: */
do i=0 to 4194304-1;
xs = ieor(xs, isll(xs, 13));
xs = ieor(xs, isrl(xs, 17));
xs = ieor(xs, isll(xs, 5));
cng = bin(69069)*cng+bin(13579);
Q(i)=cng+xs;
end;
/* Then generate 10^9 b32MWC()s */
do i=1 to 1000000000; x=b32MWC(); end;
if decimal(x) ^= 2769813733 then
put skip list ('b32MWC: x is not the expected value 2769813733', ' x=', x);
/* followed by 10^9 KISSes: */
do i=1 to 1000000000; x=KISS(); end;
if decimal(x)^=3545999299 then
put skip list ('KISS: x is not the expected value 3545999299', ' x=', x);
Finish_time = secs();
put skip edit (Finish_time - Start_time, ' secs') (F(10,2), A);
end rng;
______________
"geo" <gmarsaglia@gmail.com> wrote in message news:603ebe15-a32f-4fbb-ba44-6c73f7919a33@t35g2000yqj.googlegroups.com...
|
| I offer here a MWC (Multiply-With-Carry) RNG whose
| period I cannot determine---nor is it likely that
| anyone else can---but that unknown period is almost
| certainly greater than 10^40000000, i.e, 10^(40million).
| Versions for either 32- and 64-bit integers are given,
| as well as suggestions for similar versions with
| periods as great or even greater magnitude.
|
| MWC RNGs are ordinarily based on primes of the form
| p=ab^r-1, with b=2^32 or 2^64 and a<b. They generate,
| 32- or 64-bits at a time, and in reverse order, the
| binary expansions of rationals j/p, 0<j<p and the
| period of that binary expansion is the order of 2 mod p.
| The key part of the MWC process is: given 32-bit a,x,c,
| extract the top and bottom 32 bits of a 64 bit a*x+c,
| or, for given 64-bit a,x,c, extract the top and bottom
| 64 bits of a 128-bit a*x+c.
|
| I have put various computers to work establishing the
| order of 2 for MWC primes of the form ab^r-1, b=2^32,
| or CMWC (Complimentary-Multiply-With-Carry), p=ab^r+1.
| For example, it took 24 days on an early laptop to
| establish that the order of 2 for p=54767*2^1337287+1
| is (p-1)/2^3, and just last week I put a Samsung 64-bit
| laptop to finding the order of 2 for the 11th largest
| known prime, p=27653*2^9167433+1; in CMWC form,
| p=14158336*b^286482+1. I expect it may take months.
|
| Such searches are motivated by more than a Mt. Everest
| syndrome, as extremely-long-period RNGs not only seem to
| perform very well on tests of randomness, but have uses
| for cryptography. For example, suppose, from a random
| set of seeds, you have observed 100000 successive 32-bit
| numbers from a CMWC RNG based on p=14158336*b^286482+1.
| You will then be somewhere in an immense loop that has
| over 10^(2million) appearances of exactly that same
| sequence of 100000, leaving you virtually no chance of
| finding your particular location and thus able to
| predict forthcoming values.
|
| Of course, if you can observe 286482 successive values,
| then, after a little modular arithmetic to determine
| the carry, you are OK. It is ensuring a large exponent
| for b that makes the RNG desirable---the larger the better.
|
| But in seeking large exponents for b and record periods,
| why bother with an Everest when there are ranges with
| peaks far beyond surveyor's skills but certain
| to be several times as tall?
|
| Such ranges come from considering, for very large r,
| composites of the form m=ab^r-1, rather than primes of
| that form. The example for this posting: the composite
| m=(2^28-1)*2^(2^27)-1=(2^28-1)*b^(2^22)-1.
| m has no prime divisors less than 21000000, and order(2,m)
| is the lcm of the orders mod the prime divisors of m.
| We cannot---and may never be able to--factor such an
| m>10^40403570, but for primes q, the average value of
| order(2,q)/q is around .572. Thus, even if m has, say,
| 40 prime factors, since (.572)^40>10^(-10), we can only
| expect to "lose" around ten 10's, and perhaps ten more
| through common factors for the lcm, providing an estimate
| around 10^40403550 for the order of 2 mod m.
|
| That the order of 2 mod m is less than 10^(40million)
| seems extremely unlikely; we would have to "lose" over
| 400thousand 10's, when something in the range
| of three to sixty 10's seems more likely.
|
| So users should feel confident that the periods of the
| following MWC RNGs far exceed 10^(40million). Based on
| m=(2^28-1)b^(2^22)-1=(2^28-1)b^4194304-1 with b=2^32,
| m=(2^28-1)B^(2^21)-1=(2^28-1)B^2091752-1 with B=2^64,
| they are simple and very fast.
| You would have to observe more than 134 million
| successive bits produced be these MWC RNGs
| in order to get close to cracking them.
|
| Note that choosing a=2^i-1 in m=ab^r-1 makes finding
| the top and bottom 32 bits of a purported 64 bit
| a*x+c feasible using only 32-bit arithmetic,
| or, within 64-bit arithmetic, finding the top
| and bottom 64 bits of a purported 128 bit a*x+c.
|
| Here are C versions using, for 32-bit integers,
| an unsigned long array Q[41943104], and for 64-bits,
| an unsigned long long array Q[2091752].
| Try them and see if you get the same results I do.
|
| --------------------32-bit MWC version-------------------
|
| #include <stdio.h>
| static unsigned long Q[4194304],carry=0;
|
| unsigned long b32MWC(void)
| {unsigned long t,x; static int j=4194303;
| j=(j+1)&4194303;
| x=Q[j]; t=(x<<28)+carry;
| carry=(x>>4)-(t<x);
| return (Q[j]=t-x);
| }
|
| #define CNG ( cng=69069*cng+13579 )
| #define XS ( xs^=(xs<<13), xs^=(xs>>17), xs^=(xs<<5) )
| #define KISS ( b32MWC()+CNG+XS )
|
| int main(void)
| {unsigned long int i,x,cng=123456789,xs=362436069;
| /* First seed Q[] with CNG+XS: */
| for(i=0;i<4194304;i++) Q[i]=CNG+XS;
| /* Then generate 10^9 b32MWC()s */
| for(i=0;i<1000000000;i++) x=b32MWC();
| printf("Does x=2769813733 ?\n x=%lu\n",x);
| /* followed by 10^9 KISSes: */
| for(i=0;i<1000000000;i++) x=KISS;
| printf("Does x=3545999299 ?\n x=%lu\n",x);
| return 0;
| }
|
| ---------------- 64-bit MWC version ---------------------
|
| #include <stdio.h>
| static unsigned long long Q[2097152],carry=0;
|
| unsigned long long B64MWC(void)
| {unsigned long long t,x; static int j=2097151;
| j=(j+1)&2097151;
| x=Q[j]; t=(x<<28)+carry;
| carry=(x>>36)-(t<x);
| return (Q[j]=t-x);
| }
|
| #define CNG ( cng=6906969069LL*cng+13579 )
| #define XS ( xs^=(xs<<13), xs^=(xs>>17), xs^=(xs<<43) )
| #define KISS ( B64MWC()+CNG+XS )
|
| int main(void)
| {unsigned long long i,x,cng=123456789987654321LL,
| xs=362436069362436069LL;
| /* First seed Q[] with CNG+XS: */
| for(i=0;i<2097152;i++) Q[i]=CNG+XS;
| /* Then generate 10^9 B64MWC()s */
| for(i=0;i<1000000000;i++) x=B64MWC();
| printf("Does x=13596816608992115578 ?\n x=%llu\n",x);
| /* followed by 10^9 KISSes: */
| for(i=0;i<1000000000;i++) x=KISS;
| printf("Does x=5033346742750153761 ?\n x=%llu\n",x);
| return 0;
| }
|
| ---------------------------------------------------------
|
| Note that I advocate the KISS, (Keep-It-Simple-Stupid),
| principle. Even though the MWC RNGs perform very well on
| tests of randomness, combining with Congruential (CNG)
| and Xorshift (XS) probably provides extra insurance
| at the cost of a few simple instructions, (and making
| the resulting KISS even harder to crack).
|
| Also note that for unsigned integers such as in Fortran,
| the C instruction -(t<x) means: subtract 1 if t<s,
| and can be implemented for signed integers as t'<s',
| where the ' means: flip the sign bit.
|
| Full random seeding of the carry and the Q[] array
| calls for more than 17 megabytes of random bits.
| For ordinary randomness testing, filling, as above,
| the Q[] array with Congruential+Xorshift may serve, but
| for stronger crypto uses, many more than the random 64-
| or 128-bits that CNG and XS require are needed.
| Random seeding of the carry and the Q[] array amounts to
| producing, either 32- or 64-bits at a time, the reverse-
| order binary expansion of j/m for some 0<=j<=m. Indeed,
| for m=ab^r-1, given the carry and the Q[] array, that j is
|
| j=carry+a*sum(Q[i]*b^i),i=0..r-1.
|
| There are phi(m) such j's with period order(2,m).
| Since m has no prime factors less than 21000000, with
| probability close to 1 a random seeding of Q[] should
| produce a j that is relatively prime to m and thus
| produce the full period---unknown but almost certainly
| far in excess of 10^(40million). Even if random
| seeding of carry and Q[] provides a j with gcd(j,m)>1,
| the period is still likely to exceed 10^(40million).
| There are two exceptions:
| All Q[i]=0, carry=0 yields j=0 and the binary expansion
| 0/m=.00000000000000000...
| All Q[i]=b-1=2^32-1 and carry=a-1=(2^28-2) yields j=m
| and the binary expansion
| m/m=.11111111111111111...
|
| Almost any choice of m=(2^i-1)*2^(2^27)-1 is likely
| to provide an immense order for 2 mod m.
| Choices i=16,18,22,28 resulted in m's that have no
| factors<10^7. Memory seems to be the key restriction
| for the resulting Q[] array. With enough memory,
| one might seek i's for which m=(2^i-1)*2^(2^28)-1
| has no small factors, or even m=(2^i-1)*2^(2^29)-1,
| boosting the unknown and unknowable periods of
| MWC RNGs to ranges beyond 10^(40million).
|
| George Marsaglia
^ permalink raw reply [flat|nested] 2+ messages in thread* Re: RNGs with periods exceeding 10^(40million).
[not found] <603ebe15-a32f-4fbb-ba44-6c73f7919a33@t35g2000yqj.googlegroups.com>
2011-01-17 16:10 ` RNGs with periods exceeding 10^(40million) robin
@ 2011-01-18 10:09 ` robin
1 sibling, 0 replies; 2+ messages in thread
From: robin @ 2011-01-18 10:09 UTC (permalink / raw)
"geo" <gmarsaglia@gmail.com> wrote in message news:603ebe15-a32f-4fbb-ba44-6c73f7919a33@t35g2000yqj.googlegroups.com...
|
| I offer here a MWC (Multiply-With-Carry) RNG whose
| period I cannot determine---nor is it likely that
| anyone else can---but that unknown period is almost
| certainly greater than 10^40000000, i.e, 10^(40million).
Here is the 64-bit version using PL/I :--
/* 64-bit version in PL/I using unsigned arithmetic 18/1/2011 */
(nosize):
rng64: procedure options (main);
declare (Q(0:2097151), carry initial (0)) unsigned fixed binary (64)
static;
B64MWC: procedure () returns (unsigned fixed binary (64));
declare (t,x) unsigned fixed binary (64);
declare j fixed binary (31) static initial (2097151);
j=iand(j+1, 2097151);
x=Q(j); t=isll(x, 28)+carry;
carry=isrl(x, 36)-(t<x);
Q(j)=t-x;
return (Q(j));
end B64MWC;
KISS: procedure() returns (unsigned fixed binary (64));
declare t unsigned fixed binary (64);
t = B64MWC();
cng=unsigned(6906969069, 64)*cng+unsigned(13579);
xs = ieor(xs, isll(xs, 13));
xs = ieor(xs, isrl(xs, 17));
xs = ieor(xs, isll(xs, 43));
return (t + cng + xs);
end KISS;
declare (i,x,cng static initial (123456789987654321),
xs static initial (362436069362436069))
unsigned fixed binary (64);
declare (start_time, finish_time) float (18);
start_time = secs();
/* First seed Q with CNG+XS: */
do i=0 to 2097151;
xs = ieor(xs, isll(xs, 13));
xs = ieor(xs, isrl(xs, 17));
xs = ieor(xs, isll(xs, 43));
cng=unsigned(6906969069, 64)*cng+unsigned(13579);
Q(i)=CNG+XS;
end;
/* Then generate 10^9 B64MWC()s */
do i=1 to 1000000000; x=B64MWC(); end;
put skip list (x);
if decimal(x) ^= 13596816608992115578 then
put list ('is not the expected value 13596816608992115578');
/* followed by 10^9 KISSes: */
do i=1 to 1000000000; x=KISS(); end;
put skip list (x);
if decimal(x) ^= 5033346742750153761 then
put list ('is not the expected value 5033346742750153761');
Finish_time = secs();
put skip edit (Finish_time - Start_time, ' secs') (F(10,2), A);
end rng64;
^ permalink raw reply [flat|nested] 2+ messages in thread