ヘンゼルの補題を使って法がp**kの合同方程式を解く

### [例題 1] 以下の合同方程式を解け

```math
\begin{align}
&f(x) = x^3+8x-5のとき\\
&f(x) \equiv 0 \pmod{5^5}　\\
&となるxを求めよ
\end{align}
```

このような問題を解くのに役に立つのが[ヘンゼルの補題(Wikipedia)](https://ja.wikipedia.org/wiki/%E3%83%98%E3%83%B3%E3%82%BC%E3%83%AB%E3%81%AE%E8%A3%9C%E9%A1%8C)です。このWikiの説明はかなり難解ですが、ぶっちゃけて言うと。

* $f(x) \equiv 0 \pmod{p}\ (pは素数)$の解が見つかれば$\pmod{p^k}$での解を逐次的に計算するアルゴリズムがある（条件あり）ということです。

今回の投稿では$\pmod{p^k}$の答えから$\pmod{p^{k+1}}$を求める「ヘンゼル持ち上げ」(hensel lift)のアルゴリズムをpythonで実装して、例題を解いてみたいと思います。

### ヘンゼル持ち上げ(hensel lift)

```math 
\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}
```
```math
\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}
```
```math
\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$は一度計算すればよいの先に求めておきます。

```python
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}$となっていることも確認できました。

```python 
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】以下の合同方程式を解け

```math
\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つあります。

```python 
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]$になります。

```python 
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）