Patent Application
20120290632
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
r
1
≡n mod(d), r1+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 realizing the arithmetic
v
j
:=r
j
/d, j=1,2,3, . . . .
2. Description of the Related Art
An odd prime modulus d=p≧3 has the primitive root multiplier z that generates the cyclic sequence {z, z2, z3, . . . , zp−1≡1}. In the sense of congruence modulo p, members of this sequence sweep over all integers coprime with p in a period T=p−1≈d, and present a realization of uniform distribution. Yet, an odd d=p≧3 gives q=(p−1)/2 as an integer, and the cyclic sequence of any primitive root z includes zT/2=zq≡−1 mod (p) inevitably. Thus, the last half of the sequence {zq+1≡−z, zq+2≡−z2, . . . , z2q≡−zq} is essentially a repetition of the first half. The use of a primitive root cyclic sequence for independent random numbers should be restricted to its half period T/2≈d/2.
The preceding patent application Ser. No. 12/379,964 proposed the method to generate uniform and independent random numbers comprising the steps of:
q
1:=(p1−1)/2, q2:=(p2−1)/2
z≡z
1 mod(p1), z≡Z2 mod(p2)
v
j
=r
j
/d, 0<vj<1, 1≦j≦T=2q1q2,
(z1q1q2,z2q1q2)≡((−1)q2,(−1)q1)≈(1,−1)≠−(1,1)mod(p1,p2).
Thus, the sequence {vj|1≦j≦T=2q1q2≈d/2} of patent application Ser. No. 12/379,964 may be used wholly for uniform and independent random numbers. If n is a member of the cyclic sequence with n≡zk mod (d), then {nz≡zk+1, nz2≡zk+2, . . . } is a mere translation of the cyclic sequence. If n is coprime with d but is not a member of the cyclic sequence generated by z, then {nz, nz2, nz3, . . . , nzT≡n} with T=2q1q2 shares the same period, but has no common elements, with the cyclic sequence. In terms of groups this is the coset of the cyclic subgroup generated by z. In the way of construction of patent application Ser. No. 12/379,964, members of cyclic and coset sequences form the division of integers coprime with d.
The present invention adopts the same setting as patent application Ser. No. 12/379,964 for the modulus d, but replaces the primitive root z1 with any non-primitive integer y1 with the half-full order q1 modulo p1, and adopts the multiplier z′ determined by congruence relations
z′≡y
1 mod(p1), z′≡z2 mod(p2).
Sequences {(z′)j mod (d=p1p2)|j=1, 2, . . . } and {(y1j,z2j) mod (p1,p2)|j=1, 2, . . . } again share the same period. As ord(y1)=q1 mod (p1) is odd, ord(z2)=2q2 mod (p2) is even and as q1 and q2 are coprime, the period is T=2q1q2. And, the integer y1 with an odd order never gives (y1)j≡−1 mod (p1) and the sequence {(z′)j mod (d=p1p2)|1≦j≦T} is devoid of −1. The whole period of the coset sequence of the present invention,
{rj≡n(z′)j mod(d)|0<rj<d=p1p2, 1≦j≦T=2q1q2}
may again be used for uniform and independent random numbers.
The design of generators of uniform and independent random numbers necessarily calls for some tests of the resulting sequence. Spectral tests utilize the following fact. Consider a modulus d and a multiplier z coprime with d. Consecutive s-tuples of the cyclic sequence generated by z modulo d, as well as consecutive s-tuples of any coset sequence, generate points in the s-dimensional Euclidean space E, that notably occupy a portion (but not all) of points of a lattice determined by d and z. A little more surprising fact is that any consecutive s-tuples taken from any sequence of integers are again approximated by points of some lattice. This general statement is of little help in evaluating the performance of an arbitrary integer sequence in giving independent and uniform random numbers; as the presence of any lattice cannot usually be recognized in plotting s-tuples of such integer sequence in Es, integer sequences almost invariably give points in E, that occupy only a negligibly small portion of the noted lattice. Yet, there are conspicuous exceptions. A prime modulus d=p and its primitive root multiplier z generate a cyclic sequence that gives s-tuples forming p−1 points in Es, where the noted lattice gives p lattice points. In the setting of the present and the preceding inventions with the modulus d=p1p2 formed by odd primes p1 and p2, multipliers z′ as well as z constructed in noted special ways generate cyclic sequences whose consecutive s tuples occupy 2q1q2≈d/2 of points that are nearly the half of d lattice points, with similar 2q1q2 points from the coset sequence occupying almost all of the remaining half. In these cases the configuration of the lattice gives significant indications of the mode of distribution of consecutive random numbers, hence on the mode of independence, given by multiplicative congruential methods. These structures, and in particular the new perception of the following symmetries (or asymmetries) motivated for the present invention.
Lemma 1. Let an odd prime p1≧3 give the odd integer q1=(p1−1)/2. Then, any integer y1 with the half full order ord(y1)=q1 modulo p1 is necessarily the minus of a primitive root z1 of p1,y1≡−z1 mod (p1). Conversely, any primitive root z1′ of p1 is the minus of an integer y1′ with the half full order modulo p1, z1′ mod (p1).
Lemma 2 Let an odd prime p2≧3 give the even integer q2=(p2−1)/2, and let z2 be a primitive root of p2. Then the integer z2″≡−z2 mod (p2) is also a primitive root of p2.
Corollary 3. Let odd primes p1≧3 and p2>3 give mutually coprime integers
q
1:=(p1−1)/2, q2:=(p2−1)/2
with odd q1 and even q2. Call an integer z to be in class I modulo d=p1pe if it is coprime with d and determined by congruence relations
z≡z
1 mod(p1), z≡z2 mod(p2),
with a primitive root z1 of p1 and a primitive root z2 of p2. Also call an integer z′ to be in class II modulo d if it is coprime with d and determined by congruence relations
z′≡y
1′ mod(p1), z′≡z2′ mod(p2),
with an integer y1′ that has the half-full order ord(y1′)=q1 mod (p1) and with a primitive root z2′ of p2. Then a class II multiplier is the minus of a class I multiplier modulo d, and a class I multiplier is the minus of a class II multiplier modulo d. In this correspondence a class I multiplier and a class II multiplier share identical valuations in spectral tests, and the spectral tests need be performed on either one of multipliers.
(Proof of Lemma 1) Put j=ord(−y1) mod (p1). By definition j is the smallest positive integer that gives (−y1)j≡1 mod (p1). There exist two alternative possibilities:
(−1)j=−1, (y1)j≡−1 mod(p1), (a)
or
(−1)j=1, (y1)j≡1 mod(p1). (b)
The case (a) requires that j is odd and, since the order of y1 is q1 modulo p1, also that j is an odd multiple of q1/2. This latter requirement cannot be met by an integer j because q1 is odd by assumption. The remaining (b) requires that j is even and j is a multiple of q1, thus selecting even multiples of the odd q1 for j. Hence ord(−y1)=2q1 mod (p1) holds true, and z1=−y1 is a primitive root of p1, proving y1≡−z1 mod (p1) for some primitive root z1 of p1. Conversely, take any primitive root z1′ of p1, put y1′≡−z1′, and let j=ord(y1′) mod (p1). There should hold 1≡(−z1)j mod (p1), giving the following two possibilities:
(−1)j=−1, (z1′)j≡−1 mod(p1), (c)
or
(−1)j=1, (z1′)j≡1 mod(p1). (d)
The case (c) requires an odd integer j that is also an odd multiple of the half order q1 of the primitive root z1′. Hence any odd multiples of odd q1 is admitted by (c) for j. The second half of the case (d) admits any integral multiples of ord(z1′)=2q1, and these automatically let the first half hold. The smallest positive of these gives j=q1=ord(y1′). Thus, any primitive root z1′ of p1 has the expression z1′=−y1′ by an integer y1′ with the half-full order q1.
(Proof of Lemma 2) Let the prime p2 have a primitive root z2 with the order ord(z2)=p2−1=2q2; by assumption q2 is even. Let j be the order of y2≡−z2 mod (p2). The integer j fulfills one of the following:
(−1)j=1, z2j≡−1 mod(p2), (a)
or
(−1)j=1, z2j≡1 mod(p2). (b)
The case (a) requires that j is odd and is an odd multiple of q2, the half order of the primitive root z2. But q2 is even by assumption, and (a) can never be satisfied. The case (b) requires an even j that is also an integral multiple of the even ord(z2)=2q2. As the smallest positive of such, the order j of y2≡−z2 mod (p2) is 2q2, and is again a primitive root of p2.
(Proof of Corollary 3) Let z be a class I multiplier modulo d corresponding the vector (z1, z2) with a primitive root z1 of p1 and a primitive root z2 of p2. Lemmas 1 and 2 give
(z1,z2)≡(−y1′,−z2′)≡−(y1′,z2′)mod(p1,p2),
where y1′ has the half-full order q1 modulo d and z2′ is a primitive root modulo p2. Let the class II multiplier z′ correspond to the vector (y1′, z2′). The correspondence (group isomorphism) ensured by Chinese remainder theorem translates the above relation of vectors as follows:
z≡−z′ mod (d).
The converse statement z′≡−z mod (d) is manifest. Arithmetically, the s dimensional spectral test for a multiplier z and a modulus d takes the vector (j1, j2, . . . , js) of integers fulfilling
j
1
j
2
z′+ . . . +j
s(z′)s−≡0 mod(d)
with {jk=±1, ±2, . . . |1≦k≦s}, calculates its Euclidean length
λ=∥(j1,j2, . . . , js)∥=(j12+j22+ . . . +js2)1/2,
searches for the smallest positive value λmin of such lengths, and gives the valuation d/λmin for the multiplier z. Thus, the test evaluates the class I multiplier z and its corresponding class II multiplier z′≡−z identically. The converse will be obvious.
The aim of the present invention is to provide a method to generate uniform and independent random numbers with long periods and certified statistical quality that may be realized on computers and in their operating systems with the smallest computational cost.
The invention claims a method of obtaining uniform and independent random numbers comprising
{rj≡nzj mod(d)|1≦j≦T=2q1q2},
by solving recursively the congruence relation
r
1
≡n mod(d), rj+1≡zrj mod(d), 0<rj<d, 1≦j≦T=2q1q2,
and gives the output sequence {v1, v2, . . . , vT} for uniform and independent random numbers of period T by realizing the arithmetic,
v
j
r
j
/d, 0<vj<1, 1≦j≦T.
The heart of the invention is the same as patent application Ser. No. 12/379,964 in selecting said odd primes p1,p2 that give an odd integer q1=(p1−1)/2 and an even integer q2=(p2−1)/2 where q1 and q2 are coprime, but differs in constructing the multiplier z with a non-primitive integer y1 with ord(y1)=q1 mod (p1).
Merits obtained by this new method are:
The present invention aims to realize designs of generators for uniform and independent random numbers on computers with sufficiently long periods. The invention applies the multiplicative congruential method with moduluses formed by products of two large odd primes, and gives the largest conceivable extension of the method of patent application Ser. No. 12/379,964 with the same facilities for the design of output performances.
As usual with multiplicative congruential methods the preparation of the initial value or seed n should be left to users. The principle is to choose n so as to be coprime with the modulus d=p1p2 formed by two distinct primes p1,p2. The matter may be simplified for users by asking them to choose two integers n1,n2 satisfying 0<n1<p1 and 0<n2<p2. Programs for random number generation will readily solve congruence relations
n≡n
1 mod(p1), n≡n2 mod(p2),
and provide the seed n to start the program.
The greatest merit of the composite modulus d=p1p2 formed by two odd primes arises with spectral tests. Take the Mersenne prime modulus d=p=231−1 adopted by Fishman and Moore in 1986 in their monumental work. This modulus is certainly too small for simulations of today. Yet, it still typifies fundamental aspects of problems and guides us with the clear overview. Fishman and Moore set the criterion that the passable performance should be less than 1.25. First of all, experiments (say, performed on our small computers of today at least to some partial extent) will readily convince us that this criterion 1.25 is an invaluable finding of Fishman and Moore. Second, experiments will also convince us that exhaustive spectral tests, sweeping over all primitive roots as initiated by Fishman and Moore (1986) on this modulus, are inevitable in finding good primitive roots. The statement is all true with moduluses and multipliers of the present invention. Third, experiments will show that the computing time of 2nd to 6th sweeping spectral tests increases very rapidly with the magnitude of the modulus. Compared to the complete set of 2nd-6th spectral tests of a single primitive root, tests sweeping all primitive roots for the modulus d=231−1 take about 4×104 or more of computing time. The reduction of the difficulty enabled by adopting two primes p1,p2 of the order of magnitudes of d1/2 is magnificent, as to be reported.
As regards the quality of combined multipliers, experiments suggest that the performance of the multiplier z corresponding to ±(z1,z2) may be expected to be no worse than the product of respective performances of z1 and z2. Thus, if z1 and z2 are chosen with the performances smaller than 1.25, i.e. if z1, z2 both construct lattices that are respectively within 125% departure from their geometrically ideal configurations, then we may expect the departure of the lattice formed by z within 1.252≈1.50 or 150% from its geometrical ideal, at least for some, but of course not all, pair ±(z1, z2) of multipliers. This result might seem not so good as best primitive root multipliers for a prime modulus. Imagine, however, the case d≈1022≈273. The selection of good multipliers by exhaustive tests for this magnitude of a prime modulus d=p will not yet be possible; we lack means to find good primitive roots for p. But exhaustive tests on primes with the order of magnitude of 237 will now be possible, and we may think of tests of combinations of selected primitive roots in the composite modulus d=p1p2≈274. The noted experience will also help us to adopt the computing scheme that discards those pairs with performances not less than 1.50, say; then the speed of computation will be remarkably increased in sweeping all selected combinations. Since the obtained uniform and independent random numbers have a very large precision, of the order of d−1≈2−74, and are to be substituted into real double precision variables with its recognizable smallest units of about 2−53, each of the noted units contains approximately 221 random numbers in a period. Besides, the performance of 1.50 to be assured in all 2 to 6 dimensional distribution might be said satisfactory, in particular in circumstances where bad performances such as 7 or even 20 are not rare to arise. See, however, discussions to follow. This exposition will also indicate that composite multipliers formed with three or more primitive roots will not be adequate as design.
Some accounts on geometrical aspects related to spectral tests will be in order at this end. Let Es denote the s dimensional Euclidean space. Spectral tests are based on the following fact. Take a consecutive s-tuple of the cyclic sequence, or more generally of the coset sequence as a point in Es with coordinates
(nzjnzj+1, . . . , nzj+s−1)≡nzj(1,z, . . . , zs−1).
These coordinates and their equivalents modulo d are notably linear combinations of the following with integer coefficients:
Certainly, nzje1 gives the noted s-tuple of the coset sequence with an integer nzj, and d translations along 2nd, 3rd, . . . , and sth coordinate axes are realized by adding suitable linear combinations of e2, e3, . . . , es with integer coefficients. As to the 1st axis, integer multiples of the following suffice to give any d translations:
de1ze2−z2e3− . . . −zs−les=(d,0,0, . . . , 0).
These structures are stated that points of consecutive s-tuples of the cyclic or coset sequences are in the lattice spanned by {e1, e2, . . . , es} of Es. Consider the hypercube Cd of sides [0, d) located at the origin of Es and ask how many lattice points are in the hypercube Cd. The answer comes from the coordinates of je1 with an integer j. Its first coordinate j can give d non-equivalent values 0, 1, . . . , d−1 in the hypercube Cd. All other coordinates are retracted back to give respectively unique values in Cs by adding a suitable linear combination of {e2, e3, . . . , es} with integer coefficients. Thus, the noted lattice has d distinct points in Cd. Similar constructions prove that any hypercube of sides d issuing from (k1d, k2d, . . . , ksd) with integers {k1, k2, . . . , ks} contains d lattice points. Spectral tests of s-dimension evaluate the geometrical distribution of these d points in Cd. A valuation closer to its lower limit 1 indicates that d lattice points are nearer to the geometrically ideal packing in Es; lattice points are then distributed in more homogeneous and isotropic manner, suggesting that consecutive s random numbers (which have coordinates of coset sequences divided by d and should be considered in the unit hypercube of Es) will have better distribution in view of the aimed independence and uniformity. We said suggesting, because random numbers of a primitive root multiplier for a prime modulus, of patent application Ser. No. 12/379,964 or of the present invention occupy about only half of these lattice points in their periods, and we are merely hoping that they appear to choose their seats randomly and in an unpredictable sequential manner. This is a point to be supplemented by statistical tests as performed in the far-reaching work of Fishman and Moore (1986), but their conclusions suggest that good performances in spectral tests will imply little problematic behavior in statistical tests.
We should close descriptions of the present invention with comments on a problem that still remain dim. Consider the 2 dimensional Euclidean space E2 and its unit square with sides of length 1. Take a simple image of the double precision real numbers in the sense of fixed point, with the smallest unit, say 2−53. The unit square is then divided into small cells (call them microcells) which are totally (253)2=2106 in number. If the modulus is chosen d≈274 for the multiplicative congruence random number sequence, then said lattice for s=2 has only d≈274 lattice points in the unit square, which are very small compared to the number of microcells of double precision. More generally in the s-dimensional space, the unit hypercube contains 253s microcells, but the noted multiplicative congruential sequence of random numbers can be distributed only on the half of lattice points, and the number of lattice points is fixed to d≈274 irrespective of the dimension s. We should stress that this is not a problem to be overcome by increasing the period of random numbers. However large this period may be chosen, properties of the random number sequence should be considered within the upper limit T of numbers usable in a simulation, and the heart of the problem is what kind of statistical properties this length T sequence would show as regards uniformity and independence. Our simulations might well be taken to use T≈274 or so of them, and the problem is the distribution of this amount of random numbers. The present invention and patent application Ser. No. 12/379,964 ensure that their output sequences give this order of magnitude of lattice points (seats) in the unit hypercube in Es and they may be tested of homogeneity and isotropy of their distribution. If this distribution of seats are nearly homogeneous and isotropic, then consecutive s-tuples of random number outputs for length T will have less reasons to be denied statistically of their uniform and independent distribution. It should of course be questioned as other problems of statistics whether or not the output s-tuples of truly uniform and independent random numbers can plausibly be on such homogeneous and isotropic distribution of seats in this period T≈273, or whether the particular order of appearance of random numbers generated by the proposed invention on these seats gives the disguise of uniform and independent distribution. Inventors feel that we cannot resort to an optimism, without tests, of expecting a random number sequence, say with its gigantic period >>T and theoretical assurance of beautiful equidistribution in its one period, to be able to have its (negligible) portion of this length T≈274 retaining the literal equidistribution. Experiences of sweeping spectral tests motivate us strongly to choose the security and the clarity ensured by them. Yet, questions will remain to be analyzed and tested in the future.
