# Patent

## Title
Method of generating random numbers

## Abstract

A method of obtaining uniform and independent random numbers is given (a) comprising two distinct odd primes p1, p2 that give mutually coprime integers q1=(p1−1)/2 and q2=(p2−1)/2 with different parity to form the modulus d=p1p2; (b) comprising primitive roots z1, z2 of primes p1, p2, respectively, giving congruence relations z≡zj mod(pj) for j=1, 2 that determine the multiplier z; and (c) comprising the initial value n coprime with d=p1p2. The method generates the coset sequence n<z>={r1=n, r2, r3, . . . } of period T=2q1q2 recursively by rj+1=zrj mod(d) for j=1, 2, . . . in the reduced residue class group Z*d, giving {v1=r1/d, v2=r2/d, . . . } for output.

## Background
BACKGROUND OF THE INVENTION
1. Field of the Invention
The present invention relates to methods of obtaining uniform and independent random numbers on computers. More specifically, said methods relate to the multiplicative congruential generator comprising

    
    
        1. a positive integer d called modulus,
        2. a positive integer z (0<z<d) called multiplier that is coprime with d, namely that shares no common prime factor with d, and
        3. a positive integer n (0<n<d) called initial value or seed coprime with d, generating a sequence of integers {r1, r2, r3, . . . } by recursive congruence relations

r1=n, rj+1≡zrj mod(d), 0<rj<d, j=1, 2, 3, . . . ,

and giving an output sequence {v1, v2, v3, . . . } of rational numbers in the interval (0,1) by

vj:=rj/d, j=1,2,3, . . . ,

In contrast to the prior art that takes the modulus d for an odd prime p or a power 2i for a large index i, the invention selects those d comprising two odd distinct primes as the most appropriate, based on the survey of conceivable candidates of composite d. With other key specifications for the choice of said odd primes and for due restrictions on the design of the multiplier z, the invention will render forces to equip any computers and any operating systems feasibly with generators of uniform and independent random numbers with long periods and tested statistical high quality without improbable symmetry restrictions.

    
    


2. Description of the Related Art
Random number sequences on computers are required to have reproducibility and transportability for debugging simulation programs. Said reproducibility is the nature that the identical sequence of random numbers may be generated by demands of users. Said transportability refers to the generation of one and the same sequence of random numbers on transplantation to any other computers or computer languages.
Since the vast amount of random numbers needed for simulations never allows their storage in memory, the sole possible means left is some computational procedure to generate random numbers sequentially. Said requirements of reproducibility and transportability restrict the procedure to the arithmetic of integers which by nature is devoid of truncation and round-off errors on any computers and in any computer languages. Hence a uniform random number sequence on computers should exclusively be realized by some arithmetic as a sequence {x1, x2, . . . } of integers in a range [0, z) for a large but finite integer z, and should be transformed as {u1=x1/z, u2=x2/z, . . . } to the output sequence of rational numbers in the interval [0, 1) by division at every use.
For arbitrary positive integers z and T, any sequence {x1, x2, . . . , xT} of zero or positive integers in [0, z) may be concatenated to form a base-z numeral sequence,

x=0·x1x2 . . . xTx1x2 . . . xT . . . =x1/z+x2/z2+x3/z3+ . . .

with period T. Being excluded of uninteresting cases that said integers {x1, x2, . . . , xT} are all zero or all z:=z−1, said x is a rational number in (0, 1) and has the expression

x=(x1zT-1+x2zT-2+ . . . +xT)/(zT−1)=n/d, 

where and hereafter n/d stands for an irreducible fraction with 0<n<d. Since d divides zT−1, d and z are coprime. Any length-T sequence {x1, x2, . . . , xT} of integers in [0, z), which may be an output of any random number generator on a computer, or a length-T sequence taken from natural or physical uniform random numbers discretized with the precision 1/z, thus has an obvious correspondence to a period of the numeral sequence to base z, representing an irreducible fraction n/d with z and n coprime to d.

The division process of the irreducible fraction n/d with 0<n<d to the base z is expressed in the identity and the inequality

rj+1=zrj−xjd, 0<rj<d, j=1, 2, 3, . . . ,

where {r1=n, r2, r3, . . . } are remainders of division that never vanish since n and z are both coprime with d. Upon division by zd these give a crucial estimate found by Nakazawa (2004),

0<rj/d−xj/z=rj/d−uj=(rj+1/d)(1/z)<1/z, uj=xj/z, j=1, 2, 3, . . . .

Thus, numbers {u1=x1/z, u2=x2/z, u3=x3/z, . . . } in [0, 1), irrespective of whether they are created by some computational methods or obtained from any physical random numbers, are approximated from above within a uniform error bound 1/z, by numbers of the sequence {v1=r1/d, v2=r2/d, v3=r3/d, . . . } formed by remainders of the base-z division process of an irreducible fraction n/d, if the denominator d is allowed to be any positive integer to which the base z is coprime.

    
    
        Nakazawa (2004)/H. Nakazawa: Coset representations of uniform and independent random number sequences I—Cosets associated with periodic sequences. (URL) http://www10.plala.or.jp/h-nkzw (04eprsq1.pdf/Jul. 16-Aug. 12, 2004).
    
    


The term rj in the sequence {r1=n, r2, r3, . . . } of remainders of base-z division n/d is congruent to nzj−1 modulo d for j=1, 2, . . . , with the value of rj taken in (0, d). If the modulus d is an odd prime or a large power 2i of 2, the sequence {v1=r1/d, v2=r2/d, . . . } is the multiplicative congruential random numbers of the prior art generated by the multiplier z, for the initial value (or seed) n and in the modulus d. We analyze the most general form of the modulus d characterized by its factorization into prime factors

d=p1i1p2i2 . . . psis, ij≧1, 1≧j≧s, 

in order to find and exploit the optimal technological setting.

Let d≧2 be any positive integer, and denote Z*d for the set of integers coprime to d. Equivalence or identity of elements of Z*d is defined by the congruence modulo d; constituents of Z*d may thus be represented by integers in the interval (0, d), which will be replaced with [0, d) for later conveniences. The set Z*d forms a group with respect to the multiplication modulo d, and is called the reduced residue class group modulo d. The order, or the number of distinct elements, of the group Z*d is defined to be Euler's function and denoted as φ(d). For a prime d=p, Z*p consists of {1, 2, . . . , p−1}, giving φ(p)=p−1. For a power pi (i≧2) of a prime p the group Z*pi may be taken to consist of integers in the interval [0, pi) with the exclusion of multiples of p,

0=0·p,1·p,2·p,3·p, . . . ,(pi-1−1)·p, 

totaling to pi-1 in number. These give the following fundamental formula for any prime p:

φ(pi)=pi−pi-1=pi(1−1/p), i≧1,

which is always even with the sole exception of φ(21)=1. For any z∈Z*d define the cyclic sequence <z>, i.e. the cyclic (sub)group generated by z denoted in sequential form

<z>:≡{1=z0,z1,z2, . . . } mod(d).

With all constituents taken in the interval [0, d), <z> is identical to the multiplicative congruential sequence generated by the multiplier z with the modulus d and the initial value n=1. For the initial value n∈Z*d the multiplicative congruential sequence is the coset sequence

n<z>≡{n,nz,nz2, . . . }≡{r1,r2,r3, . . . } mod(d),

again with all constituents taken in the interval (0, d). Periods of <z> and n<z> are one and the same, to be denoted T. Lagrange's theorem proves that T is a divisor of the order φ(d) of Z*d. If the group Z*d is itself cyclic and z is its generator, then Z*d is exhausted by the sequence <z> with the full period T=φ(d). Any coset n<z> is then identical with <z> as a periodic sequence except for the starting element n. The group Z*d is known to be cyclic only for following cases of the modulus d:

    
    
        1. d=2, 4; cases d=2i, i≧3 are excluded.
        2. d=pi with an odd prime p and any integer i=1, 2, . . . .
        3. d=2pi with an odd prime p and any integer i=1, 2, . . . .

Case 1 is not interesting and discarded. Case 3 will find its natural treatments in the decomposition of the modulus d into prime factors.

    
    


The group Z*d with composite modulus d formed by coprime factors acquires a significant structure regulated by Chinese remainder theorem. Fix an arbitrary integer s≧2 and take pairwise coprime integers d1, d2, . . . , ds. Define d:=d1d2 . . . ds, and

ej:=d/dj, 1≦j≦s. 

Since ej consists of dk other than dj, dj and ej are coprime. Euclidean algorithm gives integers Dj and Ej that yield the identity

GCD(dj,ej)=1=djDj+ejEj, 1≦j≦s. 

Define Uj:=ejEj determined by d1, d2, . . . , ds. Since ej consists of dk's other than dj, the above identity gives, with δjk for Kronecker's delta,

Uj≡δjk mod(dk), 1≦j, k≦s. 

Here and hereafter 1≦j, k≦s will be used for short of 1≦j≦k and 1≦k≦s. Take arbitrary integers {zj|1≦j≦s}. The set of congruence relations

z≡zj mod(dj), 1≦j≦s, 

have a solution z=z1U1+ . . . +zsUs. If z′ is another solution of the above congruence relations, then z−z′≡0 mod(dj) holds for 1≦j≦s, implying z−z′ is a common multiple of d1, d2, . . . , ds. Hence z−z′ is divisible by their least common multiple d, proving that the solution of simultaneous congruence relations is unique modulo d with the form

z≡z1U1+z2U2+ . . . +zsUs mod(d=d1d2 . . . ds)

for integers U1, U2, . . . , Us determined by d1, d2, . . . , ds. Since this is a congruence relation modulo d, integers U1, U2, . . . , Us may be replaced with any set of integers equivalent to them modulo d. This is a form of the Chinese remainder theorem that holds true for any number s=2, 3, . . . of coprime factors.

Let again d1, d2, . . . , ds be pairwise coprime positive integers for any s≧2, and set d=d1d2 . . . ds. Any element z∈Z*d is coprime to any of d1, d2, . . . , ds, so that there is an element zj∈Z*dj giving z≡zj mod(dj) for each j, 1≦j≦s. Chinese remainder theorem thus proves that this element z taken arbitrarily in Z*d has the expression

z≡z1U1+z2U2+ . . . +zsUs mod(d=d1d2 . . . ds),

or that any z∈Z*d corresponds to a unique set (z1, z2, . . . , zs) which is an element of the product set Z*d1×Z*d2× . . . ×Z*ds. Conversely, any set of integers

(z1,z2, . . . ,zs)∈Z*d1×Z*d2× . . . ×Z*d,

corresponds to simultaneous congruence relations

z≡zk∈Z*dk mod(dk), 1≦k≦s, 

which determine an element z∈Z*11d2 . . . ds in the form z≡z1U1+z2U2+ . . . +zsUs mod(d) uniquely modulo d=d1d2 . . . ds by Chinese remainder theorem. These establish the one to one correspondence between all elements of Z*d1d2. . . ds and all elements of the product set Z*d1×Z*d2× . . . ×Z*ds. Counting elements in two sets, the multiplicative formula of Euler's function for pairwise coprime integers d1, d2, . . . , ds is obtained:

φ(d1d2 . . . ds)=φ(d1)φ(d2) . . . φ(ds).

As already noted, a composite modulus d for the group Z*d must not include any power 2i (i≧2) of 2 in order for <z> and n<z> in Z*d to realize uncorrelated sequences of random numbers. See Nakazawa and Nakazawa (2008a,b) for this discovery; Detailed Description of the Invention also gives related inferences with FIGS. 7A and 7B. Therefore, the most general form of a modulus d relevant to uniform and independent random number problems is summarized in the decompositon to prime factors as

d=d0d1d2 . . . ds, s≧1, dk=pkik, ik≧1, 1≦k≦s. 

Here p1, p2, . . . , ps are distinct odd primes, and d0 is either 20=1 or 21=2. The reduced residue class group Z*2={1} is a trivial group and φ(2)=1. The multiplicative formula of Euler's function gives the order of the group Z*d for the above noted modulus d as

φ(d)=φ(p1i1p2i2 . . . psis)={p1i1-1(p1−1)}{p2i2-1(p2−1)} . . . {psis-1(ps−1)},

for all cases of our interest.

    
    
        Nakazawa and Nakazawa (2008a)/H. Nakazawa and N. Nakazawa: Designs of uniform and independent random numbers with long period and high precision—Control of the sequential geometry through product group structures and lattice configurations—. (URL) http://www10.plala.orjp/h-nkzw (03090403.pdf/Mar. 9-Apr. 2, 2008).
        Nakazawa and Nakazawa (2008b)/H. Nakazawa and N. Nakazawa: Designs of uniform and independent random numbers with long period and high precision—Control of the sequential geometry through product group structures and lattice configurations—(URL) http://www10.plala.or.jp/h-nkzw (3978erv.pdf/Mar. 9-Jul. 8, 2008).
    
    


Two or more odd prime factors in the composite modulus d introduce into the group Z*d new structures that are absent in moduluses formed by a single prime or by its power. Take for s≧2 the modulus d=d1d2 . . . ds with pairwise coprime positive integers d1, d2, . . . , ds, each of which is assumed to give a non-trivial group Z*dk for 1≦k≦s. Any elements z, n∈Z*d determine corresponding elements in Z*dk by

z≡zk mod(dk), n≡nk mod(dk), 1≦k≦s. 

These give simultaneous congruence relations nzj−1≡nkzkj−1 mod(dk) for 1≦k≦s, which have the solution by Chinese remainder theorem in the following form:

nzj−1≡n1z1j−1U1+ . . . +nszsj−1Us mod(d=d1d2 . . . ds), j=1,2,3, . . . ,

As noted, integers U1, U2, . . . , Us are determined by {dk|1≦k≦s} and independent of {nk, zk|1≦k≦s}. The sequences <z> or its coset n<z> in Z*d1d2 . . . ds is thus a shuffling, in the sense of random number problems, of component cyclic sequences <zk> or cosets nk<zk> for 1≦k≦s, respectively. Note that this shuffling, though misstated to be linear in Nakazawa and Nakazawa (2008a) based on the appearance of the r.h.s. of nzj noted above, is nonlinear in its essence. This is a congruence relation and the value of the l.h.s. jumps nonlinearly as the r.h.s. crosses any multiple of d=d1d2 . . . ds. See Nakazawa and Nakazawa (2008b). Features of the generator with composite modulus d=d1d2 . . . ds thus stem from those in component groups {Z*dk|1≦k≦s}. Take typically the period of the coset sequence n<z>, including the case n=1 of cyclic sequence <z>. The expression of nzj−1 noted above, which is unique modulo d, proves that the period T of any coset n<z> is realized if and only if all component cosets {nk<zk>|1≦k≦s} resume their respective initial values simultaneously. Let nk<zk> have the period Tk for 1≦k≦s. The period T (or the order) of <z> and n<z> is thus given by the least common multiple,

T=LCM(T1,T2, . . . ,Ts).

Resume the most general composite modulus d of our concern with the form

d=d0d1d2 . . . ds, d0=1 or 2, dk=pkik, ik≧1, 1≦k≦s, s≧2

where p1, p2, . . . , ps are distinct odd primes. We now specify the largest period attainable with this modulus. The trivial group Z*2={1} for d0=21 always has the period T0=1 and, together with the case d0=20=1, Z*d0 is irrelevant in problems of the period. Lagrange's theorem states that the period Tk of nk<zk> for any elements zk, nk∈Z*pkik is a divisor of

(the order of the group Z*pkik)=φ(pkik)=2pkik-1qk, qk:=(pk−1)/2, 1<k≦s; 

note that pk is odd and pk−1 is even. Therefore, the period T of the coset sequence n<z> in the group Z*d satisfies

T=LCM(T0,T1,T2, . . . , Ts)≦LCM(1,2q1p1i1-1,2q2p2i2-1, . . . ,2qspsis-1).

Since the group Z*pkik is cyclic with generating elements for 1≦k≦s, the summary of related facts given below in Theorem A is now obvious

    
    
        Theorem A Let the modulus d be given with s≧2 and with distinct odd primes p1, p2, . . . , ps as

 d:=d0d, d:=d1d2 . . . ds, dk=pkik, ik≧1, 1≦k≦s, 

where d0 is either 20 or 21 and gives φ( d)=φ(d). Define integers

qk:=(pk−1)/2≧1, 1≦k≦s, 

with φ(dk)=φ(pkik)=2qkpkik-1 for 1≦k≦s. For any elements z, n in the group Z* d the period T, which is common to the cyclic sequence <z> and its coset sequence n<z>, is a divisor of φ(d), and fulfills inequalities

T≦2p1i1-1p2i2-1 . . . psis-1LCM(q1,q2, . . . ,qs)≦φ(d)/2s-1.

The largest possible period is φ(d)/2s-1. It is realized when d0=20=1 and the following set of conditions (I) and (II0) are satisfied, or when d0=21=2 and the set of conditions (I) and (II1) are satisfied:

        (I) Distinct odd primes p1, p2, . . . , ps are chosen so as to give pairwise coprime integers q1, q2, . . . , qs.
        (II0) For d0=20=1 with d=d the multiplier z is defined by generating elements zk of the component group Z*dk=Z*pkik for 1≦k≦s as

z≡zk mod(dk), 1≦k≦s 

        in the form z≡z1U1+z2U2+ . . . +zsUs mod(d).
        (II1) For d0=21=2 with d=2d the multiplier z is likewise defined by generating elements zk of the component group Z*dk=Z*pkik for 1≦k≦s as

z≡1 mod(2), z≡zk mod(dk), 1≦k≦s 

        in the form z≡1·d+(d+1)(z1U1+z2U2+ . . . +zsUs)mod( d=2d).

For the multiplier z realizing the largest period, the group Z* d is divided precisely into 2s-1 portions by the cyclic sequence <z> and its cosets.

    
    


The shuffling arising with moduluses formed by two or more odd prime factors induces effects reflecting structures of component groups, as is the case with periods noted in Theorem A. However, some drastic changes may also arise with the shuffling. The plot of the present invention is to make essential use of such a mechanism, which has its root in an obvious property of cyclic sequences and their cosets in reduced residue class groups.

    
    
        Lemma B Let d>4 be any integer, and z and n in (0, d) be any elements of the group Z*d which is not assumed to be cyclic. Denote the coset sequence n<z> as,

n<z>={r1,r2, . . . }, <rj<d, r1=n, rj≡nzj−1 mod(d), j=1,2, . . . ,

If the cyclic sequence <z> includes −1≡d−1 mod(d) as an element, the order T of n<z>, the number of elements or the period, is even, and any halves of the coset sequence are completely correlated as follows with q=T/2,

C+:={r1,r2, . . . ,rq}≡{n,nz, . . . ,nzq−1} mod(d),

C−:={rq+1,rq+2, . . . ,r2q}≡{nzq,nzq+1, . . . ,nz2q−1}≡{−n,−nz, . . . ,−nzq−1}≡−C+mod(d).

        Corollary to Lemma B Let d>4 give a cyclic group Z*d with a generator z. Then an arbitrary integer n in Z*d gives a coset sequence n<z> sweeping all elements of Z*d, but its first and second halves C± are completely correlated. Therefore, n<z> may be used, if at all, only for the half period T/2=φ(d)/2 for a sequence of independent random numbers.

(Proof) Let z have the order or the period T giving zT≡1 mod(d). By the assumption of −1∈<z> there is an integer q in 1<q<T≦φ(d) that gives −1≡zq, z2q≡1 mod(d). Thus, 2q=T, T is even, and the rest of Lemma B is obvious. By −1≡d−1 mod(d) the group Z*d invariably includes −1. Though a cyclic sequence <z′> in Z*d may be a proper subgroup and −1 need not necessarily be included in it, the assumption of Corollary gives −1∈Z*d=<z>, proving that Lemma B holds true with <z>. ▪

Primitive root generators have frequently been warned that only small portions of periods should be used for random numbers. Irrespective of whether a group Z*d is cyclic or not, its arbitrary cyclic sequence <z> with −1=zq mod(d) should be warned similarly. We may proceed a step further to the original sequences of integers corresponding to cosets C±, namely sequences corresponding to irreducible fractions n/d and (d−n)/d=1−n/d, respectively. Numeral expressions of these fractions to base z,

n/d=0·x1x2 . . . xqxq+1xq+2 . . . x2q . . . ,

1−n/d=0·xq+1xq+2 . . . x2qx1x2 . . . xq . . . ,

reveal simple relations xj+xq+j=z−1= z for all j=1, 2, . . . , and show how the correlation between C+ and C− of the cyclic sequence <z> including −1 reflects the original structure in the sequence {x1, x2, . . . } of integers with improbable symmetry between its first and second halves.

    
    


Take for a while the simplest composite modulus d=p1p2 formed by two odd distinct primes p1 and p2. The largest order of a multiplier z, or the largest period of <z> on this modulus d=p1p2 was shown to be T=φ(d)/22−1=φ(d)/2, which is realized if primes p1, p2 are so chosen as to give coprime integers

q1:=(p1−1)/2, q2:=(p2−1)/2;

and if the multiplier z∈Z*p1p2 is constructed with primitive roots z1, z2 of moduluses p1, p2, respectively. Note that the coset sequence n1<z1> for any n1∈Z*p1 has the excellence that the sequence sweeps almost all numbers in [0, p1), realizing a nearly perfect uniformity in one period. The matter is the same with n2<z2>. However, there is a spell of Corollary to Lemma B that each of n1<z1> or n2<z2> in its solitary use should be restricted to halves of their respective periods, thus spoiling their noted excellences in uniformity. It is to be admired that the shuffling that works in the composite modulus d=p1p2,

nzj−1≡n1z1j−1U1+n2z2j−1U2 mod(d=p1p2), j=1, 2, . . . ,

breaks this spell on component sequences. First, they are allowed to repeat their periods any number of times, and problems of (repeated appearances of) their correlated first and second halves are dismissed. Second, the shuffling gives chances to construct the coset n<z> that is free from the similar spell. The following Theorem C, which is the heart of the present invention, clarifies this magical work of Chinese remainder theorem. This is a generalization of the result obtained in Nakazawa and Nakazawa (2008a,b).

    
    
        Theorem C Let the modulus d be given with s≧2 and with powers of distinct odd primes p1,{grave over (p)}2, . . . , ps as

 d=d0d, d=d1d2 . . . ds, dk=pkik, ik≧1, 1≦k≦s, 

where d0 is either 20 or 21. As in Theorem A define integers

qk:=(pk−1)/2, 1≦k≦s, 

assuming that primes p1, p2, . . . , ps are so given that q1, q2, . . . , qs are pairwise coprime. Take likewise a generator zk in the component group Z*dk for 1≦k≦s and construct a multiplier z∈Z*d with congruence relations z≡zk mod(dk) for 1≦k≦s to obtain

zj≡z1jU1+z2jU2+ . . . +zsjUs mod(d), j=1,2, . . . ,

realizing the largest period T=φ(d)/2s-1 for <z>.

    
    
    (I) If integers q1, q2, . . . , qs consist of one even and all other odd integers, the cyclic subgroup <z> does not include −1.
    (II) If the remainder of the dichotomy holds true, i.e. if q1, q2, . . . , qs are all odd, then the cyclic sequence <z> includes −1.

(Proof) Assume the case d0=1 with d=d. Take a component cyclic group Z*dk (1≦k≦s) with dk=pkik and its arbitrary generator zk. This group contains −1≡dk−1 mod(dk) and has the order φ(dk)=pkik−1(pk−1)=2qkpkik−1. Therefore, zkj≡−1 mod(dk) arises when and only when j is an odd multiple of tk:=φ(dk)/2=pkik−1 qk. Since zj≡1 mod(d) occurs for the first time in j>0 at j=T=φ(d)/2s-1, we have zq≡−1 mod(d) at q=T/2 with

q=T/2=φ(d)/2s={φ(d1)/2}{φ(d2)/2} . . . {φ(ds)/2}=t1t2 . . . ts.

There is also an expression

zq≡z1qU1+z2qU2+ . . . +zsqUs mod(d),

which means that zq≡−1 mod(d) can arise only when zkq≡−1 mod(dk) holds for all 1≦k≦s. Therefore, −1∈<z> is true if all of factors

tk=φ(dk)/2=pkik-1 qk, 1≦k≦s, 

are odd, and false if one or more of them are even. Since all of primes p1, p2, . . . , ps are odd, the former corresponds to the case (II). The latter can arises only as the case (I), because q1, q2, . . . qs are restricted to be pairwise coprime. As regards the case d0=2 with d=2d, the largest period multiplier takes the form

zj≡1·d+(d+1){z1jU1+z2jU2+ . . . +zsjUs} mod( d=2d);

here zk, Uk for 1≦k≦s are those given in the preceding case d0=1. Note that d+1 is even, and zj≡1≡−1 mod(2) holds for any j. Thus zj≡−1 mod( d=2d) realizes if and only if some j gives zkj≡−1 mod(dk) for every 1≦k≦s, which was proved above to be the case in (II), and not the case in (I). ▪

The conclusion of Theorem C reveals that an appropriate choice of the composite modulus d, namely
    
        among pairwise coprime integers q1, q2, . . . , qs one is even and all others are odd,

breaks the spell that a largest period multiplier z is forced to generate a sequence <z> whose every pair of halves are stringently tied by symmetry. It will be notable that this is accomplished magically by using component sequences <zk> under the same spell constructed on a generator zk of the group Z*pkik. Further implications of this significant comprehension will be seen below in intuitive, geometrical terms.

    
    


Spectral tests have been the most powerful tool to assess the performance of a multiplier z of multiplicative congruential generators, typically with the modulus d formed by an odd prime d=p, or by a power of 2, d=2i. They remain to be an indispensable and powerful weapon to test multipliers of composite modulus generators, hence to test any random number sequences of finite length for their uniformity and independence, if they are computable at all. The inclusion of composite modulus cases, however, necessitates many reflections on structures to grasp correctly what has so far been passed unnoticed. The modulus in this paragraph is any positive integer d>4. The multiplier z and the initial value n are any elements in the group Z*d. The group Z*d need not be cyclic, nor z need be a generator. Let l=2, 3, . . . be an arbitrary integer, and let the coset sequence n<z> have the period T. Denote

n<z>≡{rj≡nzj−1 mod(d)|0<rj<d, rj+T=rj, j=1,2,3, . . . },

which includes the case n=1 of cyclic (subgroup) sequence <z>. We take sets of l consecutive numbers of n<z> as T points forming a sequence

{Pj(l)=(rj,rj+1, . . . ,rj+l-1)|j=1,2, . . . ,T}, 

in the hyper cube Cd in the 1-dimensional Euclidean space El with sides [0, d). As regards the geometry of these points, there are some crucial comprehensions that are indispensable in proceeding to full applications of spectral tests. We list their essence below.

    
    
        1. The points {Pj(l)|j=1, 2, . . . } taken from n<z> in El are in a lattice (that is, occupy a portion, not necessarily all, of lattice points of a lattice) determined by d, z and l. The lattice used here agrees in definition with the one for single odd prime modulus cases in the prior art; see Detailed Description of the Invention for the definition of the lattice and its dual.
        2. Spectral tests of our concern examine the configuration of lattice points in the hypercube Cd of El, not the distribution of points {Pj(l)|j=1, 2, . . . } taking their places in lattice points or seats.
        3. The total number of seats in the hypercube Cd is d, irrespective of the choice of z∈Z*d or of the dimension l=2, 3, . . . .

In more precise terms, we need to take all points in El that have coordinates congruent modulo d to points {Pj(l)|j=1, 2, . . . , T} noted above. This is equivalent to regard El as an l-dimensional torus with period d in all coordinates. See Nakazawa and Nakazawa (2008a,b) for plain but unabridged expositions including proofs and visual information. The total number of elements in the group Z*d is φ(d)<d. And the total number (or the period) T of elements in the coset n<z> is a divisor of φ(d) by Lagrange's theorem and T≦φ(d) is always the case. Therefore, points in El, generated by the l-tuples of <z> or n<z> in a period T leave always vacant or unoccupied lattice points in the hypercube Cd. These complications cannot be examined by spectral tests. Whatever d and z∈Z*d may be, the lattice is a definite existence, and a spectral test in El computes the largest distance λd(l)(z) between neighboring (l−1) dimensional hyperplanes that contain l−1 lattice vectors with linear independence. If z, z′ are multipliers for the same modulus d and λd(l)(z)<λd(l)(z′) holds true, the configuration of seats in El arising from z is valuated better than that given by z′. See below for more intuitive meanings of this criterion. The valuation of spectral tests has distances from the desired information on the distribution of points formed by <z> or n<z> in El. It should be stressed that the excellent performance in the spectral test is an indispensable qualification for any usable multiplier z, but not sufficient. We should be prepared for the existence of blind spots of spectral tests, and should proceed with their valuations carefully, helped by other information such as group structures of cyclic sequences or by various technical know-how that is to be accumulated in practice.

    
    


In order to have a clear overview on spectral tests, resume the simplest case of an odd prime modulus d=p and its primitive root multiplier z. The number of lattice points is d=p, while the period T of <z> is the order of Z*p given by T=φ(p)=p−1. Points in El,

{Pj(l)|j=1,2, . . . ,T}, 

generated by l-tuples of the sequence <z> for any dimension l=2, 3, . . . , therefore, fill the lattice points almost fully, leaving only the origin vacant. Thus, the spectral test of a primitive root multiplier of an odd prime modulus gives an almost precise valuation of the geometrical distribution of points generated by one period of <z> in El. FIGS. 2 to 4 depict visible cases of l=2 in the plane E2; the modulus is d=p=251; the multipliers are primitive roots z=162, 61 and 127, respectively; periodic sequences <z> and n<z> differ only in the choice of starting values and correspond to the same plot of T=250 points. The right of each plot shows points with normalized coordinates

{(rj/d,rj+1/d)|j=1,2, . . . ,T}

within a square that indicates the range −0.05≦x, y≦1.05 in E2 or xy-plane. Left plots show points with appropriately normalized coordinates

{(xj/z,xj+1/z)|j=1,2, . . . ,T}

corresponding to quotients that arise in the division process of the irreducible fraction n/d=1/251 to bases z=162, 61 and 127, respectively for FIGS. 2, 3 and 4. See that closeness of point configurations on the left and the right of these figures are in the order of FIG. 2, FIG. 4 and FIG. 3, in proportion to respective precision 1/z=1/162, 1/127 and 1/61. All of these figures contain the same number T=p−1=250 of generated points, or d=251 lattice points including the origin. There are wide differences in distributions of these points. If the largest distance between neighboring lattice lines (i.e. neighboring lines including the same l−1=1 linearly independent, or nonzero lattice vector in the present cases in E2) is larger, then the square can contain lesser number of those lattice lines, so that the density of points on the relevant lattice lines increases in inverse proportion. The lattice with the smallest value of its largest distance between neighboring lattice lines has the smallest density of points on lattice lines, realizing the most even or unbiased distribution of points. This will be seen clearly by comparing FIGS. 2-4. Our intuition that FIG. 4 is the worst and FIG. 2 is the best, as the distribution of points representing 2-tuples of numbers appearing uniformly and independently in sequence, is thus made quantitative in the spectral test as the largest distance λd(2)(z) of lattice lines for respective multiplier z; the smaller is λd(2)(z), the better is the multiplier z. Among all lattices in 2-dimension which are periodic with period 1 along coordinate axes and with d lattice points in the unit square, the triangular lattice realizes the smallest (hence the best) value of the largest distance between neighboring parallel lattice lines. However, in spectral tests with l-consecutive numbers generated by <z>, lattice vectors can have only rational components; see explicit forms of lattice vectors in terms of d, z and l given in Detailed Description of the Invention. Since the triangular lattice requires lattice vectors with irrational components, λd(2)(z) is invariably larger than this ideal case of the triangular lattice. The matter is the same with any other dimension l. The criterion, that the configuration is better in E2-spectral test for the lattice closer to the triangular lattice, will be utilized below as a quick and useful way of comprehension.

An epoch was marked by Fishman and Moore (1986) who initiated exhaustive spectral tests that examine all primitive roots for the prime modulus d=p=231−1. They adopted the criterion that a primitive root multiplier z is passable if λd(l)(z) is within the bound a=125% of geometrically ideal cases of the lattice (similar to the triangular lattice in 2-dimension) in all of dimensions 2≦l≦6. We note the smallest value λd(l), realized by the geometrically ideal case of lattices, for dimensions l=2, 3, . . . , 6 below citing Fishman and Moore (1986), but with a normalization that defines the lattice with period d in coordinates of El, which is adapted for the representation of lattice vectors with integer coordinates to be given later.













dimension l













2
3
4
5
6
















 λd(l)
2−1/231/4d1/2
2−1/6d2/3
2−1/4d3/4
2−3/10d4/5
2−1/231/12d5/6









Fishman and Moore (1986) note that, when the upper bound a for passable primitive root multiplier z was taken more loosely as a=133.3%, intractably many primitive roots passed the test, so that they tightened to a=125%, reducing the total number of successful primitive roots jumpingly to 414, viz. 7.74×10−7%≅1/(1.29×106) of the total number φ(231−2)=5.346×108 of primitive roots for the prime p=231−1. Later for another heavily used modulus of the type d=2i Fishman (1990) reported exhaustive spectral tests on relevant multipliers z≡5 mod(8) for d=232, and also gave tests on the modulus d=248 examining (effectively) 227 of multipliers among 245 total candidates, viz. 1/(2.6×105) of relevant generators. Remembering with dearness that Cray 2 in 1991 recorded 2 giga flops and thinking over supercomputers of today, we may contend that exhaustive spectral tests for the modulus d=248 are now realizable to obtain the period T=246. Yet this period will not suffice for simulations on supercomputers today.

    
    
        Fishman and Moore (1986)/G. S. Fishman and L. R. Moore: An exhaustive analysis of multiplicative congruential random number generators with modulus 231−1. SLAM Journal on Scientific and Statistical Computing, Vol. 7 (1986), pp. 24-45.
        Fishman (1990)/G. S. Fishman: Multiplicative congruential random number generators with modulus 2β: An exhaustive analysis for β=32 and a partial analysis for β=48. Mathematics of Computation, Vol. 54 (1990), pp. 331-334.
    
    


In preceding paragraphs features of spectral tests in prior art were seen by cases in which the set of points, generated by cyclic sequences and their cosets, occupy almost all of lattice points. A significant insight on the geometry of generated points, not the geometry of lattice points, are obtained from the simplest case of composite modulus d=p1p2 formed by two distinct odd primes p1 and p2. Structures of cyclic and coset sequences <z> and n<z> for this modulus have been analyzed in Theorems A, C. Visual information will give the quickest comprehension of implications. FIG. 5A depicts consecutive 2-tuples taken from one period of the cyclic sequence <z>, emitted by a generator embodying a simple case of the present invention, obeying prescriptions of Theorem A for the case s=2. The generator has the modulus d=2867=47×61. Note that p1=47 and p2=61 give

q1=(p1−1)/2=23; q2=(p2−1)/2=30,

which are coprime and realize the case (I) of Theorem C. Each square represents the range −0.05≦x, y≦1.05 again. The multiplier z is

z=678≡20U1+7U2 mod(2867)

where z1=20 and z2=7 are primitive roots of primes p1=47 and p2=61, respectively, both with high scores in their respective 2-dimensional spectral tests. The group Z*2867 has the order φ(2867)=φ(47)φ(61)=46×60=2760, and the period of <z> and its cosets are the largest T=2q1q2=1380=2760/2. Points on the right square represent (rj/d, rj+1/d), while those on the left square show (xj/z, xj+1/z) obtained from the numeral sequence of the irreducible fraction 1/d=1/2867=0.{dot over (x)}1x2 . . . {dot over (x)}T to base z=678. Since 1/z is sufficiently small, points on the right and the left are close to each other. FIG. 5B shows FIG. 5A with added points; the right is supplemented with points from the coset—<z>, and the left is added of points from a period of base-678 numeral sequence of the irreducible fraction −1/2867 or 2866/2867. Vacancies are seen in FIG. 5B to total to p1+p2−1=107 in number including the origin. Plots of FIG. 5B suggest a neat distribution of d=2867 lattice points in E2 on (nearly half portions of) which points of <z> and its cosets are nested. It is again tempting to use the sequence of <z> and its coset that together show the neat uniform distribution of numbers of Z*2867. But the successive appearance these sequences is the occurrence of random number sequence corresponding to 1/2867 and then the one to 2866/2867, which is highly improbable. Any concatenated sequence, n<z> and then −n<z>, should be rejected as implausible. Though the vacant seats in FIG. 5B are more conspicuous than in FIG. 2 with the sole vacancy at the origin, they are distributed uniformly so to say, and are not easy to be identified in the point sets of FIG. 5A. There is also the assurance of the shuffling by Chinese remainder theorem that the elements of the group Z*2867 are distributed in the interval [0, 2867) without accumulation. As regards the use of the sequence <z> or −<z> alone in this case, no statistical dubiousness is found in the present status of understanding.

FIG. 6 shows on the right the cyclic sequence <z> for d=2537=43×59 with z=289≡34U1+13U2 mod(2537). Here z1=34 and z2=13 are excellent primitive roots of p1=43 and p2=59, respectively. Since integers

q1=(p1−1)/2=21, q2=(p2−1)/2=29

are coprime, <z> has the largest period T=2q1q2=1218, and the group Z*2537 is divided into two by the cyclic subgroup <z> and its coset. However, both of q1 and q2 are odd, and the setting falls in the case of Theorem C (II) realizing −1∈<z>. Thus, only one of these halves of sequences <z> and n<z> may be used as uncorrelated random numbers. This inclusion of −1 in <z> manifests itself in the right plot of FIG. 6 as an obvious symmetry about the center M of the unit square; such a symmetry is absent in FIG. 5A. In fact, the following general inference holds true:

    
    
        Lemma D Let there be given a modulus d>4 which may be composite. Take any z, n∈Z*d and an arbitrary integer l≧1. In the l-dimensional Euclidean space El denote Sl for the set of points with coordinates congruent modulo d to consecutive l-tuples of integers in the coset sequence n<z> in Z*d,

Sl:={Pj(l)≡(nzj,nzj+1, . . . ,nzj+l-1)mod(d) j=0,1,2, . . . }.

Let S′l and Sl be point sets in El symmetric to Sl with respect to the point Md=(d/2, d/2, . . . , d/2) and to the origin O, respectively.

        (I) If −1∈<z> is true, there holds Sl=S′l. Namely, the point set formed by l consecutive numbers of n<z> defined modulo d is symmetric about the point Md, irrespective of the choice of n∈Z*d including the case n=1.
        (II) If −1 ∉<z> is the case, then Sl∩S′l=φ holds true. Namely, the point set formed by l consecutive numbers of n<z> defined modulo d is asymmetrical about the point Md, irrespective of the choice of n∈Z*d including the case n=1.

(Proof) We prove the case l=1; this will give indications of the proof for l≧2 given in Nakazawa and Nakazawa (2008a,b), where the case l=1 was not discussed. Integers of n<z> form a point set S on the line E1, to be called x-axis in this proof. Since n<z> is defined modulo d the set S is invariant under d-translations. Let S′ and S denote point sets on the x-axis which are symmetric to S about points x=d/2 and x=0, respectively; these are also d-translation invariant, and S may be denoted S=−S. The following is obvious:

S′={−x+d|x∈S}={−x|x∈S}= S. 

If −1∈<z> holds true, then—<z> is a sequential shift of <z>, so that S=S holds as sets. Thus we obtain S′= S=S. If −1 ∉<z> is the case, then <z> is a proper subgroup of Z*d, and disjoint as a set with the coset—<z>, implying S′∩S= S∩S=φ. ▪FIGS. 2-4 show primitive root multipliers for odd prime moduluses, and cyclic sequences <z> therein all include −1. Accordingly, right plots show point sets that are symmetric w.r.t. the center of the square. Similarly, FIG. 6 shows the composite modulus case with −1 included in <z>. The right plot is again symmetric w.r.t. the center of the square. In contrast, FIG. 5A stands for the case that <z> does not include −1, and its right plot shows no conceivable symmetry. Cases (I) and (II) in Lemma D give the dichotomy, so that we may state:

        Corollary to Lemma D Let z be a multiplier taken in the group Z*d of an arbitrary modulus d>4, and let Sl be the set of points in El formed by the consecutive l-tuples of a coset sequence n<z> for an arbitrary element n in Z*d. Inclusion of −1 in <z> is equivalent to the symmetry of Sl w.r.t. the center Md of the hyper cube of sides [0, d) for any one of dimensions l=1, 2, . . . .

Inclusion of −1 in <z>, or the symmetry of point sets formed by l-tuples of integers of n<z> in El, is in fact a more serious symptom than apparent symmetry of formed point sets; the second half of generated points take symmetric positions reproducing the order of appearance of their correspondents in the first half.

    
    


The present invention chooses settings of cyclic sequence <z> in the reduced residue class group Z*d for a composite modulus d, where −1 is excluded from <z> by design. An example was shown in FIG. 5A. As another we show in FIG. 1 the result for the setting

d=2867=47×61, z=1298≡29U1+17U2 mod(2867),

where z1=29 and z2=17 are excellent primitive roots for p1=47 and p2=61. The right plot of two-dimensional spectral test shows nearly triangular lattice. Comparisons with FIG. 5A will impress us that the multiplier z=678 for FIG. 5A will also be excellent in the two-dimensional spectral test. This is in fact the case. However, there can be a large difference in higher dimensions. The following table gives results of their spectral tests for l=2, 3 as comparisons to λd(l) given by geometrically ideal cases, indicating how can differences be for choices of the multiplier z.














2-nd degree
3-rd degree


multiplier z
score
score







z = 1298 ≡ 29U1 + 17U2 mod (2867)
104.96%
138.27%


z = 678 ≡ 20U1 + 7U2 mod (2867)
105.67%
196.36%









We saw what to be seen. We make use of a group Z*d with a composite modulus d in order to exploit the noted improvements given by the shuffling, and also to exclude −1 from the cyclic sequence <z>. The group Z*d for any modulus d has integers in [0, d) as its elements, with vacancies that may also be said to be distributed evenly therein. We choose parameters so as for these vacancies to be small in number, and sequences <z> and n<z> will sweep as large a portion of Z*d as possible, for uniformity and for the strength of spectral tests. Lagrange's theorem restrains that the (sub)group <z> may occupy only an integral portion of Z*d. And a composite modulus group cannot be cyclic. The largest portion that <z> can occupy in Z*d is thus the half, and the structure of the period shown in Theorem A rules that the modulus should then comprise s=2 odd distinct primes. The factor d0 of 1 or 2 in the modulus has no effect on the period and its participation in the shuffling is merely as a constant factor for d0=2:

nzj≡d+(d+1)(n1z1j+n2z2j+ . . . +nszsj)mod(2d),

so that we discard d0 altogether. We are thus left with the form d=d1d2, with d1=p1i1 and d2=p2i2 for odd primes p1 and p2.

Component groups Z*pkik for k=1, 2 have generating elements with periods

Tk=φ(pkik)=pkik(1−1/pk)≅pkik, k=1,2.

For the shuffling in <z> and its coset to be effective, we require T1≅T2≅t. Then the period of <z> is T=T1T2≅t2. Now the computability of spectral tests puts a practical limit on T or on t≅T1/2. Component moduluses dk=pkik (k=1, 2) should thus be determined by the requirement

t≅T1/2≅pkik, pk≅t1/ik, k=1,2.

Therefore, the ratio V of the number of voids, i.e. integers in [0, d) left unoccupied by the elements of Z*d, to the total number d of lattice points is

V={d−φ(d)}/d=1−(1−1/p1)(1−1/p2)≅1/p1+1/p2≅t−1/i1+t−1/i2.

This factor for a fixed t is minimized at i1=i2=1. The best setting for our aim is thus provided by moduluses of the form d=p1p2 comprising two distinct odd primes of similar magnitudes, on the premise that we exploit the largest period setting of Theorem A.

Summarizingly, we have first to determine the computable largest period T=t2. Then we choose distinct odd primes, p1, p2≅t; the choice should obey prescriptions of Theorem A and Theorem C that integers q1=(p1−1)/2 and q2=(p2−1)/2 are coprime and have different parity, one even and the other odd. The multiplier z is then formed by primitive roots z1, z2 of primes p1, p2, respectively, by congruence relations z≡z1 mod(p1) and z≡z2 mod(p2).
BRIEF SUMMARY OF THE INVENTION
The aim of the present invention is to provide a method to obtain uniform and independent random numbers with long periods and certified statistical quality, which may be realized on computers and their operating systems with the smallest computational cost.
The invention claims a method of obtaining uniform and independent random numbers comprising:

    
    
        1. two odd primes p1 and p2 which are prescribed to give mutually coprime integers

q1:=(p1−1)/2, q2:=(p2−1)/2

        with different parity, i.e. one of {q1, q2} is odd and the other is even,
        2. primitive roots z1 and z2 of prime moduluses p1 and p2, respectively,
        3. the modulus d=p1p2,
        4. the multiplier z coprime with d=p1p2, and determined modulo d uniquely by congruence relations

z≡z1 mod(p1), z≡z2 mod(p2),

        in the interval (0, d), and
        5. an integer n coprime with d chosen arbitrarily in the interval (0, d).

Said method constructs a coset sequence of integers for the initial value or seed n,

n<z>={r1,r2, . . . ,rT}

of the cyclic sequence <z> of said multiplier z, which has the order T=φ(d)/2=2q1q2 in the reduced residue class group Z*d modulo d=p1p2, by solving the congruence relations

r1=n, rj+1≡zrj mod(d), 0<rj<d, 1≦j≦T=φ(d)/2=2q1q2,

and gives the output sequence {v1, v2, . . . , vT} for uniform and independent random numbers of period T by realizing the arithmetic,

vj=rj/d, 0<vj<1, 1≦j≦T. 

    
    


The heart of the invention is in selecting said primes p1, p2 that give coprime integers q1, q2 with different parity. This prescription expels −1 from the cyclic sequence <z> generated by said multiplier z, enabling the use of the whole sequence, obtained from <z> or its coset by normalization, for uniform and independent random numbers without recognizable correlations over the whole period.



BRIEF DESCRIPTION OF THE DRAWINGS
FIG. 1 shows the 2-dimensional spectral test of a generator designed by the present invention; the right plots consecutive 2-tuples of outputs as points in the square −0.05≦x, y≦1.05; the generator comprises the composite modulus d=2867=47×61, the multiplier z=1298≡29U1+17U2 mod(2867) and the initial value n=1 with the period T=1380; the left plot shows the corresponding consecutive two tuples of quotients, in the division process of the irreducible fraction n/d to base z, which are divided by z to be normalize to the interval [0,1).
FIG. 2 shows the same setting as FIG. 1 by a generator of the prior art; prime modulus d=p=251, primitive root z=162, and T=p−1=250 for the period.
FIG. 3 shows plots of FIG. 2 taken on a different primitive root multiplier z=61.
FIG. 4 shows plots of FIG. 2 taken on still another primitive root multiplier z=127.
FIG. 5A shows a result of the present invention in the same setting as FIG. 1; the generator comprises the composite modulus d=2867=47×61, the multiplier z=678=−20U1+7U2 mod(2867) and the initial value n=1 with the period T=1380.
FIG. 5B shows plots of FIG. 5A with added points; the right is supplemented by points from the coset—<z>, and the left is added of points from the base-z numeral sequence of the irreducible fraction −1/2867 or rather 2866/2867.
FIG. 6 gives plots for the setting out of the present invention for composite modulus d=2537=43×59, the multiplier z=678≡12U1+13U2 mod(2537), and with n=1, T=1218.
FIG. 7A shows another failure occurring with the composite modulus d=2ip with odd prime p; d=21248=83×28, z=7461≡74U1+37U2 mod(21248); the right shows 2-tuples from <z> normalized by 1/d, and the left shows 2-tuples from the numeral sequence to base 7461, normalized by 1/7461, of the fraction 1/21248.
FIG. 7B shows plots of FIG. 7A with added points; the right is supplemented by points from the coset (74−1U1+U2) <z>, and the left is added of corresponding points from the base-z numeral sequence of the irreducible fraction 10753/21248.
FIG. 8 shows an inadequate setting that can arise in the present invention; the composite modulus d=3599=59×61 is formed by twin primes which are too close for shuffling to be effective; z=89≡13U1+44U2 mod(3599) and T=1740.



DETAILED DESCRIPTION OF THE INVENTION
In the presented invention, the design of a generator of uniform and independent random numbers starts by choosing the period T, which will need some trials in order to find the largest magnitude of T so as for resultant spectral tests to be computable.
Given the period T, two odd, distinct primes p1, p2, for the formation of the modulus d=p1p2, should be chosen with similar magnitude, typically of the order of T1/2, so as to minimize the ratio to d of the number of voids, viz. the number of integers in the interval [0, d) left unoccupied by the elements of the group Z*p1p2 (i.e. multiples of p1 and p2 in the interval [0, d) including 0),

{d−φ(d)}/d=1−(1−1/p1)(1−1/p2)=1/p1+1/p2−1,

under the fixed period T=φ(d)/2≅p1p2/2. There exist further indispensable restrictions in choosing primes p1, p2 Namely, integers

q1=(p1−1)/2, q2=(p2−1)/2

should be coprime, and should have different parity. Caution should also be paid to prevent the excessive closeness of p1 and p2. FIG. 8 will give an intuitive comprehension in this regard. Here the choice is twin primes p1=59 and p2=61 for the modulus d=3599=59×61, together with the multiplier z=89≡13U1+44U2 mod(3599) formed with excellent primitive roots z1=13 for p1=59 and z2=44 of p2=61. Component cyclic sequences <13> mod(59) and <44> mod(61) have too close periods T1=58 and T2=60, and the shuffling mechanism between them does not work well. Or the shuffling (or rather mixing) in this case introduces new correlations similar to beat phenomena in sound. The right of FIG. 8 shows poor performance of consecutive 2-tuples taken from <z>, together with the left plot of corresponding consecutive 2-tuples of numerals (divided by z for normalization) arising in the division process of 1/d=1/3599 to base z. By the way, q1=29 and q2=30 are different in parity, but this virtue is made invalid by the closeness of periods.

A successful case of the invention p1=47 and p2=61 is shown in FIG. 1 with the multiplier z=1298≡29U1+17U2 mod(2867). Though this is practically a toy model, the example certainly gives some measure of how large the difference |p1−p2| should be, or how (p2/p1)1/2 should be apart from 1.
As the procedure to construct the generator, the established prior art then follows to choose best primitive root multipliers, {z1, z′1, . . . } of the prime p1 and {z2, z′2, . . . } of the prime p2, by exhaustive spectral tests of 6-th degree. The list shown below,

 λd(2)=2−1/231/4d1/2=0.930605d1/2,

 λd(3)=2−1/6d2/3=0.890899d2/3,

 λd(4)=2−1/4d3/4=0.840896d3/4,

 λd(5)=2−3/10d4/5=0.812252d4/5,

 λd(6)=2−1/231/12d5/6=0.774899d5/6,

reproduces the smallest values λd(l) realized by geometrically ideal lattices, typically the triangular lattice in the plane, with a normalization adapted for the computation with lattice and dual lattice vectors of integer components to be given shortly. The largest distance between neighboring hyperplanes containing l−1 linearly independent lattice vectors can have only larger values than listed above, in respective Euclidean space El of l-dimension. Once these exhaustive spectral tests are done, say in accordance with the criterion adopted by Fishman and Moore (1986), we need to find the multiplier z∈Z*p1p2 as the solution, unique modulo d=p1p2, of congruence relations

z≡z1 mod(p1), z≡z2 mod(p2), 0<z<d=p1p2.

The procedure of solution given by Chinese remainder theorem may be briefed as follows.

    
    
        1. Find integers A, B that give the identity 1=p1A+p2B by Euclidean algorithm, and put U1≡p2B mod(p1p2), U2≡p1A mod(p1p2).
        2. Compute z by the congruence relation z≡z1U1+z2U2 mod(d=p1p2) in the interval (0, d).
    
    


The final work is to perform the spectral test, up to the appropriate degree l, on the multiplier z in the group Z*d with the composite modulus d=p1p2. The procedure is the same as in prior art, but some exposition will be in order. The aim is to obtain the largest distance λd(l)(z) between the neighboring lattice hyperplanes determined by d, z and l. Given a possibly composite modulus d>0, a multiplier z∈Z*d and l=2, 3, . . . for the dimension, the spectral test examines the lattice which may be regarded to be spanned by the following bases or basis vectors,

e1=(1,z,z2, . . . ,zl-1),

e2=(0,d,0, . . . ,0),

e3=(0,0,d, . . . ,0),

. . .

el=(0,0,0, . . . ,d).

The lattice itself is the set of vectors formed by all of linear combinations of these basis vectors with integer coefficients. Since these are manifestly linearly independent and any integral multiples of e2, e3, . . . , el may be added to e1 without destroying this linear independence, the first basis vector e1 may be taken in the sense of modulo d, though we proceed as above. Let Cl(d) denote the hypercube of El whose coordinates are all in the interval [0, d). For any lattice vector to be in Cl (d) the first coordinate requires that the coefficient c1 of e1 should be among 0, 1, . . . , d−1. Upon assignment of one of these d values to c1, the coefficient cj of ej for 2≦j≦l is determined uniquely by the requirement that c1zj−1+cjd is in the interval [0, d). There exist thus d lattice points (or seats for the l-tuples of numbers from n<z>) in the hypercube Cl(d), irrespective of l. It may also be noted that the hyper parallelepiped, so to say, spanned by basis vectors has the volume (corresponding to the determinant obtained by arranging basis vectors in rows) dl-1, which is another way of saying that there are d lattice points in Cl(d) of volume dl. For a given set of numbers d, z∈Z*d and l=2, 3, . . . , the aim of the l-dimensional spectral test is to search the largest distance between neighboring (l−1)-dimensional hyperplanes that contain the same set of l−1 linearly independent lattice vectors. The problem is transformed tractably by the use of dual lattice (or reciprocal lattice of solid state physics) spanned by its basis vectors

f1=(d,0,0, . . . ,0),

f2=(−z,1,0, . . . ,0),

f3=(−z2,0,1, . . . ,0),

. . .

fl=(−zl-1,0, . . . ,1),

which give characteristic inner products with lattice basis vectors as follows:

(ej,fk)=dδjk, 1≦j, k≦l. 

The setting reveals that the problem is identical with the spectral test for a primitive root multiplier z for an odd prime modulus p, if only we replace p with d and translate the distribution of points from <z> to that of their seats. We thus need to find the length |f|min of the shortest vector with non-vanishing length in the dual lattice, which gives λd(l)(z)=d/|f|min. In short, the spectral test of the present generalized setting is completely the same with the prior art for primitive root multipliers of odd prime moduluses.

If z turns out to be unsatisfactory in the closeness of its λd(l)(z) to the theoretical optimum λd(l), different choices of component multipliers should be tried and the test need be repeated. If the result is passable, we have a design of the uniform and independent random number generator with the modulus d=p1p2 and the multiplier z with the period T=2q1q2.
The computability of spectral test is a difficult problem. Some rough estimates were discussed in Addendum of Nakazawa and Nakazawa (2008b). The core of spectral tests is in the preparatory procedure to find a small length L of the dual vector formed by linear combinations of integral multiples of dual basis vectors using, typically, Schmidt orthogonalization. Once this is done giving a small enough L, the spectral test needs to take only dual lattice vector f=(a1, a2, . . . , al) in Cartesian coordinates fulfilling the condition

a1+a2z+a3z2+ . . . +alzl-1≡0 mod(d), |aj|≦L, 1≦j≦l. 

Namely, the test sifts all vectors filling the above restriction out of (2L+1)l/2 candidates, computes their lengths, sorts out the shortest |f|min of such lengths, and gives finally the valuation λd(l)(z)=d/|f|min. The computational load of this whole process will be hard to be measured quantitatively. However, we know for certain that at present an exhaustive spectral test of d≅248 will be possible. If the present invention is to be applied to this case, then odd primes p1, p2 of the order of 225 in magnitude, say, should be chosen and exhaustive spectral tests on their primitive roots should be performed. These are easily realized now, and we need to perform several spectral tests on the combined modulus d=p1p2, which is certainly computable again. The computational economy to be realized by the present invention will be manifest, together with the extensibility of attainable period T. See Addendum of Nakazawa and Nakazawa (2008b) for detailed arguments. Fishman and Moore (1986) notes that the best primitive root showed the score of 120.21% as the average of dimensions l=2, 3, . . . , 6. This will give us an indication how far we should or could proceed.

A composite modulus of the type d=2id′ with large index i and an odd d′ should not be used. Though this was found, and noted in detail, in Sec. 7B of Nakazawa and Nakazawa (2008a,b), some visual information will be in order at this place. FIG. 7A takes the case d=21248=83×256, z=7461=74U1+37U2 mod(21248) with the period T=2624=41×64. The right of FIG. 7A depicts 2-tuples taken from <z> for a period T. The left shows 2-tuples from the numeral sequence of the irreducible fraction 1/21248 to base z=7461. Patterns, or figures of correlations, are conspicuous. Nevertheless, the right of FIG. 7B, showing the plot of FIG. 7A supplemented with 2-tuples from a coset (among three) of <z> shows uncriticizably neat and regular occupation of (half of) the group Z*21248. The spectral test for this z and the modulus d can never detect this concealed correlation appearing in the occupation of seats by 2-tuples from <z>.
The key prescription of the invention is that one of said {q1, q2} is even and the other is odd. This ensures generically that −1≡d−1 mod(d) is expelled from the cyclic sequence <z>, that the whole of the sequences <z> and n<z> may be used for random numbers, and that for any l=1, 2, 3, . . . these sequences in the group Z*p1p2 yield l-tuples of consecutive integers which, regarded as points in El, form a set devoid of improbable symmetry with respect to the center Md=(d/2, d/2, . . . , d/2) of the hypercube in El of sides [0, d]. In the contrary case, shown for example in FIG. 6, the cyclic subgroup <z> may be used for independent random numbers only for the half of the period, and the inferences of spectral tests on the geometrical configuration of d seats in the cube of El considerably lose their power.
In closing, we reflect on our path through the jungle we have waded; it started from the setting of finite sequences of integers, periodic sequences, irreducible fractions, and their replacement by cyclic sequences in reduced residue class groups. We arrived at the group Z*p1p2 for distinct odd primes p1, p2 forming the modulus, its element z for the multiplier giving the coset sequence n<z>, and the normalized coset sequence {vj} with period T=φ(p1p2)/2,

{vj=rj/d|r1=n, rj+1=zrj mod(d), 0<rj<d, 1≦j≦T}. 

Note that {vj=rj/d} is not the aimed random numbers. The output should be the original sequence {uj=xj/z|1≦j≦T} corresponding to the integer sequence {x1, x2, . . . , xT} from which we have started. Yet, it is at once to realize this final link. The identity of the division process for n/d in numerals to base z readily gives, with the use of periodicity vj+T=vj for any j,

uj:=xj/z=(zrj−rj+1)/(dz)=rj/d−(rj+1/d)/z=vj−vj+1/z, 1≦j≦T. 

which form a set of linear transformations from {vj} to {uj}. Here we remind ourselves that vj=rj/d has so far been regarded to be defined as an irreducible fraction without any approximation, while the above noted transformations realize uniform random numbers {uj} that were taken as the discretization of natural random numbers with the smallest unit 1/z. We come here also to realize that the original sequence of integers {xj} was taken in the interval [0, z) without any restrictions on the magnitude of z; it may thus happen that the irreducible fraction n/d representing the concatenated sequence of {xj} may have the base z larger than d, as with the case of 2/7 with its decimal sequence 0.285714 for z=10. Even if z>d is the case, the closeness of {uj=xj/z} to {vj=rj/d} within 1/z holds true. But the relation 1/z<1/d implies that the sequence {uj} traces properties of the normalized coset sequence {vj}. Thus, 0 can never arise in {uj}, or duplicated values can never occur in the original sequence {xj}. If the given sequence of integers {xj} in the interval [0, z) happens to give an irreducible fraction n/d with d smaller than z, the sequence should be discarded as inadequate for uniform and independent random numbers. We should tighten starting assumptions: The upper bound z of integers x1, x2, . . . , xT should fulfill 0<z<d in the result of the representation by an irreducible fraction x=n/d. Since this is a natural and usual premise set at starting the search for a good sequence {uj} by spectral tests of multiplicative congruential generators, there arises no practical problem. Usual practices tell us that a good multiplier z<d has the magnitude d1/2 or larger, which may be interpreted that its multiplication on the preceding integer rj should give rj+1=zrj>d to be pulled back modulo d with sufficient frequency. This gives an empirical rule that good random number sequence should be realized with z fulfilling 1>>1/z≅1/d1/2>>1/d. Now that vj=rj/d is in effect the discretization in unit of 1/d, the use of these {vj} in the transformation uj=vj−vj+1/z is excessive in precision. It is then all natural to discard the term −vj+1/z. In order to avoid confusion we denote {v′j} for the representation of {vj} in the smallest unit ∈, obtained by truncating and rounding off vj=rj/d. The choice of ∈≅1/z makes the term −v′j+1/z irrelevant, and we may use {v′j} itself as the approximation of {uj} within small errors of the order of 2/z, say. This is not only a simple and economical way, but also realizes the range [0, 1] for {u′j:=v′j} by 1/d≦vj=rj/d≦1−1/d. It notably admits the appearance of duplicated values in one period T of {u′j=v′j} reproducing a natural property of random numbers in the original sequence {uj}, but not in {vj}. Last but not least, the valuations of spectral tests on {vj} are inherited also by {u′j=v′j}, just as by the exact {uj}.

We thus recommend in the present invention the use of said {v′j}, with the smallest unit of values taken as ∈≅1/z or smaller, as the most economical and natural method of generating uniform and independent random numbers. We note that all of FIGS. 1-8 show points generated by {uj=vj−vj+1/z} and {vj} with units far smaller than 1/d. Plots corresponding to {u′j=v′j} will have point configurations which are within the range of 3/z from those in the right plots. Comparisons of right and left plots with differences of 1/z will convince us that the direct use of {v′j} for ∈≅1/z will suffice for a practical use. Even a larger choice of E may well be taken, if the structure of the generator is hoped to be more concealed to users; but an excessively large E may weaken the power of spectral test valuations.
The present invention is believed to have found a remedy of faults in primitive root multiplicative congruential generators and devices to realize more economical spectral tests, together with general theoretical prospects on uniform and independent random number problems. We advocate ways suggested in Nakazawa and Nakazawa (2008b) to extend the period of random number sequences with tested statistical properties. Structures found might well be said rich and significant, but many problems are certainly left unclarified. Learning many pitfalls met in the history of random number theory, we hope in earnest that this survey is free from any such pitfalls. Contributions are waited for from people in various fields, in presenting excellent multipliers they find for their own rights along with the progress of computers, as well as in giving still new methods realizing excellent uniform and independent random numbers.

## Claims

1. A system for debugging a simulation program on a computer, the system comprising:
a processor for calculating a random number sequence using the steps of:
taking two distinct odd primes p1, p2 in such a way that integers


q1:=(p1−1)/2, q2:=(p2−1)/2



are mutually coprime with one of q1, q2 is odd and the other is even;
taking the integer d=p1p2 as the product of said odd primes p1 and p2;
taking a primitive root z1 modulo said odd prime p1 and a primitive root z2 modulo said odd prime p2;
constructing the sequence of integers {r1, r2, r3, . . . } by realizing recursive congruence relations


r1≡n mod(d), rk+1≡zrk mod(d), 1≦rk<d, k=1,2, . . .



for an arbitrarily given initial value n coprime with said d, where z is determined modulo d uniquely by congruence relations


z≡z1 mod(p1), z≡z2 mod(p2);

obtaining the sequence {v1, v2, . . . } by realizing the arithmetic


vk=rk/d, k=1,2,3, . . . ; and


outputting {v1, v2, v3 . . . } as the aimed sequence for uniform and independent random numbers,
wherein the multiplicative congruential random number generator (d; z) defined by said modulus d and said multiplier z shows 2nd degree and 3rd degree spectral test performances, which is the same as that two dimensional vectors of the form v=(j1,j2) with integer coordinates fulfilling


j1+zj2≡0 mod(d)


and with non-zero length ∥v∥=(j12+j22)1/2>0 have their shortest vector vmin with the property


∥vmin∥>21/2d1/2/(1.25×31/4)



and as that three dimensional vectors of the form v=(j1, j2, j3) with integer coordinates fulfilling


j1+zj2+z2j3≡0 mod(d)


and with non-zero length ∥v∥=(j12+j22+j32)1/2>0 have their shortest vector vmin, with the property


∥vmin∥>21/6d1/3/1.25; and


a memory that stores a simulation program suitable for operation on a computer,
wherein the calculated random number is reproducible and transportable allowing for debugging, via the processor, the simulation program stored in memory, and saving the debugged simulation in the memory.

2. The system of claim 1, further comprising a simulatable element that is simulated using the saved debugged simulation program stored in memory.

3. A method of debugging a simulation program on a computer using a reproducible and transportable random number sequence, the method comprising:
identifying a simulation program stored in memory that needs to be debugged;
debugging the simulation program via a processor using the processor to generate at least one reproducible and transportable random number sequence by:
taking two distinct odd primes p1, p2 in such a way that integers


q1:=(p1−1)/2, q2:=(p2−1)/2



are mutually coprime with one of q1, q2 is odd and the other is even;
taking the integer d=p1p2 as the product of said odd primes p1 and p2;
taking a primitive root z1 modulo said odd prime p1 and a primitive root z2 modulo said odd prime p2;
constructing the sequence of integers {r1, r2, r3, . . . } by realizing recursive congruence relations


r1≡n mod(d), rk+1≡zrk mod(d), 1≦rk<d, k=1,2, . . .



for an arbitrarily given initial value n coprime with said d, where z is determined modulo d uniquely by congruence relations


z≡z1 mod(p1), z≡z2 mod(p2);

obtaining the sequence {v1, v2, . . . } by realizing the arithmetic


vk=rk/d, k=1,2,3, . . . ; and


outputting {v1, v2, v3 . . . } as the aimed sequence for uniform and independent random numbers,
wherein the multiplicative congruential random number generator (d; z) defined by said modulus d and said multiplier z shows 2nd degree and 3rd degree spectral test performances, which is the same as that two dimensional vectors of the form v=(j1,j2) with integer coordinates fulfilling


j1+zj2≡0 mod(d)


and with non-zero length ∥v∥=(j12+j22)1/2>0 have their shortest vector vmin with the property


∥vmin∥>21/2d1/2/(1.25×31/4)



and as that three dimensional vectors of the form v=(j1, j2, j3) with integer coordinates fulfilling


j1+zj2+z2j3≡0 mod(d)


and with non-zero length ∥v∥=(j12+j22+j32)1/2>0 have their shortest vector vmin, with the property


∥vmin∥>21/6d1/3/1.25; and


storing the debugged simulation program in memory.

4. The method of claim 3, comprising simulating an element that is simulated using the saved debugged simulation program stored in memory.