[例題 1] 以下の合同方程式を解け
\begin{align}
&f(x) = x^3+8x-5のとき\\
&f(x) \equiv 0 \pmod{5^5} \\
&となるxを求めよ
\end{align}
このような問題を解くのに役に立つのがヘンゼルの補題(Wikipedia)です。このWikiの説明はかなり難解ですが、ぶっちゃけて言うと。
- $f(x) \equiv 0 \pmod{p}\ (pは素数)$の解が見つかれば$\pmod{p^k}$での解を逐次的に計算するアルゴリズムがある（条件あり）ということです。
今回の投稿では$\pmod{p^k}$の答えから$\pmod{p^{k+1}}$を求める「ヘンゼル持ち上げ」(hensel lift)のアルゴリズムをpythonで実装して、例題を解いてみたいと思います。
ヘンゼル持ち上げ(hensel lift)
\begin{align}
&f(x) \equiv 0 \pmod{p}\ (pは素数) の解をx_1とする \\
&k番目の式、f(x_k) \equiv 0 \pmod{p^{k}} が成り立つとき \\
&x_{k+1} \equiv x_k \pmod{p} となるので（ヘンゼルの補題より）\\
&x_{k+1} = x_k+p^ktと置ける \cdots (1)\\
&0 \equiv f(x_{k+1}) = f(x_k+p^kt) \\
&右辺をテイラー展開すると \\
&= f(x_k)+p^ktf'(x_k) + (p^kt)^2f''(x_k)/2!+... \\
&この３項目以降はp^{k+1}で割ると0になるので\\
& f(x_k)+p^ktf'(x_k) \equiv 0 \pmod {p{k+1}} \\
\end{align}
\begin{align}
&f(x_k) はp^kで割れるので、両辺をp^kで割る \\
&\frac{f(x_k)}{p^k}+tf'(x_k) \equiv 0 \pmod {p} \\
&x_k \equiv x_{k-1} ... \equiv x_1 \pmod {p} より f'(x_k) \equiv f'(x_1) \pmod {p}\\
&ここで a= modinv(f'(x1)) \pmod pと置くとｔが求まり \\
&t \equiv -a \frac{f(x_k)}{p^k} \\
&これを(1)に代入して \\
&x_{k+1} = x_k + p^kt = x_k - af(x_k) \pmod p^{k+1} \\
&以下の漸化式が求まります\\
\end{align}
\begin{array}{|c|} \hline \\
【ヘンゼル持ち上げの漸化式】 \\
\large \ \ \ \ x_{k+1} = x_k - af(x_k) \pmod p^{k+1} \ \ \ \ \ \\
ただし \ a = modinv(f'(x1)) \pmod p \\
\\ \hline
\end{array}
ヘンゼル持ち上げの実装（その１）
$modinv$の所をsympyのmod_inverseを使うと$x_k$から$x_{k+1}$を求める関数hl_stepはそのままの実装になります。$a=modinv(f'(x)) \pmod p$は一度計算すればよいの先に求めておきます。
from sympy import mod_inverse
def hl_step(x, a, k):
return (x - a*f(x))%p**k
p, x, x1, f(x), f1(x) = ...
a = mod_inverse(f1(x1),p)
例題では、$x_1=0$なので$k=5$までループして答えの$x_5 = 1985$が得られます。各々の$f(x_k) \equiv 0 \pmod {p^k}$となっていることも確認できました。
def f(x): return x**3 + 8*x - 5
def f1(x): return 3*x**2 + 8
def hl_step(x, a, k):
return (x - a*f(x))%p**k
p, x, x1 = 5, 0, 0
a = mod_inverse(f1(x1),p)
for k in range(2,5+1):
x = hl_step(x, a, k)
print(f"k={k}, x={x}, f({x})%{p}**{k}={f(x)%p**k}")
# k=2, x=10, f(10)%5**2=0
# k=3, x=110, f(110)%5**3=0
# k=4, x=110, f(110)%5**4=0
# k=5, x=1985, f(1985)%5**5=0
ヘンゼル持ち上げの実装（複数解対応）
【例題 1】は$x_1$の解が１つでしたが、三次方程式なら最大3つあるので、そのような複数解に対応するようなプログラムを実装します。
【例題 2】以下の合同方程式を解け
\begin{align}
&f(x) = x^3+3x-7のとき\\
&f(x) \equiv 0 \pmod{7^5} \\
&となるxをすべて求めよ
\end{align}
まず$x_1 = f(x) \equiv 0 \pmod{7}$の解ですが、以下のように$[0, 2, 5]$の3つあります。
def f(x): return x**3 + 3*x - 7
p = 7
x1list = [i for i in range(p) if f(i) % p == 0]
print(f"p = {p}, x1={x1list}")
# p = 7, x1=[0, 2, 5]
従って$x_1$のリストの各々に対してhl_stepを再帰的に呼びだすように実装しました。
【例題 2】の解は$[3318, 3061, 11275]$になります。
from sympy import mod_inverse
def hensel(f, f1, x1list, p, k):
#--- hensel lift step ---
def hl_step(x, a, k):
return (x - a*f(x))%p**k
def _hl(p, k): # return [(x1n, a),(x2n,a),...]
if k == 1: return [(i, mod_inverse(f1(i),p)) for i in x1list]
xl = _hl(p,k-1)
return [(hl_step(x, a, k), a) for (x, a) in xl ]
return [x for (x, a) in _hl(p,k)]
def f(x): return x**3 + 3*x - 7
def f1(x): return 3*x + 3
p = 7
x1list = [i for i in range(p) if f(i) % p == 0]
print(f"p = {p}, x1={x1list}")
for k in range(2,5+1):
hl = hensel(f, f1, x1list, p, k)
print(k, hl)
# p = 7, x1=[0, 2, 5]
# 2 [35, 23, 33]
# 3 [231, 156, 5]
# 4 [917, 1374, 2140]
# 5 [3318, 3061, 11275]
この考え方はEuler Project: Problem 457: A polynomial modulo the square of a primeと解くのに役に立ちます。
（開発環境：Google Colab）