0. はじめに
2020年9月7日にAtCoder公式のアルゴリズム集 AtCoder Library (ACL)が公開されました。
私はACLに収録されているアルゴリズムのほとんどが初見だったのでいい機会だと思い、アルゴリズムの勉強からPythonでの実装までを行いました。
この記事では Modint をみていきます。
対象としている読者
- Modintってなに？という方。
- ACLのコードを見てみたけど何をしているのかわからない方。
- C++はわからないのでPythonで読み進めたい方。
参考にしたもの
Pythonのドキュメントです。
https://docs.python.org/ja/3/reference/datamodel.html
1. Modintとは
Modintとは自動であまりをとる構造体です。
1.1. Modintを使うとこんなに綺麗
例として $(a \times b + c - d) / e$ を $m$ で割ったあまりを求める場合を見てみます。
通常の場合
一般的にはオーバーフローや負数の剰余に注意して
(((((a * b) % m + c) % m - d + m) % m) * modular_inverse(e, m) % m
# modular_inverse(e, m)はmを法としたeの逆元
などのように書きます。
Modintを使う場合
aをModint型にキャストしておくと
a = Modint(a, mod=m)
(a * b + c - d) / e
と自然な形で書けます。
1.2. メリット・デメリット
Modintを使うメリットとしては
- コードの可読性が上がる
- オーバーフローや負数の剰余などに気を遣う必要がなくなる
などが考えられます。
一方デメリットとしては
- 言語によっては非常に遅くなる
- Modint自体にバグがあると見つけるのは非常に困難
などが考えられます。
1.3. 完成までに必要な知識
コンテストなどで使えるModintを作るためには
- 「Modの世界での演算」に関する数学の知識
- 最低限機能するModintを作るプログラミングの知識
- 実用的なModintを作る言語の知識
が必要です。
３番目の、実用的なModintを作る知識というのは
- いかにして高速化するか（低速化を防ぐか）
- (任意のmodに対応するか等)どの程度自由度を持たせるか
などのもので、これらは言語に強く依存するのと同時に個々人のこだわりが出るところだと思うので、本記事では扱いません。
本記事では「最低限機能するModintを作る」ことを目標とし、そのために必要な１番目、２番目の知識について説明します。
そして最後にPythonでの実装例を載せます。
2. Modの世界での演算
Modの世界での演算方法についてみていきます。
2.1. 基本事項の確認
整数 $a, b$ を自然数 $m$ で割ったあまりが等しいとき、「$m$ を法として $a$ と $b$ は合同である」といい、
a \equiv b \pmod{m}
と書きます。
例えば、3を法とするとき
\begin{aligned}
4 &= 1 \cdot 3 + 1\\
13 &= 4 \cdot 3 + 1
\end{aligned}
なので
4 \equiv 13 \pmod{3}
です。
また、本記事では「$a$ を $m$ で割ったあまり」を
a \% m
と書くことにします。つまり
\begin{aligned}
5 \% 3 = 2\\
13 \% 3 = 1
\end{aligned}
となります。
より具体的には、負の無限大への丸め込みである床関数($\lfloor x \rfloor$)を用いて
a \% m := a - m \left\lfloor \frac{a}{m} \right\rfloor
と定めます。これにより、整数 $a$ が負の場合でも自然数 $m$ で割った余りは常に $0 \leq a%m < m$ となります。
例）
\begin{aligned}
-4 \% 3 &= (-4) - 3 \left\lfloor \frac{-4}{3} \right\rfloor\\[2ex]
&= (-4) - 3(-2)\\
&= 2
\end{aligned}
2.2. 加法
加法（足し算）について以下の関係が成り立つことを示します。
(a + b) \% m = ((a \% m) + (b \% m)) \% m
(証明)
あまりをとる操作"%"の定義に沿って計算します。
(左辺) = (a + b) - m \left\lfloor \frac{a + b}{m} \right\rfloor
一方、
\begin{aligned}
(右辺) &= \left(a - m \left\lfloor\frac{a}{m} \right\rfloor + b - m \left\lfloor \frac{b}{m} \right\rfloor\right) \% m\\[3ex]
&= \left\{a + b - m \left(\left\lfloor \frac{a}{m} \right\rfloor + \left\lfloor \frac{b}{m} \right\rfloor\right)\right\} \% m\\[3ex]
&= a + b - m \left(\left\lfloor \frac{a}{m} \right\rfloor + \left\lfloor \frac{b}{m} \right\rfloor\right) - m \left\lfloor \frac{a + b - m \left(\left\lfloor \frac{a}{m} \right\rfloor + \left\lfloor \frac{b}{m} \right\rfloor\right)}{m}\right\rfloor
\end{aligned}
ここで、
\frac{m \left(\left\lfloor \frac{a}{m} \right\rfloor + \left\lfloor \frac{b}{m} \right\rfloor\right)}{m}
は整数なので
\begin{aligned}
(右辺) &= a + b - m \left(\left\lfloor \frac{a}{m} \right\rfloor + \left\lfloor \frac{b}{m} \right\rfloor\right) \\[3ex]
&\;\;\;\;- m \left\lfloor \frac{a + b}{m} \right\rfloor + m \left(\left\lfloor \frac{a}{m} \right\rfloor + \left\lfloor \frac{b}{m} \right\rfloor\right)\\[3ex]
&= (a + b) - m \left\lfloor \frac{a + b}{m} \right\rfloor
\end{aligned}
となります。したがって
(a + b) \% m = ((a \% m) + (b \% m)) \% m
が成り立ちます。
(証明終)
この関係式はどういう意味でしょうか。
左辺は
- $a+b$を計算する
- その結果を $m$ で割ったあまりにする
という順番です。これは、「$a+b$ を $m$ で割ったあまり」を求める自然な計算です。
一方、右辺は
- $a, b$ をそれぞれ $m$ で割ったあまりにする
- それらを足し合わせる
- その結果を $m$ で割ったあまりにする
となっています。
右辺はわざわざ余計な手間をかけているように見えますが、実用的な意味はあります。いま、$a% m$ や $b% m$ は $m$ 未満の数なので右辺の計算過程ではせいぜい $2m$ 程度にしかならないのです。つまり、$2m$ が「扱える整数の最大値」を超えなければオーバーフローせずに計算できます。
具体例として符号付１バイト整数型(最大値127)を使って計算する場合をみてみます。
$a = 99, b = 88, m = 55$ のとき、左辺のように足し算をしてからあまりをとると
a + b = 187 > 127
なのでオーバーフローしてしまいます。
しかし右辺のように計算することで
\begin{aligned}
(a \% m + b \% m) \% m &= (44 + 33) \% 55\\
&= 77 \% 55\\
&= 22
\end{aligned}
と常に符号付１バイト整数型の範囲内で計算できます。
コンテストでよく現れる $m = 1000000007$ は２倍しても 符号付４バイト整数の最大値 $2147483647$ を超えない数になっています。
また、
\begin{aligned}
(a + b) \% m = ((a \% m) + b) \% m\\
(a + b) \% m = (a + (b \% m)) \% m
\end{aligned}
なども先ほど示した証明と同様の手順で示すことができます。これは加法の場合は最後にあまりをとりさえすれば途中で好きな時にあまりをとって良いということです。
まとめると、Modの世界での加法について以下のことがわかりました。
- あまりをとる操作はいつ何度でも行って良い
- 毎度あまりをとって計算することで計算過程での最大値を $2m$ に抑え、オーバーフローが防げる。
2.3. 減法
減法（引き算）についても加法と同様に以下の関係が成り立ちます。
(a - b) \% m = ((a \% m) - (b \% m)) \% m
証明は加法と同様なので割愛します。
加法と同様の考察からModの世界での減法については以下のことがわかります。
- あまりをとる操作はいつ何度でも行って良い
- 毎度あまりをとって計算することで計算過程での絶対値の最大値を $m$ に抑え、オーバーフローが防げる。
2.4. 乗法
乗法（掛け算）については以下の関係が成り立ちます。
(ab) \% m = ((a \% m) \cdot (b \% m))\% m
(証明)
"％"の定義から
(左辺) = ab - m\left\lfloor\frac{ab}{m}\right\rfloor
一方、
\begin{aligned}
(右辺) &= \left\{\left(a - m\left\lfloor\frac{a}{m}\right\rfloor\right)\cdot\left(b - m\left\lfloor\frac{b}{m}\right\rfloor\right)\right\}\% m\\[3ex]
&= \left\{ab - m\left(a\left\lfloor\frac{b}{m}\right\rfloor + b\left\lfloor\frac{a}{m}\right\rfloor\right) + m^2\left\lfloor\frac{a}{m}\right\rfloor\cdot\left\lfloor\frac{b}{m}\right\rfloor\right\}\% m\\[3ex]
&= ab - m\left(a\left\lfloor\frac{b}{m}\right\rfloor + b\left\lfloor\frac{a}{m}\right\rfloor\right) + m^2\left\lfloor\frac{a}{m}\right\rfloor\cdot\left\lfloor\frac{b}{m}\right\rfloor\\[3ex]
&\;\;\;\;-m\left\lfloor\frac{ab - m\left(a\left\lfloor\frac{b}{m}\right\rfloor + b\left\lfloor\frac{a}{m}\right\rfloor\right) + m^2\left\lfloor\frac{a}{m}\right\rfloor\cdot\left\lfloor\frac{b}{m}\right\rfloor}{m}\right\rfloor
\end{aligned}
ここで、
\frac{- m\left(a\left\lfloor\frac{b}{m}\right\rfloor + b\left\lfloor\frac{a}{m}\right\rfloor\right) + m^2\left\lfloor\frac{a}{m}\right\rfloor\cdot\left\lfloor\frac{b}{m}\right\rfloor}{m}
は整数なので床関数の外に出すことが出来、
\begin{aligned}
(右辺) = ab - m\left\lfloor\frac{ab}{m}\right\rfloor
\end{aligned}
となります。したがって、
(ab) \% m = ((a \% m) \cdot (b \% m))\% m
が成り立ちます。
(証明終)
乗法も加法、減法と同様の関係が成り立つことがわかりました。
また、加法の場合には計算過程の最大値は $2m$ でしたが、乗法の場合には $m^2$ となります。
例えば、$m=1000000007$ のとき右辺のように計算することで符号付８バイト整数で扱えるようになります。
以上よりModの世界での乗法について以下がわかりました。
- あまりをとる操作はいつ何度でも行って良い
- 毎度あまりをとって計算することで計算過程での絶対値の最大値を $m^2$ に抑えられる
2.5. 除法
除法（割り算）はこれまでとは少し変わります。
ここまで加法、減法、乗法はModを考えない通常の世界での演算と同様に自然に計算できました。
例えば加法の場合、通常の世界で
4 + 7 = 11
であることと同様にModの世界でも
4 + 7 \equiv 11 \pmod{3}
のようにできました。
しかし、除法の場合はどうでしょうか。
4 \div 7 \equiv ? \pmod{3}
通常の世界と同様に考えると
4 \div 7 = 0.5714\dots
なので 「$0.5714\dots$ を $3$ で割ったあまり」となってしまいよくわかりません。
実は「割り算とは何か」ということを考えると Modの世界での除法も通常の世界と同様に考えられることがわかります。
2.5.1 通常の世界での割り算
いま、$a$ を整数、$b$ を $0$ でない整数とします。すると $a$ を $b$ で割ることができます。
割り算は分数と掛け算に直すことができ
a \div b = a \times \frac{1}{b}
でした。
つまり言葉で書くと
「$b$ で割る」は「$\frac{1}{b}$ を掛ける」と等しい
ということです。では、$\frac{1}{b}$ とは何でしょうか。もちろん「$1$ を $b$ で割った数」ではあるのですがここでは次のように見ます。
$\frac{1}{b}$ は「$b$ に掛けると $1$ になる数」
このような数を(乗法についての) $b$ の逆数、$b$ の逆元などと呼びます。本記事では「乗法についての」は省略し、単に「 $b$ の逆元」と呼ぶことにし、 $b^{-1}$ と書くことにします。数式で書くなら $b^{-1}$ は
b\cdot b^{-1} = 1
を満たす数ということです。
以上より、
「$b$ で割る」は「$b$ の逆元($b^{-1}$)を掛ける」と等しい
ということがわかりました。
まとめると
- 割る = 逆元を掛ける
- 通常の世界では $0$ でない数 $x$ の逆元は分数 $\frac{1}{x}$
となります。
2.5.2. Modの世界での割り算
Modの世界でも通常の世界と同じように考えます。
つまり、「$m$ を法とした世界で $b$ で割る」というのは「$m$ を法とした世界における $b$ の逆元($b^{-1}$)を掛ける」という意味になります。そしてこの逆元は
b \cdot b^{-1} \equiv 1 \pmod{m}
を満たす数です。このよう逆元を $m$ を法とした $b$ の逆元と言います。
具体例を見てみます。
$a = 4, b = 7, m = 3$ のとき
a \div b \pmod{m}
を考えます。
7 \cdot 4 \equiv 28 \equiv 1 \pmod{3}
なので $3$ を法とした $7$ の逆元は $4$ です。
したがって
\begin{aligned}
a \div b &\equiv a \cdot b^{-1}\\
&\equiv 4 \cdot 4\\
&\equiv 16 \\
&\equiv 1 \pmod{3}
\end{aligned}
となります。
2.5.3. 逆元の存在条件と求め方
通常の世界で $0$ の逆元が存在しないのと同様にModの世界でも逆元が常に存在するとは限りません。
例えば、$3$ を法とした $6$ の逆元 $x$ は
6x \equiv 1 \pmod{3}
を満たす数ですが $6x$ は $3$ を法として常に $0$ と合同なのでこの式を満たす $x$ は存在しません。
$m$ を法とした $b$ の逆元が存在するのは $m$ と $b$ が互いに素のときのみになります。
そして逆元を求める方法として
- フェルマーの小定理
- 拡張ユークリッドの互除法
の２つがあります。
フェルマーの小定理とは次のようなものです。
$p$ を素数とし、$a$ を $p$ の倍数でない整数（$a$ と $p$ は互いに素）とするときに、
a^{p-1} \equiv 1 \pmod{p}
が成り立つ。
この合同式の両辺に $p$ を法とした $a$ の逆元 $a^{-1}$ をかけると
a^{p-2} \equiv a^{-1} \pmod{p}
となります。すなわち逆元は $a^{p-2}$ で与えられます。
例えば、$p = 1000000007$ を法とした $6$ の逆元は
6^{1000000005} \pmod{p}
を計算すれば良いことになります。これは繰り返し二乗法を用いて高速に計算できます。
ただし、この方法は法が素数の場合にしか使えないことに注意が必要です。(フェルマーの小定理が素数に関する定理なので当然です。)
もっとも、$1000000007$ や $998244353$ など多くの場合は与えられる法が素数だと思います。
一方、拡張ユークリッドの互除法は一次不定方程式
ax + by = \gcd(a, b)
を解くアルゴリズムです。（$a, b$ は $0$ でない整数、$\gcd(a, b)$ は $a$ と $b$ の最大公約数）
$m$ を法とした $a$ の逆元 $x$ は
a x \equiv 1 \pmod{m}
を満たす数でした。この合同式は整数 $y$ を用いて
ax + my = 1
という一次不定方程式に書き直せます。逆元が存在するとき、$a$ と $m$ は互いに素なのでこの方程式は
ax + my = \gcd(a, m)
となります。したがってこれを拡張ユークリッドの互除法で解くことで逆元を求めることができます。
繰り返し二乗法や拡張ユークリッドの互除法についてはinternal_math編①で説明していますのでそちらをご覧ください。
3. Modintに必要なプログラミングの要素
3.1. クラス
Modintは
- 四則演算などの演算ができる
- 自動であまりをとる
などの機能を持ったオブジェクトです。このような何らかの機能や属性を持ったオブジェクトを作るために多くの言語ではクラスが導入されています。よってやるべきことはModintというクラスを作成し、そこに必要な要素を詰め込むという作業です。
これは整数型、浮動小数点型、文字列型といった様々な型に加えて新たにModint型という型を作る作業とも言えます。
3.2. 演算子のオーバーロード
Modの世界での演算をModint型に組み込んでも例えば $a+b$ や $a \times b$ をするために
a.add(b), a.mul(b)
などと書かなければならないようではModintの恩恵があまりありません。やはり
a + b, a * b
と書けるようにしたいです。そのためには演算子の意味を書き換える必要があります。これを演算子のオーバーロードと言います。（演算子のオーバーロードがユーザーに提供されていない言語もあります）
例えば、通常 "+" という演算子は数の足し算を意味しますが、文字列型においては文字列の連結という意味になるなど、演算子のオーバーロードはすでに多く使われています。
極端な例として天邪鬼型という型を次のように作ったとします。(コードはイメージで特定の言語の文法に則ったものではありません)
class Amanojaku{
def constructor(int: x){
int _x = x
}
// 演算子のオーバーロード
op{+}(lhs, rhs) := lhs._x - rhs._x
op{-}(lhs, rhs) := lhs._x + rhs._x
// 右辺の + や - は lhs._xやrhs._xがint型なことから
// int型の演算における +,- (つまり通常の足し算、引き算としての +,-)の意味になる。
}
すると、天邪鬼型にキャストされた数 $a, b$ 同士の演算は
a = Amanojaku(5)
b = Amanojaku(3)
print(a + b) --> 2
print(a - b) --> 8
のようになります。
このように演算子のオーバーロードによって演算子の意味を「通常の世界での演算」から「Modの世界での演算」に変えることでModintを作ることができます。
4. Pythonでの実装
それでは実装していきます。
まず、Modintクラスを作成します。Modintの計算ではmodを全体で共有する必要があるのでクラス変数としてmodという変数を持つことにします。クラス変数はクラス名.変数名でアクセスできます。また、この値は計算途中で変更されてはいけないのでhas_been_setというフラグを持つことにします。
値は_vというインスタンス変数に格納され、このとき負かmod以上であればmodで割ったあまりに置き換えます。
class Modint:
mod = 0
has_been_set = False
def __init__(self, v=0, m=None):
if m != None:
assert m >= 1
assert not Modint.has_been_set
Modint.mod = m
Modint.has_been_set = True
assert Modint.has_been_set
self._v = v if 0 <= v < Modint.mod else v % Modint.mod
Pythonでは演算子のオーバーロードは特殊メソッドという形で提供されます。特殊メソッドは前後にアンダースコア(_)が２つずつ付いたものでコンストラクタも特殊メソッドの１つです。
数の演算に関する演算子と対応する特殊メソッドを表にまとめると次のようになります。
| 演算子 | 特殊メソッド |
|---|---|
+ |
__add__ |
- |
__sub__ |
* |
__mul__ |
/ |
__truediv__ |
// |
__floordiv__ |
** |
__pow__ |
これらをModの世界での演算となるように定義していきます。
例えば a + bをいうコードがあったとき、a.__add__(b)が呼ばれます。なのでbの型をチェックする必要があります。bがModint型であればb._vの値を取り出して演算し、bがint型であればそのまま演算します。そしてその結果の値を新たなModint型として返します。
Modの世界で割り算は整数同士の演算の結果整数が得られるので//を使うことにしました。逆元は組み込み関数pow()で求めることができます。Python3.8以降では第２引数を-1とすれば拡張ユークリッドの互除法でmod-2とすればフェルマーの小定理を使った方法になります。それ以外のバージョンの場合はmod-2とするか拡張ユークリッドの互除法を実装してください。実装方法はinternal_math編①にあります。フェルマーの小定理を使う方法はmodが素数に限定されますので拡張ユークリッドの互除法を使うことをお勧めします。
def __add__(self, other):
if isinstance(other, Modint):
res = self._v + other._v
if res > Modint.mod: res -= Modint.mod
else:
res = self._v + other
return Modint(res)
def __sub__(self, other):
if isinstance(other, Modint):
res = self._v - other._v
if res < 0: res += Modint.mod
else:
res = self._v - other
return Modint(res)
def __mul__(self, other):
if isinstance(other, Modint):
return Modint(self._v * other._v)
else:
return Modint(self._v * other)
def __floordiv__(self, other):
if isinstance(other, Modint): other = other._v
inv = pow(other, -1, Modint.mod)
return Modint(self._v * inv)
def __pow__(self, other):
assert isinstance(other, int) and other >= 0
return Modint(pow(self._v, other, Modint.mod))
さらに代入演算子も使えると便利なのでこちらもオーバーロードします。
代入演算子は +=, *=といったもので、対応する特殊メソッドは通常の演算子の特殊メソッドに接頭辞iをつけたもの(__iadd__など)です。
代入演算はselfとotherの演算の結果新たなものが生まれるのではなくselfをotherによって変更するものなので、返り値はself自身となります。
def __iadd__(self, other):
if isinstance(other, Modint):
self._v += other._v
if self._v >= Modint.mod: self._v -= Modint.mod
else:
self._v += other
if self._v < 0 or self._v >= Modint.mod: self._v %= Modint.mod
return self
def __isub__(self, other):
if isinstance(other, Modint):
self._v -= other._v
if self._v < 0: self._v += Modint.mod
else:
self._v -= other
if self._v < 0 or self._v >= Modint.mod: self._v %= Modint.mod
return self
def __imul__(self, other):
if isinstance(other, Modint):
self._v *= other._v
else:
self._v *= other
if self._v < 0 or self._v >= Modint.mod: self._v %= Modint.mod
return self
def __ifloordiv__(self, other):
if isinstance(other, Modint): other = other._v
inv = pow(other, -1, Modint.mod)
self._v *= inv
if self._v > Modint.mod: self._v %= Modint.mod
return self
def __ipow__(self, other):
assert isinstance(other, int) and other >= 0
self._v = pow(self._v, other, Modint.mod)
return self
ここまでで演算子の左側がModint型の場合は計算できるようになりました。しかし、左側がint型、右側がModint型の場合には計算できません。
例えばint型の変数aとModint型の変数bの足し算a+bはa.__add__(b)が呼ばれますが、int型の__add__メソッドにはModint型の変数との"+"演算が定義されていないのでエラーとなります。しかし、すぐにエラーを吐くのではなく「演算子の右側の変数(Modint型)には演算が定義されているかもしれない」と考え、右側の型の中身を参照します。このとき参照されるのが各特殊メソッドに接頭辞rをつけたもの(__radd__など)です。これを定義しておくことで、(int) + (Modint)ができるようになります。
このメソッドでは非可換の演算の実装に注意する必要があります。
例えばint型の変数aとModint型の変数bの引き算a-bを考えます。このときint型には(int) - (Modint)が定義されていないので、b.__rsub__(a)が参照されます。そのため、__rsub__の定義ではself = b, other = aであることに注意してother - selfにする必要があります。
また、__rpow__は(int) ** (Modint)という状況がないと思うので、実装していません。
def __radd__(self, other):
return Modint(self._v + other)
def __rsub__(self, other):
return Modint(other - self._v)
def __rmul__(self, other):
return Modint(self._v * other)
def __rfloordiv__(self, other):
inv = pow(self._v, -1, Modint.mod)
return Modint(other * inv)
ここまでで算術演算子の定義が終了したので、次に比較演算子の定義をします。
定義するのは==, !=でそれぞれmodを法として
\begin{aligned}
a &\equiv b\\
a &\not \equiv b
\end{aligned}
であることとします。
これらはint型とModint型が混在していても機能します。
def __eq__(self, other):
if isinstance(other, Modint):
return self._v == other._v
else:
if other < 0 or other >= Modint.mod:
other %= Modint.mod
return self._v == other
def __ne__(self, other):
if isinstance(other, Modint):
return self._v != other._v
else:
if other < 0 or other >= Modint.mod:
other %= Modint.mod
return self._v != other
最後にその他の特殊メソッドを定義します。
__str__はprint()によって出力する際に呼ばれます。これがないと<Modint object at 0x10aa43310>のように表示されてしまいます。
__repr__はオブジェクトを表す正式な文字列を返します。数値型の場合には数値そのものの文字列を返すのでそれに倣い実装しました。
__int__はint()で呼ばれます。これによってint型へのキャストができます。
def __str__(self):
return str(self._v)
def __repr__(self):
return str(self._v)
def __int__(self):
return self._v
5. まとめ
全てをまとめたものを載せておきます。
class Modint:
mod = 0
has_been_set = False
def __init__(self, v=0, m=None):
if m != None:
assert m >= 1
assert not Modint.has_been_set
Modint.mod = m
Modint.has_been_set = True
assert Modint.has_been_set
self._v = v if 0 <= v < Modint.mod else v % Modint.mod
def __add__(self, other):
if isinstance(other, Modint):
res = self._v + other._v
if res > Modint.mod: res -= Modint.mod
else:
res = self._v + other
return Modint(res)
def __sub__(self, other):
if isinstance(other, Modint):
res = self._v - other._v
if res < 0: res += Modint.mod
else:
res = self._v - other
return Modint(res)
def __mul__(self, other):
if isinstance(other, Modint):
return Modint(self._v * other._v)
else:
return Modint(self._v * other)
def __floordiv__(self, other):
if isinstance(other, Modint): other = other._v
inv = pow(other, -1, Modint.mod)
return Modint(self._v * inv)
def __pow__(self, other):
assert isinstance(other, int) and other >= 0
return Modint(pow(self._v, other, Modint.mod))
def __radd__(self, other):
return Modint(self._v + other)
def __rsub__(self, other):
return Modint(other - self._v)
def __rmul__(self, other):
return Modint(self._v * other)
def __rfloordiv__(self, other):
inv = pow(self._v, -1, Modint.mod)
return Modint(other * inv)
def __iadd__(self, other):
if isinstance(other, Modint):
self._v += other._v
if self._v >= Modint.mod: self._v -= Modint.mod
else:
self._v += other
if self._v < 0 or self._v >= Modint.mod: self._v %= Modint.mod
return self
def __isub__(self, other):
if isinstance(other, Modint):
self._v -= other._v
if self._v < 0: self._v += Modint.mod
else:
self._v -= other
if self._v < 0 or self._v >= Modint.mod: self._v %= Modint.mod
return self
def __imul__(self, other):
if isinstance(other, Modint):
self._v *= other._v
else:
self._v *= other
if self._v < 0 or self._v >= Modint.mod: self._v %= Modint.mod
return self
def __ifloordiv__(self, other):
if isinstance(other, Modint): other = other._v
inv = pow(other, -1, Modint.mod)
self._v *= inv
if self._v > Modint.mod: self._v %= Modint.mod
return self
def __ipow__(self, other):
assert isinstance(other, int) and other >= 0
self._v = pow(self._v, other, Modint.mod)
return self
def __eq__(self, other):
if isinstance(other, Modint):
return self._v == other._v
else:
if other < 0 or other >= Modint.mod:
other %= Modint.mod
return self._v == other
def __ne__(self, other):
if isinstance(other, Modint):
return self._v != other._v
else:
if other < 0 or other >= Modint.mod:
other %= Modint.mod
return self._v != other
def __str__(self):
return str(self._v)
def __repr__(self):
return str(self._v)
def __int__(self):
return self._v
6. 使用例
$mod = 13$ において $a = 16, b = 7, c = 12, d = 8, e = 9$ として
(a + b - c) * d \div e
を計算してみます。
まず、modを設定する必要があります。
これは最初に計算される変数をModint型にキャストする際に設定しても良いですし、modだけを先に与えても構いません。
mod = 13
a = 16
b = 7
c = 12
d = 8
e = 9
#modの設定
Modint(m=mod)
#最初に計算する部分をModint型にキャスト
a = Modint(a)
#a = Modint(a, m=mod) #ここでmodを設定しても良い
result = (a + b - c) * d // e
print(result) # 4
print(type(result)) # <class '__main__.Modint'>
7. 速度について
実際にModintが使える問題を解いてみました。
使った問題はEducational DP Contest H-Grid 1です。この問題はDPを使った数え上げで1000000007で割ったあまりを出力します。
Modintを使うと以下のようになります。
"""
Modintを貼る(長いので省略)
"""
def main():
mod = 1000000007
h, w = map(int, input().split())
wall = [[i == '#' for i in input()] for _ in range(h)]
dp = [[0] * w for _ in range(h)]
dp[0][0] = Modint(1, mod) # Modint型にキャスト
# dp[0][0] = 1 # Modintを使わない場合
for i in range(h):
for j in range(w):
if wall[i][j]: continue
if i > 0: dp[i][j] += dp[i-1][j]
if j > 0: dp[i][j] += dp[i][j-1]
# dp[i][j] %= mod # Modintを使わない場合
print(dp[h-1][w-1])
if __name__ == '__main__':
main()
Modintを使う場合と使わない場合でそれぞれ提出した結果、実行時間は以下のようになりました。
| 使う | 使わない |
|---|---|
| $1870\rm{ms}$ | $524\rm{ms}$ |
本記事で作成したModintでは非常に遅くなります。
8. PythonでModintを実装するメリットは...?
個人的な意見ですが、PythonではModintを整備する優先度はかなり低いと思います。
なぜなら、Pythonの場合メリットが小さく、デメリットは大きいと思うからです。
１章で述べたModintのメリットの内、オーバーフローや負数の剰余については
- オーバーフロー $\rightarrow$ 多倍長整数
- 負数の剰余 $\rightarrow$ 正の数で割ったあまりは正の数になる
と解決されていますし、拡張ユークリッドの互除法を用いた逆元の計算や繰り返し二乗法を用いた高速な累乗計算も組み込み関数に実装されています。
一方、前章で述べた通り（本記事の実装では）実行速度は非常に低下します。もともと遅い言語であるPythonにとってこのデメリットは非常に大きく感じます。
もちろん、コードはスッキリするので問題の本質に集中できますし、デバッグもしやすくなる思います。また、$O(1)$ で解ける問題など、実行時間を気にしなくて良い問題などでは活躍すると思います。
(2020.12.28 追記)
多倍長整数によって解決されるのはオーバーフローによる未定義動作の心配です。
巨大な数の演算は非常に時間がかかり、TLEの原因になります。これを回避するためには計算過程においてある程度の頻度で（$\fallingdotseq$演算のたびに）あまりをとる必要があります。
9. おわりに
あまりを自動でとる構造体Modintを作成しました。非常に便利ですがPythonで実用的なものを作るのは大変そうです。
説明の間違いやバグ、アドバイス等ありましたらお知らせください。