Description
给定
n
,
k
,
p
n,k,p
n,k,p,保证
k
≡
1
(
m
o
d
p
)
k\equiv 1\pmod p
k≡1(modp)
求
∑
i
=
0
⌊
n
k
⌋
(
n
k
i
)
F
(
k
i
)
(
m
o
d
p
)
\sum_{i=0}^{\lfloor\frac{n}{k}\rfloor}{\binom{n}{ki}F(ki)}\pmod {p}
i=0∑⌊kn⌋(kin)F(ki)(modp)
其中F(i)是斐波那契数列第i项
n
≤
1
0
18
,
k
≤
20000
,
p
≤
1
0
9
n\le10^{18},k\le20000,p\le{10^9}
n≤1018,k≤20000,p≤109保证p为质数
Solution
很妙的姿势,感觉这个不太好出题啊,做法貌似很套路
我们直接枚举ki,那么柿子就是
∑
i
=
0
n
[
k
∣
i
]
(
n
i
)
F
(
i
)
(
m
o
d
p
)
\sum_{i=0}^{n}{[k|i]\binom{n}{i}F(i)}\pmod {p}
i=0∑n[k∣i](in)F(i)(modp)
然后单位根反演有个结论,就是
[
n
∣
k
]
=
1
n
∑
i
=
0
n
−
1
w
n
k
i
[n|k]=\frac{1}{n}\sum_{i=0}^{n-1}{{w_n}^{ki}}
[n∣k]=n1i=0∑n−1wnki
这是一个公比鬼畜的等比数列求和,证明的话我们分情况算一下就行了
往上面代就有
∑
i
=
0
n
1
k
∑
j
=
0
k
−
1
w
k
i
j
(
n
i
)
F
(
i
)
\sum_{i=0}^{n}\frac{1}{k}\sum_{j=0}^{k-1}{{w_k}^{ij}\binom{n}{i}F(i)}
i=0∑nk1j=0∑k−1wkij(in)F(i)
交换求和符号然后交换字母,那么就是
1
k
∑
i
=
0
k
−
1
∑
j
=
0
n
w
k
i
j
F
(
i
)
(
n
j
)
\frac{1}{k}\sum_{i=0}^{k-1}\sum_{j=0}^n{{w_{k}}^{ij}F(i)\binom{n}{j}}
k1i=0∑k−1j=0∑nwkijF(i)(jn)
根据斐波那契数列的矩阵推法,我们设
A
=
[
1
1
1
0
]
A= \left[ \begin{matrix} 1 & 1 \\ 1 & 0 \\ \end{matrix} \right]
A=[1110]
那么
F
(
i
)
=
A
i
[
0
,
0
]
F(i)={A^i}[0,0]
F(i)=Ai[0,0]。这里矩阵可能长成别的样子,同样也是对的
于是那个就可以二项式弄回去了,也就是
1
k
∑
i
=
0
k
−
1
(
w
k
i
A
+
I
)
n
\frac{1}{k}\sum_{i=0}^{k-1}{\left({w_k}^{i}A+I\right)}^{n}
k1i=0∑k−1(wkiA+I)n
这里
I
I
I是单位矩阵,
A
∗
I
=
A
A*I=A
A∗I=A
注意到这里的答案是膜意义下的,那么 w k = g p − 1 k w_k=g^{\frac{p-1}{k}} wk=gkp−1这里g是p的原根,这也是为什么题目保证了p膜k余1
于是我们求一个原根,做快速幂和矩阵快速幂就可以得到答案了
Code
#include <stdio.h>
#include <string.h>
#include <algorithm>
#include <vector>
#define rep(i,st,ed) for (int i=st;i<=ed;++i)
#define fill(x,t) memset(x,t,sizeof(x))
typedef long long LL;
int v[233333],MOD;
void upd(LL &x,LL v) {
x+=v; (x>=MOD)?(x-=MOD):0;
}
struct Mat {
LL r[2][2];
Mat() {fill(r,0);}
LL* operator [](int x) {
return r[x];
}
Mat operator *(LL x) {
Mat C=*this;
C[0][0]=C[0][0]*x%MOD;
C[0][1]=C[0][1]*x%MOD;
C[1][0]=C[1][0]*x%MOD;
C[1][1]=C[1][1]*x%MOD;
return C;
}
Mat operator +(Mat B) {
Mat C=*this;
upd(C[0][0],B[0][0]);
upd(C[0][1],B[0][1]);
upd(C[1][0],B[1][0]);
upd(C[1][1],B[1][1]);
return C;
}
Mat operator *(Mat B) {
Mat A=*this,C;
C[0][0]=(A[0][0]*B[0][0]+A[0][1]*B[1][0])%MOD;
C[0][1]=(A[0][0]*B[0][1]+A[0][1]*B[1][1])%MOD;
C[1][0]=(A[1][0]*B[0][0]+A[1][1]*B[1][0])%MOD;
C[1][1]=(A[1][0]*B[0][1]+A[1][1]*B[1][1])%MOD;
return C;
}
} A,I,B,S;
LL read() {
LL x=0,v=1; char ch=getchar();
for (;ch<'0'||ch>'9';v=(ch=='-')?(-1):v,ch=getchar());
for (;ch<='9'&&ch>='0';x=x*10LL+ch-'0',ch=getchar());
return x*v;
}
LL ksm(LL x,LL dep) {
LL res=1;
for (;dep;dep>>=1,x=x*x%MOD) {
(dep&1)?(res=res*x%MOD):0;
}
return res;
}
Mat ksm(Mat A,LL dep) {
Mat C=I;
for (;dep;dep>>=1,A=A*A) {
if (dep&1) C=C*A;
}
return C;
}
int get_g(int MOD) {
int tmp=MOD-1; v[0]=0;
for (int i=2;i*i<=tmp;++i) {
if (tmp%i==0) {
while (tmp%i==0) tmp/=i;
v[++v[0]]=i;
}
}
if (tmp!=1) v[++v[0]]=tmp;
rep(g,2,MOD) {
bool flag=false;
rep(i,1,v[0]) {
if (ksm(g,(MOD-1)/v[i])==1) {
flag=true; break;
}
}
if (!flag) return g;
}
}
int main(void) {
A[0][0]=A[0][1]=A[1][0]=I[0][0]=I[1][1]=1;
int T; scanf("%d",&T);
for (;T--;) {
LL n,k; scanf("%lld%lld%d",&n,&k,&MOD);
LL g=get_g(MOD),wn=ksm(g,(MOD-1)/k),w=1;
LL ans=0;
rep(i,0,k-1) {
S=ksm(A*w+I,n);
upd(ans,S[0][0]);
w=w*wn%MOD;
}
ans=ans*ksm(k,MOD-2)%MOD;
printf("%lld\n", ans);
}
return 0;
}