Sympyを使ってn次合同方程式を解く

### 【例題】N=500のとき以下のf(x)が911で割り切れるようなxの合計を求めよ

```math 
\Large f(x) = x^{13} + 1  \ \ \ \  (1 \le x \le N)
```

### 素朴にループで解く

もちろんN=500程度であればこのような単純ループで解くことが出来ますがNが大きくなると$f(x)$の値が大きくなり効率がよくありません。そこでタイトルのようにn次合同方程式を解くことを考えます。

```python 
ans = [x for x in range(1,500+1) if (x**13+1)%911 == 0]
print(f"Answer: {sum(ans)} ({ans})")
# Answer: 1606 ([11, 14, 330, 334, 420, 497])
```

### n次合同方程式を解く

はっきり言ってn次合同方程式を解くのは簡単ではありませんが、こちらの動画を何度も見てなんとかやってみました。
> [【初等整数論44】11x≡5 (mod 13) や x^3≡5 (mod 13) をIndの力で解く！【数学　整数論】(みつのきチャンネル)](https://www.youtube.com/watch?v=epatfpO3gak)

```math 
\begin{align}
&まずf(x)が911で割り切れるので \\
&x^{13} + 1  \equiv 0 \pmod{911} \\
&移項して正の数に変えます \\
&x^{13} \equiv -1 = 910 \pmod{911} \cdots (1)\\
\end{align}
```
次に必要となる911の原始根の一つ**17**をsympyの**primitive_root**で求めます。
また後で必要になる$Ind_{17}910$を**discrete_log**を使って求めておきます。

```python 
from sympy.ntheory import primitive_root, discrete_log
p = 911
g = primitive_root(p)
dl = discrete_log(p,p-1,g)
print(f"p = {p}, primitive_root ={g}, Ind({g},{p-1})={dl}")
# p = 911, primitive_root =17, Ind(17,910)=455
```

```math 
\begin{align}
&原始根が17なので式(1)の両辺のInd_{17}を取ります　\\
&このときmodはp-1の910になります　\\
&Ind_{17}x^{13} \equiv  Ind_{17}910 \pmod{910}  \\
&13が前に出て、Ind_{17}910=455なので\\
&13 \times Ind_{17}x \equiv  455 \pmod{910} \\
&両辺を13で割ります（modも含めて）ここで解が13個あることが分かります。\\
&Ind_{17}x \equiv  35 \pmod{70} \\
&以下の13個です\\
&Ind_{17}x = (35, 105, 175, 245, 315, 385, 455, 525, 595, 665, 735, 805, 875) \\
&これらをaとするとx=17^a \pmod {911}なので　\\
&x = (11, 14, 330, 334, 420, 497, 715, 757, 783, 790, 846, 881, 910) \\
&となります\\
\end{align}
```
従って例題の答えはこの$x$の500以下のもの$(11, 14, 330, 334, 420, 497)$となりループでの解答と一致しました。

### n次合同方程式を解く汎用関数rootn_mod

このアルゴリズムを汎用化して$x^m \equiv a \pmod p$を解く関数**rootn_mod**を作りました。これを使って例題を解くと当然同じ答えになります。

```python 
from sympy.ntheory import primitive_root, discrete_log
from sympy import mod_inverse
import math

def rootn_mod(a,m,p):     # Return x array for x**m = a (mod p)
  g = primitive_root(p)       # primitive root of p
  Inda = discrete_log(p,a,g)  # Ind (a,g) mod p
  d = math.gcd(m, p-1)        # d: number of solutions
  if Inda%d != 0: return []   # No solution
  if Inda == 0: return [1]    # x=1
  m //= d; Inda //=d; p1 = (p-1)//d   # divid both sides by d
  Inda = (Inda*mod_inverse(m, p1))%p1    # ax = b -> x = b * mod_inverse(a)
  return sorted([pow(g,a,p) for a in range(Inda,p,p1)])

m, p, N = 13, 911, 500
ans = [a for a in rootn_mod(p-1,m,p) if a <= N]
print(f"m={m}, p={p}, answer={ans}")
# m=13, p=911, answer=1606([11, 14, 330, 334, 420, 497])
```
さらにこのアルゴリズムが有効である証として$N=10^{10}$として例題を解いてみます。

### 【例題 2】N=10**10として例題 1を解け

これだけNが大きくなるとループでは不可能になりますが、n次合同方程式を解く形では以下のようにして解くことができます。

```math 
\begin{align}
&例えば例題 1の解の一つが11のとき解はmod\ 911で同じなら解なので \\
&x = 11 + 911 \times y　(x \le N)\\
&も解になります。\\
\end{align}
```
これらの解の和は等差数列の和の公式から求まるので、その関数**get_sum**は

```python 
def get_sum(a, b, d):  # return a+(a+d)+(a+2d)+,,,  < b
    if a > b: return 0
    n = (b - a)//d + 1
    return n*(2*a+(n-1)*d)//2
```
となります。これを既に求まった解すべてに適用すればよいので

```python 
m, p, N = 13, 911, 10**10
print(sum([get_sum(a, N, p) for a in rootn_mod(p-1,m,p)]))
# 713501648457738379
```
と簡単に求まりました。

### Sympyのnthroot_mod

ちなみにSympyには**rootn_mod**と同等の関数**nthroot_mod**が既に存在してそれは使えば同様に簡単に書く事が出来ました。
ただn次合同方程式の解き方を勉強するためにアルゴリズムを理解することは重要だと考えたのでこの記事を書いてみました。

```python
from sympy.ntheory import nthroot_mod
m, p, N = 13, 911, 10**10
print(sum([get_sum(a, N, p) for a in nthroot_mod(p-1,m,p, True)])) 
```


（開発環境：Google Colab）

このアルゴリズムはProject Euler, Problem 421: Prime factors of n15+1を解くのに役に立ちます。