BZOJ 3328: PYXFIB
题意
给定(n,p,k(1le nle 10^{18},1le kle 20000,1le ple 10^9,p is prime,k|(p-1)))，求
[sum_{i=0}^{lfloorfrac{n}{k}
floor}inom{n}{ik}F_{ik}
]
其中(F_{i})代表斐波那契数列的第(i)项。
首先分析一波题目条件，有个很奇怪的条件(k|varphi(p))
如果对NTT或者单位根有所了解，我们可以猜到这个条件可以取出单位根来。
[w_k=g^{frac{p-1}{k}}
]
然而如果没发现，推一波式子之后也可以发现需要使用单位根。
不过可以先考虑如何处理fib，一个简单的想法是构造矩阵乘法，比如
[A=egin{bmatrix}1&1\1& 0end{bmatrix}
]
可以得到
[A^i=egin{bmatrix}F_i&F_{i-1}\dots&dotsend{bmatrix}
]
然后我们直接把式子转换成求矩阵
[egin{aligned}
&sum_{i=0}^{lfloorfrac{n}{k}
floor}inom{n}{ik}A^{ik}\
=&sum_{i=0}^ninom{n}{i}A^i[k|i]\
=&sum_{i=0}^ninom{n}{i}A^ifrac{1}{k}sum_{j=0}^{k-1}w_k^{ij}\
=&frac{1}{k}sum_{j=0}^{k-1}sum_{i=0}^ninom{n}{i}A^i(w_k^j)^i\
=&frac{1}{k}sum_{j=0}^{k-1}(w_k^jA+I)^n
end{aligned}
]
复杂度(O(Tklog n))
Code:
#include <cstdio>
#include <cctype>
#define ll long long
const int SIZE=1<<21;
char ibuf[SIZE],*iS,*iT;
//#define gc() (iS==iT?(iT=(iS=ibuf)+fread(ibuf,1,SIZE,stdin),iS==iT?EOF:*iS++):*iS++)
#define gc() getchar()
template <class T>
void read(T &x)
{
x=0;char c=gc();
while(!isdigit(c)) c=gc();
while(isdigit(c)) x=x*10+c-'0',c=gc();
}
ll n;
int k,p;
inline int add(int a,int b){return a+b>=p?a+b-p:a+b;}
#define mul(a,b) (1ll*(a)*(b)%p)
inline int qp(int d,int k)
{
int f=1;
while(k)
{
if(k&1) f=mul(f,d);
d=mul(d,d);
k>>=1;
}
return f;
}
struct Matrix
{
int a,b,c,d;
Matrix friend operator *(Matrix a,Matrix b)
{
Matrix c;
c.a=add(mul(a.a,b.a),mul(a.b,b.c));
c.b=add(mul(a.a,b.b),mul(a.b,b.d));
c.c=add(mul(a.c,b.a),mul(a.d,b.c));
c.d=add(mul(a.c,b.b),mul(a.d,b.d));
return c;
}
}A;
inline Matrix mqp(Matrix d,ll k)
{
Matrix f=d;--k;
while(k)
{
if(k&1) f=f*d;
d=d*d;
k>>=1;
}
return f;
}
int getg(int p)
{
int s[30]={},phi=p-1;
for(int i=2;i*i<=phi;i++)
{
if(phi%i==0)
{
s[++s[0]]=i;
while(phi%i==0) phi/=i;
}
}
if(phi!=1) s[++s[0]]=phi;
for(int i=2;;i++)
{
int flag=1;
for(int j=1;j<=s[0];j++)
if(qp(i,(p-1)/s[j])==1)
{
flag=0;
break;
}
if(flag) return i;
}
}
int main()
{
int T;read(T);
while(T--)
{
read(n),read(k),read(p);
int w=qp(getg(p),(p-1)/k),ans=0;
for(int j=0;j<k;j++)
{
A.a=A.b=A.c=qp(w,j),A.d=1;
A.a=add(A.a,1);
A=mqp(A,n);
ans=add(ans,A.a);
}
printf("%lld
",mul(ans,qp(k,p-2)));
}
return 0;
}
2019.5.8