Re: RNGs with periods exceeding 10^(40million).

Re: RNGs with periods exceeding 10^(40million).

[robin]

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

Re: RNGs with periods exceeding 10^(40million).

[robin]

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