Sol
设
[A=egin{bmatrix}1 & 1 \ 1 & 0end{bmatrix}
]
那么要求的相当于是
[sum_{i=0}^{n}[k|i]inom{n}{i}A^i
]
求出其中的 (A_{0,0}) 即可
引入单位根(单位根反演?)
[[n mid k]=frac{1}{n}sum_{i=0}^{n-1} omega_{n}^{ik}
]
证明：
若 (n|k)，那么根据单位根的消去引理可以得到就是 (1)
否则，等比数列求和，发现分子为 (0)
所以带入单位根
[sum_{i=0}^{n}[k|i]inom{n}{i}A^i=frac{1}{k}sum_{i=0}^{n}inom{n}{i}A^isum_{j=0}^{k-1}w_k^{ij}=frac{1}{k}sum_{i=0}^{k-1}sum_{j=0}^{n}inom{n}{j}A^jw_k^{ij}
]
由二项式定理得到
[frac{1}{k}sum_{i=0}^{k-1}(Aw_k^i+I)^n
]
(I) 为单位矩阵
这道题 (k=p-1)
所以直接求出 (p) 的原根 (g)，令 (w_k^i=g^{frac{p-1}{k}i}) 即可
只要矩乘+快速幂就好了
# include <bits/stdc++.h>
using namespace std;
typedef long long ll;
int k, test, mod, g, pr[233333], cnt;
ll n;
struct Matrix {
ll a[2][2];
inline Matrix() {
a[0][0] = a[0][1] = a[1][0] = a[1][1] = 0;
}
inline Matrix operator *(Matrix b) const {
register Matrix c;
c.a[0][0] = (a[0][0] * b.a[0][0] + a[0][1] * b.a[1][0]) % mod;
c.a[0][1] = (a[0][0] * b.a[0][1] + a[0][1] * b.a[1][1]) % mod;
c.a[1][0] = (a[1][0] * b.a[0][0] + a[1][1] * b.a[1][0]) % mod;
c.a[1][1] = (a[1][1] * b.a[1][1] + a[1][0] * b.a[0][1]) % mod;
return c;
}
inline Matrix operator +(Matrix b) const {
register Matrix c;
c.a[0][0] = a[0][0] + b.a[0][0] >= mod ? a[0][0] + b.a[0][0] - mod : a[0][0] + b.a[0][0];
c.a[0][1] = a[0][1] + b.a[0][1] >= mod ? a[0][1] + b.a[0][1] - mod : a[0][1] + b.a[0][1];
c.a[1][0] = a[1][0] + b.a[1][0] >= mod ? a[1][0] + b.a[1][0] - mod : a[1][0] + b.a[1][0];
c.a[1][1] = a[1][1] + b.a[1][1] >= mod ? a[1][1] + b.a[1][1] - mod : a[1][1] + b.a[1][1];
return c;
}
inline Matrix operator *(int b) const {
register Matrix c;
c.a[0][0] = a[0][0] * b % mod, c.a[0][1] = a[0][1] * b % mod;
c.a[1][0] = a[1][0] * b % mod, c.a[1][1] = a[1][1] * b % mod;
return c;
}
} ANS, I, A;
inline int Pow(ll x, int y) {
register ll ret = 1;
for (; y; y >>= 1, x = x * x % mod)
if (y & 1) ret = ret * x % mod;
return ret;
}
inline void Getroot() {
register int i, phi = mod - 1, x = phi, j;
for (cnt = 0, i = 2; i * i <= x; ++i)
if (x % i == 0) {
pr[++cnt] = i;
while (x % i == 0) x /= i;
}
if (x > 1) pr[++cnt] = x;
for (i = 2; i <= phi; ++i) {
for (x = 0, j = 1; !x && j <= cnt; ++j)
if (Pow(i, phi / pr[j]) == 1) x = 1;
if (x) continue;
g = i;
return;
}
}
inline Matrix MatrixPow(Matrix X, ll y) {
register Matrix RET = I;
for (; y; y >>= 1, X = X * X)
if (y & 1) RET = RET * X;
return RET;
}
int main() {
register int i, w, wn;
I.a[0][0] = I.a[1][1] = A.a[0][0] = A.a[0][1] = A.a[1][0] = 1;
scanf("%d", &test);
while (test) {
scanf("%lld%d%d", &n, &k, &mod), --test;
ANS = Matrix(), Getroot(), wn = 1, w = Pow(g, (mod - 1) / k);
for (i = 0; i < k; ++i) ANS = ANS + MatrixPow(A * wn + I, n), wn = (ll)wn * w % mod;
ANS = ANS * Pow(k, mod - 2);
printf("%lld
", ANS.a[0][0]);
}
return 0;
}