えびちゃんの日記

えびちゃん(競プロ)の日記です。

mod についてふんわりとしか理解していない競プロ er のみなさん

ちゃんと理解していない人が大半だと思っています。あるいは、法が素数の場合に過学習して、法が合成数になった途端に不安になる人も多いでしょう*1

競プロ er の三大得意分野として「計算量削減」「mod なんとかの数え上げ」「なんか」があるはずなのに、big-O 記法や mod の性質やなんか*2に関して適当な理解をしている人が多いように見受けられるのはやや残念なところです。 水色から青程度であってもちゃんとやっていない人が多い気がします。AtCoder Beginner Contest においてパターンマッチで問題を解くマシーンになる上では必要ない(身につかない)部分かもしれないので仕方がない気もします。

今回は、適当な理解で適当にやっているライバルに差をつけましょうという感じの回です。適当な導入をしました。

「これを読むことで、解けるようになる問題がばりばり増えて、短期的にレートに寄与する」みたいな記事は目指していません。雰囲気でなんとなくやっている部分を減らし、必要に応じて自分で証明できるようになってもらうのを目指しています。

数式をがちゃがちゃして証明しているパートは、冗長になりそうな部分は折り畳んだりしていますが、必要に応じて展開して読んでもらえればよい気がします。 数式が苦手な人も、落ちついて調べたり考えたりしながら読んでもらえるとうれしいです。

数式慣れしていない人向けメモ

$\Z$ と書いて整数全体からなる集合 $\{\dots, -1, 0, 1, 2, \dots\}$ を表します。 「$x\in\Z$ とします」と言えば「$x$ を何らかの整数とします」のような意味になります。より直訳風に言えば「整数全体からなる集合に属する何かを取ってきて、それを $x$ と呼びます」のような感じです。

$\N^+$ は正整数全体からなる集合 $\{1, 2, \dots\}$ を表します。$\Q$ は有理数全体からなる集合 $\{\dots, -1/3, \dots, 0, \dots, 1/2, \dots, 1, \dots\}$ です。

$A \triangleq B$ は、「$A$ を $B$ と定義しますよ($A$ と書いたら $B$ のことですよ)」のような意味になります。たとえば先の話は $\Z \triangleq \{\dots, -1, 0, 1, 2\dots\}$ と書けます。$\triangleq$ の代わりに $\coloneqq$ や $\stackrel{\text{def}}{=}$ のように書く人もいます。

他にも、集合の記法などでわからないところがあれば、適宜調べてみるとよい気がします。 herb.h.kobe-u.ac.jp

定義

任意の正整数 $m\in\N^+$ と整数 $a\in\Z$ に対し、整数 $q$ および $0\le r\lt m$ が一意に定まり $a = qm+r$ と書けます。この $q$ を (quotient)、$r$ を あまり (remainder) と言います。

$m\in\N^+$ および $a, b\in\Z$ に対して、$a$ が $b$ と($m$として)合同 である ($a$ is congruent to $b$ modulo $m$) とは、$a-b$ が $m$ の倍数であることを言い、$a\equiv b\pmod{m}$ と書きます。これは上記における $r=0$ のことであり、$a-b = qm$ なる整数 $q$ が存在することと言い換えられます。 文脈から $m$ が明らかな場合は「$m$ を法として」や「$(\operatorname{mod}\, m)$」などは省略することも多いです*3

ここで、$\equiv$ は 同値関係 (equivalence relation) になっています。

同値関係の説明と、同値関係になっていることの証明

一般に、集合 $S$ 上の関係 $\sim$ が同値関係であるとは、

  • 任意の $a\in S$ に対して $a\sim a$(反射律 (reflexivity))
  • 任意の $a, b\in S$ に対して $a\sim b \iff b\sim a$(対称律 (symmetry))
  • 任意の $a, b, c\in S$ に対して $a\sim b \wedge b\sim c \implies a\sim c$(推移律 (transitivity))

を満たすということです*4。ここでは集合 $\Z$ 上の関係 $\equiv$ が同値関係となることを示したいので、$\equiv$ の定義を使って上記三つを示せばよいです。

(反射律)$a-a = 0 = 0\cdot m$ であり、$0\in\Z$ なので $a\equiv a$ が成り立つ。

(対称律)$\implies$ のみ示せば十分。$a\equiv b$ より、ある $q\in\Z$ に対して $a-b=qm$ と書ける。$b-a=(-q)\cdot m$ であり $-q\in\Z$ なので $b\equiv a$ が成り立つ。

(推移律)$a\equiv b$ と $b\equiv c$ より、整数 $q$, $q'$ が存在して $a-b=qm$, $b-c=q'm$ と書ける。$a-c = (a-b)+(b-c) = qm+q'm = (q+q')\cdot m$ であり、$q+q'\in\Z$ なので $a\equiv c$ が成り立つ。

反射律・対称律・推移律が成り立つことが示せたので、$\equiv$ は同値関係である。証明終わり。

ほぼ自明に感じられる命題でも、「なんとなくそんな気がするから」「成り立たないと直感に反するから」のような理由づけではなく、定義に従って証明することが大事です。そういう癖をつけておかないと、いざちゃんと証明する必要がある状況でも「なんとなく成り立ってほしいから成り立つと思う」とするしかできなくなってしまいます。

一般に、ある同値関係 $\sim$ に対して、$a$ と同値になるものの集合 $\{x\mid x\sim a\}$ を $a$ の 同値類 (equivalence class) と言い、$[a]$ のように書きます。この $a$(同値類を作るための基準として選んできた要素)のことを 代表元 (representative) と呼びます*5

さて、$\equiv$ による同値類を考えます。式で書くと下記のようになります。

$$ \begin{aligned} [a] &\triangleq \{x\in\Z\mid x\equiv a\pmod{m}\}. \end{aligned} $$

たとえば、$[0] = \{\dots, -m, 0, m, 2m, \dots\}$ や $[1] = \{\dots, -m+1, 1, m+1, 2m+1, \dots\}$ などとなります。$m=1$ であれば $[0] = [1]$ になっていることに注意しておくとよいかもしれません。あるいは、$[0] = [m] = [2m]$ などが成り立ちます。

さて、これも一般に、集合 $S$ 上の同値関係 $\sim$ によって作られる同値類の集合を $S/{\sim}$ と書き、$S$ の $\sim$ による 商集合 (quotient set) と呼びます。

ここでは $(\Z/{\equiv}) = \{[0], [1], \dots, [m-1]\}$ となります。 $\Z/{\equiv}$ のことを $\Z/m\Z$ や $\Z_m$ などと書いたりします。法が素数 $p$ のときは $\mathbb{F}_p$ と書いたりします。なぜ素数のときに限るのかという話の詳細は一旦ここでは触れません。$\mathbb{F}$ は field(たい)の頭文字で、法が素数のときには $\Z/m\Z$ が体になるからだと思います。

また、法 $m$ による $a$ の同値類ということを明示するために、以下では $[a]$ の代わりに $[a]_{(m)}$ と書く場合がありますが、一般的な記法ではないかもしれません。 いずれにせよ、$a\in\Z$ と $[a]_{(m)} = \{\dots, -m+a, a, m+a, \dots\}\in\Z_m$ を区別して考えておいた方がよいと思われる文脈では、そのように明示したいという意図です(特定の整数を持ってきて議論しているのか、$a$ と合同な整数であれば何でもよいという前提で議論しているのか、など)。慣れてきたら $[a]_{(m)}$ の意味合いで $a$ と書いてしまってもよいでしょう。単に区別したいだけなら $[a]$ でなくても $\overline{a}$ や $\underline{a}$ などでも構いません。

競プロ er が勝手に勘違いするかもしれないので補足ですが、代表元は $0\le a\lt m$ で選ばなければいけないという数学的な決まりはありません(普段は特に $0\le a\lt m$ 以外から選ぶ理由がないというだけ。競プロで言えばジャッジ側の都合に合わせているというくらい)。たとえば big-O 記法で $O(2n+3)$ のように書くのが記法としては全く問題ないことと、話としては似ているかもしれません。

また、$a = qm+r$ ($0\le r\lt m$) のとき、二項演算 $\bmod\colon \Z\times\N^+\to\Z$ を $a\bmod m = r$ で定義します。$a \equiv r \pmod{m}$ と混同しないように気をつけましょう。たとえば $5 \equiv 3 \pmod{2}$ は真ですが、$5\bmod 2 = 3$ は偽です。

(追記)文献によっては、$a \bmod m$ が同値類 $[a]_{(m)}$ を返す演算として定義されているものもあるようです。

演算

法 $m$ を固定して、商集合 $(\Z/{\equiv}) = \{[0], [1], \dots, [m-1]\}$ の各要素間における演算を定義していきます。

いわゆる modint(a) * modint(b) みたいなのに相当するパートです。 $(a \times b)\bmod m$ の代わりに $( (a\bmod m) \times (b\bmod m))\bmod m$ のように計算することの正当性に関連します。

また、以下で「$[a]$ を計算する」「$[a]$ を前計算で求めておく」などと言ったときは、必ずしも $0\le a\lt m$ の代表元を使うことを意味していません。適切に実装すれば $0\le a\lt m$ の範囲から取れるようなものだけ書いているつもりではありますが、たとえば $-m\le a\lt 0$ の範囲などの方が都合がよい人はそうしても問題なさそうです。とはいえ、わざわざ $10^{10^{100}}\cdot m\le a\lt (10^{10^{100}}+1)\cdot m$ のような範囲から持ってくるのは、計算量やパソコンがこわれるのでやめてほしいです。

加減算と乗算

さて、加減算と乗算を以下のように定義します。複号同順です。

  • $[a] \pm [b] \triangleq [a\pm b]$,
  • $[a] \times [b] \triangleq [a\times b]$.

この定義が代表元の取り方に依存していないことを確かめましょう。

なぜ確かめたいのかという話

たとえば、$[1] + [2] \neq [m+1] + [3m+2]$ だと困るので、そんなことないですよーというのを示したいということです。$[1] = [m+1] = \{\dots, -m+1, 1, m+1, \dots\}$ なのに、$1$ を使うのか $m+1$ を使うのかで結果が変わってしまっては、妥当な定義ではなさそうです。代表元を取ってきて演算を定義するような状況ではいつもこういうことに気を払う必要があります。

確かめる

$a, a'\in [a]$ と $b, b'\in [b]$ に対し、$[a]\pm [b] = [a']\pm [b']$ と $[a]\times [b] = [a']\times [b']$ を示したいです。

まず $\pm$ の方から示します。 $a, a'\in[a]$ と $b, b'\in[b]$ に対し、$x\in[a]\pm[b] \implies x\in[a']\pm[b']$ を示せば十分です(対称性から $x\in[a]\pm[b] \impliedby x\in[a']\pm[b']$ もわかるので)。 定義から $$ [a]\pm [b] = [a\pm b] = \{x\in\Z\mid x\equiv a\pm b\pmod{m}\} $$ なので、$x\in[a]\pm [b]$ ならば $x\equiv a\pm b\pmod{m}$ です。 このとき $x\in[a'+b']$ が成り立つことを示すのが目標です。すなわち、$x\equiv a'\pm b'\pmod{m}$ を目指します。

$a, a'\in[a]$ から $a\equiv a'$ であることに注意しておきます。すなわち、$q_a\in\Z$ が存在して $a-a'=q_a\cdot m$ と書けます。同様に $b-b' = q_b\cdot m$ とします。 また、$x-(a\pm b) = qm$ なる $q\in\Z$ が存在します。さて、 $$ \begin{aligned} x-(a'\pm b') &= x-( (a-q_a\cdot m) \pm (b-q_b\cdot m)) \\ &= x-(a\pm b) + (q_a\pm q_b)\cdot m \\ &= qm + (q_a\pm q_b)\cdot m \\ &= (q+q_a\pm q_b)\cdot m \end{aligned} $$ であり、$q+q_a\pm q_b\in\Z$ なので、$x\equiv a\pm b\implies x\equiv a'\pm b'$ が示せました。

同様に $\times$ の方も示します。やるべきことは同じで、$x\equiv(a\times b) \implies x\equiv(a'\times b')$ を示します。$q_a$, $q_b$, $q$ を同様に定義します。 $$ \begin{aligned} x-(a'\times b') &= x-( (a-q_a\cdot m)\times (b-q_b\cdot m)) \\ &= x-(a\times b - aq_b \cdot m - bq_a\cdot m + q_aq_b\cdot m^2) \\ &= (x-a\times b) + (aq_b + bq_a - q_aq_bm)\cdot m \\ &= (q+aq_b + bq_a - q_aq_bm)\cdot m \end{aligned} $$ で、$q+aq_b + bq_a - q_aq_bm\in\Z$ なので $x\equiv a\times b \implies x\equiv a'\times b'$ が示せました。

以上より、$[a]\pm[b] = [a+b]$ や $[a]\times [b] = [a\times b]$ が破綻していない定義であることが確かめられました。

代表元同士の計算で済むようになったので、下記のような性質は整数の性質からそのまま従うことになります。

  • $([a]+[b])+[c] = [a]+([b]+[c])$
  • $[a]+[b] = [b]+[a]$
  • $([a]\times[b])\times[c] = [a]\times([b]\times[c])$
  • $[a]\times[b] = [b]\times[a]$
  • $[a]\pm[0] = [a]$
  • $[a]\times[0] = [0]$
  • $[a]\times[1] = [a]$
  • $[a]\times([b]+[c]) = [a]\times[b]+[a]\times[c]$

冪乗

$[a]$ に関する冪乗を、(整数に関する冪乗と同様に)下記のように定義します。 $$ [a]^k \triangleq \underbrace{[a]\times[a]\times\cdots\times[a]}_{k\text{ many}}. $$ 再帰的に書けば、$[a]^0 = [1]$ および $[a]^k = [a]^{k-1}\times[a]$ ($k\gt 0$) です。

さて、結合法則から、$[a]^{2k} = ([a]^k)\times([a]^k)$ および $[a]^{2k+1} = ([a]^k)\times([a]^k)\times[a]$ が成り立ちます。 よって、半分のサイズを解けば定数時間で元の問題も解けることになり、$[a]^k$ は $O(\log(k))$ 時間で求められることがわかります*6。俗に言う「繰り返し二乗法」です。

念のため補足ですが、ここまでの計算においては法が素数であることは用いていないため、法が合成数や $1$ だった場合でも適用することが可能です。 代表元を $0\le a\lt m$ で選ぶ必要がある状況においては、$m=1$ のケースで $[0]^0$ を $[1]$ ではなく $[0]$ と答える必要があることに注意しましょう。

なお、ここで定義したのは $[a]_{(m)}^k$ であって $[a]_{(m)}^{[k]_{(m)}}$ ではないことに注意してください。$ab\bmod m = ( (a\bmod m) \times (b\bmod m))\bmod m$ などからの類推から、なんでも $m$ で割ったあまりを使おうとする人がよくいますが、そういうのはだめです。 実際、$3^2 \bmod 5 = 4$ ですが $3^7 \bmod 5 = 2$ です。$3^6 \bmod 5 = 3^{10} \bmod 5 = 3^{14} \bmod 5 = \dots = 4$ は成り立ちます。 理由は後述します。

代表元を任意に取ってきて $[a]_{(m)}^{[k]_{(m)}} \triangleq [a^k]$ で定義しようとすると破綻するわけです。先ほどの加減算のところでやったような確認が必要であることがわかりましたね。

除算

さて、四則演算のうち $[a]\div[b]$ に相当するものはないのかと思うのはある程度妥当でしょう。 一方で、$[a]\div[b] \triangleq [a\div b]$ という定義は妥当ではなさそうです。たとえば $m=4$ として $[3]\div[2] = [3/2] = [7/10] = [7]\div[10]$ のような式は役に立たなさそうです。そもそも $[3/2]$ や $[7/10]$ のような同値類は定義されていないです($3/2, 7/10\notin\Z$ なので)。

それよりも、「除算に成り立ってほしい等式、な〜んだ?」という方針で考えるのがよさそうです。 有理数で考えると、$a\in\Q$ かつ $b\in\Q\setminus\{0\}$ に対して $(a\div b)\times b = a$ が成り立ちます。乗算の逆演算ということです。これは成り立っていてほしいです。あるいは、$a\div b = a\times b^{-1}$ なる要素 $b^{-1}\in\Q$ が存在してほしいです。 $b\times b^{-1}=1$ が成り立ちます。このように、ある元 $x$ に掛けると $1$ になる元を、$x$ の乗法逆元と呼びます(演算が明らかな場合は単に逆元と呼んだりします)。$b^{-1}$ は $b$ の乗法逆元です。$1$ のように「掛けても相手の値が変わらない元」のことを乗法単位元と呼びます(これも文脈によっては単に単位元と呼びます)。

そこで、まず $[b]^{-1}$ を「$[c]$ であって $[c]\times[b] = [1]$ を満たすもの」と定義し、$[a]\div[b] \triangleq [a]\times[b]^{-1}$ と定義します。ただし、$[b]^{-1}$ が存在しない場合は $[a]\div[b]$ は定義されないものとします。どの法における乗法逆元なのかを明示して、$[b]^{-1}_{(m)}$ のように書いたりもすることにします。

具体例を挙げておくと、$m = 9$ のとき、$[5]\times [2] = [10] = [1]$ であり、$[5]^{-1}_{(9)} = [2]$ です。一方、$m = 10$ のとき、$[5]\times [x] = 1$ なる $x$ は存在しないため、$[5]^{-1}_{(10)}$ は定義されません。

さて、こういう定義のされ方で気にするべきこととしては、

  • そうした $c$ が存在するのはどういうときか?
  • そうした $c$ は一意か?
    • $[c]$, $[c']$ に対して $[c]\times[b] = [c']\times[b] = [1]$ となるとき $[c]=[c']$ か?
    • あるいは、$c$, $c'$ に対して $[c]\times[b] = [c']\times[b] = [1]$ となるとき $c\equiv c'$ か?
  • そうした $c$ はどのように求めるのか?
    • 計算量がバカみたいにかかると嫌である。
    • たとえば $0$, $1$, ..., $m-1$ から全探索というのは非常に嫌である*7
      • 探索範囲がこれでいい理由は $\times$ が代表元の計算に帰着できることから*8

があると思います。

存在性

$[c]\times[b] = [1]$ というのは、$c\times b \equiv 1 \pmod{m}$ のことです。 定義から、ある $q\in\Z$ が存在し $cb - 1 = qm$ が成り立ちます。 移項して $cb - qm = 1$ です。

さて、これを満たす $(c, q)$ のペア(実際には $c$ だけわかれば十分ですが)を求めたいです。「与えられた $a$, $b$ に対して $ax-by=1$ を満たす $(x, y)$ を求めよ」という形式の方が見慣れているでしょうか。以下では後者の文字を用いることにし、後々になってから適宜置き換えることにしましょう。

$ax-by$ の形で表せる最小の正整数を考える

$S = \{ax-by\mid x, y\in\Z, ax-by\gt 0\}$ とおきます。 $a=b=0$ のとき $S=\emptyset$ です。以下、$a$ と $b$ の少なくとも一方は $0$ でないとします。

$x, y \in\{-1, +1\}$ から適切なものを選べば $|a|+|b|\gt 0$ を作ることができるため、$S\ne\emptyset$ です。 よって最小元が存在するため、$\min(S) = d = as-bt$ とおきます($S$ の最小元を $d$ とおき、それを達成する $(x, y)$ を $(s, t)$ とおくということです)。

$d$ による商とあまりを考え、$a = qd+r$ ($0\le r\lt d$) とおくと、 $$ \begin{aligned} r &= a - qd \\ &= a - q(as - bt) \\ &= a(1-sq) - (-bq)t \end{aligned} $$ となり、$r$ も $ax-by$ の形で書けることになります。$r \lt d$ より、$r\in S$(つまり $r\gt 0$)だと最小性 $\min(S)=d$ に矛盾します。よって $r = 0$ とわかります。 あまりが $0$ になったので、$d$ は $a$ の約数であることがわかります。同様にして $b$ の約数でもあることがわかります。

ここまでで、$d$ は $a$ と $b$ の約数であることがわかりました。 さて、$c$ を $a$ と $b$ の任意の約数とし、$a = ca'$, $b = cb'$ とおくと、 $$ \begin{aligned} d &= as + bt \\ &= ca's + cb'u \\ &= c(a's + b't) \end{aligned} $$ です。$d \gt 0$ なので $d\ge c$ が従い、$d$ は $a$ と $b$ の最大公約数であることがわかりました。 すなわち、$\min(S) = \gcd(a, b)$ です。

上記より、$ax-by$ の形で表せる最小の正整数は $\gcd(a, b)$ であることがわかりました。 よって、$\gcd(a, b) \le 1$ である必要があります。 いま、$a=b=0$ ではないことから $\gcd(a, b)$ は正なので、すなわち $\gcd(a, b) = 1$ のことですね。

文字を置き換え直せば、$\gcd(b, m) = 1$ が必要であることがわかりました(既知の文字が $b$, $m$ で、未知の文字が $c$, $q$ だったことに注意)。 別の言い方をすると、$b$ と $m$ が互いに素であることが必要です。 これが十分でもあることは、後で実際に構成することで示します。

一意性

$c$, $c'$ に対して $[b]\times [c] = [b]\times[c'] = [1]$ とします。 このとき $[c] = [c']$ を示したいです。

交換法則や結合法則などは事前に示したので、 $$ \begin{aligned} [b] \times [c] &= [b] \times [c'] \\ ([c] \times [b]) \times [c] &= ([c] \times [b]) \times [c'] \\ [1] \times [c] &= [1] \times [c'] \\ [c] &= [c'] \end{aligned} $$ となり、示されました。

$[c]$ が一意であることがわかったので、$[b]^{-1} = [c]$ を満たすような $c$ を $b^{-1}$ や $b^{-1}_{(m)}$ と書くことにします。 $[c]$ は一意なものの $c$ 自体の選び方には任意性があり($-m+c$ や $3m+c$ などもよいので)、どれを選んでも構いませんが、ここでは便宜上 $0\le b^{-1}\lt m$ となるように選ぶことにしておきます。

構成

さて、$ax+by=\gcd(a, b)$ なる $(x, y)$ を構成します。先ほどと $by$ の符号が変わっていますが本筋としては影響ありません。 互除法に少し手を加えたアルゴリズムによって求めることができます。

fn ecl(a: i64, b: i64) -> (i64, i64) {
    let (mut x0, mut y0, mut r0) = (1, 0, a);
    let (mut x1, mut y1, mut r1) = (0, 1, b);
    while r1 != 0 {
        let q = r0.div_euclid(r1);  // this means floor(r0 / r1)
        let (x2, y2, r2) = (x0 - q * x1, y0 - q * y1, r0 - q * r1);

        // (*)

        x0 = std::mem::replace(&mut x1, x2);  // this means x0, x1 = x1, x2
        y0 = std::mem::replace(&mut y1, y2);
        r0 = std::mem::replace(&mut r1, r2);
    }
    // at this point, r0 == gcd(a, b)
    (x0, y0)
}

正当性に関して

a * x0 + b * y0 == r0, a * x1 + b * y1 == r1 がループ不変条件として成り立ちます。ループ開始前では確かに成り立っています。 下記では、(*) の時点での x0 などの値を用いて証明をしていきます。

(本来は、添字を使い回さずに $(x_i, y_i, r_i)$ のような形で議論するのがよさそうな気もします。)

上記の仮定の下で、a * x2 + b * y2 == r2 が成り立つことを示します。

$$ \begin{aligned} ax_2 + b y_2 &= a(x_0 - qx_1) + b(y_0 - qy_1) \\ &= (a x_0 + by_0) - q(ax_1 + by_1) \\ &= r_0 - q r_1 = r_2. \end{aligned} $$

よって、帰納的にループ不変条件が示されたことになります。

あとは、いつかしらは $r_2 = 0$ になることと、その時点で $r_1 = \gcd(a, b)$ になることを示せば十分です(ループ不変条件から返り値が所望の条件を満たすことが従うため)。 $r$ に関して注目すると、ただの互除法に他ならないので、互除法の正当性と計算量解析をすればよいです。 競プロ er は互除法の証明も適当にやっているはずなので、ちゃんと考え直します。

停止性を示すだけなら $r_2 = (r_0 \bmod r_1)$ から $0\le r_2 \lt r_1$ であることさえ言えばよさそうです。 計算量解析もしたいのでもう少し厳しく評価します。一般に、$0 \lt x \le y$ に対して $(y\bmod x) \lt y/2$ が成り立ちます。(添字の入れ替わりなどに注意すると)ループ二回ごとに r0 が半分未満になることがわかるので、$O(\log(\min(\{a, b\})))$ 回程度のループで停止します。 この性質については付録で証明します。

あとは正当性です。 最後のループの (*) の時点で $r_1$ が $a$ と $b$ の公約数であることさえ示せば、以前示した $\min(S) = \gcd(a, b)$ から、$r_1 = \gcd(a, b)$ は従いそうです(最大公約数でないと $S$ における最小性に矛盾するため)。

$r_0 = q r_1 + r_2$ のとき $\gcd(r_0, r_1) = \gcd(r_1, r_2)$ であることを示します。 $c$ を $r_0$ と $r_1$ の公約数とし、$r_0 = cr_0'$, $r_1 = cr_1'$ とします。このとき、$r_2 = r_0 - q r_1 = c(r_0' - q r_1')$ より、$c$ は $r_2$ の約数でもあります。 同様に、$d$ を $r_1$ と $r_2$ の公約数とすると、$d$ は $r_0$ の約数になります。 よって、$r_0$ と $r_1$ の公約数全体の集合は、$r_1$ と $r_2$ の公約数全体の集合と等しいです。

ところで、$r_0$ と $r_1$ の初期値は $a$ と $b$ であり、$a$ と $b$ の公約数全体の集合は、(上記の性質が推移的に成り立つので、)各ループにおいて $r_1$ と $r_2$ の公約数全体の集合と等しいです。 また、最後のループにおいて、$r_1$ と $r_2$ の公約数全体の集合は、($r_2 = 0$ なので)$r_1$ の約数全体の集合に等しいです。 以上より、$r_1$ は $a$ と $b$ の公約数となります。

よって、手続きの終了時において $ax_0 + by_0 = \gcd(a, b)$ となることが示されました。

さて、これにより、$\gcd(b, m) = 1$ であれば、$cb - qm = 1$ なる $c$ を $O(\log(b))$ 時間で求められるようになりました。 $\gcd(b, m)$ も同時に求められるので、存在性の判定についても行えます。

note: 互除法で使っている行列 $\left(\begin{smallmatrix} x_0 & y_0 & r_0 \\ x_1 & y_1 & r_1\end{smallmatrix}\right)$ は Stern–Brocot 木などとも関係が深いです。近いうちにまた話します。

別の求め方

競プロ界隈では、互除法に基づく逆元の計算はやや敬遠され、Fermat の小定理に基づく方法が初心者に好まれがちな印象があります。 予想される理由としては、

  • 互除法による方法は実装の理解が難しい
  • Fermat の小定理さえ fact として認めれば実装は比較的簡単
  • こと AtCoder においては法が素数であるケースが多く、うまくいくケースが多い
    • もちろん素数であれば常に逆元が存在するわけではないのだが、それすら無視できる場合が多い
  • Fermat の小定理に基づく方法しか知らない初心者が他の初心者にそれを教えている

などでしょうか。

というわけで、Fermat の小定理の話をします。 法が素数 $p$ のとき、$a$ と $p$ が互いに素(すなわち $a\not\equiv 0 \pmod{p}$)であれば $a^{p-1} \equiv 1 \pmod{p}$ であるという定理です。 $(a^{p-2})\cdot a \equiv 1\pmod{p}$ であることから、$[a]^{-1} = [a]^{p-2}$ として求められるというわけです。 Fermat の小定理は Euler の定理の特殊ケースと見なすことができるため、Euler の定理の証明のみを付録に載せます。

法が合成数のときにも Fermat の小定理を用いようとする人もいますが、前提条件をきちんと確認する癖をつける方がよいでしょう。

速度の観点の話など

互除法に基づく方法では、割り算する値が変わるのでコンパイラの最適化の恩恵を受けにくい。 一方、Fermat の小定理に基づく方法においては同じ値によるあまりのみ求めればよいため、(特に $998244353$ のように予めコンパイル時に値が決まっているような場合は)定数除算最適化の恩恵を受けやすく、定数倍が軽くなることが期待される。

とはいえ、Fermat の小定理に基づく方法では毎回 $\Theta(\log(m))$ 回のループが発生するのに対し、互除法のループ回数は値によっては数回で終わることも多く、最悪時も $\Theta(\log(m))$ ではなく $\Theta(\log(b))$ で抑えられるため、速度面に関しては状況によるのではなかろうかという感じがする。

注意

法が素数であっても $[a]$ の逆元 $[a]^{-1}_{(p)}$ が存在するとは限らないことに注意してください。 たとえば $[a] = [0] = [p]$ などのケースです。 素数の法 $p$ において除算をしたくなった際、「$[a]$ の逆元を考える。ただし $a\ne 0$ とする」と書いても $a=p$ のケースなどは除外できていなさそうなので気をつけましょう。「$[a]\ne [0]$ とする」などが正しそうです。

法を $p$ とした数え上げの問題において、「ゼロ除算は起きないから平気」などと油断して $p$ で除算しないように気をつけましょう。 人が「法が $998244353$ だったら全方位木 DP が逆元を仮定して簡単にできた」と言っているときは十中八九 EDPC V の話をしているのですが、この問題においては割りたい値が常に $M$ と互いに素となるとは限らないため、そうした考察は正しくないです。 たとえば、頂点 $1$ を黒く塗る方法が $M$ 通り ($1\le M\le 10^{18}$) となるような木を構築してみましょうという問題を考えてみるとよいかもしれません。

ところで、$m=1$ であれば $\gcd(0, 1) = 1$ ですから、$[0]^{-1}_{(1)}$ は定義できて $[0]$ ですね。$[0]_{(1)} \times [0]^{-1}_{(1)} = [1]_{(1)} = [0]_{(1)}$ が成り立ちます。$[0]_{(1)} = \{\dots, -1, 0, 1, \dots\} = \Z$ です。

基本編

基本となる定理や操作などです。必ずしも「理解などが簡単」という意味合いで「基本」と言っているわけではないので、わからなくても焦らず手を動かしたりしてみてください。

中国剰余定理

いわゆる CRT (Chinese remainderリメインダ theorem) と呼ばれているものです。

中国の算術書『孫子算経』に

今物が有るが、その数はわからない。三つずつにして物を数えると、二余る。五で割ると、三余る。七で割ると、二余る。物はいくつあるか?

答え:二十三。

という問題が載っているらしいです。見慣れた表記にすると次のような感じです。

$$ \begin{cases} x \equiv 2 \pmod{3}, \\ x \equiv 3 \pmod{5}, \\ x \equiv 2 \pmod{7}. \end{cases} $$

実際、$23 = 7\times 3 + 2 = 4\times 5 + 3 = 3\times 7 + 2$ で正しそうです。 しかし、他にも $x=128$ などでも条件を満たしそうです。 そこで、気になることを問題の形式で書くと次のようになります。

ある整数 $x$ があり、$x \equiv r_0 \pmod{m_0}$ かつ $x\equiv r_1 \pmod{m_1}$ を満たすという。 $x$ としてありうる整数全体の集合 $X$ を求めよ。また、それが空集合になる条件も求めよ。

さらに、条件式が 3 つ以上のケースについてはどうか。

世間的には中国剰余定理と言った場合には、$\gcd(m_0, m_1)$ の仮定をして、「$0\le x\lt m_0m_1$ を満たす $x$ であって上記の条件を満たすものが存在する」という主張を指しそうです。 競プロ界隈においては、特にその仮定はせずに、また存在性の判定だけではなく解の構成も行うアルゴリズムを指して「CRT」と呼んでいるケースが多そうです*9

がんばって求めます

$X$ が空でないとして $x\in X$ を固定すると、 ある整数 $q_0$, $q_1$ が存在して $x = q_0m_0 + r_0 = q_1m_1 + r_1$ です。 ここから $q_0m_0 - (r_1 - r_0) = q_1m_1$ であり、$q_0m_0 \equiv r_1-r_0 \pmod{m_1}$ が得られます。

ここで両辺を $m_0$ で割って $q_0$ に関する式にしたいですが、$m_0$ と $m_1$ が互いに素とは限らないので困ってしまいます。 そこで、互いに素となる問題に帰着させることを考えます。

$g = \gcd(m_0, m_1)$ とおき、$m_0 = gm_0'$, $m_1 = gm_1'$ とします。 $q_0gm_0' + r_0 = q_1gm_1' + r_1$ となります。 このとき、$r_1 - r_0 = g(q_0m_0' - q_1m_1')$ なので、$r_1-r_0$ が $g$ の倍数である必要があります。 そうでないとき $X = \emptyset$ です。

以下、$r_1 - r_0$ は $g$ の倍数とします。整数 $r$ を $r = (r_1-r_0)/g$ で定義しましょう。 $q_0gm_0' - gr = q_1gm_1'$ ですから、$g\gt 0$ で割って $q_0m_0' - r = q_1m_1'$ を得ます。 すなわち、$q_0m_0' \equiv r \pmod{m_1'}$ であり、$q_0 \equiv (m_0')^{-1}\cdot r\pmod{m_1'}$ です。

以上より、ある整数 $q$ が存在して $q_0 = (m_0')^{-1}_{(m_1')}\cdot r + qm_1'$ と書け、 $$ \begin{aligned} x &= q_0m_0 + r_0 \\ &= ( (m_0')^{-1}_{(m_1')}\cdot r + qm_1')\cdot m_0 + r_0 \end{aligned} $$ となります。

逆に、任意の $q\in\Z$ に対して、この $x$ が問題文の条件を満たすことを示します。$x\equiv r_0 \pmod{m_0}$ は明らかなので、$x\equiv r_1\pmod{m_1}$ を示します。

$(m_0')^{-1}_{(m_1')}\cdot m_0'$ は、ある整数 $i$ に対して $im_1' + 1$ と書けることに注意しておきます。 $$ \begin{aligned} x &= ( (m_0')^{-1}_{(m_1')}\cdot r + qm_1')\cdot m_0 + r_0 \\ &= (m_0')^{-1}_{(m_1')}\cdot r\cdot gm_0' + qgm_0'm_1' + r_0 \\ &= (im_1'+1)\cdot gr + qm_0'm_1 + r_0 \\ % &= im_1 + gr + qm_0'm_1 + r_0 \\ &= (gr + r_0) + (ir+qm_0')\cdot m_1 \\ &= r_1 + (ir+qm_0')\cdot m_1. \end{aligned} $$ $ir+qm_0'\in\Z$ なので、$x\equiv r_1 \pmod{m_1}$ が確かめられました。

よって、任意の $q\in\Z$ に対して $$ x = ( (m_0')^{-1}_{(m_1')}\cdot r + qm_1')\cdot m_0 + r_0 $$ であり、$m_0m_1' = m_0m_1/g$ に気をつけて $$ x \equiv (m_0')^{-1}_{(m_1')}\cdot m_0\cdot r + r_0 \pmod{m_0m_1/g} $$ となります。

まとめると、$g = \gcd(m_0, m_1)$ として以下が成り立ちます。 $$ X = \begin{cases} \emptyset, & \text{if }r_0 \not\equiv r_1\pmod{g}; \\ [(m_0')^{-1}_{(m_1')}\cdot m_0\cdot r + r_0]_{(m_0m_1/g)}, & \text{if }r_0 \equiv r_1\pmod{g}. \end{cases} $$

$X$ の代表元 $x$ を $0\le x\lt m_0m_1/g$ から取ってくるアルゴリズムに関しては、ACL math の解説記事などをご覧ください。

rsk0315.hatenablog.com

条件式三つ以上の場合に関しては、条件式二つの場合の答えが $x \equiv r \pmod{m}$ という形で得られたため、逐次的に求められます。 計算量は、条件式の個数を $n$ として $O(n+\sum_i \log(m_i))$ 時間ですが、$m$ がすぐにオーバーフローしそうなので困ります。

Garner’s algorithm

条件を満たす最小の値 $x$ と与えられた整数 $\mu$ に対する $(x\bmod\mu)$ を(オーバーフローさせないように気をつけながら)求めるアルゴリズムも存在し、$O(n^2 + \sum_i \log(m_i))$ 時間です。Garner’s algorithm などと呼ばれていそうです。文脈によっては、先の $O(n+\sum_i \log(m_i))$ の方を指してそう呼んでいることもありそうです。

kopricky.github.io

Garner’s algorithm の紹介

条件式で与えられる $m_i$ たちは、いずれも互いに素であるとします。 $$ \begin{aligned} x &= s_0 + s_1\cdot m_0 + s_2\cdot (m_0m_1) + \dots + s_{n-1}\cdot (m_0m_1\dots m_{n-2}) \\ &= \sum_{i=0}^{n-1} s_i\cdot \prod_{j=0}^{i-1} m_j \end{aligned} $$ の形で管理することを考えます。ただし、$0\le s_i\lt m_i$ です。

条件式を順に処理していき、$x \equiv r_i\pmod{m_i}$ を処理すると $s_i$ が求まります。 $i = k$ での処理は以下の通りです(この時点で $s_0$, ..., $s_{k-1}$ は既知です)。 $$ \begin{aligned} x \equiv \sum_{i=0}^k s_i\cdot \prod_{j=0}^{i-1} m_j &\equiv r_k \pmod{m_k} \\ r_k - \sum_{i=0}^{k-1} s_i\cdot \prod_{j=0}^{i-1} m_j &\equiv s_k\cdot \prod_{j=0}^{k-1} m_j \\ \left(r_k - \sum_{i=0}^{k-1} s_i\cdot \prod_{j=0}^{i-1} m_j\right)\cdot \left(\prod_{j=0}^{k-1} m_j\right)^{-1} &\equiv s_k. \end{aligned} $$ ここで、$\prod_j m_j$ などは $m_k$ で割ったあまりとして計算することで、オーバーフローを回避できます。 一方で、$m_k$ で割ったあまりとして計算する関係で、これらの GCD などを計算することが(陽に素因数分解をするなどを経由しないと)困難です。 そのため、入力の時点でこれらが互いに素であることを仮定しています。

さて、最終的に $\sum_i s_i\cdot \prod_j m_j$ の形が得られたら、これを $\mu$ で割ったあまりとして計算していけばよいです。 適切に実装することで、計算過程の値は $(\max(M)-1)^2$ 以下で抑えられます。 ここで、$M = \{m_0, m_1, \dots, m_{n-1}, \mu\}$ です。 計算量を $O(\log(\max(M)))$ 倍悪化させれば、乗算を加算の繰り返しとして計算することで $2\cdot(\max(M)-1)$ 以下で抑えることも可能そうです。

計算量は、$i$ 番目の条件式の処理で $O(n+\log(m_i))$ 時間かかるため、全体では $O(n^2+\sum_i \log(m_i))$ 時間となります。

法の変換

$a \equiv b \pmod {m}$ という式が得られているときに $a'\equiv b' \pmod{m'}$ という別の式を導けたらうれしいなぁという話です。 たとえば、二つのタイムスタンプがあって、24 時間表記でそれらの hour が等しければ、12 時間表記にしても(午前・午後は無視して)hour は等しくなるということがわかると思います。逆に、12 時間表記で hour が等しくても 24 時間表記に直したときに等しいとは限らないですね(5 時と 17 時だった場合など)。

また、$\Z_{24}$ から $\Z_{12}$ を構成する方法としては、午前・午後の同一視の他に、「1 時間を 120 分ということにする」というものも考えられそうです。 あるいは文字盤が 2 時間刻みで書かれている時計をイメージする方が直感的でしょうか。

表にしてみます。

$\Z_{24}$ $\Z_{12}$ (A) $\Z_{12}$ (B)
$[0]$ $[0]$ $[0]$
$[1]$ $[1]$ N/A
$[2]$ $[2]$ $[1]$
$[3]$ $[3]$ N/A
... ... ...
$[10]$ $[10]$ $[5]$
$[11]$ $[11]$ N/A
$[12]$ $[0]$ $[6]$
$[13]$ $[1]$ N/A
... ... ...
$[22]$ $[10]$ $[11]$
$[23]$ $[11]$ N/A

(A) の方では $[10]_{(12)} = [10]_{(24)}\cup[22]_{(24)}$ を用いているのに対し、(B) の方では $[10]_{(24)} = [5]_{(12)}$ を用いて対応づけています。 さて、だいたいイメージはできたかもしれないですが、以下のような操作をすることが可能そうです。

  • $a\equiv b \pmod{km} \implies a\equiv b \pmod{m}$ ... (A)
  • $ka\equiv kb \pmod{km} \iff a\equiv b\pmod{m}$ ... (B)

法が $0$ になることはないので、$k\ne 0$ です。証明は、定義に立ち返ってやれば単純そうです。

(A) の証明

ある整数 $q$ に対して $a-b=q\cdot km$ が成り立つ。 このとき、$a-b = (qk)\cdot m$ も成り立ち、$qk\in\Z$ なので $a\equiv b\pmod{m}$ となる。

(B) の証明

($\implies$) ある整数 $q$ に対して $ka-kb=q\cdot km$ が成り立つ。 このとき、$k\ne 0$ で割ることで $a-b = q\cdot m$ も成り立ち、$q\in\Z$ なので $a\equiv b\pmod{m}$ となる。

($\impliedby$) ある整数 $q$ に対して $a-b=q\cdot m$ が成り立つ。 このとき、$k\ne 0$ を掛けることで $ka-kb=(kq)\cdot m$ が成り立ち、$kq\in\Z$ なので $ka\equiv kb\pmod{km}$ となる。

ところで、$\gcd(l, m) = 1$ のとき $la\equiv lb\pmod{m} \iff a\equiv b\pmod{m}$ が成り立ちます($l^{-1}$ を掛ければよい)。 よって、(B) と合わせて $\gcd(kl, km) = k$ のとき $kla\equiv klb\pmod{km} \iff a\equiv b\pmod{m}$ となります。

あるいは、$ka\equiv kb\pmod{m} \iff a\equiv b\pmod{m/{\gcd(m, k)}}$ という表現もできそうです。

関連しそうな話

$\gcd(m, m') = 1$ のとき、$x \equiv a \pmod{m}$ から $x \equiv a' \pmod{m'}$ のような式を導くのは不可能そうです。 なぜなら、CRT から、$x \equiv i \pmod{m'}$ を満たす $x$ が任意の $i$ に対して構成できてしまうためです。

ところで、$[0]_{(4)} = [0]_{(12)} \cup[4]_{(12)} \cup[8]_{(12)}$ のことを考えていたところ、以下のことに気づきました。

応用編

頻出テクニックの紹介です。

階乗・組合せ

$[0!] = [1]$ と $[i!] = [(i-1)!] \times [i]$ を用いて前計算しておくことで、$i! \bmod m$ を $O(1)$ 時間で答えられます。 また、$[n!]$ から $[n!]^{-1}$ を求めれば、$[(i-1)!]^{-1} = [i!]^{-1} \times [i]$ によって降順に求めていけるため、階乗の逆元の列挙も $O(n+\log(m))$ 時間で行えます。$[n!]$ の逆元を求めているので、$\gcd(n!, m) = 1$ が必要そうです。

さらに、$[i]^{-1} = [(i-1)!]\times[i!]^{-1}$ を用いて逆元の列挙も(前述の前計算を含めて)$O(n+\log(m))$ 時間でできます。 実際には、後述の方法で逆元の列挙を $O(n)$ 時間で行えるため、階乗や階乗の逆元の列挙も $O(n)$ 時間になります。

また、数え上げ頻出の組合せについても、 $$ \binom{n}{k} = {}_n\mathrm{C}_k = \frac{n!}{k!\,(n-k)!} $$ が成り立つので、階乗と階乗の逆元を前計算しておくことで、クエリあたり $O(1)$ 時間で求められます。

上記の方法は有名ですが、$\gcd(n!, m) = 1$ を仮定していました。$10^6! \bmod 998244353$ などを求める際にはこれで困らないですが、そうでない状況も考えておくといいかもしれません。

たとえば、$\binom{19}{5} \bmod{7}$ を求めてみましょう。$\binom{19}{5} = 11628 \equiv 1 \pmod{7}$ です。 一方、前計算した値では $19! \equiv 0 \pmod{7}$ と $5!^{-1} \equiv 1$ で、$14! \equiv 0 \pmod{7}$ で $14!^{-1}$ は存在しません。これは困りました。 あるいは、$\binom{13}{6} \bmod{8} = 1716 \bmod{8} = 4$ などはどうでしょうか。こちらもつらそうです。

法が合成数だったり法が十分に大きくなかったりする状況にあまり出くわさないかもしれないですが、上記の方法がいつでも使えるわけではないことには注意しておいた方がよさそうです。合成数であっても $\gcd(n!, m) = 1$ であればよいので、たとえば $\binom{40}{12} \bmod (101\cdot 103\cdot 107\cdot 109) = 5586853480 \bmod 121330189 = 5664786$ などは上記の方法で求められそうです。

調べたところいろいろな記事があったので、詳しい解説はそちらに任せます。中国剰余定理が重要になってくるようです。 法がうれしくないものの場合でも $n$ の上限が小さい場合は Pascal の三角形を利用する方法もありますね。

また、素数 $p$ に対して $(p-1)! \equiv -1 \pmod{p}$ が成り立つという話もあり、Wilson の定理として知られています。逆も成り立つようです。

有理数の表現

確率や期待値などを求める問題で、「答えは既約分数 $p/q$ ($\gcd(q, m) = 1$) の形で書くことができるため、$pq^{-1}\bmod m$ を出力せよ」という形式のものは頻出です。

ABC 298 E から引用すると、

この問題で求める確率は必ず有理数になることが証明できます。 また、この問題の制約下では、求める確率を既約分数 $\frac{y}{x}$ で表したときに $x$ が $998244353$ で割り切れないことが保証されます。 このとき $xz\equiv y \pmod{998244353}$ を満たすような $0$ 以上 $998244352$ 以下の整数 $z$ が一意に定まります。この $z$ を答えてください。

とのことです。$998244353$ が素数であることを明示しないと、「$x$ が $998244353$ で割り切れない」だけでは不十分な気もしますが、各々素数判定くらいやってくれ(あるいは有名事実なのでわかるでしょう)という感じでしょうか*10

逆元の存在性の議論を無視している人たちは「$\gcd(q, m) = 1$ の条件はなんで必要なんですか」「小数は使わずどんな問題でも $pq^{-1}\bmod m$ でやらせたらいいんじゃないですか」と思いがちな気がしますが、上記で除算に関する議論をし終えたので、読者の皆さんはもう大丈夫でしょう。

また、$x = pq^{-1} \bmod m$ のみを与えられた状態で $(p, q)$ を求める方法もあります。 適当にやれば無限個出てきてしまうというのはそうなんですが、$p$, $q$ の範囲を適切に絞ることで一意になることが知られています。

rsk0315.hatenablog.com

Re: 除算

考察の末、答えが $x/y\in\Z$ になるとわかり、これを $m$ で割ったあまりを求めたくなったとします。 たとえば、等比数列の和 $(a^k-1)/(a-1)$ のような例がわかりやすいでしょうか。

しかし、$\gcd(y, m) = 1$ とは限らず困ってしまいました。どうしましょう。

等比数列の和であれば、そういうモノイド積 $(s_l, t_l)\circ (s_r, t_r) \triangleq (s_l\cdot t_r+s_r, t_l\cdot t_r)$ によって $$ \underbrace{(1, a)\circ\dots\circ(1, a)}_{k\text{ many}} = \left(\sum_{i=0}^{k-1}a^i, a^k\right) $$ が求まります(ダブリングで高速化。単位元は $(0, 1)$)が、ここではそれ以外の方法を考えます。

先ほど $ka\equiv kb\pmod{km} \iff a\equiv b\pmod{m}$ を示したので、まず $z = x \bmod (ym)$ ($0\le z\lt ym)$ を求め、(仮定から $z$ は $y$ の倍数なので)$z/y = (x/y \bmod m)$ となります(ちょっと行間があるように感じるかもしれません)。計算途中の値が $(ym-1)^2$ 程度まで大きくなりうるので、状況によっては 64-bit 整数でもオーバーフローに気をつけましょう。 そもそも $ym$ がめちゃくちゃ大きくなるときは大変かもしれません。

$x/(y+km)\bmod m$ ($0\le y\lt m$, $\gcd(y, m) \gt 1$, $k\gg m$) を求めたいときはどうしましょう。 そもそも $x/(y+km)$ が整数だとしても $x/y$ が整数になるとは限らないですし、適当な方針では無理そうです($x = y\cdot (y+km)\cdot x'$ の場合は?)。

Re: Re: 除算

$x/y \bmod{m}$ の計算をするためには $O(\log(y))$ 時間かかってしまうのでした。

ここで、分子と分母を分けてペアで持つことにし、これに対する演算を考えます。 $\gcd(y_0, m) = \gcd(y_1, m) = 1$ を前提とします。除算については $\gcd(x_1, m) = 1$ とします。 $$ \begin{aligned} ([x_0], [y_0]) \pm ([x_1], [y_1]) &\triangleq ([x_0]\times[y_1]\pm[x_1]\times[y_0], [y_0]\times[y_1]) \\ ([x_0], [y_0]) \times ([x_1], [y_1]) &\triangleq ([x_0]\times[x_1], [y_0]\times[y_1]) \\ ([x_0], [y_0]) \div ([x_1], [y_1]) &\triangleq ([x_0]\times[y_1], [y_0]\times[x_1]) \end{aligned} $$ この形式で管理することで、最終的な答え $(x, y)$ に対して一度だけ $x/y\bmod{m}$ を計算すればよくなりそうです。

互いに素の制約に関して

ある $q$ に対して $zy = x + qm$ なる整数 $z$ が求める答えです。$g = \gcd(y, m)$ とし $y = gy'$, $m = gm'$ とすると $x = g(zy' - qm)$ なので $x$ も $g$ の倍数です。$x = gx'$ とします。

そこで、$z' \equiv x'/y'\pmod{m'}$ がわかるわけですが、$m'\lt m$ のときは実際の $(z\bmod m)$ がわかりません。たとえば、$(gz' \bmod m)$ や $(g(z'+1)\bmod m)$ などが解としてありえます。

具体例として、$18/2 \bmod 10$ を求めるとします。$(18 \bmod 10, 2\bmod 10) = (8, 2)$ を得たとき、答えが $[8/2]_{(10/2)} = [4]_{(5)}$ に属することはわかりますが、$[4]_{(5)} = [4]_{(10)} \cup[9]_{(10)}$ なので、解が定まらなさそうです。

別の例を出すと、$\frac{2}{5}\times\frac{3}{2}\times\frac{5}{3}\times\frac{x}{1} = \frac{30x}{30}$ を $30$ を法として求めたいとき、得られているペアは $(0, 0)$ ですから、困ってしまいそうです。

法を $g^2m'$ として計算し直す方法も考えられそうです? しかし、DP などの計算過程で、分子・分母が $m^{10^{10^{100}}}$ などの倍数に相当する状態になったりしていたらもう不可能そうです。

素数冪に分解して CRT で復元することを考えればゼロ除算にのみ気をつければいいですが、これだと結局 CRT の計算量の方が重くなってしまいそうです。

たとえば、素数の法 $p\gt n$ に対して $a(x) = \sum_{i=0}^{n-1} a_i x^i$ ($a_0\ne 0$, $0\le a_i\lt p$ ($0\le i\lt n$)) が与えられて、次の条件を満たす $b(x) = \sum_{i=0}^{n-1} b_i x^i$ を求めよという問題を考えます。

  • $c(x) = a(x)\cdot b(x) = \sum_{i=0}^{2n-2} c_i x^i$ とするとき、以下が成り立つ
    • $c_0 \equiv 1\pmod{p}$
    • $c_1 \equiv \dots \equiv c_{n-1} \equiv 0\pmod{p}$
    • ($n$ 次以上の係数については不問)

これは、 - $a_0b_0 = c_0$ - $c_i = \sum_{j=0}^i a_j b_{i-j}$ ($0\le i\lt n$)

となっていることから、$b_0 = c_0 a_0^{-1}$($a_0\equiv 0$ に注意)と $b_i = -a_0^{-1} \sum_{j=1} a_j b_{i-j}$ によって、$\Theta(n^2+n\log(p))$ 時間で求められそうです。

しかし、計算過程を分子と分母のペアで持ちつつ計算することで、逆元の計算を遅延させることができます。

続き

最終的に $b_i = (x_i, y_i)$ ($0\le i\lt n$) が得られますが、そこからまず $\prod_{i=0}^{n-1} y_i$ を求めます。 その後、前述の階乗の逆元を $O(n+\log(m))$ で求める方法と同様にして、$(\prod_{i=0}^{n-1} y_i)^{-1}$ を求めてから、$i$ の降順に $y_i^{-1}$ が $O(1)$ 時間で求められます。 これにより、全体で $\Theta(n^2+\log(p))$ 時間にできます。

🍑🍝🦃

実は、奇素数 $p = q\cdot 2^e+1$ に対して $n\lt 2^e$ であれば、$O(n\log(n)+\log(p))$ 時間で上記の $b(x)$ を求めることができます。(安直にやると $n\lt 2^{e-1}$ が必要そうですが、工夫すると $n\lt 2^e$ でできそうな気がします。嘘だったら訂正します。) 長さ $n$ の多項式の積の計算に必要な計算量を $M(n)$ として $O(M(n))$ 時間になる気がします*11

そうでない素数に対しても、いい感じの素数を寄せ合わせて CRT を用いることなどで高速に積を計算できそうですが、そうした素数を探す部分が高速にできないような気がしています。

今回の話では詳しく触れないものの、Fourierフーリエ 変換を $\Z_p$ で高速に行うアルゴリズムが存在し、それを用いていろいろなことをするとできます。 Newtonニュートン 法などを用います。

有名なものとしては、$998244353 = 119\cdot 2^{23}+1$ などがあります。998244353もも パスタ ターキー*12

$1000000007 = 500000003\cdot 2^1+1$ 笑

発展編

ここは発展編であり、詳しい説明はしません。紹介だけするので興味を持った方は各自調べてください。

逆元の列挙

与えられた $n$ と法 $m$ に対して、$1\le i\le n\le m$ なる各 $i$ ($\gcd(m, i) = 1$) に対して $[i]^{-1}_{(m)}$ を求めます(存在しない $i$ についてはその旨を報告)。$O(n)$ 時間です。

よく知られた方法は

dp[1] = 1;
for i in 2..n {
    dp[i] = m - (m / i) * dp[m % i] % m;
}

のようなものですが、法が合成数のときはうまくいくとは限らないことに注意してください。たとえば、$[3]^{-1}_{(8)}$ を求める際に $[8\bmod 3]^{-1}_{(8)} = [2]^{-1}_{(8)}$ が必要になりますが、これは存在せず、うまくいきません。

www.youtube.com

$m$素数でない場合は線形篩を使う方法が必要になります。

37zigen.com

平方根

与えられた $a$ に対して $[x]^2 = [a]$ なる $x$ を求める。たとえば、法を $12$ とするとき、$[3]^2 = [9]^2 = [9]$ となるが、$[x]^2 = [2]$ となる $x$ は存在しない。

$[1]^2 = [m-1]^2$ であることと $\Z_m$ の要素数が有限であることから、$m\gt 2$ では少なくとも一つの要素については平方根が存在しないことがすぐにわかる。 なお、平方根は高々二つとは限らず、$12$ を法とするときは $[1]^2 = [5]^2 = [7]^2 = [11]^2 = [1]$ となる。

関連キーワード:Tonelli–Shanks algorithm, Cipolla’s algorithm, Legendre symbol

離散対数 (discrete logarithm)

$\gdef\dlog{\operatorname{dlog}}$

与えられた $a$, $b$ と法 $m$ に対して、$[a]^x = [b]$ なる $x$ の最小値を $\dlog_{[a]}([b])$ と書く。 これを求める。ただし、存在しない場合はその旨を報告する。

関連キーワード:baby-step giant-step

テトレーション (tetration)

${}^k a = \underbrace{a^{(a^{(\cdots ^a)})}}_{k\text{ many}}$ と定義する。このとき、${}^k a \bmod m$ を求める。

なお、$a_1^{(a_2^{(\cdots ^{a_k})})} \bmod m$ を求める問題が解ければ十分である。

また、十分大きい $k$ に対して ${}^k a \equiv {}^{k+1} a$ となるので、この値を ${}^{\infty} a \bmod m$ と書くことにする。これも求めよ。 なお、この $k$ は $\phi^{\star}(m) \in O(\log(m))$ で抑えられる*13

note: $\gcd(a, m) \gt 1$ の場合は平気?

hack: ${}^3 2 = 16$ だが、${}^3 2 \bmod 32 = 0$ とするコードを書く人が多いので注意。

関連キーワード:Japanese Exponentation*14

trap.jp

原始根 (primitive root)

整数 $m$ に対して、$m$ 以下の正整数のうち $m$ と互いに素な整数からなる($\phi(m)$ 要素の)集合を $\Phi(m)$ とおきます。 $$ \Phi(m) \triangleq \{a\mid 1\le a\le m, \gcd(a, m) = 1\} = \{a_1, a_2, \dots, a_{\phi(m)}\}. $$ 特に、$m$素数のときは $\Phi(m) = \{1, 2, \dots, m-1\}$ となります。

$1\le a\lt m$ であって、$\{a^i\mid 1\le i\le m-1\} = \Phi(m)$ となるものを、法 $m$ における原始根と言います。 原始根は、存在すれば $\phi(\phi(m))$ 個あるらしいです。

$m$素数の場合の求め方などについては ACL math の解説記事に書きました。

rsk0315.hatenablog.com

その他トピック

  • 定数除算最適化
    • 除算命令(あまりも含む)は重い。
    • コンパイル時定数であれば最適化される場合が多いはず。
      • グローバルに置くのであれば const と書く必要があるかも。
    • Barrett reduction や Montgomery multiplication など。
  • $n! \bmod p$ を $O(\sqrt{p}\log(p))$ 時間で求める
    • むずかしい。
  • Library Checker の Math
    • 興味深いものがあるかも。むずかしい。
  • 思い出したら追記予定

おまけ

おきもち

「この性質は法が素数のときしか使えない」「この性質は法と互いに素じゃないと使えない」のようなものを丸暗記して使おうとするのは、たとえばデータ構造の各操作の計算量を(表とかを見て)丸暗記しようとするのに似ていそうです。 それよりも、証明を追ったり実装を理解したりする経験を通して直観を得るのがよさそうです。データ構造の計算量に関しても、内部を理解していればいちいち計算量を丸暗記する必要はないでしょう。

おわり

がんばってたくさん書きました。えびちゃんも mod についてふんわりとしか理解していない競プロ er のみなさんの一員であることを自覚しながら記事を書きました。苦しかったです。


付録

証明をしました。

証明

Claim. $b\gt 0$ に対して $a - \floor{a/b}\cdot b = (a\bmod b)$ が成り立つ。

Proof. $0\le a - \floor{a/b}\cdot b \lt b$ を示す。 床関数の性質から、 $$ \begin{aligned} \floor{a/b} &\le a/b \lt \floor{a/b+1} \\ \floor{a/b}\cdot b & \le a \lt \floor{a/b+1}\cdot b \\ 0 & \le a - \floor{a/b}\cdot b \lt b \end{aligned} $$ が成り立つ。 よって、$q = \floor{a/b}$ とすると $a = qb + (a-qb)$ かつ $0\le a-qb \lt b$ であり、$a-qb$ は $a\bmod b$ に他ならない。

Claim. $0\lt x \le y$ に対して $(y\bmod x) \lt y/2$ が成り立つ。

Proof. $y = qx+r$ とおく ($0\le r\lt x$)。このとき $q\ge 1$ であることに注意する。 $$ \begin{aligned} y &= qx+r \\ &\ge x + r \\ &\gt r + r \\ &= 2r. \end{aligned} $$ よって $y \gt 2r$ であり、$r\lt y/2$ が従う。

Claim. $\gcd(a, m) = 1$ とする。 $a^{\phi(m)} \equiv 1\pmod{m}$ が成り立つ(Euler の定理)。特に、$m$素数のときは $\phi(m) = m-1$ であり $a^{m-1} \equiv 1 \pmod{m}$ となる(Fermat の小定理)。

Proof. $m$ と互いに素な $m$ 以下の正整数は $\phi(m)$ 個あり、これらからなる集合を $B = \{b_1, b_2, \dots, b_{\phi(m)}\}$ とおく。 また、$B' = \{ab\bmod m\mid b\in B\} = \{(ab_1 \bmod m), \dots, (ab_{\phi(m)}\bmod m)\}$ とする。

前提より、$\gcd(a, m) = 1$ なので、各 $b\in B$ に対して $\gcd(ab, m) = 1$ である。 また、$a^{-1}$ が存在し、$B = \{a^{-1}b'\bmod m\mid b'\in B'\}$ と書けるので、$|B| = |B'|$ とわかる*15。 すなわち、$B'$ は $m$ と互いに素な $m$ 以下の正整数 $\phi(m)$ 個からなる集合であり、$B = B'$ である。

$\prod_{b\in B} = P$ とすると、$\prod_{b'\in B'} \equiv a^{\phi(m)}P$ となる。$B = B'$ よりこれらの積は等しく、$\gcd(P, m) = 1$ より $P^{-1}$ が存在するため、$a^{\phi(m)} \equiv 1 \pmod{m}$ が従う。

note: $\gcd(a, m) \ne 1$ の場合、$B = B'$ とならないことから破綻する。たとえば、$\phi(12) = 4$ のとき、互いに素なケースでは $5^4 = 625 \equiv 1 \pmod{12}$ のようになるが、互いに素でないケースでは $3^4 = 81 \equiv 9$ や $6^4 = 1296 \equiv 0$ となる。

note: Carmichael の定理というのもある。適宜調べてみよ。

Claim. 素数 $p$ に対し、$(p-1)!\equiv -1\pmod{p}$ である。

Proof. $f(x) = x^{p-1} - 1$ と $g(x) = (x-1)(x-2)\dots(x-(p-1))$ を考える。 Fermat の小定理より、$f(1) \equiv f(2) \equiv \dots \equiv f(p-1) \equiv 0$ となる。 $f(x)$ を因数分解することで(最高次の係数が $1$ であることと次数が $p-1$ であることから)$f(x) = f(x-1)(x-2)\dots(x-(p-1))$ を得る。 よって $f(x) = g(x)$ とわかる。

以上より、$g(0) = (-1)^{p-1} (p-1)!$ かつ $f(0) = 0^{p-1}-1$ より、($-1\neq 0$ から Fermat の小定理が使えて $(-1)^{p-1}=1$ なので)$(p-1)!\equiv -1\pmod{p}$ となる。

note: $f$ の次数個の根を見つけて因数分解して比較する部分は、$p$ が素数であり $\Z_p = \mathbb{F}_p$ が体になることが必要な気がします ($ab = 0 \implies a=0\vee b=0$)。 たとえば、法を $12$ として $f(x) = x^2-1$ を考えると $f(1) \equiv f(5) \equiv 0$ となりますが、$g(x) = (x-1)(x-5)$ とすると $f(0) \equiv -1 \not\equiv 5 \equiv g(0)$ となり破綻しています。

素数冪 $p^k$ に対して、$f(x) = x^{\phi(p^k)} = x^{p^{k-1}\cdot(p-1)}$ を考えることで、$\prod_{\gcd(i, p)=1} i$ を求めようとしたのですが、うまくいかなさそうでした。逆元が一意に定まることからペアを作って $1$ にするタイプの証明をするしかないでしょうか。 とはいえ Euler の定理の証明のところで逆元を使ったので、実質同じようなことをしようとしている気もします。

*1:ところで法が $1$ というケースもあります。$(0^0 \bmod 1) = 1$ としてしまう嘘ライブラリを持っている人も多いのではないでしょうか。

*2:各々思いついたものを入れてください。アルゴリズムの正当性の証明とかでしょうか。

*3:えびちゃんは $(\operatorname{mod}\, m)$ と書く代わりに $a\stackrel{m}{\equiv} b$ と書いて済ますこともありましたが、あまり一般的ではない気もします。

*4:ある集合 $S$ 上の関係 $\sim$ というのは、$x, y\in S$ に対して「$x\sim y$ ですよ」または「$x\sim y$ ではないですよ」というのが定まる演算のことだと思ってよさそうです。ここでは $(\Z, \equiv)$ の話をしていますが、$(\R, \le)$ なども関係です(同値関係ではない)。

*5:集合 $\{x\mid x\sim a\}$ の中で $a$ が普遍的に特別ですよーという意味ではなく、自分が代表としてそれを選んだよーという程度の意味です。別の $a'\in\{x\mid x\sim a\}$ を持ってきて $\{x\mid x\sim a'\}$ を作ったら、その文脈では $a'$ が代表元となります。

*6:半分サイズの問題を実際に二回解いてしまうと計算量は落ちないので注意しましょう。

*7:ですが、なんちゃらの定理やらなんちゃら法などを盲目的に使って求めていた人は、このようにして求めることが可能であることを認識していなかったのでは? それはあまり好ましくないのでは? とも思います。

*8:「$0$ は除いてもよいのでは?」「$m=1$ のケースにも一応気を払いました」

*9:それが妥当かの議論は置いておきます。

*10:ところで、ユーザがクリックしないと見えない場所に「答えは何々を満たすことが証明できます」というのが書いてあるのは、ちょっと不親切に感じられる気もします。わかる人は出力形式をそう指示されただけでわかる気もしますが。

*11:ここで、積のための計算量を $O(M(n) )$ ではなく $M(n)$ の書いていることに違和感を持つ人もいるのかもしれません。計算量は何でも $O(\bullet)$ のような形で書く必要があるというのは誤解な気がします。

*12:文脈は ARC 153 A で、意味はありません。

*13:一般に、$f$ を $n$ に繰り返し適用して $1$ になる回数を $f^{\star}(n)$ と書く。$\underbrace{f(f(\dots f(n)\dots) )}_{f^{\star}(n)\text{ applications}} = 1$。

*14:正しくは Exponentiation だが問題名に typo があるので注意。

*15:重複する要素ができてしまっていないことに気をつけている。