【例題】N=500のとき以下のf(x)が911で割り切れるようなxの合計を求めよ
\Large f(x) = x^{13} + 1 \ \ \ \ (1 \le x \le N)
素朴にループで解く
もちろんN=500程度であればこのような単純ループで解くことが出来ますがNが大きくなると$f(x)$の値が大きくなり効率がよくありません。そこでタイトルのようにn次合同方程式を解くことを考えます。
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の力で解く！【数学 整数論】(みつのきチャンネル)
\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を使って求めておきます。
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
\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を作りました。これを使って例題を解くと当然同じ答えになります。
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次合同方程式を解く形では以下のようにして解くことができます。
\begin{align}
&例えば例題 1の解の一つが11のとき解はmod\ 911で同じなら解なので \\
&x = 11 + 911 \times y (x \le N)\\
&も解になります。\\
\end{align}
これらの解の和は等差数列の和の公式から求まるので、その関数get_sumは
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
となります。これを既に求まった解すべてに適用すればよいので
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次合同方程式の解き方を勉強するためにアルゴリズムを理解することは重要だと考えたのでこの記事を書いてみました。
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を解くのに役に立ちます。