BZOJ 3328: PYXFIB 解题报告
BZOJ 3328: PYXFIB
题意
给定\(n,p,k(1\le n\le 10^{18},1\le k\le 20000,1\le p\le 10^9,p \ is \ prime,k|(p-1))\)，求
\[\sum_{i=0}^{\lfloor\frac{n}{k}\rfloor}\binom{n}{ik}F_{ik}
\]
其中\(F_{i}\)代表斐波那契数列的第\(i\)项。
首先分析一波题目条件，有个很奇怪的条件\(k|\varphi(p)\)
如果对NTT或者单位根有所了解，我们可以猜到这个条件可以取出单位根来。
\[w_k=g^{\frac{p-1}{k}}
\]
然而如果没发现，推一波式子之后也可以发现需要使用单位根。
不过可以先考虑如何处理fib，一个简单的想法是构造矩阵乘法，比如
\[A=\begin{bmatrix}1&1\\1& 0\end{bmatrix}
\]
可以得到
\[A^i=\begin{bmatrix}F_i&F_{i-1}\\\dots&\dots\end{bmatrix}
\]
然后我们直接把式子转换成求矩阵
\[\begin{aligned}
&\sum_{i=0}^{\lfloor\frac{n}{k}\rfloor}\binom{n}{ik}A^{ik}\\
=&\sum_{i=0}^n\binom{n}{i}A^i[k|i]\\
=&\sum_{i=0}^n\binom{n}{i}A^i\frac{1}{k}\sum_{j=0}^{k-1}w_k^{ij}\\
=&\frac{1}{k}\sum_{j=0}^{k-1}\sum_{i=0}^n\binom{n}{i}A^i(w_k^j)^i\\
=&\frac{1}{k}\sum_{j=0}^{k-1}(w_k^jA+I)^n
\end{aligned}
\]
复杂度\(O(Tk\log 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\n",mul(ans,qp(k,p-2)));
}
return 0;
}
2019.5.8
