PNP    P(N1)0P \equiv NP ~~ \rightarrow ~~ P(N - 1) \equiv 0

と式変形をします. すると,N1N-1の範囲は[1,M2][-1, M-2]になりますが,1M1-1 \equiv M-1 であり,NNは整数であるため,N1N-1の範囲を[0,M1][0, M-1]としても問題ありません.よって,N=N1N' = N-1とすると,

0P,N<M0 \leq P, N' \lt M

PN0PN' \equiv 0

を満たすP,NP, N'の個数を求めてあげればよいです.

ここで, PP の値を固定した時に条件を満たす NN'の個数を考えます.これは,PPMMの最大公約数をggと置くと,NN'(M/g)(M / g)の倍数の時に与えられた式を満たすため,条件を満たすNN'の個数は M/(M/g)=gM / (M / g) = g 個になります.よって,本問題は以下のように言い換えることができます.

  • 正整数MMが与えられる.gcd(a,b)\operatorname{gcd} (a, b)a,ba, bの最大公約数とするとき, k=0M1gcd(k,M)\sum_{k=0}^{M-1} \operatorname{gcd}(k, M) を求めてください

以降,言い換えた問題を考えます.

M=a×bM = a \times b, gcd(a,b)=1\operatorname{gcd}(a, b) = 1 となるような(a,b)(a, b)を考えると,

gcd(i,M)=gcd(i,a)×gcd(i,b)\operatorname{gcd}(i, M) = \operatorname{gcd}(i, a) \times \operatorname{gcd}(i, b)

が成り立ちます.この変換を用いると,

k=0M1gcd(k,M)=k=0M1(gcd(k,a)×gcd(k,b))\sum_{k=0}^{M-1} \operatorname{gcd}(k, M) = \sum_{k=0}^{M-1} \left(\operatorname{gcd}(k, a) \times \operatorname{gcd}(k, b)\right)

となります.ここで,gcd(k,a)=gcd(k,a%k)\operatorname{gcd}(k, a) = \operatorname{gcd}(k, a \% k) と,中国の剰余定理より「互いに素である (a,b)(a, b)a×b=Ma \times b = M を満たすとき,0以上MM未満の正整数の中に,各 0i<a0 \leq i \lt a, 0j<b0 \leq j \lt b について aa で割ると ii 余り,bb で割ると jj 余る数が丁度1つずつ存在する」ことを利用すると,

k=0M1gcd(k,a)×gcd(k,b) \sum_{k=0}^{M-1} \operatorname{gcd}(k, a) \times \operatorname{gcd}(k, b)

=(gcd(0,a)+gcd(1,a)++gcd(a1,a))×(gcd(0,b)+gcd(1,b)++gcd(b1,b)) = \left(\operatorname{gcd}(0, a) + \operatorname{gcd}(1, a) + \cdots + \operatorname{gcd}(a - 1, a)\right) \times \left(\operatorname{gcd}(0, b) + \operatorname{gcd}(1, b) + \cdots + \operatorname{gcd}(b - 1, b) \right)

=k=0a1gcd(k,a)×k=0b1gcd(k,b)= \sum_{k=0}^{a - 1} \operatorname{gcd}(k, a) \times \sum_{k=0}^{b - 1} \operatorname{gcd}(k, b)

を解けば良いことになります.よって,MMが素数pip_i及び正整数eie_iを用いて

M=pieiM = \prod p_i ^ {e_i}

と表されるとき,

(k=0piei1gcd(k,piei)) \prod \left( \sum_{k=0}^{p_i^{e^i}-1} \operatorname{gcd}(k, p_i^{e_i}) \right)

が答えになります.

後は,各pi,eip_i, e_i について,k=0piei1gcd(k,piei)\sum_{k=0}^{p_i^{e^i}-1} \operatorname{gcd}(k, p_i^{e_i})の値を求めてあげればよいです. この答えについては,[0,piei)[0, p_i^{e_i})の範囲内にpieip_i^{e_i}の倍数が1個,piei1p_i^{e_i-1}の倍数がpi{p_i}個,piei2p_i^{e_i-2}の倍数がpi2{p_i^2}個 ... ずつあることを利用すると,最大公約数の総和は

piei+k=0ei1pik×(pieikpieik1)p_i^{e_i} + \sum_{k=0}^{e_i-1} p_i^{k} \times (p_i^{e_i-k} - p_i^{e_i-k - 1})

=piei+ei×(pieipiei1)= p_i^{e_i} + e_i \times (p_i^{e_i} - p_i^{e_i-1})

=piei1×(pi+(pi1)×ei)= p_i^{e_i - 1} \times (p_i + (p_i - 1) \times e_i)

式で計算できます. よって,各素因数ごとに登場回数を調べ,最後の式で計算される値を掛け合わせることによって,本問題を解くことが出来ます.

計算量は,各素因数ごとに登場回数を数える箇所に依存し,C++ の map 等を使って求めるとO(klogk)O(k \log k)に,unordered_map を使ったり入力がソートされた状態で与えられることを利用するなどして登場回数を数えると O(k)O(k) になります.

想定解 (C++ with ACL)

想定解 (Python)