64bit数の素数判定 今回の記事はZennでも公開しています https://zenn.dev/mizar/articles/791698ea860581 64bit数の素数判定 20220927 Mizar/みざー 1

自己紹介 Mizar/みざー : アカデミック関係者ではありません VRChatは 2020年5月より 2022年春頃まではコンピュータ将棋の開発に参加 最近は AtCoder にて 競技プログラミング始めました 最近の参加成績は現状の一つ上のレーティング水 色（その上はレーティング400毎に青・黄色・橙・ 赤）くらいで伸び悩みそうかなという位 AtCoder では Rust言語 をメイン言語として参加 64bit数の素数判定 20220927 Mizar/みざー 2

高等学校学習指導要領（平成30年告示）解説 数学編 理数編 平成30年7月 文部 科学省 に一部追記 数学Ａ (２単位) 1. 図形の性質 平面図形 (三角形の性質・円の性質・作図) 空間図形 2. 場合の数と確率 場合の数 (数え上げの原則・順列・組合せ) 確率 (確率とその基本的な法則(＊余事象，排反，期待値)・独立な試行 と確率・条件付き確率) 3. 数学と人間の活動 (※改定前: 整数の性質) 数量や図形と人間の活動 (＊整数の約数や倍数，ユークリッドの互除法 （最大公約数を求める・一次不定方程式の整数解），二進法，平面や 空間における点の位置) 遊びの中の数学 64bit数の素数判定 20220927 三つの内容のすべてを履修するときには 3 単位程度の単位数を必要とす るが，標準単位数は 2 単位であり，生徒の特性や学校の実態，単位数等 に応じて内容を適宜選択して履修させることとしている。 “ “ Mizar/みざー 3

今回の内容 2022年8月に競技プログラミング関係のサイト (yukicoder) で目に留まった素数判定問 題を掘り下げてみたのが発端 : 同サイトにて私の提出したコードも見ることが出来ます 2022年8月31日にはこのスライドの主な内容は Zenn にて公開済み 「64bit数」の「とにかく速い」「素数判定」処理の実装が今回の第一目標 「多倍長整数」の素数判定や、「素因数分解」の手法についても触れることはある が、今回そちらはメインではない それでも詰めていくと、様々な沼が広がっているという事を紹介 暗号処理目的の実装にはあまり向かない所も多々ある 処理タイミングや消費電力量の変化など、装置の物理的な挙動を観察する事で情報 の抜き取りを図る、いわゆるサイドチャネル攻撃に対する耐性が低い手法も含む 暗号実装への攻撃に対する攻撃手法はサイドチャネル攻撃には限らないが、これも 今回の主眼ではない 64bit数の素数判定 20220927 Mizar/みざー 4

素数判定問題 No.3030 ミラー・ラビン素数判定法のテスト https://yukicoder.me/problems/no/3030 実行時間制限: 1ケース 9.973秒 / メモリ制限: 509MB / ソースコードのサイズ制限: 64kiB 問題文 与えられた 個の正整数 それぞれについて、それが素数かどうかを判 定してください。 は正整数。 (64ビット符号つき整数に収まります) 出力 正整数のそれぞれについて、その値と素数か否か（素数なら 1 そうでないなら 0 ）を半角 空白区切りで、入力順を維持したまま出力し、改行してください。 64bit数の素数判定 20220927 Mizar/みざー 5

試し割りで解けるか？ (922京3372兆0368億5477万5808) (1844京6744兆0737億9551万1616) は素数ではない、 は素数、 以上の偶数は素数ではない。それ以外のケース( の奇数)な ら、 の奇数 で全て割り切れなければ素数と判定できるが… … 1万個の整数にそれぞれ約15億回試し割りする必要がある計算 仮に1秒で10億回試し割りできるとしても(除算は加減乗算より計算が重いので、実際にはより悪 く見積もった方が良い)、10秒以内には解けない 判定する整数1万個 × 除算15億回 / (楽観で)秒速10億回 = 15000秒? の全ての素数( : 約1.5億個)を表にするのが困難(実行時間10秒制 限・ソースコード64kiB制限・メモリ509MB制限)、その素数表を1から生成する時間もない、 の素数表のコストを無視しても試し割り法ではやっぱり時間が足りない 素数表生成 + 判定する整数1万個 × 除算1.5億回 / (楽観で)秒速10億回 = 素数表生成 + 1500秒? 64bit数の素数判定 20220927 Mizar/みざー 6

このスライドにおける数学記号の補足説明(1) : モジュロ(modulo)・剰余演算、割り算における余りの値を取得する演算記号。 ja.wikipedia.org - 剰余演算 ja.wikipedia.org - 除法 ja.wikipedia.org - 整数の合同 : 両辺の値が において合同、つまり両辺の値をそれぞれ で割った余り が等しい事を示す : 両辺の値が において合同ではない、つまり両辺の値をそれぞれ で割 った余りが等しくない事を示す : それぞれの演算子（加算記号・減算記号・乗算記号）の左 右にある値の入力およびその結果出力をそれぞれ、 で割った余りの値を使う事を 示す 64bit数の素数判定 20220927 Mizar/みざー 7

このスライドにおける数学記号の補足説明(2) : 最大公約数 (greatest common divisor)。 greatest common measure と言う事も。 ja.wikipedia.org - 最大公約数 : 最小公倍数 (least common multiple)。 ja.wikipedia.org - 最小公倍数 : 最小値 (minimum)。 : 最大値 (maximum)。 ja.wikipedia.org - 最大と最小 : 絶対値 (absolute value)。 ja.wikipedia.org - 絶対値 : 平方根 (square root)。 ja.wikipedia.org - 平方根 : べき乗、冪乗 (exponentiation)。 ja.wikipedia.org - 冪乗 64bit数の素数判定 20220927 Mizar/みざー 8

このスライドにおける数学記号の補足説明(3) : 対数 (logarithm)。 ja.wikipedia.org - 対数 : 自然対数 (natural logarithm)。 ja.wikipedia.org - 自然対数 : 床関数 (floor function)。実数に対してそれ以下の最大の整数を返す関数。正の数 の場合は小数点以下切り捨ての操作と等しい。 ja.wikipedia.org - 床関数と天井関数 : 近似 (approximation)。ほぼ等しい。 ja.wikipedia.org - 近似 : ランダウの記号 (Landau symbol)。このスライドでは計算量の問題を大まかに述 べるのに使っている。 ja.wikipedia.org - ランダウの記号 64bit数の素数判定 20220927 Mizar/みざー 9

このスライドにおける数学記号の補足説明(4) : ルジャンドル記号 (Legendre symbol)・ヤコビ記号 (Jacobi symbol)・クロネッカ ー記号 (Kronecker symbol)。 のように縦に並べた形式が取られる事も多いが、こ のスライドでは横に並べる形式で表記する。 ja.wikipedia.org - ルジャンドル記号 が 奇素数の場合に定義された記号。 ja.wikipedia.org - ヤコビ記号 ルジャンドル記号を拡張して が 正の奇数の場合に定義された記号。 en.wikipedia.org - Kronecker_symbol ヤコビ記号を拡張して が任意の整数の場合に定義された記号。線形代数などで 使われる「クロネッカーのデルタ」とは別物。 64bit数の素数判定 20220927 Mizar/みざー 10

(剰余) の基本的な性質(1) 同値律 反射律: 任意の整数 に関して ; 対称律: 任意の整数の対 に関して ; 推移律: 任意の整数の組 に関して かつ ならば . 例えば 64bit数の素数判定 20220927 Mizar/みざー 11

(剰余) の基本的な性質(2) 代数学的性質 かつ ならば ならば 任意の整数 に対して 任意の整数 に対して 任意の整数 に対して 正整数 に対して 64bit数の素数判定 20220927 Mizar/みざー 12

(剰余) の基本的な性質(3) 単位的環 整数 および 任意の において 加法の結合性: 加法の可換性: 加法単位元: ならば 加法逆元: となる が存在する 乗法の結合性: 乗法単位元: ならば 左右分配性: 64bit数の素数判定 20220927 Mizar/みざー 13

RSA暗号 桁数が大きい合成数の素因数分解が現実的な計算量では困難であることを安全性の前提とした公開鍵暗号 の一つ。1977年に発明された。ja.wikipedia.org - RSA暗号 SSL公開鍵(WebのHTTPS通信など)、SSH公開鍵などで現在も使われている所が多いが、RSA暗号からは楕 円曲線暗号(ECDSA, ED25519, ...)などのより強力な暗号方式への移行が徐々に進んでいる。例えばRSA 1024bit鍵は2011年以降使用禁止扱いであり、RSA 2048bit鍵は2031年以降使用禁止扱いとなる見込み。 nict.go.jp - 暗号の安全性評価 en.wikipedia.org - RSA_Factoring_Challenge: RSA素因数分解チャレンジ（抜粋） RSA_Number 10進 桁数 2進 桁数 素因数分解 された日 RSA_Number 10進 桁数 2進 桁数 素因数分解 された日 RSA-100 100 330 1991-04-01 RSA-260 260 862 未踏 RSA-640 193 640 2005-11-02 RSA-896 270 896 未踏 RSA-768 232 768 2009-12-12 RSA-1024 309 1024 未踏 RSA-240 240 795 2019-12-02 RSA-1536 463 1536 未踏 RSA-250 250 829 2020-02-28 RSA-2048 617 2048 未踏 64bit数の素数判定 20220927 Mizar/みざー 14

2020-02-28 に解かれた RSA-250 (10進250桁・2進829桁) の素因数分解: RSA-250 = 214032465024074496126442307283933356300861471514475501779775492088141802344714013664 334551909580467961099285187247091458768739626192155736304745477052080511905649310668 7691590019759405693457452230589325976697471681738069364894699871578494975937497937 RSA-250 = 641352894770715802787901901705773890848250147429434472081168596320245323446302386235 98752668347708737661925585694639798853367 × 333720275949781565562260106053551142279407603447675546667845209870238417292100370802 57448673296881877565718986258036932062711 おそらくまだ素因数分解されていない RSA-2048 (10進617桁・2進2048桁) の合成数: RSA-2048 = 25195908475657893494027183240048398571429282126204032027777137836043662020707595556264018 52588078440691829064124951508218929855914917618450280848912007284499268739280728777673597 14183472702618963750149718246911650776133798590957000973304597488084284017974291006424586 91817195118746121515172654632282216869987549182422433637259085141865462043576798423387184 77444792073993423658482382428119816381501067481045166037730605620161967625613384414360383 39044149526344321901146575444541784240209246165157233507787077498171257724679629263863563 73289912154831438167899885040445364023527381951378636564391212010397122822120720357 64bit数の素数判定 20220927 Mizar/みざー 15

素数は幾つあるのか 素数計数関数 は、実数 に対しそれ以下の素数の個数を 表す関数。近似式は や が知られてい るが、後者の方が収束が優れている。 以下の素数は 個 以下の素数は 個 以下の素数は 個 以下の素数は 個 以下の素数は 個 以下の素数は 個 以下の素数は 個 (2022-02-28) 以下の素数は 個 OEIS6880: Number of primes < 10^n. OEIS7053: Number of primes <= 2^n. 64bit数の素数判定 20220927 Mizar/みざー 16

64bit数の素数判定 20220927 Mizar/みざー 17

64bit数の素数判定 20220927 Mizar/みざー 18

64bit数の素数判定 20220927 Mizar/みざー 19

64bit数の素数判定 20220927 Mizar/みざー 20

素数の存在する割合 64bit数の素数判定 20220927 素数計数関数の近似式 からも見て取 れる通り、桁数が2倍になる毎に素数の存在する割合 はおよそ半分になる。 右の表では、その割合の逆数 の値を示す。 における の値は を級数展開したもの ( Maxima expintegral_li (z) function ) を用いて、 における の値は を用いて計算している。 Mizar/みざー 21

64bit数の素数判定 20220927 Mizar/みざー 22

64bit数の素数判定 20220927 Mizar/みざー 23

Ubuntu 22.04 以降には primecount という素数計数関数の計算をするパッケージが収録さ れています。(他にもWindows版などもあります) https://github.com/kimwalisch/primecount/ $ sudo apt install primecount $ primecount 10000 1229 $ primecount 2**32 203280221 $ primecount 1e10 455052511 $ primecount 2**64 425656284035217743 $ primecount 1e20 2220819602560918840 64bit数の素数判定 20220927 Mizar/みざー 24

さまざまな素数判定法 決定的素数判定法 試し割り法・エラトステネスの篩 ( ) 遅い Miller素数判定法 APR(Adleman–Pomerance–Rumely)素数判定法 ( ) 高速だが多項式時間ではない 楕円曲線素数判定法 ( ) 2022年6月、 の素数判定の記録 AKS(Agrawal–Kayal–Saxena)素数判定法 多項式時間で決定的なアルゴリズムとして2002年発表 確率的素数判定法 (今回はこちら、 未満の範囲で決定的に判定する) Miller-Rabin 素数判定法 Miller素数判定法を確率的判定法にしたもの Lucas 確率的素数判定法 Baillie–PSW 素数判定法 (Miller-Rabin + Lucas の複合) 特殊な形の数に固有の素数判定法 など メルセンヌ数( の形の自然数、 は自然数 )に対する素数判定法、Lucas–Lehmer test プロス数( の形の自然数、 は自然数、 は奇数かつ )に対する素数判定法 の形の自然数 に対する素数判定法、 は自然数、 は奇数かつ Lucas–Lehmer–Riesel test 64bit数の素数判定 20220927 Mizar/みざー 25

フェルマーの小定理 を素数とし、 を整数とすると、 が成立する。 を素数とし、 を と互いに素な整数とすると が成立する。 オイラーの定理 を と互いに素な整数とすると、 が成立する。 は 未満の と互いに素な自然数の個数を表す。(オイラー関数) が素数の時 となり、フェルマーの小定理に一致する。 カーマイケルの定理 を と互いに素な整数とすると、 が成立する。 なら、 が奇素数 を用いて と書けるなら、 が と素因数分解できるなら、 はカーマイケル関数。 は最小公倍数。 64bit数の素数判定 20220927 Mizar/みざー 26

高校数学の美しい物語 : フェルマーの小定理の証明と例題 https://manabitimes.jp/math/680 (補題) と が互いに素なとき を で割った余りは全て異なる。 (補題の証明) と を で割った余りが同じだと仮定すると、 は の倍数となるはずだが、 かつ と は互いに素なのでこれは矛盾。 と が互いに素なとき、 は全て 素数 の倍数でないので、補 題より を で割った余りを並べると から までが全て 1回ずつ登場する。 よって、 ここで、 は と互いに素なので両辺を で割って、 を得る。 64bit数の素数判定 20220927 Mizar/みざー 27

英語版Wikipedia : フェルマーの小定 理の証明たち https://en.wikipedia.org/wiki/Proof s_of_Fermat's_little_theorem 数学的背景をほとんど必要としな い、組み合わせの観点からの証明。 通りの文字の組み合わせを、 種 類のシンボルからなる 長さ の輪状 のネックレス に当てはめて考えた 時、 同一シンボル 個 のみで構成さ れる 個の文字列を取り除くと、残 りは全て 個 ずつのグループを作れ る と説明した もの。 つまり、 となる。 64bit数の素数判定 20220927 Mizar/みざー 28

Fermat 素数判定法・Carmichael 数 フェルマー Fermat の小定理の対偶「 と互いに素な整数 が を満たせば は合成数である」 ただしこの判定法のみでは全ての合成数を判定するには不十分で、 と互いに素な整数 がすべて を満たしてしまう カーマイケル Carmichael 数(OEIS2997) も存在する。 561, 1105, 1729, 2465, 2821, 6601, 8911, 10585, 15841, 29341, 41041, 46657, 52633, 62745, ... カーマイケル数 は3つ以上の素数の積で構成されている奇数で、その全ての素因子 に対して が を割り切るという特徴を持っている。10333229505 個の素因数を持つ 295486761787桁 のカーマイケル数も知られている(W.R. Alford; et al. 2012 arXiv, 2014 doi)。 , 64bit数の素数判定 20220927 Mizar/みざー 29

Miller 素数判定法 (決定的素数判定) の奇数に対して以下のように判定する。 1. を のべき乗で割り、 の形式 (ただし は奇数) にする。 2. 以下を の範囲の全ての base 値 について繰り返す: かつ の範囲の全ての について なら、 composite(合成数) を返す。 3. 2.で composite(合成数) と判定されなかった場合、 prime(素数) を返す。 拡張リーマン予想が真であることを前提とした決定的素数判定法のため、根拠がやや弱 い。証明されていない予想に依存しない、他の決定的素数判定法もある の値は だと 、 だと 。でも、 の時に 回弱もこの通りに調べる必要があるかと言うと…。 64bit数の素数判定 20220927 Mizar/みざー 30

Miller-Rabin 素数判定法(1) (確率的素数判定) の奇数に対して以下のように判定する。 は判定の正確度のパラメータ。 1. を のべき乗で割り、 の形式 (ただし は奇数) にする。 2. 以下を 回繰り返す: の範囲から base 値 をランダムに選ぶ。 かつ の範囲の全ての について なら composite(合成数) を返す。 3. 2.で composite(合成数) と判定されなかった場合、 probably prime(おそらく素数) を返す。 (n が合成数なのに誤って probably prime(おそらく素数) を返す確率は 未満) これは確率的な素数判定法だが、 の選び方を工夫すると…。 64bit数の素数判定 20220927 Mizar/みざー 31

Miller-Rabin 素数判定法(2) (範囲限定:決定論的変形) の範囲では7回のループで判定可能な事が判明している。(Jim Sinclair, 2011年) 1. を のべき乗で割り、 の形式 (ただし は奇数) にする。 2. base 値 の値をそれぞれ として以下 を行う: かつ かつ の範囲の全ての について なら composite(合成数) を返す。 3. 2.で composite(合成数) と判定されなかった場合、 probably prime(おそらく素数) を返す。 (少 なくとも の範囲で probably prime(おそらく素数) を返した場合、 は必ず素数) の範囲であれば、 は の3つで判定可能。 (Gerhard Jaeschke, 1993年) となる場合があるため、 の条件が追加されている点に注意。 64bit数の素数判定 20220927 Mizar/みざー 32

The best known SPRP bases sets https://miller-rabin.appspot.com/ count date best solution bases discoverer 1 2013- 03-03 Steve Worley 2 2013- 07-26 Wojciech Izykowski, Marcin Panasiuk 3 2013- 05-27 Steve Worley 4 2013- 05-23 Steve Worley 5 2013- 05-23 Steve Worley 6 2013- 05-23 Steve Worley 7 2011- 04-20 at least Jim Sinclair 64bit数の素数判定 20220927 Mizar/みざー 33

冪剰余 べきじょうよ の計算・左向きバイナリ法(1) 例えば を計算してみましょう…? 流石に384京回も掛け算は出来ない → まずは を2進数に分解 64bit数の素数判定 20220927 Mizar/みざー 34

冪剰余 べきじょうよ の計算・左向きバイナリ法(2) を別途求めれば冪剰余を効率的に計算できる 64bit数の素数判定 20220927 Mizar/みざー 35

Miller-Rabin 素数判定法 の実装 では、 Miller-Rabin 素数判定法 を実装してみましょう PyPy3、C++17(gcc)、Rust それぞれでの実装例 63bit数1万個の素数判定の実行時間の例は以下の通り: https://yukicoder.me/submissions/792430 : PyPy3: 394ms https://yukicoder.me/submissions/794259 : C++17(gcc): 211ms https://yukicoder.me/submissions/794257 : Rust: 198ms 64bit数の素数判定 20220927 Mizar/みざー 36

# Python : nは64bitで表す事ができる自然数とし、1以下の整数もしくは合成数ならFalse、素数ならTrueを返す def miller_rabin(n: int): if n == 2: return True if n < 2 or (n & 1) == 0: return False n1 = n - 1 d = n1 s = 0 while (d & 1) == 0: d //= 2 s += 1 for a in [2,325,9375,28178,450775,9780504,1795265022]: if a % n == 0: continue t = pow(a, d, n) # = pow(a, d) % n if t == 1 or t == n1: continue for _ in range(s - 1): t = pow(t, 2, n) # = pow(t, 2) % n if t == n1: break # else節を実行せず次の for a in ... ループへ else: # breakでループを抜けず、普通に for _ in range() ... ループが終了した時 # for文に付くelse節の説明: # https://docs.python.org/ja/3/reference/compound_stmts.html#the-for-statement return False return True 64bit数の素数判定 20220927 Mizar/みざー 37

// C++ : miller_rabin(n) : 1以下の整数もしくは合成数ならfalse、素数ならtrueを返す #include #include uint64_t modmul(uint64_t a, uint64_t b, uint64_t n) { return (uint64_t)(((__uint128_t)a) * ((__uint128_t)b) % ((__uint128_t)n)); } uint64_t modpow(uint64_t a, uint64_t b, uint64_t n) { uint64_t t = ((b & 1) == 0) ? 1 : a; for (b >>= 1; b != 0; b >>= 1) { a = modmul(a, a, n); if ((b & 1) == 1) { t = modmul(t, a, n); } } return t; } const uint64_t bases[] = {2,325,9375,28178,450775,9780504,1795265022}; bool miller_rabin(uint64_t n) { if (n == 2) { return true; } if (n < 2 || (n & 1) == 0) { return false; } uint64_t n1 = n - 1, d = n - 1; uint32_t s = 0; for (; (d & 1) == 0; d >>= 1) { s += 1; } for (const auto& base : bases) { uint64_t a = base; if (a >= n) { a %= n; if (a == 0) { continue; } } uint64_t t = modpow(a, d, n); if (t == 1) { continue; } for (uint32_t j = 1; t != n1; ++j) { if (j >= s) { return false; } t = modmul(t, t, n); } } return true; } 64bit数の素数判定 20220927 Mizar/みざー 38

// Rust : miller_rabin(n) : 1以下の整数もしくは合成数ならfalse、素数ならtrueを返す pub fn modmul(a: u64, b: u64, n: u64) -> u64 { ((a as u128) * (b as u128) % (n as u128)) as u64 } pub fn modpow(mut b: u64, mut p: u64, n: u64) -> u64 { let mut r = if (p & 1) == 0 { 1 } else { b }; loop { p >>= 1; if p == 0 { return r; } b = modmul(b, b, n); if p & 1 != 0 { r = modmul(r, b, n); } } } pub fn miller_rabin(n: u64) -> bool { if n == 2 { return true; } if n < 2 || n & 1 == 0 { return false; } let n1 = n - 1; let s = n1.trailing_zeros(); let d = n1 >> s; [2,325,9375,28178,450775,9780504,1795265022].iter().all(|&base| { let a = if base < n { base } else { base % n }; if a == 0 { return true; } let mut t = modpow(a, d, n); if t == 1 || t == n1 { return true; } for _ in 1..s { t = modmul(t, t, n); if t == n1 { return true; } } false }) } 64bit数の素数判定 20220927 Mizar/みざー 39

ご清聴ありがとうございました 前編終了です 64bit数の素数判定 20220927 Mizar/みざー 40

mod(剰余)演算 整数の合同 mod n での 加減乗算・合同 冪剰余の計算 ax mod n ミラー・ラビン 素数判定法 C++/Rust/Python プログラミング modでの乗法逆元 ×a mod n の結果が1になる値 ニュートン・ラフソン法 逆数・平⽅根の 反復的な近似法 最⼤公約数(gcd) 互いに素 ユークリッドの互除法 最⼤公約数(gcd)の求め⽅ 拡張ユークリッド互除法 ax+by=gcd(a,b) の x, y, gcd(a,b) の求め⽅ モンゴメリ剰余乗算 ヤコビ記号 リュカ数列 リュカ素数判定法 強いリュカ擬素数テスト BPSW素数判定法 モンゴメリ剰余乗算+ ミラー・ラビン 素数判定法 モンゴメリ剰余乗算+ BPSW素数判定法 前編の範囲 平⽅剰余 ニュートン・ラフソン法 による 整数平⽅根 の求め⽅ 拡張ユークリッド互除法 による 任意の mod n での 乗法逆元 の求め⽅ ニュートン・ラフソン法 による mod 264, mod 1016 など での 乗法逆元 の求め⽅ mod 8, mod 4095 での平⽅数判定 mod 264 での乗法逆元を 使った倍数判定法 最⼤公約数(gcd) 互いに素 mod(剰余)演算 整数の合同 mod n での 加減乗算・合同 冪剰余の計算 ax mod n C++/Rust/Python プログラミング 平⽅剰余 ヤコビ記号 試し割り法による 素数判定/素因数分解 素数計数関数と その近似式 後編予告 前編ではMiller-Rabin素数判定法を 実装するまでを話しました 後編では、これをさらに10倍速で実 行する取り組みを話します 10倍速の鍵は 「モンゴメリ剰余乗算」 「Baillie-PSW素数判定法」 次回もお楽しみに Zennでも記事を公開中です: 64bit数の素数判定 20220927 Mizar/みざー 41

宣伝 2022-09-24(土) 21:00～22:40 に、 AtCoder にて初心者向けの競技プログ ラミングコンテスト、 トヨタ自動車プログラミングコンテス ト2022(AtCoder Beginner Contest 270) が行われます。 VRChat競技プログラミング部では、 毎週 AtCoder Beginner Contest 終了後 の 23:00頃より VRChat 内で ABC感想 会 を実施しています。 興味のある方は是非 Discord や VRChat 感想会 に参加してみてください。 64bit数の素数判定 20220927 Mizar/みざー 42

もっとspeedupしてみよう 題意は既に満たしているが、 C++言語 や Rust言語 でさらに 10倍速 を目指す 128bit剰余算は遅いので多用を避けたい → モンゴメリ剰余乗算を使う 128bit剰余算1回の代わりに、おおむね 64bit乗算2回+64bit加減算2回+上位64bit切 り出し2回+下位64bit切り出し2回 を使う Pythonではモンゴメリ剰余乗算での高速化は難しいのでこれは行わない 例えばCPythonは整数型に 進数の内部表現を用いており、モンゴメリ剰余乗算 を使うには実行コストが大きい C++ では128bit整数だけではなく、GCCのビルトイン関数( __builtin_uaddll_overflow , __builtin_usubll_overflow , __builtin_clzll , __builtin_ctzll )も使う Rustでは整数型で overflowing_add , overflowing_sub , leading_zeros , trailing_zeros などが使える Miller-Rabin 素数判定法より速い方法は無いか → Baillie–PSW 素数判定法がある 64bit数の素数判定 20220927 Mizar/みざー 43

mod(剰余)演算 整数の合同 mod n での 加減乗算・合同 冪剰余の計算 ax mod n ミラー・ラビン 素数判定法 C++/Rust/Python プログラミング modでの乗法逆元 ×a mod n の結果が1になる値 ニュートン・ラフソン法 逆数・平⽅根の 反復的な近似法 最⼤公約数(gcd) 互いに素 ユークリッドの互除法 最⼤公約数(gcd)の求め⽅ 拡張ユークリッド互除法 ax+by=gcd(a,b) の x, y, gcd(a,b) の求め⽅ モンゴメリ剰余乗算 ヤコビ記号 リュカ数列 リュカ素数判定法 強いリュカ擬素数テスト BPSW素数判定法 モンゴメリ剰余乗算+ ミラー・ラビン 素数判定法 モンゴメリ剰余乗算+ BPSW素数判定法 前編の範囲 平⽅剰余 ニュートン・ラフソン法 による 整数平⽅根 の求め⽅ 拡張ユークリッド互除法 による 任意の mod n での 乗法逆元 の求め⽅ ニュートン・ラフソン法 による mod 264, mod 1016 など での 乗法逆元 の求め⽅ mod 8, mod 4095 での平⽅数判定 mod 264 での乗法逆元を 使った倍数判定法 最⼤公約数(gcd) 互いに素 mod(剰余)演算 整数の合同 mod n での 加減乗算・合同 冪剰余の計算 ax mod n C++/Rust/Python プログラミング 平⽅剰余 ヤコビ記号 試し割り法による 素数判定/素因数分解 素数計数関数と その近似式 次は: 剰余の乗法逆元（モジュラ逆数） ニュートン・ラフソン法 ニュートン・ラフソン法による 累乗数を法 とする乗法逆元 での乗法逆元を用いた倍 数判定 64bit数の素数判定 20220927 Mizar/みざー 44

剰余の乗法逆元（モジュラ逆数） のように、剰余において乗算すると になるような数の事を 乗法逆元 と呼びます。 において、乗法逆元が存在するのは と 互いに素な整数の時 です。 ここでは、このような乗法逆元を計算する手法を2つ紹介しておきます。 などの累乗数を法とする場合の乗法逆元の求め方 (ニュートン・ラフソン法) 素数などを法とする場合の乗法逆元の求め方 (拡張ユークリッドの互除法) 日本語版Wikipediaのモンゴメリ乗算の項 には が の冪の時の と なる の効率的な求め方も載っているがここでは省略。（実質的な乗算は不要だが、必要 なbit数と同じ回数の繰り返しが必要なため、今回のターゲットにおいては性能が劣る） 64bit数の素数判定 20220927 Mizar/みざー 45

Newton-Raphson法 アイザック・ニュートンとジョゼフ・ラフソンに由来。方程式を数値計算 によって反復的に解く求根アルゴリズム。例えば を実数における関数で ある として となる の値を求めたい時、 に近い適当な初期値 を取り の漸化式で、 に収束する数列を求める手法。右図は で を求める例 逆数 平方根 逆数平方根: 64bit数の素数判定 20220927 Mizar/みざー 46

Newton-Raphson法 の収束 初期値や関数の性質などの条件が適切であれば、 Newton-Raphson法 による漸化式は二次収束する（真の解に二乗の速度で近づく）。 平方根 64bit数の素数判定 20220927 Mizar/みざー 47

Newton-Raphson法 の高次への拡張 Newton-Raphson法 は一階微分を用いて二次収束する（1回の反復で精度の桁数を2倍にする） 漸化式を用いたが、これを拡張して更に高次への収束を用いる方法もある。ここでは のマクローリン展開を利用した方法を示す。 逆数: 64次収束: 次収束: 平方根: 7次収束: 逆数平方根: 7次収束: 他にも、高次収束漸化式の求め方には幾つか方法がある。Halley法、Ford-Pennline法、etc. 64bit数の素数判定 20220927 Mizar/みざー 48

高次収束漸化式の実行例 平方根: 7次収束: 64bit数の素数判定 20220927 Mizar/みざー 49

における奇数 の乗法逆元(1) 与えられた奇数 に対し、整数 を掛けて で割った余りが になる、整数 の 値を求める問題 ニュートン・ラフソン法と同じ漸化式を用いる。 64bit数の素数判定 20220927 // C++ #include #include uint64_t invmod_u64(uint64_t n) { assert(n % 2 == 1); uint64_t ni = n; for(int i = 0; i < 5; ++i) { ni = ni * (2 - n * ni); } assert(n * ni == 1); return ni; } // Rust fn invmod_u64(n: u64) -> u64 { assert_eq!(n % 2, 1); let mut ni = n; for _ in 0..5 { ni = ni.wrapping_mul( 2u64.wrapping_sub(n.wrapping_mul(ni)) ); } assert_eq!(n.wrapping_mul(ni), 1); ni } Mizar/みざー 50

における奇数 の乗法逆元(2) 64bit数の素数判定 20220927 アルゴリズム は正の奇数 証明 が、 を満たす時、 を と定義するなら、 を満たす。 とおくと、 奇数 を とおくと、 は偶数であるから とした時 Mizar/みざー 51

における奇数 の乗法逆元(3) 64bit数の素数判定 20220927 アルゴリズム は正の奇数 証明 が、 を満たす時、 を と定義するなら、 を満たす。 とおくと、 奇数 を とおくと、 は偶数であるから とした時 Mizar/みざー 52

における整数 の乗法逆元 64bit数の素数判定 20220927 アルゴリズム の場合で触れた証明を踏まえると、 進数の 場合でも同様に計算できる。有効桁数が1ループ毎に2 倍,3倍になる事の応用として、初期のループでは計算 する桁数を抑え、徐々に桁数を増やして計算する例と しても示す。 例題: , と なる を求める と は互いに素 なので における の乗法逆元 は存在する 計算する桁数: 初期値:10進1桁 → 1ループ目:10進3桁 → 2ループ目:10進9桁 Mizar/みざー 53

における整数 の乗法逆元(計算例) 64bit数の素数判定 20220927 アルゴリズム 例 , となる を求める 検算: Mizar/みざー 54

, となる を求める (続き) 検算: 64bit数の素数判定 20220927 Mizar/みざー 55

での乗法逆元を使った倍数判定(1) における の乗法逆元を として、整数 に関して が の倍数なら、 であ る。 が の倍数ではないなら、 である。 この事から の場合、奇数 の における乗法逆元 と、 を事前に準備しておけば、剰余算の代わりに乗算と比較によってこのように実装できる。 なら は の倍数 なら は の倍数ではない 64bit数の素数判定 20220927 Mizar/みざー 56

での乗法逆元を使った倍数判定(2) // C++ #include #include bool is_multiple_of_3_divmod(uint64_t n) { return n % 3 == 0; } bool is_multiple_of_3_mulinvmod(uint64_t n) { // C++: n が 3 の倍数なら true、上の関数と同等 return (uint64_t)(n * 0xAAAAAAAAAAAAAAABULL) <= (uint64_t)0x5555555555555555ULL; } // Rust fn is_multiple_of_3_divmod(n: u64) -> bool { n % 3 == 0 } fn is_multiple_of_3_mulinvmod(n: u64) -> bool { // Rust: n が 3 の倍数なら true、上の関数と同等 // 0xAAAA_AAAA_AAAA_AAAB : mod 2**64 における 3 の乗法逆元 // 0x5555_5555_5555_5555 : floor((2**64 - 1) / 3) n.wrapping_mul(0xAAAA_AAAA_AAAA_AAABu64) <= 0x5555_5555_5555_5555u64 } 64bit数の素数判定 20220927 Mizar/みざー 57

での乗法逆元 を使った倍数判定(3) いずれも同じアセンブリ出力になり、 コンパイラも実は同様の最適化をして いる事が分かる。 C++: https://godbolt.org/z/M8z1W1hqM Rust: https://rust.godbolt.org/z/cvzKrvYTo この例では除数が定数なのでコンパイ ル時の最適化が効いているが、変数の 場合は必ずしもこのような最適化が行 われるとは限らない。 64bit数の素数判定 20220927 Mizar/みざー 58

での乗法逆元を使った倍数判定(4) 通常の除算→乗法逆元を使った試し割り法での、64bit数の素因数分解の実行時間例: 666,196μs → 1,838,621μs: 事前準備 以下の素数( 個)読み込み(256MiB)+乗法逆元計算+ 計算 4μs → 3μs: 3: 3 6μs → 4μs: 11111111111111: 11 239 4649 909091 5μs → 3μs: 111111111111111: 3 31 37 41 271 2906161 3μs → 2μs: 1111111111111111: 11 17 73 101 137 5882353 1,157μs → 314μs: 11111111111111111: 2071723 5363222357 33μs → 8μs: 111111111111111111: 3 3 7 11 13 19 37 52579 333667 380,952μs → 79,259μs: 1111111111111111111: 1111111111111111111 508,707μs → 105,694μs: 6111111111111111111: 3 2037037037037037037 1,153,151μs → 220,903μs: 11297479344276717331: 3263519417 3461747243 1,491,976μs → 281,662μs: 18446744073709551557: 18446744073709551557 64bit数の素数判定 20220927 Mizar/みざー 59

for &n in &[ 3u64,11111111111111,111111111111111,1111111111111111,11111111111111111,111111111111111111, 1111111111111111111,6111111111111111111,11297479344276717331,18446744073709551557 ] { let mut m = n; let mut pv = vec![]; if m == 0 { continue; } while (m & 1) == 0 { pv.push(2); m >>= 1; } // 2の素因数は別途処理 let mut sqrtm = sqrt_newton(m); // mの平方根 をニュートン法で計算 for i in 0..(primes.len()) { // 32bit素数リスト から let p = primes[i]; // 素数p let (pinv, pdiv) = primes_id[i]; // 法2**64での乗法逆元, floor((2**64-1)/p) if (p as u64) > sqrtm { break; } // mの平方根 より 素数p が大きければ終了 let mut f = false; while m.wrapping_mul(pinv) <= pdiv { // m が 素数pの倍数なら (判定は乗法逆元,floor((2**64-1)/p)の値で) m /= p as u64; // 素数p で除算 pv.push(p as u64); // 素因数リストに追加 f = true; // 素因数が見つかったフラグ } if f { sqrtm = sqrt_newton(m); } // 素因数を見つけたら mの平方根 を再計算 } if m > 1 { pv.push(m); } // 余った素因数をリストに追加 println!("{}: {}", n, pv.iter().cloned().map(|p| p.to_string()).collect::>().join(" ")); } 64bit数の素数判定 20220927 Mizar/みざー 60

mod(剰余)演算 整数の合同 mod n での 加減乗算・合同 冪剰余の計算 ax mod n ミラー・ラビン 素数判定法 C++/Rust/Python プログラミング modでの乗法逆元 ×a mod n の結果が1になる値 ニュートン・ラフソン法 逆数・平⽅根の 反復的な近似法 最⼤公約数(gcd) 互いに素 ユークリッドの互除法 最⼤公約数(gcd)の求め⽅ 拡張ユークリッド互除法 ax+by=gcd(a,b) の x, y, gcd(a,b) の求め⽅ モンゴメリ剰余乗算 ヤコビ記号 リュカ数列 リュカ素数判定法 強いリュカ擬素数テスト BPSW素数判定法 モンゴメリ剰余乗算+ ミラー・ラビン 素数判定法 モンゴメリ剰余乗算+ BPSW素数判定法 前編の範囲 平⽅剰余 ニュートン・ラフソン法 による 整数平⽅根 の求め⽅ 拡張ユークリッド互除法 による 任意の mod n での 乗法逆元 の求め⽅ ニュートン・ラフソン法 による mod 264, mod 1016 など での 乗法逆元 の求め⽅ mod 8, mod 4095 での平⽅数判定 mod 264 での乗法逆元を 使った倍数判定法 最⼤公約数(gcd) 互いに素 mod(剰余)演算 整数の合同 mod n での 加減乗算・合同 冪剰余の計算 ax mod n C++/Rust/Python プログラミング 平⽅剰余 ヤコビ記号 試し割り法による 素数判定/素因数分解 素数計数関数と その近似式 ここまで: 剰余の乗法逆元（モジュラ逆数） ニュートン・ラフソン法 ニュートン・ラフソン法による 累乗数を法 とする乗法逆元 での乗法逆元を用いた倍 数判定 次は: ユークリッドの互除法 拡張ユークリッド互除法 拡張ユークリッド互除法による 任意の での乗法逆元 64bit数の素数判定 20220927 Mizar/みざー 61

0 の約数? 「整数 が の約数である」とは、「ある整数 をとると が成立する」という定 義がありますが、ここに「 」の条件が付く説と外れる説があるようです。 の約数 (a) : すべての整数 (上の定義で「 」の条件が外れる場合) の約数 (b) : すべての でない 整数 (上の定義で「 」の条件が付く場合) の約数 : ( は含まれない) の約数 : ( は含まれない) の約数 : ( は含まれない) の約数 : ( は含まれない) の約数に が含まれない、としてしまうと を含む数に対する最大公約数を考える時に定 義が大変になるので、ここでは次から、 の約数に は含まれるもの として説明します。 64bit数の素数判定 20220927 Mizar/みざー 62

0 同士のときに最大値とは限らない最大公約数 と の公約数 は無限に存在しますが、最大公約数 ( ) は と定義されるようです。 数が最大な公約数ではなく、それらが持つ約数の集合が最大になる公約数を指す事に近い と考えることができます。 空集合 (2つの数における最大公約数の定義の例) 2つの整数 について、 以上の整数 が次 の2つの条件を満たす時、 は と の最大公約数という。 1. は と の 公約数 である。 2. と の全ての公約数が の約数である。 64bit数の素数判定 20220927 Mizar/みざー 63

ユークリッドの互除法 を負でない整数として、最大公約数 ( ) の値は (~ is odd = ～は奇数) 順番を入れ替えても最大公約数は変わらない 片方の値をもう片方から足したり引いたり、剰余を取ったりしても変わらない 共通しない素因数を取り除いても、最大公約数は変わらない のような性質を繰り返し使い、係数の値を減らして最大公約数の値を求めるアルゴリズムがユークリ ッドの互除法と呼ばれている。 64bit数の素数判定 20220927 Mizar/みざー 64

// Rust: 剰余算を利用したユークリッドの互除法の単純な実装 pub fn gcd_eucrid(mut x: u64, mut y: u64) -> u64 { while y != 0 { // y が 0 になるまで繰り返し // x を y で割った剰余に置き換える、結果、 x < y になる x %= y; // x と y の値を入れ替え、 x > y にする std::mem::swap(&mut x, &mut y); } // y の値が 0 になった時の x の値が、 gcd(x,y) の結果なのでそれを返す x } 64bit数の素数判定 20220927 Mizar/みざー 65

// Rust: (参考) Binary GCD の実装例 (拡張ユークリッド互除法ではない) pub fn gcd_binary(x: u64, y: u64) -> u64 { // 内部処理用、再帰関数だがスタックを消費しない形で(ループ処理と見なして)上手くコンパイル可能 fn gcd_binary_intl(x: u64, y: u64) -> u64 { // 常にここに入力される x,y はそれぞれ奇数として扱う debug_assert_eq!(x & y & 1, 1); if x == y { return x; } // x と y が一致した時点で終了 let (s, t) = if x > y { (x - y, y) } else { (y - x, x) }; // 減算で生じた偶数の 2 の素因数は公約数たり得ないのでシフト演算で潰し、桁数を減らす gcd_binary_intl(s >> s.trailing_zeros(), t) } if x == 0 { return y; } // x, y が 0 だった場合は上手く行かないので別途処理 if y == 0 { return x; } // x != 0 なので TLZ(x) は x が持つ 2 の素因数の数 let u = x.trailing_zeros(); // y != 0 なので TLZ(y) は y が持つ 2 の素因数の数 let v = y.trailing_zeros(); // 2の素因数はいったん全てシフト演算で潰して、最後の出力時に公約数分だけ戻す gcd_binary_intl(x >> u, y >> v) << u.min(v) } 64bit数の素数判定 20220927 Mizar/みざー 66

拡張ユークリッドの互除法(1) でない整数 に対し、 (ベズーの等式 Bézout's identity ) を満たす 整数 : が存在する。以下は拡張ユークリッドの互除法で の値を同時に求める例: 64bit数の素数判定 20220927 Mizar/みざー 67

拡張ユークリッドの互除法(2) の式に対し、 に乗法逆元を求めたい値、 に 法 の値を設定する と、 と が互いに素 であれば、拡張ユークリッドの互除法で の乗法逆 元 を求める事ができる。 64bit数の素数判定 20220927 Mizar/みざー 68

拡張ユークリッドの互除法(3) // Rust 拡張ユークリッドの互除法の単純な実装(再帰) pub fn extgcd(a: i64, b: i64) -> (i64, i64, i64) { /* return */ match b { 0 => /* if b == 0 then */ (a, 1, 0), _ => { // if b != 0 then let (d, y, x) = extgcd(b, a % b); (d, x, (y - ((a / b) * x))) }, }} (拡張)ユークリッドの互除法の実装には、例えば他にも以下のようなバリエーションがある。 Binary法 / 2進GCD (剰余算の代わりに減算と偶数検知とシフト演算を使う、短桁整数向け) Lehmer法 (上位桁に絞って除算を繰り返す事で全桁除算の回数を削減、多倍長整数向け) 64bit数の素数判定 20220927 Mizar/みざー 69

mod(剰余)演算 整数の合同 mod n での 加減乗算・合同 冪剰余の計算 ax mod n ミラー・ラビン 素数判定法 C++/Rust/Python プログラミング modでの乗法逆元 ×a mod n の結果が1になる値 ニュートン・ラフソン法 逆数・平⽅根の 反復的な近似法 最⼤公約数(gcd) 互いに素 ユークリッドの互除法 最⼤公約数(gcd)の求め⽅ 拡張ユークリッド互除法 ax+by=gcd(a,b) の x, y, gcd(a,b) の求め⽅ モンゴメリ剰余乗算 ヤコビ記号 リュカ数列 リュカ素数判定法 強いリュカ擬素数テスト BPSW素数判定法 モンゴメリ剰余乗算+ ミラー・ラビン 素数判定法 モンゴメリ剰余乗算+ BPSW素数判定法 前編の範囲 平⽅剰余 ニュートン・ラフソン法 による 整数平⽅根 の求め⽅ 拡張ユークリッド互除法 による 任意の mod n での 乗法逆元 の求め⽅ ニュートン・ラフソン法 による mod 264, mod 1016 など での 乗法逆元 の求め⽅ mod 8, mod 4095 での平⽅数判定 mod 264 での乗法逆元を 使った倍数判定法 最⼤公約数(gcd) 互いに素 mod(剰余)演算 整数の合同 mod n での 加減乗算・合同 冪剰余の計算 ax mod n C++/Rust/Python プログラミング 平⽅剰余 ヤコビ記号 試し割り法による 素数判定/素因数分解 素数計数関数と その近似式 ここまで: ユークリッドの互除法 拡張ユークリッド互除法 拡張ユークリッド互除法による 任意の での乗法逆元 次は: での平方数判 定 ニュートン・ラフソン法による整 数平方根 64bit数の素数判定 20220927 Mizar/みざー 70

平方数の判定(1) Miller-Rabin素数判定法では関係ないが、この後紹介する素数判定法では、与えられた数が 平方数だと問題が生じるものがある。また、平方数は全て ( や も含め) 素数ではない。 そのため、平方数の判定方法についても考察する。 を取った時、平方数が取り得る値は しかない。逆に、 を取った値が である場合、その数は平方数ではない事が分かる。 特に、整数 に対して は偶数であるから、 より、 奇数の平方数に を取ると となる。 64bit数の素数判定 20220927 Mizar/みざー 71

平方数の判定(2) を取った場合は 32個中7個(約22%)が平方数かもしれない値、 を取った場合は 4095個中336個(約8.2%)が平方数かもしれない値となり、 それでも平方数かもしれない数について、最後にニュートン法で平方根を計算している。 // 擬平方数判定 mod32 7/32 0.2188 (true: 平方数かもしれない, false: 平方数ではない) fn issq_mod32(x: u64) -> bool { (0x02030213u32 >> ((x as u32) & 31)) & 1 == 1 } // 平方数判定 (true: 平方数である, false: 平方数ではない) fn issq(x: u64) -> bool { issq_mod32(x) && // 擬平方数判定 mod32 7/32 0.2188 issq_mod4095(x) && // 擬平方数判定 mod4095 336/4095 0.0821 issq_newton(x) // ニュートン法による整数平方根 } 64bit数の素数判定 20220927 Mizar/みざー 72

平方数の判定(3) // 擬平方数判定 mod4095 336/4095 0.0821 (true: 平方数かもしれない, false: 平方数ではない) fn issq_mod4095(x: u64) -> bool { const SQTABLE_MOD4095: [u64; 64] = [ 0x2001002010213,0x4200001008028001,0x20000010004,0x80200082010,0x1800008200044029, 0x120080000010,0x2200000080410400,0x8100041000200800,0x800004000020100,0x402000400082201, 0x9004000040,0x800002000880,0x18002000012000,0x801208,0x26100000804010, 0x80000080000002,0x108040040101045,0x20c00004000102,0x400000100c0010,0x1300000040208, 0x804000020010000,0x1008402002400080,0x201001000200040,0x4402000000806000,0x10402000000, 0x1040008001200801,0x4080000000020400,0x10083080000002,0x8220140000040000,0x800084020100000, 0x80010400010000,0x1200020108008060,0x180000000,0x400002400000018,0x4241000200, 0x100800000000,0x10201008400483,0xc008000208201000,0x800420000100,0x2010002000410, 0x28041000000,0x4010080000024,0x400480010010080,0x200040028000008,0x100810084020, 0x20c0401000080000,0x1000240000220000,0x4000020800,0x410000000480000,0x8004008000804201, 0x806020000104000,0x2080002000211000,0x1001008001000,0x20000010024000,0x480200002040000, 0x48200044008000,0x100000000010080,0x80090400042,0x41040200800200,0x4000020100110, 0x2000400082200010,0x1008200000000040,0x2004800002,0x2002010000080 ]; let p = (x % 4095) as usize; (SQTABLE_MOD4095[p >> 6] >> (p & 63)) & 1 == 1 } 64bit数の素数判定 20220927 Mizar/みざー 73

平方数の判定(4) // ニュートン法による整数平方根 fn sqrt_newton(x: u64) -> u64 { if x <= 1 { return x; } // x == 0 の時に上手く働かないので x <= 1 の時はショートカット let k = 32 - ((x - 1).leading_zeros() >> 1); let mut s = (1 as u64) << k; // s = 2**k let mut t = (s + (x >> k)) >> 1; // t = (s + x/s)/2 // whileループ回数=除算回数は u64:最大6回 u32:最大5回 を想定 // s > floor(sqrt(x)) -> 単調減少する : floor(sqrt(x)) <= t < s // s == floor(sqrt(x)) -> 減少しない : s == floor(sqrt(x)) <= t <= floor(sqrt(x)) + 1 while t < s { s = t; t = (s + (x / s)) >> 1; } s } // 平方数判定 (true: 平方数である, false: 平方数ではない) fn issq_newton(x: u64) -> bool { let sqrt = sqrt_newton(x); sqrt * sqrt == x } 64bit数の素数判定 20220927 Mizar/みざー 74

とする。 であるとする。 (補題) は自然数であるから、 (1) のとき、 である。(この条件の間は単調減少する保証) より であり、 を使って とすると、 (2) のとき、 である。( 到達時の挙動の保証) より、 は自然数であるから および より、 64bit数の素数判定 20220927 Mizar/みざー 75

mod(剰余)演算 整数の合同 mod n での 加減乗算・合同 冪剰余の計算 ax mod n ミラー・ラビン 素数判定法 C++/Rust/Python プログラミング modでの乗法逆元 ×a mod n の結果が1になる値 ニュートン・ラフソン法 逆数・平⽅根の 反復的な近似法 最⼤公約数(gcd) 互いに素 ユークリッドの互除法 最⼤公約数(gcd)の求め⽅ 拡張ユークリッド互除法 ax+by=gcd(a,b) の x, y, gcd(a,b) の求め⽅ モンゴメリ剰余乗算 ヤコビ記号 リュカ数列 リュカ素数判定法 強いリュカ擬素数テスト BPSW素数判定法 モンゴメリ剰余乗算+ ミラー・ラビン 素数判定法 モンゴメリ剰余乗算+ BPSW素数判定法 前編の範囲 平⽅剰余 ニュートン・ラフソン法 による 整数平⽅根 の求め⽅ 拡張ユークリッド互除法 による 任意の mod n での 乗法逆元 の求め⽅ ニュートン・ラフソン法 による mod 264, mod 1016 など での 乗法逆元 の求め⽅ mod 8, mod 4095 での平⽅数判定 mod 264 での乗法逆元を 使った倍数判定法 最⼤公約数(gcd) 互いに素 mod(剰余)演算 整数の合同 mod n での 加減乗算・合同 冪剰余の計算 ax mod n C++/Rust/Python プログラミング 平⽅剰余 ヤコビ記号 試し割り法による 素数判定/素因数分解 素数計数関数と その近似式 ここまで: での平方数判 定 ニュートン・ラフソン法による整 数平方根 次は: モンゴメリ剰余乗算 64bit数の素数判定 20220927 Mizar/みざー 76

モンゴメリ剰余乗算(1) 概要 互いに素(最大公約数が )な自然数 を考える。 ここでは、 正の奇数 として考える。 そして、加減乗算・比較演算をする際、それぞれの変数にあらかじめ という変換 (以下、モンゴメリ変換・モンゴメリ表現)をして、以下のように処理を置き換えていく。 の値での除算・剰余算が効率的に行える場合、この モンゴメリ剰余乗算 の左辺が実は効率的に計算可能(モンゴメリリダクション: , )というのが肝になっている。 64bit数の素数判定 20220927 Mizar/みざー 77

モンゴメリ剰余乗算(2) モンゴメリリダクション であり、 は 、 は 、 は を満たすとする。 モンゴメリリダクション は以下の通り。 は床関数。 証明: より、 は で割り切れる。 より、 である。 より、 であるため、 手続きが返す値は より小さい。 64bit数の素数判定 20220927 Mizar/みざー 78

モンゴメリ剰余乗算(3) リダクション処理の変種 今回紹介した上の方法以外にも、モンゴメリリダクション は 幾つかのバリエーションが存在する。 ※ 下位ワードの処理を上手くコスト削減できる訳でなかったり、そもそも が の累乗数ではなく、床関数付き除算を効率的に行える値ではなかったりなど、 状況に応じて以下のような実装を行う場合もある。 ※ 64bit数の素数判定 20220927 Mizar/みざー 79

モンゴメリ剰余乗算(4) モンゴメリ変換・剰余乗算 あらかじめ を用意しておき、 と の乗算剰余 を以下のようにして求める。 , ただし、 の値をまた直ぐに次の演算に用いる場合、 の処理は不要となる。 乗算剰余を1回だけ求め、すぐにモンゴメリ表現から元に戻す場合は以下のように計算できる。 べき剰余 を求めたい場合、まず を計算し、 べき乗の計 算は を繰り返し適用する。べき剰余の計算においても、通常のべき 乗の計算と同様、バイナリ法( を計算し、 の 進数表現を基にそれらから必要な ものを掛け合わせる)などの効率化手法が利用できる。 64bit数の素数判定 20220927 Mizar/みざー 80

モンゴメリ剰余乗算(5) 10進数でモンゴメリ剰余乗算 とする。 と は互いに素であり、 を満たし、 、 を満たす値は である 。 とし、 の手法で を計算してみる。 64bit数の素数判定 20220927 Mizar/みざー 81

モンゴメリ剰余乗算(6) 10進数での計算 よって、 64bit数の素数判定 20220927 Mizar/みざー 82

モンゴメリ剰余乗算(7) まとめ 互いに素(最大公約数が )な自然数 を考え、かつ 正の奇数 のよう に、 による除算・剰余算が容易な値を定める。 モンゴメリ表現 (モンゴメリ変換 をした値) を用いると、加減乗算と による除 算・剰余算 によって における剰余乗算を計算できる。 モンゴメリ剰余乗算は、同じ法 を繰り返し使う場合に有用。 の変更が多いと、 その分 や などの定数を求め直すコストが重くなる。 モンゴメリ剰余乗算を用いて演算効率を上げるためには、モンゴメリ表現への変換（モンゴメリ 変換）・逆変換（乗算に伴わないモンゴメリリダクション）は、特にループ処理の内側では極力 行わない。モンゴメリ表現のままで演算や比較ができるよう、必要な定数はループのなるべく外 側であらかじめ求めておく。 モンゴメリ表現された変数値とモンゴメリ表現されていない定数値の比較がループ内で必要 な時も、あらかじめモンゴメリ変換した定数値を準備しておくと、ループ内でモンゴメリ表 現された変数値をモンゴメリリダクションで逆変換する必要が無くなる 例: 64bit数の素数判定 20220927 Mizar/みざー 83

AtCoder Regular Contest 148 (2022-09-11) "F - 998244353 → 1000000007" https://atcoder.jp/contests/arc148/tasks/arc148_f 問題文 この問題は output-only です。 符号無し 64 bit 整数の加算・乗算・ を除数とする modulo 演算ができるプログラミング言語があります。この 言語を用いて における乗算を行うプログラムを作成してください。 厳密に説明すると、 以上 以下の整数 が与えられたときに を計算するプログラム を、以下の 仕様・形式 に従って作成してください。 プログラムの仕様 このプログラムでは、英大文字で表される A,B,…,Z の 個の 変数 を扱うことが出来る。 各変数は 以上 未満の整数値 (以下 符号無し 64 bit 整数 と表記) を保持することが出来る。 プログラムの実行開始時点で、 には整数 が、 には整数 が、それ以外の変数には が代入されている。 プログラムの実行終了時点で変数 に が保持されている必要がある。 プログラムの形式 プログラムの 行目にはプログラムの命令数を表す整数 が書かれる。 プログラムの 行目から 行目には 個の命令が書かれる。命令は上から下に順次実行される。 命令は次の 3 つのいずれかである。 add x y z : に を代入する。ここでは は変数、 は変数または符号無し 64 bit 整数である。 mul x y z : に を代入する。ここでは は変数、 は変数または符号無し 64 bit 整数である。 rem x y : に を代入する。ここでは は変数、 は変数または符号無し 64 bit 整数である。 64bit数の素数判定 20220927 Mizar/みざー 84

AtCoder ARC148F "998244353 → 1000000007" 問題文: https://atcoder.jp/contests/arc148/tasks/arc148_f 私が考えてみた版の18行解の解説記事: https://zenn.dev/mizar/articles/84f360c00a30af snuke氏の解説記事: さらに短い17行解やcodegolf解への言及など https://atcoder.jp/contests/arc148/editorial/4831 64bit数の素数判定 20220927 Mizar/みざー 85

mod(剰余)演算 整数の合同 mod n での 加減乗算・合同 冪剰余の計算 ax mod n ミラー・ラビン 素数判定法 C++/Rust/Python プログラミング modでの乗法逆元 ×a mod n の結果が1になる値 ニュートン・ラフソン法 逆数・平⽅根の 反復的な近似法 最⼤公約数(gcd) 互いに素 ユークリッドの互除法 最⼤公約数(gcd)の求め⽅ 拡張ユークリッド互除法 ax+by=gcd(a,b) の x, y, gcd(a,b) の求め⽅ モンゴメリ剰余乗算 ヤコビ記号 リュカ数列 リュカ素数判定法 強いリュカ擬素数テスト BPSW素数判定法 モンゴメリ剰余乗算+ ミラー・ラビン 素数判定法 モンゴメリ剰余乗算+ BPSW素数判定法 前編の範囲 平⽅剰余 ニュートン・ラフソン法 による 整数平⽅根 の求め⽅ 拡張ユークリッド互除法 による 任意の mod n での 乗法逆元 の求め⽅ ニュートン・ラフソン法 による mod 264, mod 1016 など での 乗法逆元 の求め⽅ mod 8, mod 4095 での平⽅数判定 mod 264 での乗法逆元を 使った倍数判定法 最⼤公約数(gcd) 互いに素 mod(剰余)演算 整数の合同 mod n での 加減乗算・合同 冪剰余の計算 ax mod n C++/Rust/Python プログラミング 平⽅剰余 ヤコビ記号 試し割り法による 素数判定/素因数分解 素数計数関数と その近似式 ここまで: モンゴメリ剰余乗算 次は: Baillie-PSW素数判定法 平方剰余・ヤコビ記号 リュカ擬素数テスト リュカ数列 64bit数の素数判定 20220927 Mizar/みざー 86

BPSW素数判定法(1) 概要 Baillie-PSW素数判定法は の奇数に対して以下のように判定する。 1. base 値 を として Miller-Rabin素数判定を行う。失敗すれば composite(合成数) と出力して終了する。 1.のみでは、以下の合成数は擬素数となる。(OEIS:A001262 Strong pseudoprimes to base 2.) 2. の中で最初に となる を見つける。 はヤコビ記号。 が平方数の 場合 は永遠に見つからないため、適当なステップ数の経過時に平方数か確認し、平方数なら composite(合成数) と出力して終了する。 3. とし、 である事を確認する。 のリュカ数列を用い て 強いリュカ擬素数(strong Lucas pseudoprime)テスト を行う。失敗すれば composite(合成数) と出力して終 了する。 2.3.のみでは、以下の合成数は擬素数となる。(OEIS:A217255 Strong Lucas pseudoprimes.) 4. 1.2.3.の何れでも終了しなかった場合、 probably prime(おそらく素数) と出力して終了する。これは少なくとも の範囲では正しく素数と合成数を判別可能である。この判定法をすり抜けて擬素数となる合成数は おそらく無数に存在すると推測されるが、今のところその具体例は発見されていない。 64bit数の素数判定 20220927 Mizar/みざー 87

BPSW素数判定法(2) クロネッカー記号 Kronecker symbol ・ヤコビ記号 Jacobi symbol 任意の整数 と に対し、クロネッカー記号 は の奇素因数 に対応する ルジャンドル記号 Legendre symbol の積と、追加因数 によって定義される。ヤコビ記号はそのうち、 が正の奇数 の範囲における定義である。 と が素因数分解 されるとして ルジャンドル記号 は 任意の整数 と全ての奇素数 に対して次のように定義される。 平方剰余 平方非剰余 クロネッカー記号では、追加で以下の定義を行う。ヤコビ記号ではこのうち のみ使われる。 64bit数の素数判定 20220927 Mizar/みざー 88

BPSW素数判定法(3) ヤコビ記号の性質 ヤコビ記号 ( が正の奇数) の場合の性質を示す。必ずしもクロネッカー記号 ( が整数) には当てはまらない。 の時、 と が奇数で互いに素な正整数である場合、 の場合、 は を法として平方非剰余である。 が を法として平方剰余であり、 の場合、 である。 の場合、 は を法として平方剰余であるかもしれないし、ないかもしれない。奇素数 において の場合、 は を法として平方剰余である。 64bit数の素数判定 20220927 Mizar/みざー 89

BPSW素数判定法(4) ヤコビ記号の計算 ユークリッドの互除法に似た方法を使い、 程度の処理回数で計算できる。 // Rust : ヤコビ記号 (a/n) pub fn jacobi(a: i64, mut n: u64) -> i32 { debug_assert!((n & 1) == 1 && a > i64::MIN); // n は正の奇数に限る // step 1 : n を法として剰余を取り、 a を減らす let (mut a, mut t) = ((a.abs() as u64) % n, if a < 0 && (n & 3) == 3 { -1 } else { 1 }); // step 3 : a が 1 のとき、結果は 1。 a と n が互いに素でなければ、結果は 0。 while a != 0 { // step 2 : 偶数の a を抽出する while (a & 1) == 0 { a >>= 1; if (n & 7) == 3 || (n & 7) == 5 { t = -t; } } // step 4 : 記号を反転し、 n を法として a を減らす if (a & n & 3) == 3 { t = -t; } std::mem::swap(&mut a, &mut n); a %= n; } if n == 1 { t } else { 0 } } 64bit数の素数判定 20220927 Mizar/みざー 90

BPSW素数判定法(5) リュカ擬素数テスト 自然数 と 整数 に対して とおき、 をそれぞれ対応 するリュカ数列とする。 を である正の奇数とし、 をヤコビ記号とする。また、 とおく。( は奇数) （リュカ擬素数テスト）もし が素数であれば、 が成り立つ。また、 成り立たない場合 は素数ではない。 （強いリュカ擬素数テスト）もし が素数かつ である時、ある において、 または のいずれかを満たす。 64bit数の素数判定 20220927 Mizar/みざー 91

BPSW素数判定法(6) リュカ数列 のリュカ数列は以下のように で の 進 数表現を最上位bitから見ながら を素早く計算する漸化式がある。 の分子が奇数で割り切れない場合、右辺のように を足して偶数にしてから除算する。 64bit数の素数判定 20220927 Mizar/みざー 92

モンゴメリ剰余乗算・BPSW素数判定法の実装結果 Miller-Rabin素数判定法+128bit剰余乗算 https://yukicoder.me/submissions/792430 : PyPy3: 394ms https://yukicoder.me/submissions/794259 : C++17(gcc): 211ms https://yukicoder.me/submissions/794257 : Rust: 198ms Miller-Rabin素数判定法+モンゴメリ剰余乗算 https://yukicoder.me/submissions/793551 : C++17(gcc): 31ms https://yukicoder.me/submissions/793560 : Rust: 28ms BPSW素数判定法+モンゴメリ剰余乗算 https://yukicoder.me/submissions/793997 : C++17(gcc): 23ms https://yukicoder.me/submissions/794215 : Rust: 15ms モンゴメリ剰余乗算実装前の10倍速をこれで実現できた。C言語・同様のアルゴリズムで 14ms を計 測した提出例も存在したので、 C++言語でももう少し最適化の余地はあったかもしれない。 64bit数の素数判定 20220927 Mizar/みざー 93

ご清聴ありがとうございました 後編終了です 64bit数の素数判定 20220927 Mizar/みざー 94

宣伝 2022-09-30(金) 23:00- VRChatエンジニア作業飲み会 LightningTalkイベント ～3行のPythonコードで 1024bit の素因数分解ができてしまう?～ 「Shorのアルゴリズムと量子コンピュータ」(23:40頃-登壇予定) import math;G=math.gcd;a=g=2;n=int(input()) while g<3:print("?",a);t=pow(a,int(input())//2,n);a+=1;g=max(G(n,t-1),G(n,a)) print("!",g,n//g) 64bit数の素数判定 20220927 Mizar/みざー 95