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 popular 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:lihdeng@memphis.edu. Permissiontomakedigitalorhardcopiesofpartorallofthisworkforpersonalorclassroomuseis grantedwithoutfeeprovidedthatcopiesarenotmadeordistributedforprofitordirectcommercial advantageandthatcopiesshowthisnoticeonthefirstpageorinitialscreenofadisplayalong withthefullcitation.CopyrightsforcomponentsofthisworkownedbyothersthanACMmustbe honored.Abstractingwithcreditispermitted.Tocopyotherwise,torepublish,topostonservers, toredistributetolists,ortouseanycomponentofthisworkinotherworksrequirespriorspecific permissionand/orafee.PermissionsmayberequestedfromPublicationsDept.,ACM,Inc.,1515 Broadway,NewYork,NY10036USA,fax:+1(212)869-0481,orpermissions@acm.org. (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 primitive 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 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.

ACM transactions on modeling and computer simulation (January)

Similar ACM transactions on modeling and computer simulation (January)