目次
- 自己紹介
- Pythonを使った数学ライブラリーの作成をしてみよう
- 定理の説明
- どういった使い道がある？
- 実装
- 組み込み関数pow()を使用してスマートに
- 問題例
- 大きな数をコンパクトに表現してみよう
- ProjectEuler 134
実行環境
- Python __3.8__以上
バージョンが__3.7__以下だと、今回使用する組み込み関数powの仕様が違うため、__プログラムが正常に作動しません__ので注意！
自己紹介
こんにちは。sgswといいます。
今は大学生で、好きな言語はPythonです。
公式サイトはこちら
Pythonにおける、実際の活用例をアウトプットしてゆきたいと思います！
Pythonを使った数学ライブラリーの作成をしてみよう
今回のテーマは、「Pythonを使った数学での活用」です。
具体的には__中国剰余定理__とよばれる整数分野における数学の定理について、__猫でもわかる、__ように解説をしてゆきたいと思います。
厳密な説明ではないですが、厳密性よりわかりやすさを意識して説明してゆきたいと思います。
定理の説明
wikipedia先生によると、以下が中国剰余定理の説明です。
中国の剰余定理（ちゅうごくのじょうよていり、英: Chinese remainder theorem）は、中国の算術書『孫子算経』に由来する整数の剰余に関する定理である。あるいは、それを一般化した可換環論における定理でもある。中国人の剰余定理（ちゅうごくじんのじょうよていり）、孫子の定理（そんしのていり、英: Sunzi's theorem）とも呼ばれる。
『孫子算経』には、「3で割ると2余り、5で割ると3余り、7で割ると2余る数は何か」という問題とその解法が書かれている。中国の剰余定理は、この問題を他の整数についても適用できるように一般化したものである。
数学的にいうと、以下のようになるらしいです。(同じくwikipediaから)
与えられた二つの整数 m, n が互いに素ならば、任意に与えられる整数 a, b に対し、連立合同方程式
x ≡ a (mod m),
x ≡ b (mod n)
を満たす整数 x が mn を法として一意的に存在する。
もう少し拡張すると以下のようになります。(同上)
これは明らかに次のように拡張できる。すなわち：
与えられた k 個の整数 m1, m2, ..., mk がどの二つも互いに素ならば、任意に与えられる整数 a1, a2, ..., ak に対し
x ≡ a1 (mod m1),
x ≡ a2 (mod m2),
…
x ≡ ak (mod mk)
を満たす整数 x が m1m2…mk を法として一意的に存在する。
なんのこっちゃ、って感じですね。
要は、ある数を何個かの数で割った時のあまりとして持っておき情報量を失わないようにする（＝元の数を復元できる)
ということがテーマになっていそうです。
どういった使い道がある？
唐突ですが、
X = 12345
Y = 54321
Z = 6789
W = 9876
という数が与えられたとしましょう。
そして、4つの数の積XYZWを計算しなくてはいけなくなったとします。
これは、とてもめんどくさいです。なんとかして、より小さい世界の話に置き換えることはできないでしょうか...？
ここで__中国剰余定理__の出番です。
4つの数をある数たちで割ったあまりの世界で考えて、よりコンパクトに演算をしてから、元の世界に復元して結果を求めてみましょう。
今回は、
P = [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47]
という50以下の素数のリストを、あまりの世界のベースとして選択してみたいと思います。
ここで、剰余における性質
(XYZW)%p = (X%p)(Y%p)(Z%p)(W%p)
を使うと、4つの数の積をpで割ったあまりはそれぞれをpで割ったあまりの積として表すことができます。
これを利用して、
-
(Pの素因数n):(Xをnで割った時のあまり)という辞書を返す関数
-
上記の関数を利用した、XYZWをnで割った時のあまりの辞書を返す関数
を返すコードを書くと以下のようになりますね。
P = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]
X = 12345
Y = 54321
Z = 6789
W = 9876
def mod_data(M):
"""
Mについて、Pの素因数をkey,各keyで割った時のあまりを持つ辞書を返す
"""
return {p: M % p for p in P}
def mod_multiply(NUM_LIST):
"""
NUM_LISTに入っている整数の積について、Pの素因数をkey,各keyで割った時のあまりを持つ辞書を返す
"""
res = {p: 1 for p in P}
for elem in NUM_LIST:
for (p, mod) in mod_data[elem].items():
res[p] = res[p] * mod % p
return res
この結果、XYZWをPの各素因数nで割った時のあまりを得ることができました。
あとは、中国剰余定理を用いて、元々の数(=XYZW)を復元することができれば、元々のX,Y,Z,Wだけを利用した愚直な積計算をすることなく、結果を得ることができますね。
このように、中国剰余定理を利用して、大きな数での演算をする際に小さな数における演算に戻して考えることで、言語によってはオーバーフローの危険性があったり、ただでさえ時間がかかりやすい大きな数どうしの演算を避けることができます。
これは大きなメリットですね。
3.中国剰余定理の実装&問題解決例
実装
- 組み込み関数pow()を利用する
Python3.8以降ではpow()という関数に新しい引数をとることができるようになりました。
公式ドキュメントからの引用です。
For int operands base and exp, if mod is present, mod must also be of integer type and mod must be nonzero. If mod is present and exp is negative, base must be relatively prime to mod. In that case, pow(inv_base, -exp, mod) is returned, where inv_base is an inverse to base modulo mod.
>>> pow(38, -1, mod=97)
23
>>> 23 * 38 % 97 == 1
True
バージョン 3.8 で変更: For int operands, the three-argument form of pow now allows the second argument to be negative, permitting computation of modular inverses.
何ができるようになったかというと、第1引数( = a)が整数で、第3引数( = mod)が整数の時、これらが互いに共通因数を持たなければ、第2引数としてマイナスの整数の値を取ることができ、特に-1とした時inv_base = pow(a,-1,mod)という整数の値を得ることができ、この時以下の関係式が成り立ちます。
a * inv_base = 1(mod mod)
このpow関数の挙動が後々役に立ってきます。
以下では、これを用いて、先ほど述べた2整数における中国剰余定理(以下参照)
与えられた二つの整数 m, n が互いに素ならば、任意に与えられる整数 a, b に対し、連立合同方程式
x ≡ a (mod m),
x ≡ b (mod n)
を満たす整数 x が mn を法として一意的に存在する。
を解くための関数を作っていきます。以下のような課題を解決する関数を考えましょう。
固定された二つの整数 m, n（どちらも0ではない)と、任意に与えられる整数 a, b に対し、連立合同方程式
x ≡ a (mod m),
x ≡ b (mod n)
を満たす0以上の最小の整数 x を一つ求める。存在しない場合は-1を返す。
完成版だけみたい、という人は以下のコードを参照してください。
import math
def CRT(rem_pair, mod_pair):
"""
rem_pair = (r1,r2)
mod_pair = (m1,m2) とします。この時m1,m2がどちらも0ではないとします。
この時
x = m1 (mod r1)
x = m2 (mod r2)
となる整数xが存在する場合、0以上の最小の自然数を返します。
存在しない場合は-1を返します。
"""
r1, r2 = rem_pair
m1, m2 = mod_pair
assert (m1 != 0 and m2 != 0)
assert len(rem_pair) == len(mod_pair) == 2
g = math.gcd(m1, m2)
if (r2 - r1) % g != 0:
return -1
M1, M2, R = m1//g, (-m2)//g, (r2 - r1)//g
inv = pow(M1, -1, M2) * R % M2
return (m1*inv + r1) % abs(m1*m2//g)
それでは、コードの説明にうつります。
まず、2つの整数(m1,m2)があり、それぞれで割った時のあまりが(b1,b2)となるような整数tがあるものとします。この時、上のことを言い換えると、tは以下の形で2通りに表すことができます。
t = m1 * x + b1
t = m2 * y + b2
"""
x,yはtをm1,m2で割った時の商)とします。
"""
tを消去して、xとyに関する一次式にしてみましょう。
m1 * x + (-m2) * y = (b2 - b1)
唐突ですが、(m1,-m2)という二つの数の最大共通因数 g を考えます。これを求めるには、mathライブラリのgcdモジュールを使えば良いです。
すると、上の一次式の左辺は以下のように書けるはずです。
(左辺) = g * (m1//g * x + (-m2)//g * y)
"""
(m1,-m2)どちらも g を共通因数に持つため、gをくくり出すことができます。
因数分解における係数くくりだしのイメージです。
"""
これからわかることは、__左辺はgで割り切れる__ということです。したがって、右辺もgで割り切れなくてはいけません!
もし割り切れなければ、このようなxは存在できないということがわかるため、-1を返せば良いです。以下では右辺が g で割り切れる時のことを考えましょう。
両辺 g で割って係数をわかりやすい形におきなおします。
g * (m1//g * x + (-m2)//g * y) = g * (b2//g - b1//g)
つまり
m1//g * x + (-m2)//g * y = b2//g - b1//g
"""
係数を置き直す
"""
M1, M2, R = m1 // g, (-m2) // g, b2//g - b1//g
この結果もともとの一次式を以下の簡潔な形でまとめることができました。
M1 * x + M2 * y = R
ここでこのような関係式を満たす整数(x,y)のペアは先述の__pow関数を利用することで求める__ことができます！
pow関数は以下のことができるのでした。
何ができるようになったかというと、第1引数( = a)が整数で、第3引数( = mod)が整数の時、これらが互いに共通因数を持たなければ、第2引数としてマイナスの整数の値を取ることができ、特に-1とした時inv_base = pow(a,-1,mod)という整数の値を得ることができ、この時以下の関係式が成り立ちます。
a * inv_base = 1(mod mod)
今、M1,M2はあらかじめ最大共通因数 g で割っておいたので、互いに共通因数をもちません。したがってこのpow関数の引数として使える形にあります。
pow関数を用いると、一次式を満たす(x,y)のペアの一つのxは以下のように簡潔にあらわすことができます。
```python
x = pow(M1, -1, M2) * R % M2
どうしてかというと、X = pow(M1,-1,M2)は以下の条件を満たします。
M1 * X = 1 (mod M2)
つまり
M1*X + M2 * Y = 1
となるYがある
これにRをかけると...
M1 * (R * X) + M2 * (R * Y) = R
ということで、一次式の解が見つかってしまいました。
これより、条件を満たす整数tを一つ得ることができます。
inv = R * X % M2
t = m1 * inv + r1
ここまで来れたらゴールは目前です。最後はこれを元に0以上の最小のtを得ることを考えます。
二つの整数（t1,t2）が存在してどちらも条件を満たす時、明らかに、t1 - t2は(m1,m2)の両方で割り切ることができます。つまり、t1 - t2は(m1,m2)の最小公倍数で必ず割り切れることがわかります！
これは、mathライブラリのgcdモジュールを用いることで容易に計算できます。この値をLとした時、0以上で条件を満たす最小のtはこのLを用いて、以下で得られますね！
L = abs(m1 * m2 // math.gcd(m1,m2))
t = t % L
これをまとめて、完成です。
import math
def CRT(rem_pair, mod_pair):
"""
rem_pair = (r1,r2)
mod_pair = (m1,m2) とします。この時m1,m2がどちらも0ではないとします。
この時
x = m1 (mod r1)
x = m2 (mod r2)
となる整数xが存在する場合、0以上の最小の自然数を返します。
存在しない場合は-1を返します。
"""
r1, r2 = rem_pair
m1, m2 = mod_pair
assert (m1 != 0 and m2 != 0)
assert len(rem_pair) == len(mod_pair) == 2
g = math.gcd(m1, m2)
if (r2 - r1) % g != 0:
return -1
M1, M2, R = m1//g, (-m2)//g, (r2 - r1)//g
inv = pow(M1, -1, M2) * R % M2
return (m1*inv + r1) % abs(m1*m2//g)
問題例
- 大きな数をコンパクトに表現してみよう
上で作った関数を何層も積み重ねることで、一般の場合の中国剰余定理の形を作り出せます。
コードにすると以下のようになりますね。
def many_CRT(rem_pair, mod_pair):
"""
rem_pair = (r1,r2...r_t)
mod_pair = (m1,m2...r_t) とします。この時m1,m2......m_tはいずれも0ではないとします。
この時
x = m1 (mod r1)
x = m2 (mod r2)
.
.
.
x = m_t (mod r_t)
となる整数xが存在する場合、0以上の最小の自然数を返します。
存在しない場合は-1を返します。
"""
assert len(rem_pair) == len(mod_pair)
assert all(elem != 0 for elem in mod_pair)
M, LCM = CRT(rem_pair[:2], mod_pair[:2]), mod_pair[0] * \
mod_pair[1]//math.gcd(mod_pair[0], mod_pair[1])
for (r, m) in zip(rem_pair[2:], mod_pair[2:]):
if M == -1:
return -1
M = CRT((r, M), (m, LCM))
LCM = LCM * m // math.gcd(LCM, m)
return M
これを使って巨大な数を小さな数のあまりで表すことができるようになりました！
P = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]
A = 12345678910111213
B = many_CRT([A % p for p in P], P)
print(A == B) # True
これを使ってProjectEulerというサイトにある問題を解いてみましょう。
問題文の概略です。こちらから引用です。
連続する素数 p1 = 19, p2 = 23 について考える. 1219 は末尾の桁が p1 からなり p2 で割り切られる最小の数であることが確かめられる.
実際, p1 = 3, p2 = 5 を除けば, 全ての p2 > p1 なる連続する素数のペアについて, 末尾の桁が p1 からなり p2 で割り切られる数 n が存在する. S を n の最小のものであるとする.
5 <= p1<= 1000000 を満たす連続する素数のペア全てに対し ∑ S を求めよ.
この時、
末尾の桁が p1 からなり p2 で割り切られる数 n が存在する. S を n の最小のもの
とありますが、このようなnは以下の形でかける0以上の整数である、と言い換えられますね。
n = p1 (mod 10^K) ただし,Kはp1の桁数とする。
n = 0 (mod p2)
まさに中国剰余定理の形ですね！早速ライブラリーを使ってコードを作ってみましょう。
import math
def PRIME_SIEVE(M):
"""
return [0,M]の素数リストを返します
"""
isprime = [True]*(M + 1)
isprime[0] = isprime[1] = False
primes = []
for i in range(2, M + 1):
if isprime[i] == True:
primes.append(i)
for j in range(i*i, M + 1, i):
isprime[j] = False
return primes
def CRT(rem_pair, mod_pair):
"""
rem_pair = (r1,r2)
mod_pair = (m1,m2) とします。この時m1,m2がどちらも0ではないとします。
この時
x = m1 (mod r1)
x = m2 (mod r2)
となる整数xが存在する場合、0以上の最小の自然数を返します。
存在しない場合は-1を返します。
"""
r1, r2 = rem_pair
m1, m2 = mod_pair
assert (m1 != 0 and m2 != 0)
assert len(rem_pair) == len(mod_pair) == 2
g = math.gcd(m1, m2)
if (r2 - r1) % g != 0:
return -1
M1, M2, R = m1//g, (-m2)//g, (r2 - r1)//g
inv = pow(M1, -1, M2) * R % M2
return (m1*inv + r1) % abs(m1*m2//g)
def main():
ans = 0
primes = PRIME_SIEVE(1001000)
for (i, p) in enumerate(primes):
if 5 <= p <= 1000000:
p2 = primes[i + 1]
ans += CRT((p, 0), (10**len(str(p)), p2))
print("the answer is > {}".format(ans))
return
if __name__ == "__main__":
main()
中国剰余定理は偉大ですね！
4. おわりに
中国剰余定理は仰々しい名前とは裏腹に、意外と組み込み関数を用いればシンプルにかけることがお分かりいただけたのではないでしょうか？
「今までPythonを学んできたけれども、いまいち具体的な活用方法がわからない...。」
という人にきっかけを与える手助けになれたら幸いです。
それでは、よきコーディングライフを！