多階質數乘餘法之全週期(I)

行政院國家科學委員會專題研究計畫 成果報告

多階質數乘餘法之全週期(I)

研究成果報告(精簡版)

計 畫 類 別 ： 個別型 計 畫 編 號 ： NSC 95-2221-E-151-032- 執 行 期 間 ： 95 年 08 月 01 日至 96 年 07 月 31 日 執 行 單 位 ： 國立高雄應用科技大學工業工程管理系 計 畫 主 持 人 ： 唐惠欽 計畫參與人員： 碩士班研究生-兼任助理：徐茂為 大學生-兼任助理：黃士科 講師級-兼任助理：陳琪琪 處 理 方 式 ： 本計畫涉及專利或其他智慧財產權，1 年後可公開查詢

中 華 民 國 96 年 10 月 23 日

行政院國家科學委員會補助專題研究計畫

5

成 果 報 告

□期中進度報告

多階質數乘餘法之全週期(I)

計畫類別：

5

個別型計畫 □ 整合型計畫

計畫編號：NSC

95-2221-E-151-032

執行期間：

95 年 08 月 01 日起至 96 年 07 月 31

日

計畫主持人：

唐惠欽

共同主持人：

計畫參與人員：

成果報告類型(依經費核定清單規定繳交)：

5

精簡報告 □完整報告

本成果報告包括以下應繳交之附件：

□赴國外出差或研習心得報告一份

□赴大陸地區出差或研習心得報告一份

□出席國際學術會議心得報告及發表之論文各一份

□國際合作研究計畫國外研究報告書一份

處理方式：除產學合作研究計畫、提升產業技術及人才培育研究計畫、

列管計畫及下列情形者外，得立即公開查詢

5

涉及專利或其他智慧財產權，

5

一年□二年後可公開查詢

執行單位：

高雄應用科技大學工業工程與管理系

中 華 民 國 96 年 10 月 22 日

一、中英文摘要

1.1 中文摘要

關鍵詞：全週期、一般梅森質數、多階質數乘餘法、亂數 一有效率檢查多階質數乘餘法亂數產生器之全週期為質數階數、質數模數及質數 ，即為一般梅森質數。模數為梅森質數或Sophie Germain質數之多階質 數乘餘法，則具計算有效性及週期性。由三個引理可設計五個實驗以搜尋較大模數及階數 使得r為質數，結果中最大週期為當模數 且階數 ) 1 /( ) 1 ( − − = m m r k 69 231₋ = m k =4217，其週期為 ， 此多階質數乘餘法亂數產生器符合現今模擬所需之長週期。 7482 . 39352 10

1.2 英文摘要

Keywords: Full period, Generalized Mersenne Prime, Multiple recursive generator, Random

number.

An interesting special case of efficiently checking the full period conditions of a kth-order multiple recursive generator (MRG) is that of a prime order k, a prime modulus m and a prime

11 −−

m

mk which is called a Generalized Mersenne prime (GMP). A MRG with either Mersenne prime

modulus or Sophie Germain prime modulus is a way of improving on computational efficiency and full period, respectively. Computational experiments are conducted to search a prime order k, for a given prime modulus m , such that ₁1

−−

m

mk is a prime number. We find a GMP for _k =4217

and ₌231₋69. Therefore, a MRG with order

m k =4217 and modulus can

efficiently check the full period conditions and its full period length is .

69 231₋ = m 7482 . 39352 10

二、報告內容

2.1 前言

多階質數乘餘法亂數產生器之遞迴關係式為 2 2 1 1 − + − ≡ _{n n} n a X a X X +...+a_kX_n−k(modm) 其 中 k 為 階 數 ， 模 數 m 為 質 數 ， 乘 數 為 介 於-m至m之間的整數，起始狀態 為介於0至m-1之間的整數，但不可全為0。 k a a a₁, ₂,..., k X X X₁, ₂,..., 一個理想的多階質數乘餘法亂數產生器必須具有長週期、統計隨機性和電腦有效性， 本 研 究 乃 探 討 多 階 質 數 乘 餘 法 之 週 期 性 ， 此 時 我 們 將 面 臨 三 個 問 題(1) 因 數 分 解 及(2)計算 及(3)計算 ，其中 為多 ) 1 /( ) 1 ( − − = m m

r k xr (modf(x)) xr/q(mod f(x)and m) f(x)

階質數乘餘法之特徵函數且q為r之所有質因數。若 為質數，則問題(1)變 為較簡單，此乃因質數測試較因數分解簡單。且問題(3)將自動滿足，由此可知本研究將搜 尋模數m及階數k使得 為質數。 ) 1 /( ) 1 ( − − = m m r k ) 1 /( ) 1 ( − − = m m r k

2.2 研究目的

由於滿足全週期之多階質數乘餘法亂數產生器非常多，所以全週期為一理想亂數產生 器之必要條件。另一方面，由於實際應用常需長週期，所以一理想亂數產生器之模數及階 數皆需非常大，由此可知本研究之目的為搜尋模數m及階數k使得 為質數 且其週期 非常大。 ) 1 /( ) 1 ( − − = m m r k ) 1 ( k − m

2.3 文獻探討

多階質數乘餘法之重要相關書籍為Gentle(2003)、Fishman (1996)、Knuth (1997)、Law 和Kelton(2000)及Niederreiter (1992)，而重要相關文獻可參考Deng(2005)、L'Ecuyer和Couture (1997)、Tang (2006)、Tang和Kao(2004)、Sezgin(2004)。與本研究相關文獻有L’Ecuyer et al. (1993)發現階數及模數分別為( mk, )=(5, )及(7, )，後來L’Ecuyer(1999)搜尋 模數及階數使得 61 231_{− 2}31₋₁₅₁ 2−1 m _及 1 1 −− m mk _{為質數，發現中最大週期為 ，其中階} 數 及模數 。在給定一質階數k，Deng(2002)搜尋模數m使得 8614 . 134 7 64 _{51149) 1 10} 2 ( − − = 7 = k m=264−51149 m_mk₋−₁1_為質數， 結果中最大週期為階數 ，模數 且其週期 。與本研究相關有 限體及模數可參考Niederreiter(1992)及Riesel(1994)。 1511 = k m=231−55719 1014100.529

2.4 研究方法

Knuth(1997)指出滿足全週期k階質數乘餘法之充分必要條件為 (i) 對於所有(m−1)之質數q，(( 1) ) 1(mod ) 1 1 m a q m k k ≠ − − − _；

(ii) (mod ( )) ( 1) 1 (mod ) ， 其 中 且

m a x f xr ≡ − k− _k r=(mk −1)/(m−1) k x x f( )= −a₁xk−1−a₂xk−2 −_L−a_k−1x−a_k； (iii) 對於所有r之質數q， / (mod ( ) )之次方大於零。 m and x f xr q 理論上，滿足全週期的乘數向量組有φ(mk −1)/k，其中φ(.)是尤拉函數。 本研究主要工作為在給定一質數模數m，搜尋階數k，使得 為質數。 我們將利用下列三個引理，使得搜尋階數k更有效率。 ) 1 /( ) 1 ( − − = m m r k 引理1：若r =(mk −1)/(m−1)為質數，則k為質數。

引理2：若 =2p −1，p為質數且k為偶數，則 m w= m+1可整除r=(mk −1)/(m−1)。 證 明 ：k=2時 ， 因 為 wm+1=w， 所 以w 成 立 。 若r w 滿 足 k 為 偶 數 ， 則r 時 ， ，由數學歸納假設可知 2 + k 1 ) 1 ( 1 ) 1 /( ) 1 ( +2_{− − =} +1_{+ +} −1_{+ + + = + +} −1_{+ + +} m m m m m m m m m mk k k k _K k k _K =mkw+mk−1+ 1 + + m K −1_{+ + +}1 m m w k _K ，所以 +1 _{+ +} −1_{+ + +}1 m m m m w k k k _K 。□ 引理3：若 =2p −1，p為質數且k為偶數，若 m

∑

− = − −− = + + + = − = 2 0 1 11 1 ( 1) k i i i i k m m _{m m w aw} r k _K 且 ， 則 且 1 2 = − k a

∑

= − _{+ + + = −} + k i i i i k k w b m m m 0 1 _{1 ( 1)} K b₀ =bk =1 ， b1=a1 ， ， ， ， 。 2 1 − − + = i i i a a b 2 , , 2 , 1 − _K − =k k i ( 1) 1 2 0 = −

∑

− = k i i i a ( 1) 1 0 = −

∑

= k i i i b 證 明 ： 時 可 知 其 結 果 成 立 ， 若 滿 足 時 。 考 慮 時 ， 若 ，則 ， 其 中 2 = k k >2 k+2

∑

= + _{+ + + + =} k ₋ i i i i k k w a w m m m 0 1 _{1 ( 1)} K +2₊ +1_{+ + +}1₌ ( +1_{+ + + +}1)₊1 m m m m m m mk k _K k k _K 1 ) 1 ( ) 1 ( 1 ) 1 ( ) 1 ( 0 1 0 2 0 + − − − = + − − =

∑

∑

∑

= + = + = k i i i i k i i i i k i i i i w a w a w a w w =

∑

− − + = − 2 2 2 ) 1 ( k j j j j w a

∑

∑

∑

∑

+ = + = − + = − + = − + _{+ = − + − + = −} − 2 1 1 1 1 2 2 2 1 1 1 1 _{1 ( 1) ( 1) 1 ( 1)} ) 1 ( k j j j j k j j j j k j j j j k j j j j w b w a w a w a b₀ =b_k+2 =1， ， ， 1 1 a b = b_j =a_j−1+a_j−2 j=k+1 K,k, ,2 。 若 替 代 w=1 且 m= w−1=0 於 及 ，可得 及 。□

∑

− = − _{+ + + =} 2 ₋ 0 1 ₁ k _{( 1)} i i i i k w a w m m _K

∑

= − _{+ + + = −} + k i i i i k k w b m m m 0 1 _{1 ( 1)} K 2( 1) 1 0 = −

∑

− = k i i i a 1 ) 1 ( 0 = −

∑

= k i i i b

2.5 結果與討論

本 研 究 將 以Mathematica 軟 體 (1999) 來 執 行 搜 尋 質 數 模 數 m 及 階 數 k ， 使 得 為質數，其實驗共分五個部份。 ) 1 /( ) 1 ( − − = m m r k 第一實驗為當 ，找比 小的最大質數階數k，使得 為質 數，此部份因產生溢位問題，所以無法執行。 1 231₋ = m 231 r=(mk −1)/(m−1) 第二實驗為當 ，找第一個最小質數階數k，使得 為質數， 執行至今於 ，皆無法找到階數k，使得 為質數。另一方面， 當 ，測試 是否為質數所需電腦時間超過10小時，所以本實驗將 不再繼續執行。 1 231₋ = m r=(mk −1)/(m−1) 12893 2≤ k ≤ r=(mk −1)/(m−1) 10000 ≥ k r =(mk −1)/(m−1) 第三實驗為取模數為比231小的100個最大質數，即 ₌231₋1,231₋19, ,231₋2295，找 K m 3

質數階數k， 使得 為質數。結果如表一，由表一可知(1)有35 個模數，可找到質數階數k，使得 為質數。(2)最大的模數為 且 階數 使得 為質數。(3)最小的模數為 且階數 使得 為質數。(4)最大的週期是當模數為 且階數 ，此時週 期為 。 1000 1≤ k≤ r =(mk −1)/(m−1) ) 1 /( ) 1 ( − − = m m r k m=231−61 5 = k r=(m5−1)/(m−1) m=231−2295 k =53 ) 1 /( ) 1 ( 53_{− −} = m m r m=231−699 k =809 1 ) 699 2 ( 31₋ 809₋ 由 前 一 實 験 之 結 果 中 ， 可 知 最 大 的 模 數 為 及 最 大 的 週 期 之 模 數 為 ，第四實驗為針對此二模數，擴大質數階數k之範圍， 61 231₋ = m 699 231₋ = m 1≤ k≤5000，找更大的 質數階數k使得 為質數，結果為找不到階數k使得 為質 數。 ) 1 /( ) 1 ( − − = m m r k r=(mk −1)/(m−1) 第五實驗為擴大質數階數k之範圍為2≤ k ≤10000，模數限定為下列兩大類(1)比 小 的5個最大Sophie Germain質數，即m及 31 2 2−1 m _{皆為質數，其模數為 ， ，} ， 及 。(2)比 小的5個最大質數但不是Sophie Germain質數， 其 模 數 為 ， ， ， 及 。 找 更 大 的 質 數 階 數 k 使 得 為質數。結果如表二，由表二可知(1)10個模數中，只有3個模數，即 ， 及 ，可找到質數階數k，使得 為質數，其 分 別 為( , 5) ，( , 4217) 及( , 13)。 (2) 最 大 的 週 期 是 當 模 數 為 且階數 ，其週期為 。 69 231_{− 2}31₋₅₂₅ 885 231_{− 2}31_{−2085 2}31_{−2505 2}31 1 231_{− 2}31_{−19 2}31_{−61 2}31_{−85 2}31₋₉₉ ) 1 /( ) 1 ( − − = m m r k 61 231_{− 2}31_{−69 2}31_{−2505 =( −1)/( −1)} m m r k ( km, ) 61 231_{− 2}31_{−69 2}31₋₂₅₀₅ 69 231₋ = m k =4217 mk −1=(231−69)4217−1=1039352.7482 所以本研究不僅提供模數為 之4217階質數乘餘法亂數產生器以符合現今模 擬所需之長週期 ，且由實驗可知模數為Sophie Germain質數是未來可進一步研究之 主題。 69 231₋ = m 7482 . 39352 10 表一：模數為比₂31小的100個最大質數且階數_{1≤ k≤1000} 使得 為質 數 ) 1 /( ) 1 ( − − = m m r k 第i個最大質數 模數m 階數k 3 231₋61 ₅ 7 231₋105 ₁₅₁ 8 231₋151 _7、17 9 231₋159 ₁₁₃ 18 ₂31−_{469 379} 25 ₂31−_{595 769} 29 ₂31−_{699 101、809} 31 ₂31−_{711 79} 33 ₂31−_{771 31、107} 34 ₂31−_{775 23、439、617} 36 ₂31−_{789 11、29} 37 ₂31−_{829 37}

38 ₂31−_{831 31、463} 39 ₂31−_{837 157} 42 ₂31−_{909 19} 46 ₂31−_{985 5、193} 48 ₂31−_{1027 97} 49 ₂31−_{1057 263} 50 ₂31−_{1065 11、127} 51 ₂31−_{1071 31} 52 ₂31−_{1141 31} 53 ₂31−_{1147 53、641} 58 ₂31−_{1287 3} 60 ₂31−_{1305 541} 63 ₂31−_{1375 67} 66 ₂31−_{1425 467} 67 ₂31−_{1527 13} 71 ₂31−_{1585 7、71} 77 ₂31−_{1741 5} 78 ₂31−_{1747 11} 79 ₂31−_{1749 5} 82 ₂31−_{1785 43、263} 95 ₂31−_{2161 11} 97 ₂31−_{2281 23} 99 ₂31−_{2295 53}

表二：比 小的5個最大Sophie Germain質數及5個最大質數但不是Sophie Germain質數，搜

尋階數 k使得 為質數。 31 2 ] 10000 , 1 [ ∈ k r=(mk −1)/(m−1) 最大質數m 階數k Sophie Germain質數m 階數k 1 231_{− 2}31_{−69 4217} 19 231_{− 2}31₋₅₂₅ 61 231_{− 5 2}31₋₈₈₅ 85 231_{− 2}31₋₂₀₈₅ 99 231_{− 2}31_{−2505 13}

參考文獻

1. L.Y. Deng, Generalized Mersenne Prime Number and Its Application to Random Number Generation, in Monte Carlo and Quasi-Monte Carlo Methods 2002 (H. Niederreiter, ed.), Springer-Verlag, 167-180, (2004).

2. L.Y. Deng, Efficient and portable multiple recursive generators of large order, ACM

Transactions on Modeling and Computer Simulation 15, 1-13, (2005).

3. G.S. Fishman, Monte Carlo: Concepts, algorithms, and applications, Springer Series in

Operations Research, Springer-Verlag, New York, (1996).

4. J.E. Gentle, Random number generation and Monte Carlo methods, 2nd edition,

Springer-Verlag, (2003).

5. D.E. Knuth, The art of computer programming Vol. 2: Semi-numerical Algorithms, 3rd

edition, Addison-Wesley, Reading MA, (1997).

6. A.M. Law and W.D. Kelton, Simulation modeling and analysis, 3rd edition, McGraw-Hill, New York, (2000).

7. P. L'Ecuyer, Good parameters and implementations for combined multiple recursive random number generators, Operations Research 47, 159-164, (1999).

8. P. L'Ecuyer, F. Blouin and R. Couture, A search for good multiple recursive random number generators, ACM Transaction on Modeling and Computer Simulation 3, 87-98, (1993). 9. P. L'Ecuyer and R. Couture, An implementation of the lattice and spectral tests for multiple

recursive linear random number generators, INFORMS Journal on Computing 9, 206-217, (1997).

10. H. Niederreiter, Random number generation and quasi-Monte Carlo methods, SIAM CBMS-NSF Regional Conference Series in Applied Mathematics, vol. 63, SIAM, Philadelphia, (1992).

11. H.C. Tang and C. Kao, Searching for good multiple recursive random number generators via a genetic algorithm, INFORMS Journal on Computing 16, 284-290, (2004).

12. H.C. Tang, An exhaustive analysis of two-term multiple recursive random number generators with efficient multipliers, Journal of Computational and Applied Mathematics

192, 411-416, (2006).

13. H. Riesel, Prime numbers and computer methods for factorization, 2nd edition, Birkhäuser

Boston, (1994).

14. F. Sezgin, A method of systematic search for optimal multipliers in congruential random number generators, BIT Numerical Mathematics 44, 135-149, (2004).

15. Wolfram Research, Inc., Mathematica 4.1, Champaign, IL, (1999).

計畫成果自評

本研究為搜尋模數、階數及 為質數，使得有效率檢查多階質數乘餘

法亂數產生器之全週期。結果可得最大週期為模數 ，4217階質數乘餘法亂數產

生器以符合現今模擬所需之長週期 ，此一部部份已寫成”Generalized Mersenne

Prime Number and Full Period Multiple Recursive Random Number Generators”，並投稿於國際 性期刊。另一方面，由實驗可知模數為Sophie Germain質數是一較佳選擇，所以針對模數為 Sophie Germain質數，搜尋較佳光譜值之全週期線性同餘法及多階質數乘餘法亂數產生器。 ) 1 /( ) 1 ( − − = m m r k 69 231₋ = m 7482 . 39352 10