ACM transactions on modeling and computer simulation (January)

107 Pages·2005·1.144 MB·English

Efficient and Portable Multiple Recursive
Generators of Large Order
LIH-YUANDENG
UniversityofMemphis
DengandXu[2003]proposedasystemofmultiplerecursivegeneratorsofprimemodulus pand
orderk,whereallnonzerocoefficientsoftherecurrenceareequal.Thistypeofgeneratorisefficient
because only a single multiplication is required. It is common to choose p = 231−1 and some
multiplierstofurtherimprovethespeedofthegenerator.Inthiscase,somefastimplementations
areavailablewithoutusingexplicitdivisionormultiplication.Forsucha p,DengandXu[2003]
providedspecificparameters,yieldingthemaximumperiodforrecurrenceoforderk,upto120.
One problem of extending it to a larger k is the difficulty of finding a complete factorization of
pk−1.Inthisarticle,weapplyanefficienttechniquetofindksuchthatitiseasytofactorpk−1,
with p=231−1.Thelargestonefoundisk=1597.Tofindmultiplerecursivegeneratorsoflarge
orderk,weintroduceanefficientsearchalgorithmwithanearlyexitstrategyincaseofafailed
search. For k = 1597, we constructed several efficient and portable generators with the period
lengthapproximately1014903.1.
Categories and SubjectDescriptors: F.2.1 [Analysis of Algorithms and Problems Complex-
ity]:NumericalAlgorithmsandProblems—Computationsinfinitefields;G.3[Mathematics of
Computing]:ProbabilityandStatistics—Randomnumbergeneration
GeneralTerms:Algorithms,Design,Performance
AdditionalKeyWordsandPhrases:DX-kgenerator,GMP,irreduciblepolynomial,linearcongru-
entialgenerator,MRG,primitivepolynomial.
1. INTRODUCTION
Multiple Recursive Generators (MRGs) have recently become the most popu-
lar random number generators. They compute the next value iteratively from
the previous k values. The k-th order recurrence equation corresponds to a
k-th degree primitive polynomial under a prime modulus p. When k = 1, the
MRG is reduced to Lehmer’s [1951] Linear Congruential Generator (LCG).
In Section2, we discuss two problems using the conditions of checking a
primitive polynomial given in Knuth [1998, page 30]: (1) the requirement
of the complete factorization of pk − 1, which can be difficult when p and
Author’saddress:DepartmentofMathematicalSciences,TheUniversityofMemphis,Memphis,
TN38152-3240;email:[e-mail].
Permissiontomakedigitalorhardcopiesofpartorallofthisworkforpersonalorclassroomuseis
grantedwithoutfeeprovidedthatcopiesarenotmadeordistributedforprofitordirectcommercial
advantageandthatcopiesshowthisnoticeonthefirstpageorinitialscreenofadisplayalong
withthefullcitation.CopyrightsforcomponentsofthisworkownedbyothersthanACMmustbe
honored.Abstractingwithcreditispermitted.Tocopyotherwise,torepublish,topostonservers,
toredistributetolists,ortouseanycomponentofthisworkinotherworksrequirespriorspecific
permissionand/orafee.PermissionsmayberequestedfromPublicationsDept.,ACM,Inc.,1515
Broadway,NewYork,NY10036USA,fax:[phone],[e-mail].
(cid:1)C 2005ACM1049-3301/05/0100-0001$5.00
ACMTransactionsonModelingandComputerSimulation,Vol.15,No.1,January2005,Pages1–13.
2 • L.-Y.Deng
k are large, and (2) no provision for an early exit strategy for nonprimitive
polynomials.
In Section 3, we describe our effort to find MRGs of large order k with p =
231−1. Factoring a general huge number (say, more than 200 digits) is hard,
particularlywhenithastwoprimefactorsofasimilarsize.Itissimpletofactor
anintegerifitcontains,atmost,onelargeprimefactor.Usingthispropertyand
thepolynomialfactorizationofxk−1,DengandXu[2003]foundthecomplete
factorization of pk − 1 for several (but not all) values of k ≤ 120. To avoid
the difficulty of factoring R(k, p) = (pk −1)/(p−1), one can consider other
values of p so that R(k, p) is a prime number. This idea was first used in
L’Ecuyeretal.[1993]fork ≤ 7,andlaterinL’Ecuyer[1999],fork ≤ 13.Deng
[2004] formally proposed a class of prime numbers of the form R(k, p), called
Generalized Mersenne Prime (GMP). This approach is efficient because of the
well-known fact that a primality check of a huge number is easier than its
factorization. The search of GMP is successful for k up to 4001. One can fix
p=231−1,andthensearchforaprimekforwhich R(k, p)isaGMP.However,
we do not find any such GMP for k < 25000. During this search process, we
findk = 47,k = 643,andk = 1597,suchthat R(k, p)canbefactoredbecause
itcontainsexactlyonehugefactor.AttheendofSection3,wediscussanearly
exitstrategyfornonprimitivepolynomials.
In Section 4, we extend the DX-k generators as proposed by Deng and Xu
[2003]. DX-k is a system of special MRGs where all nonzero coefficients of
the recurrence are equal, and k is the order of recurrence. Motivation for the
DX-k generator is dictated by the fact that a computer multiplication is more
time consuming than a computer addition/subtraction. Requiring a common
coefficient for a DX-k generator can achieve a high efficiency with a single
multiplication. By limiting the multiplier size, one can implement a portable
generator under various situations. If we choose the multiplier of the form
±2r±2w,thegeneratorcanbecomeevenmoreefficient.Evenwithsuchalimited
formofthemultipliers,wecanfindageneratorintheextendedDX-k class.
InSection5,wetabulatesomeportableandefficientDX-k(anditsextended
class) generators with p = 231 − 1, and k = 47, k = 643, and k = 1597.
Finally, in Section 6, we apply Deng’s [2004] automatic generation method to
constructmanymaximalperiodMRGswiththesameorderfromasingleDX-
1597generator.
2. MULTIPLERECURSIVEGENERATORS
2.1 MRGandPrimitivePolynomial
Amultiplerecursivegenerator(MRG),isdefinedas
Xi =(α1Xi−1+···+αkXi−k)mod p, i ≥k (1)
foranyinitialseeds(X0,...,Xk−1),notallofthembeingzero.Herethemodulus
p is a large prime number and X can be transformed using U = X /p. To
i i i
avoid the possibility of obtaining 0 or 1, Deng and Xu [2003] recommended
ACMTransactionsonModelingandComputerSimulation,Vol.15,No.1,January2005.
EfficientandPortableMRGofLargeOrder • 3
U = (X +0.5)/p. It is well-known that the maximum period of an MRG is
i i
pk −1,whichisreachedifitscharacteristicpolynomial
f(x)= xk −α xk−1−···−α , (2)
1 k
is a primitive polynomial. A set of necessary and sufficient conditions under
which f(x) is a primitive polynomial has been given in Alanen and Knuth
[1964]andKnuth[1998,page30]:
(i)(−1)k−1α mustbeaprimitiverootmod p. (3)
k
(ii) xR =(−1)k−1α mod (f(x), p),where R =(pk −1)/(p−1). (4)
k
(iii)Foreachprimefactorq of R,thedegreeof xR/q mod (f(x), p)
ispositive. (5)
A maximum period MRG is known to have the property of equidistribution
uptok dimensions:everym-tuple(1≤m≤k)ofintegersbetween0and p−1
appearsexactlythesamenumberoftimes(pk−m)overitsentireperiod pk−1,
withtheexceptionoftheall-zerotuplewhichappearsonetimeless(pk−m−1).
See,forexample,LidlandNiederreiter[1994,Theorem7.43].
2.2 NotationsandSimpleKnownFacts
Throughoutthisarticle,weuse f(x)todenoteak-thdegreepolynomial,given
in(2),whosecoefficientsarefromZ ,afinitefieldofintegerstakingvalues0to
p
p−1.Occasionally,asin(2),wewritethecoefficientas−a,whichequals p−a
in Z . We also use Z [x] to denote a set of all polynomials with coefficients in
p p
Z . An f(x) ∈ Z [x] is a monic polynomial when the coefficient of its leading
p p
termisone.
It is easy to check whether or not an integer is a primitive root mod p via
Knuth’s[1998]condition(i),givenin(3).Aninteger Aisaprimitiverootmod
p,ifthesmallestpositiveexponente suchthat Ae =1mod pise =(p−1).A
necessary and sufficient condition for an integer Ato be a primitive root mod
p is that A(p−1)/q (cid:5)= 1mod p for any prime factor q of (p−1). The number of
primitive roots in Z is φ(p−1), where φ(x) is the Euler “totient” function of
p
x, which is the number of integers between 1 and x that are relatively prime
to x.Thereareexactlyφ(pk −1)/k primitivepolynomialsofdegreek [e.g.,see
Knuth[1998]].
Forarealnumberx,weuse(cid:6)x(cid:7)and(cid:8)x(cid:9)todenotethefloorfunctionandthe
ceilingfunction,respectively.Thatis,(cid:6)x(cid:7)isthelargestinteger≤ x,and(cid:8)x(cid:9)is
thesmallestinteger≥ x.Ofcourse,when x isaninteger, x =(cid:6)x(cid:7)=(cid:8)x(cid:9).
Next,wediscusstheproblemswithconditions(ii)and(iii)givenin(4)and(5).
2.3 TwoBottlenecks
When k and p are large, there are two limiting factors to find a primitive
polynomialusingtheconditionsgiveninKnuth[1998]:
(1) Condition(iii)canbedifficulttocheckbecauseitrequiresacompletefactor-
izationof pk−1.Giventhecurrentstateoftechnology,itisextremelyhard
ACMTransactionsonModelingandComputerSimulation,Vol.15,No.1,January2005.
4 • L.-Y.Deng
TableI. Listofkandqin
R(k,p)=q×Hfor p=231−1
k q:smallcofactor log (pk−1)
10
47 123244717 438.6
643 7717 6000.4
1597 634021777 14903.1
to factor a general number with 200 decimal digits (or more), especially if
it contains only two prime factors of the similar size. For p = 231−1 and
k ≥22,wehave pk −1>10200.Wehavefoundacompletefactorizationof
pk −1foreveryintegerk ≤22,andforseveralother(butnotall)k ≤120.
However, it appears to be hard to extend to a much larger value of k for
p=231−1.
(2) Forcondition(ii),thereisnoearlyexitstrategyavailablewhen f(x)in(2)is
notaprimitivepolynomial.Whetherornot f(x)isaprimitivepolynomial,
thetimerequiredtoverifycondition(ii)isexactlythesame.Whilewecan
useanefficientmethodtoevaluatex(pk−1)/(p−1) mod f(x),itstillrequiresa
lotofcomputingtimeforalargek.Furthermore,thechanceoffindingak-
thdegreeprimitivepolynomialislessthan1/k.Therefore,whenkislarge,
a successful search for a primitive polynomial requires lots of computing
time.
Inthenextsection,wedescribesomesolutionstotheaboveproblems.
3. FINDINGMRGSWITHLARGEkANDp=231−1
3.1 Factoringpk−1,Withp=231−1
Inthisarticle,weareconcernedmainlywiththeproblemoffindingMRGsfor
p = 231 −1, which is a popular modulus because there are some advantages
fromthecomputingefficiencyviewpoint.
Foragivenprimek,Deng[2004]proposedthatwefindaprime psuchthat
R(k, p)=(pk −1)/(p−1) (6)
is also a prime number. Such R(k, p) is called a Generalized Mersenne Prime
(GMP). With p = 231 − 1 fixed, we search for a prime number k such that
R(k, p)isaGMP.Unfortunately,suchaneffortisnotsuccessfulfork <25000.
Tothebestofourknowledge,therearenoknownmathematicalresultsabout
theexistenceofsuchaGMP R(k, p)for p=231−1.Duringthesearchprocess,
however,wedidfindseveralk’ssuchthatR(k, p)canbeeasilyfactoredbecause
itcontainsexactlyonehugefactorlargerthanafixedconstant.Inpractice,we
choosetheconstant(relatively)smallsothatwecanquicklyremovesmallprime
cofactors(say,≤1,000,000,000=(cid:3))ofR(k, p).Inparticular,weconsideronly
R(k, p) = q× H, where q ≤ (cid:3), both q and H are prime numbers. Three such
k’sarefound,andtheyarelistedinTableI.
ForthevaluesofklistedinTableI,thecompletefactorizationof pk−1isthus
knownandwecanthentrytofindsomeprimitivepolynomialsofdegreek.The
MRGs with the corresponding primitive polynomials will have the maximum
ACMTransactionsonModelingandComputerSimulation,Vol.15,No.1,January2005.
EfficientandPortableMRGofLargeOrder • 5
periodof pk−1.AlthoughMRGsoforder47havea(relativelyspeaking)shorter
period length of about 10438.6, they can be found rather quickly using Knuth
[1998]conditions(i)–(iii)asgivenin(3)–(5).ForMRGsoforderk,oneneedsto
keepanarrayofsizektostorethepreviouskvalues.Ifthememoryrequirement
fork =1597istoohigh,thenonecanconsiderMRGsoforder643withperiod
length approximately 106000.4. If maximum period, or a higher dimension of
equidistribution property, is the major consideration, we use MRGs of order
k =1597,forwhichtheperiodlengthofsuchMRGisapproximately1014903.1.
For a large k, it is not efficient to verify a k-th degree primitive polynomial
usingtheconditionsgivenbyKnuth[1998].Wedescribeamoreefficientalgo-
rithmnext.
3.2 EarlyExitStrategy
Knuth’scondition(ii),asgivenin(4),isnotveryefficientbecausetherequired
computing time is exactly the same regardless of whether f(x) is a primitive
polynomial or not. There is no early exit strategy when f(x) is not a primi-
tive polynomial. Deng [2004] proposed to replace Knuth’s condition (ii) with a
new condition that can be checked more quickly when f(x) is not a primitive
polynomial:
(ii)Initially,let g(x)= x.Fori =1,2,3,...,(cid:6)k/2(cid:7),do
(7)
(a) g(x)= g(x)p mod f(x);
(b) d(x)=gcd(f(x), g(x)−x);
(c) ifd(x)(cid:5)=1,then f(x)cannotbeaprimitivepolynomial.
This new condition is Algorithm 2.2.9 (irreducibility test) in Crandall and
Pomerance [2000, page 88] to test whether f(x) is an irreducible polynomial.
At the first iteration (i = 1), it is checking whether f(x) has any linear factor
(x−a)whichcanbedonequicklybecause
(cid:1)p−1
g(x)−x = xp−x = (x−a)mod p.
a=0
To check whether f(x) has any linear factor, we simply check whether or not
d(x) = gcd(f(x),xp − x) is a polynomial of degree ≥ 1. In general, at the i-
thiteration,itischeckingwhether f(x)hasanyirreduciblepolynomialfactor
of degree i. Again, it is based on the fact that xpi − x is the product of all
monicirreduciblepolynomialswithdegreesthatcandividei.See,forexample,
LidlandNiederreiter[1994],CrandallandPomerance[2000].Indeed,thetotal
numberofmonicirreduciblepolynomialsofexactdegreek is
(cid:2)
1
N (p)= pdµ(k/d),
k
k
d|k
where µ(n) is the Mo¨bius function: µ(n) = 0, if n is not square-free, µ(n) = 1,
if n is square-free with an even number of prime factors, µ(n) = −1, if n is
ACMTransactionsonModelingandComputerSimulation,Vol.15,No.1,January2005.
6 • L.-Y.Deng
square-free with an odd number of prime factors. In addition, N (p) ≈ pk/k
k
becauseofthelasttermwhered =kinthesummationof N (p)isdominating.
k
The search program based on this early exit strategy can be fast because
a failed search will be likely to terminate from the i-th loop when i is small.
Therefore, a lot of computing time can be saved. Indeed, using this early exit
strategy, we are able to find many primitive polynomials of degree 1597. The
correspondingMRGscanbeeasilyconstructedwiththeperiodlengthapprox-
imately 1014903.1. The next problem is to find generators in a class of efficient
andportableMRGs.
4. EFFICIENTANDPORTABLEGENERATOR
4.1 EfficientMRG
Whenk islarge,ageneralMRGmaybelessefficientbecauseitneedsseveral
multiplications, where an LCG needs only one multiplication. To improve the
efficiencyofMRGs,manyauthorsconsideredonlytwononzerocoefficients,α
j
andα (1≤ j <k),oftheMRGin(1).Forexample,see,Grube[1973],L’Ecuyer
k
andBlouin[1988],L’Ecuyeretal.[1993].
Deng and Lin [2000] proposed a Fast MRG (FMRG) which has a slightly
simplerformandrequiresasinglemultiplication.ExtendingtheideaofFMRG,
Deng and Xu [2003] introduced a class of DX-k-s generator which is a special
MRG in (1) that has s nonzero coefficients, all of them being equal, and the
nonzero coefficients indices are about k/(s−1) apart. However, if we further
limitchoices/typesofnonzerocoefficients,theremaynotexistsuchgenerators,
especiallywhenkislarge.OnesimplesolutionistoexpandtheDX-k-sclassto
theDX-k-s-t classgenerators:
(1) DX-k-1-t (α =1,α = B),1≤t <k.
t k
Xi = Xi−t +BXi−k mod p, i ≥k. (8)
(2) DX-k-2-t (α =α = B),1≤t <k.
t k
Xi = B(Xi−t +Xi−k)mod p, i ≥k. (9)
(3) DX-k-3-t (αt =α(cid:8)k/2(cid:9) =αk = B),1≤t <(cid:8)k/2(cid:9).
Xi = B(Xi−t +Xi−(cid:8)k/2(cid:9)+Xi−k)mod p, i ≥k. (10)
(4) DX-k-4-t (αt =α(cid:8)k/3(cid:9) =α(cid:8)2k/3(cid:9) =αk = B),1≤t <(cid:8)k/3(cid:9).
Xi = B(Xi−t +Xi−(cid:8)k/3(cid:9)+Xi−(cid:8)2k/3(cid:9)+Xi−k)mod p, i ≥k. (11)
Notice that there are slight changes in the indices of nonzero coefficients
from the original definition of Deng and Xu [2003]. We have introduced an
additional parameter t to slightly expand the class of generators. In addition,
weusetheceilingfunction(cid:8)x(cid:9)toreplacethefloorfunction(cid:6)x(cid:7)inEquations(10)
and (11). In Deng and Xu [2003], the values of k of main interest are k = 102
and k = 120, where both definitions yield the same result. The parameters
given in the Tables of Deng and Xu [2003] and Deng [2004] are unaffected.
Oneadvantagethismodificationisthatitmayworkbetter,evenforverysmall
value of k. For example, if k = 3 and s = 3, then (cid:8)k/2(cid:9) = 2, and the DX-3-3
ACMTransactionsonModelingandComputerSimulation,Vol.15,No.1,January2005.
EfficientandPortableMRGofLargeOrder • 7
generator(witht =1)in(10)isreducedto
Xi = B(Xi−1+Xi−2+Xi−3)mod p, i ≥3. (12)
Marsaglia[1996]proposedtheabovegeneratorwith B =210, p=232−5,and
B = 220, p = 232−209. The value of k in Equation (12) is only 3, which is too
small.
FollowingDengandXu[2003],thegeneratorsinEquations(8)–(11)arere-
ferredtoasDX-k-s-t generators.Here,sisthenumberoftermswiththesame
coefficient B, and t is the smallest index j for which α = B. This expanded
j
system of generators is still referred to as DX-k generators. Clearly, DX-k-s is
a special case with t = 1. FMRG is another special case with s = 1 and t = 1.
By allowing t to change, we can expand the search space for the desired form
of Basdiscussednext.
4.2 PortabilityIssue
We prefer a generator to be portable so that it produces exactly the same se-
quence of random numbers, regardless of which computing platform is used.
The problem is that there is no standard way of storing a 64-bit number in a
32-bit computer. According to L’Ecuyer [1997], a necessary (but not su(cid:3)fficient)
conditionfora“good”MRGisthatthesumofthesquaresofcoefficients, k α2,
i=1 i
shouldbelarge.ForDXgenerators,thismeanssand B needtobeaslargeas
possible, while maintaining the efficiency and portability property. It is com-
mon to impose certain limits on the size of B so that the exact result of the
multiplicationcanbeproduced:
√
(1) Use integer operations with B < p. This is useful if one wishes to use
onlyintegeroperationsina32-bitcomputerwiththetechniquesproposed
by Payne et al. [1969] and L’Ecuyer [1988]. The problem is that the upper
limit(≤46340,for p≤231−1)israthersmall.
(2) FollowIEEEdoubleprecisionstandard,anduseupperlimitsfor B:
B<2d, whered =20, whens=1,2; d =19, whens=3,4. (13)
Thisupperboundisabout10to20timeslargerthantheonefoundusingthe
firstapproach.Thereareportableimplementationsforvariouscomputing
platforms.
(3) Useasystem-dependent64-bitintegerdatatypewhichisavailableinmany
popularcompilersin32-bitcomputersystems.Forexample,a64-bitinteger
datatypeisavailableinMS-C,with__int64,andinGNU-C,withint64_t.
Inthiscase,weuseauniformupperlimitfor Bas
B<230, fors≤4. (14)
Without using 64-bit data types/operations, a portable implementation of
MRGscanbeusedattheexpenseofslightgeneratinginefficiency.Forthis
B,wefindthesmallestpositiveintegeresuchthat
B+e× p=C ×C , where0<C ,C <219, (15)
1 2 1 2
where the upper bound for C ,C is chosen as in Equation (13). The main
1 2
problemisthatweneedtwomultiplicationsinsteadofonemultiplication.
ACMTransactionsonModelingandComputerSimulation,Vol.15,No.1,January2005.
8 • L.-Y.Deng
TableII. ValuesofBforDX-k-sGenerators
√
k s (a)minB (b)B< p (c)B<2d (d)B<230 (e) C C
1 2
47 1 408 46193 1047811 1073741241 4 59827 161527
47 2 114 46312 1048423 1073741549 10 113749 198231
47 3 278 46316 524193 1073741553 19 194182 215653
47 4 14 46281 524254 1073741381 3 51958 144659
643 1 13082 42720 1047252 1073741462 0 29338 36599
643 2 1542 39963 1048207 1073733589 5 52321 225744
643 3 2095 42188 522057 1073741583 0 6159 174337
643 4 1679 42938 521553 1073740543 1 30559 105410
1597 1 3233 19661 1036675 1073740265 2 10891 492949
1597 2 4621 43253 1041405 1073730399 0 12999 82601
1597 3 9897 40507 520021 1073738734 0 30931 34714
1597 4 1854 44875 512675 1073741362 0 29746 36097
AnotherapproachistoconsideraspecialformofthemultiplierBforDX-k-s-t
generatorswith p=231−1:
B=δ(r)2|r|+δ(w)2|w|, |w|<|r|, (16)
where δ(x) = +1, when x ≥ 0; δ(x) = −1, when x < 0. A multiplier of this
form,calledpowers-of-twodecomposition,wasfirstsuggestedbyWu[1997]for
LCGs. It can result in a fast computation by using only shift and addition op-
erations. Another advantage of this approach is that we can find a large B,
without the limitation given in (13), while maintaining its portability. How-
ever, L’Ecuyer and Simard [1999] pointed out that LCGs with parameters of
thepowers-of-twoformhavebadstatisticalpropertiesbecausetherecurrence
does not “mix the bits” well enough. While there is no evidence (yet) that the
successivevaluesintheoutputoftheDX-k-s-t generators(especiallyforlarge
s)havesuchweakness,wewouldwelcomemoreextensivestudyorevaluation
forsuchgenerators.
Sincewepreferalargersizeof B,wemakearestrictionon|r|for B in(16)
that20≤|r|≤30.Becausethenumberofchoicesfor B=±2|r|±2|w| israther
limited, we need to expand our search space, which is the main motivation
behindtheconsiderationoftheDX-k-s-t generators.
5. LISTOFMRGSWITHLARGEORDERk
5.1 DX-k-sGeneratorsfork =47,k =643andk =1597
With the complete factorization of pk −1, as given in Table I, we then search
for DX-k class for k = 47, k = 643, and k = 1597. For each DX-k-s generator
√
with s = 1,2,3,4, we perform a complete search for 0 < B < p. For other
widerangesof B’s,acompletesearchisnotfeasible.Forsimplicity,wewillonly
√
report(a)minB,(b)maxB< p,(c)maxB<2d asin(13),and(d)maxB<230
asin(14).Sincethesearchrangefor Bislarge,wecanalwaysfindagenerator
insidetheDX-k class.TableIIgivesvariousDX-k-sgenerators.
Thevaluesunder“minB”incolumn(a)arenotrecommended,buttheycanbe
√
usefultodeterminethe“power”ofempiricaltests.Thevaluesunder“B< p”
in column (b) are useful if we use only integer operations. The values under
ACMTransactionsonModelingandComputerSimulation,Vol.15,No.1,January2005.
EfficientandPortableMRGofLargeOrder • 9
TableIII. Listingofk,s,t,(r,w)andB=δ(r)2|r|+δ(w)2|w|forDX-k-s-tGenerators
k s t (r,w) B t (r,w) B
47 1 1 (−30,4) −1073741808 1 (29,−10) 536869888
47 2 1 (−29,−9) −536871424 1 (26,−21) 65011712
47 3 1 (29,−26) 469762048 1 (−27,−19) −134742016
47 4 1 (29,−6) 536870848 1 (26,−7) 67108736
643 1 8 (−29,1) −536870910 10 (−28,11) −268433408
643 2 5 (25,13) 33562624 2 (25,11) 33556480
643 3 2 (29,10) 536871936 8 (27,−16) 134152192
643 4 7 (29,-26) 469762048 8 (26,−3) 67108856
1597 1 13 (−27,−4) −134217744 14 (−26,6) −67108800
1597 2 4 (30,−19) 1073217536 1 (29,−11) 536868864
1597 3 15 (−24,16) −16711680 35 (22,0) 4194305
1597 4 3 (29,8) 536871168 3 (28,11) 268437504
“B < 2d” in column (c) are recommended for general purpose generators. If
64-bitintegerdatatypeanditsfastoperationsareavailable,thevaluesunder
“B < 230” in column (d) are recommended because of their large multiplier
values.Incolumn(e),theconstantse,C ,C ,describedin(15),areusefulfora
1 2
portableimplementationfor B,listedincolumn(d).
One advantage of using the parameters given in column (d) is their large
multipliers, which are at least 1000 times larger than those in column (c).
Thedisadvantageisthat64-bitintegersaresystem-dependent,andthe64-bit
operationscanbeslow.If64-bitintegerdatatypeisnotavailable,wecanusea
slower,butportableimplementation:foreachBincolumn(d),wefindconstants
e,C ,C fromcolumn(e),andperformtherequiredmultiplicationasin(15).
1 2
5.2 DX-k-s-tGeneratorsWithB=±2|r|±2|w|
For more efficient generators, we consider B = ±2|r| ±2|w| within the class of
DXgenerators.Asmentioned,withrestrictedformandsizefor B,wemayhave
tosearchintheexpandedDXclass.Fork =47,wefindmanyvaluesofBinside
theDX-47-sclass.Fork =643,ork =1597,wesearchforDX-k-s-t generators
from t = 1 and beyond, until a number of generators is found. Table III gives
thetoptwoDX-k-s-t generatorsfork =47,k =643,andk =1597.
5.3 ImplementationandEmpiricalTesting
Since only minor changes were made to the definition of DX-k-s and DX-k-s-t
generators,modificationstotheoriginalCcodeareeasy.TheCcodeisavailable
attheauthor’swebsitehttp://www.cs.memphis.edu/˜dengl/dx-rng/.
L’EcuyerandTouzin[2004]conductatheoreticalanalysisofthelatticestruc-
tureofcertainMRGs.MostoftheirconclusionsareapplicabletoanyMRGswith
few,smallnonzerocoefficients,notjustspecificallyforFMRGsorDX-kgenera-
tors.TheythenperformempiricalstatisticaltestsonsomeFMRGs(aslistedin
DengandLin[2000])andsomeDX-kgenerators.Theirconclusionisconsistent
with the main motivation behind DX-k generators: taking only a few nonzero
coefficients with small values is a bad idea. For FRMGs, they study empiri-
cally several generators with k = 2,3,4 as in Deng and Lin [2000]. For DX-k
ACMTransactionsonModelingandComputerSimulation,Vol.15,No.1,January2005.
10 • L.-Y.Deng
generators, they study the cases of k = 102, s = 2, B = 23, and B = 45787,
where B and s are not the largest that we would recommend. L’Ecuyer and
Touzin[2004]reportthatbothDX-102-2generatorsstudied(evenwithB=23)
passallthegeneralpurposeteststheytried.Inthesamestudy,theyfindthat
FRMGs(hence,DX-kgenerators)arereasonablyfastcomparedwithotherwell-
known generators. Additionally, we perform some predefined batteries of sta-
tistical tests using a well-known test suite, TestU01, whose source code and
user’smanualareavailablefromhttp://www.iro.umontreal.ca/˜lecuyer/.There
are38DXgeneratorslistedinTableIofDengandXu[2003],and72DXgener-
atorslistedinTableIIandTableIIIofthisarticle.Excluding20DXgenerators
with the smallest B (under column minB), all other generators have passed
the “crush” test which is a predefined battery of stringent tests. None of the
p-values of statistical tests on these DX generators are less than 10−7. While
we do not recommend DX generators with the smallest B, all of them indeed
havealsopassedthesametest.
No random number generator can claim to be a perfect generator. One can
createaspecificsubsequencewithabadlatticestructure,oraspecializedtest
to fail a generator with a known generating algorithm. While DX-k genera-
tors have a property of equidistribution up to dimension k, one can choose a
subsequencewithabadlatticestructure.Forexample,foraDX-k-sgenerator
with s = 1, or s = 2, L’Ecuyer and Touzin [2004] show that the subsequence
Sk = {(Ui,Ui+k−1,Ui+k),i = 0,k+1,2(k+1),...}doesnothaveagoodlattice
structure. As noticed by L’Ecuyer and Touzin [2004], this subsequence S can
k
beobtainedbytakingthefirstoutputvaluefromthegenerator,skippingk−2,
taking another 3, skipping k−2, taking another 3, skipping k−2, and so on.
Theyshowthatthesubsequence,S ,failed,spectacularly,verysimplestatis-
102
tical tests for the DX-102-2 generators studied, especially for B = 23. For the
subsampling scheme corresponding to S , we can choose DX-k-s-t generators
k
witht > 1,ors > 2sothat S doesnothaveabadlatticestructure.However,
k
foranygivenDX-k-s-tgenerator,onecanfindotherspecificsubsequencewitha
badlatticestructure,iftheparametersk,s,t areknown.Abetterstrategyisto
choosealargeorderk,larges,andlargevaluesof Basproposedinthisarticle.
Asusual,oneshouldnotrelyonlyonasinglegeneratorforanimportantlarge
scalesimulationstudy.
6. CONSTRUCTIONOFMRG-k-sGENERATORS
6.1 AutomaticConstructionMethod
WerefertoMRG-k-sasaclassofmaximalperiodMRGswithsnonzerocoeffi-
cientsintheircorrespondingcharacteristicpolynomial.WhenR(k, p)in(6)isa
GMP,Deng[2004]proposedasimplemethodtoconstructMRG-k-sgenerators
fromthecharacteristicpolynomialofaDX-k-sgenerator.Evenif R(k, p)isnot
aGMP,mostofhisapproachcanstillbeuseful.
We know that an irreducible polynomial is a necessary condition, but not a
sufficient condition, for a primitive polynomial. Furthermore, if f(x) in (2) is
a primitive polynomial, then f(cx) and xk f(c/x) must also be an irreducible
ACMTransactionsonModelingandComputerSimulation,Vol.15,No.1,January2005.