はじめに
- 最近はじめた競技プログラムの某コンテスト問題の回答条件で、「期待値を求めろ。なお、解は有理数となるので(mod998244353)で出力せよ」というものがあって、解答例の期待値が3/2（＝1.5）の場合、出力を「499122178」としなければならかった。
- プログラミング以前に、回答の出力形式が意味不明だったので、一日掛けて解明した。せっかく苦労したので、これを共有したい。
modのおさらい
- プログラミング的には、5 mod 3 = 2 （5 ÷ 3の剰余は2）みたいな使い方が常識として知られていると思う。swiftやc++では、modを使わず、「5 % 3」と書くけどね。
- これを数学的な使い方で書くと、5 ≡ 8（mod3）みたいな感じ。「≡」は「合同記号」であり、発音は「ごうどう」。意味は、5mod3 = 8mod3 (= 2)ということ。つまり、自然数Aと自然数Bをそれぞれ自然数Cで割ったときの剰余が等しいときに「A ≡ B (mod C)」と表現する。このときCを「法」と呼び、法が文脈上明らかなときは、単に「A ≡ B」と書き、（mod C）は省略する。
modの公式
- （mod m）のとき A ≡ B かつ C ≡ D ならば、下記が成立する。以降、それぞれ「公式1.」「公式2.」「公式3.」と呼ぶ。
- A+C ≡ B+D （左記より、A-C ≡ B-Dも成立する)
- A×C ≡ B×D （左記より、$A^n ≡ B^n$も成立する)
- mとnが素であり、かつ、nがAとBの公約数であるとき、A/n ≡ B/n（「mとnが素である」とは、「mとnの公約数は1のみ」という意味）
- 公式2.と公式3.って実質的に同内容な気もするけど、説明しやすさから、分けて書いておく。
- 公式は、(mod 10)、つまり、自然数の1の位のみでの計算を考えれば、理解しやすいと思う。
- 例えば、公式3.について、3の倍数同士で1の位が同じである 45 ≡ 15 について考える。10と素である3で両辺を割ると、15 ≡ 5となり、正しいことが確認できる。
- もっと言えば、10と素である、3,7,11,13,･･･等で、両辺を割ることが出来れば、常に合同式が成立する。
- 例えば、7の倍数で、1の位が等しくなるものを探るため、7の倍数を順に見ていくと、{7,14,21,28,35,42,49,56,63,70,77,84,91....}となり、1倍〜10倍の1の位は、10個の自然数{7,4,1,8,5,2,9,6,3,0}が繰り返し出現する事が分かる。数列中の1の位が等しい 14 ≡ 84 について、7で割ると、2 ≡ 12 となる。
- 3の倍数についても、{3,6,9,12,15,18,21,24,27,30,33,36,39,....}であり、10個の自然数{3,6,9,2,5,8,1,4,7,0}が繰り返し出現している。9 ≡ 39 について、3で割ると 3 ≡ 13 となる。
- 公式2.について、(mod 998244353)での利用を考える。ループ処理なんかで、987654321を10回乗じるようなケースで、普通に10回乗じようとすると、途中で整数の桁が増えすぎてオーバーフローする。でも、公式2.の考え方(合同な他の数値に変えても結果は同じ)に従って、1回乗じる毎にmod計算により、桁の少ない合同な値に変換可能であり、例えば、$987654321^2 = 975461057789971041 ≡ 17678886 \pmod {998244353}$ のように桁を落とすことが出来るので、オーバーフローを回避できる。
- 下記のswiftでの解法①②は、整数値の上限制限がなければ結果は同じはずであるが、実際は制限があるため、解法①ではオーバーフローとなる。
let a = 987654321
let m = 998244353
//解法① オーバーフローする
// var ans = 1
// for _ in 0..<10 {
// ans *= a
// }
// ans = ans % m
// print(ans)
//解法② オーバーフローしない
var ans = 1
for _ in 0..<10 {
ans *= a
ans = ans % m
}
print(ans)
有理数とは
- これは、理系の人には常識だと思うけど、整数a、整数bを用いて、a/bの形で表すことが出来る数値のこと。だから、π(3.141592...)や2の平方根(1.41421356...)は有理数ではない。
有理数b/a(mod m)とは？
- まず、具体的なイメージを伝えると、1/2 ≡ 4 (mod 7) です。
- これを初めて見た人は、？？？？？...ケンカ売ってんのかっ！となると思う。
- でも、1 ≡ 8 (mod 7) はどうだろう？これは、何の問題もないでしょ？これを、法7と素である2で割ると、最初の合同式になるよ。
- 定義として言えば、b/a(mod m)とは、ax ≡ b (mod m)が成立するときのxの解を指す。
- 上記例で言えば、有理数1/2 は、2x ≡ 1 (mod 7)を満たすxの解となっている。また、1/2 ≡ 4 で示された4も同様に「2x ≡ 1 (mod 7)」のxの解の一つとなっている。よって、(mod 7)の世界では、1/2と4は合同となる。
- ここで、「はじめに」で出てきた、3/2(mod9982443539)について考えると、2x ≡ 3(mod 998244353)を満たすxを求めることとなる。法(998244353)の1/2は、499122176.5であり、これに1.5(=3/2)を足した499122178がxの解の一つとなる。
- たしかに、AtCoderでの出力例「499122178」と一致した。
- ふむ...一致はしたね。
- で、そもそも、有理数を(mod 998244353)で表示することに、どんな意味があるのかな？
有理数の回答を (mod998244353) で表示させる意味について
- 考えてみれば、解が有理数となる場合、「1/2」が正解なら、「2/4」も「3/6」も正解となる。
- でも、コンテストの実施者側としては、通常の回答受付は文字列出力のみであり、「2/4」も正答として受け付け得る場合、極端に言えば、「1000000000/2000000000」も受け付けざるを得なくなる。
- 想像するに、非常に多数の正解を受け付けるのは非効率なので、(mod998244353)で縛ることにより、「1/2」「2/4」「3/6」は受け付けず、「499122178」のみを受け付けることとした、ということかな？
- もちろん、分母の整数を最小の数にしろ、という制約を加えれば1/2のみが正解になるけど、その縛りを加えると、解答者にコードを無駄に書かせる必要が生じるから避けたいのかな？
- 次に、「998244353」という数値は何かというと、「大きな素数」である。
- 「大きな素数」だと便利な理由は、次の説明に書いてるから、もう少し付き合って欲しい。
有理数b/aを(mod 998244353)で出力するための考察
- 最初に言っておくと、c++とPythonの人は、AtCoderで使える関数が存在するので、自分でコードを書く必要ないです。残念ながらswift用の関数は無いみたいだから、自作するしかない。これに関する問題(ABC360e)が出された際のrated参加者で、swiftでのAC提出が無いのは、その辺が理由かな？
- コードを書く前に、そもそも、手作業でやる場合でも、一般的にどのように解けば良いのか？
- まず、1/a(mod m)について考える。これはaの逆元とも呼ばれ$a^{-1}$と書き、ax ≡ 1 (mod m) を満たすxの解である。
- 分かりやすい(mod 10)の場合を考える。
- $1^{-1}$ ≡ 1 、言うまでもなく、全てのmに対して、1×1 ≡ 1 (mod m)となる。
- $2^{-1}$ は存在しない。法10と2は素でないため、2x ≡ 1 を満たす解xは存在しない。感覚的にも、2にどんな数を乗じても、下一桁が1になることがないと分かるよね。
- $3^{-1}$ ≡ 7 、これは 3×7=21 より分かる。
- $4^{-1},5^{-1},6^{-1}$ は存在しない。
- $7^{-1}$ ≡ 3 、これは 7×3=21 より分かる。
- $8^{-1}$ は存在しない。
- $9^{-1}$ ≡ 9 、これは 9×9=81 より分かる。
- $10^{-1}$ は存在しない。
- $11^{-1}$ ≡ 1 、これは 11×1=11 より分かる。
- $12^{-1}$ は存在しない。
- $13^{-1}$ ≡ 7 、これは 13×7=91 より分かる。
- $14^{-1},15^{-1},16^{-1}$ は存在しない。
- $17^{-1}$ ≡ 3 、これは 17×3=91 より分かる。
- $18^{-1}$ は存在しない。
- $19^{-1}$ ≡ 9 、これは 19×9=171 より分かる。
- ....(なんだか、周期性がありそう)....
- つぎに、b/a(mod 10)を考えてみる。
- 例えば、5/3については、$5×3^{-1}$と書けるから、公式1.より、5/3 ≡ $5×3^{-1} ≡ 5×7 ≡ 35 ≡ 5$ となる。
- 実際に、3x ≡ 5(mod 10)を満たすxの解として、5は正しい。
- 次に、8/7について考えると、8/7 ≡ 8×3 ≡ 24 ≡ 4 であり、7x ≡ 8(mod 10)を満たすxとして、4は正しい。
- 上記より、b/a(mod m)は、aとmが素であるときに算出することが出来、その一つは、$b×a^{-1}$となる事が分かる。
- また、(mod 998244353)を考えれば、法が素数であるから、AtCoderで回答すべき有理数b/aについて、1 ≦ a ＜ 998244353 であれば法とaは素であり、aの逆数は存在するから、b/aは$b×a^{-1}$により算出可能。998244353が『大きな素数』である理由がここにある。
実際にコードを書いてみる(swift)
- まずは、aの逆数を算出するコードを書いてみよう。
- で、コードを書くに当たっては、拡張ユークリッド互除法を利用することになるんだが、詳しくは、wikiのモジュラ逆数の解説を見て欲しい。
- 結論から言うと、下記のベズーの等式を用いる。
ax + by = gcd(a,b)
※ gcd:最大公約数(greatest common divisor)
- 上記のベズーの等式は、「整数a,bの最大公約数gcd(a,b)に対して、上記を満たす整数(x,y)の組み合わせが存在する」というものであり、例えば、整数32,20の最大公約数4に対して、(x,y)=(2,−3)があり、検算すると、下記の通り。(x,y)の求め方は【おまけ】を見てほしい。
32 × 2 + 20 × (-3) = 64 - 60 = 4
- ここで、逆元の話に戻ると、下記の式を導出できることが分かる。
\displaylines{
n×n^{-1} ≡ 1 \pmod m \:\:\:※\:mは素数\\
n×n^{-1} + m×k ≡ 1 \pmod m \:\:\:※\:kは任意の整数\\
}
- 上記のmは素数としたので、gcd(n,m) = 1 となる。これにより、ベズーの等式が成立する$(n^{-1},k)$が存在するといえる。
n×n^{-1} + m×k = gcd(n,m)
- 実際に、【おまけ】の関数exGCDの引数に(n,m)を代入すると、$(n^{-1},k)$を得ることが出来る。(n,m)=(2,998244353)でexGCDを実行すると、$(n^{-1},k)$=(-499122176,1)になったよ。
- あれ？逆数が負の数？う〜ん．．．まぁ、ベズーの等式を満たす(x,y)の組は無数にあるからxが負の数になる事自体は間違ってないし、$n^{-1}$は、$n×n^{-1} ≡ 1 \pmod m$の答えだから、負数でも矛盾はしてないけどね。
- でも、見た目がよろしくないし、なにより、AtCoderの求める回答は、0..<998244353の範囲内だろうから、この間に収める必要があるね。
- swiftで実装すると
func invMOD(_ n:Int,_ mod:Int)->Int{
var inv = exGCD(n,mod).0
if inv < 0 {inv += mod}
return inv
}
let n = 2
let mod = 998244353
print(invMOD(n,mod)) //499122177
- 次に、目標としていた3/2(mod 998244353)を算出するけど、3×499122177(mod 998244353)より、1497366531 % 99824435 = 499122178 となる。
- AtCoder向けに(mod 998244353)を固定化して、有理数r/q(mod 998244353)を求めるコードは
func rqMOD_998244353(_ r:Int,_ q:Int)->Int{
let mod = 998244353
var inv_q = exGCD(q,mod).0
if inv_q < 0 {inv_q += mod}
return r * inv_q % mod
}
let (r,q) = (3,2)
print(rqMOD_998244353(r,q)) // 499122178
- 完成した！！
rqMOD_998244353で遊んでみる
- 下記の分数計算例について、(mod 998244353)でも成立するか、計算して確認する。
- 足し算：$1/2 + 1/3 = 5/6$
- 掛け算：$4/5 × 25/6 = 10/3$
let mod = 998244353
let add0 = rqMOD_998244353(1,2)
let add1 = rqMOD_998244353(1,3)
let addA = rqMOD_998244353(5,6)
print(add0,add1,addA) // 499122177 332748118 831870295
print(add0 + add1 % mod) // 831870295
let mlt0 = rqMOD_998244353(4,5)
let mlt1 = rqMOD_998244353(25,6)
let mltA = rqMOD_998244353(10,3)
print(mlt0,mlt1,mltA) // 399297742 166374063 332748121
print(mlt0 * mlt1 % mod) // 332748121
let zero = rqMOD_998244353(0,1)
let plus5 = rqMOD_998244353(5,1)
let minus5 = rqMOD_998244353(-5,1)
print(zero,plus5,minus5) // 0 5 -5
- 計算もバッチリ！
最後に
- ABC360の問題Eで有理数の回答を(mod 998244353)で出力しろ、とか書いてあるのを見たとき、サッパリ意味が分からなかったので、公式監修の入門書を読んだけど、有理数のmod出力なんて目次にも索引にも見当たらなくて、ガックリきた。いや、この本自体はいい本だと思うけどね。
- しゃあないので、modに関する高校数学レベルの「合同記号≡」の意味を調べることから始めたから、結構疲れたよ。
- でも、なんだかんだ苦労したけど、関数rqMOD_998244353を作るところまで漕ぎ着けたから、今週末はなかなか満足のいく週末だったね。
- swiftでAtCoderを頑張っている人が使ってくれると嬉しい。
【おまけ】ベズーの等式を満たす(x,y)を算出する関数exGCD(a,b)
ベズーの等式
a と b を 0 でない整数とし、d をそれらの最大公約数とする。このとき整数 x と y が存在して
ax + by = d
となる。さらに、
- d は ax + by と書ける最小の正の整数であり、
- ax + by の形のすべての整数は d の倍数である。
- その(x,y)の組み合わせをどのように求めるかというと、拡張ユークリッド互助法を用いることとなる。
- これは、最大公約数を求める方法であり、図で表すと、下記gif(wikiからの借り物)のようにa,bの整数の組を縦横とする長方形を作り、角から45度の斜線を引き、辺に当たる度に「斜線により生成された正方形を除いて『新しく作った長方形』でまた、斜線を引き始める」を繰り返したとき、『新しく作った長方形』が正方形になったとき、正方形の辺が最大公約数となる。
- 数値を用いて表すと、下記の流れとなる。当初数値をa0=32,b0=20として、
32(a0) = 20(b0⇒a1) × 1(商⇒s1) + 12(余り⇒b1)
20(a1) = 12(b1⇒a2) × 1(商⇒s2) + 8(余り⇒b2)
12(a2) = 8(b2⇒a3) × 1(商⇒s3) + 4(余り⇒b3)
8(a3) = 4(b3⇒a4) × 2(商⇒s4) + 0(余り⇒b4) 余り(b4)==0 -> ループ終了
- 上記よりgcd(32,20)=4(b3)が判明した。これを逆算していき、32(a0)と20(b0)を使ったベズーの公式をみたすx,yを算出する。
4(b3) = 12(a2) - 8(b2) × 1
= 12(b1) - ( 20(a1) - 12(b1) )
= 32(a0) - 20(b0) - ( 20(b0) - ( 32(a0) - 20(b0) ) )
= 32(a0) × 2 + 20(b0) × (-3)
- 上記より(x,y)は(2,-3)となった。
- このx,yを求めるコードをswiftで書くと以下の通り。関数exGCDは整数の組(a,b)を引数とし、戻り値はベズーの等式を満たす(x,y)を格納したタプルとなる。よって、ax + by が最大公約数となる。
func exGCD(_ a:Int,_ b:Int)->(Int,Int){
if (b == 0) {return (1,0)}
let temp = exGCD(b,a % b)
let q = a / b
return (temp.1 , temp.0 - q * temp.1)
}
let a = 32
let b = 20
let ans = exGCD(a,b)
print( ans.0, ans.1, a*ans.0 + b*ans.1 ) // 2 -3 4