BZOJ3328: PYXFIB
题目：http://www.lydsy.com/JudgeOnline/problem.php?id=3328
题解：关键在于只处理i%k的项，那么我们就需要用一个式子来表达这个东西。
p%k==1.会让我们想到NTT的w=power(g,(p-1)/k)。而w的性质就是w^i=1%p当且仅当i%k=0。g是p的一个原根。
所以sigma（w^i）（0<=i<k)=0
然后我们构造一个A[x]=x^(-n)*(I*x+T)^n 其中I是单位矩阵，T是fib矩阵。
然后做A[W^0],A[W^-1]……A[W^-K+1]求和左上角就是K*ans。
因为sigma（w^ij）（0<=i<k)=k[j%k==0]
构造技巧实在太科幻，不愧是业界毒瘤出的题。
代码：
1 #include<cstdio> 2 #include<cstdlib> 3 #include<cmath> 4 #include<cstring> 5 #include<algorithm> 6 #include<iostream> 7 #include<vector> 8 #include<map> 9 #include<set> 10 #include<queue> 11 #include<string> 12 #define inf 1000000000 13 #define maxn 10000+5 14 #define maxm 10000000 15 #define eps 1e-10 16 #define ll long long 17 #define pa pair<int,int> 18 #define for0(i,n) for(int i=0;i<=(n);i++) 19 #define for1(i,n) for(int i=1;i<=(n);i++) 20 #define for2(i,x,y) for(int i=(x);i<=(y);i++) 21 #define for3(i,x,y) for(int i=(x);i>=(y);i--) 22 #define for4(i,x) for(int i=head[x],y=e[i].go;i;i=e[i].next,y=e[i].go) 23 #define lch t[k].l,l,mid 24 #define rch t[k].r,mid+1,r 25 using namespace std; 26 inline ll read() 27 { 28 ll x=0,f=1;char ch=getchar(); 29 while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();} 30 while(ch>='0'&&ch<='9'){x=10*x+ch-'0';ch=getchar();} 31 return x*f; 32 } 33 struct matrix 34 { 35 int d[3][3]; 36 }; 37 matrix A,I,T,ans; 38 ll n; 39 int p,k,g,w,pr[maxn]; 40 inline void print(matrix a) 41 { 42 for1(i,2)for1(j,2)printf("%d%c",a.d[i][j],j==2?'\n':' '); 43 } 44 matrix operator *(matrix a,matrix b) 45 { 46 matrix c=A; 47 for1(i,2)for1(j,2)for1(k,2)(c.d[i][j]+=(ll)a.d[i][k]*b.d[k][j]%p)%p; 48 return c; 49 } 50 matrix operator +(matrix a,matrix b) 51 { 52 matrix c; 53 for1(i,2)for1(j,2)c.d[i][j]=(a.d[i][j]+b.d[i][j])%p; 54 return c; 55 } 56 matrix operator *(matrix a,int b) 57 { 58 matrix c; 59 for1(i,2)for1(j,2)c.d[i][j]=(ll)a.d[i][j]*b%p; 60 return c; 61 } 62 inline int power(int x,int y) 63 { 64 int t=1; 65 for(;y;y>>=1,x=(ll)x*x%p) 66 if(y&1)t=(ll)t*x%p; 67 return t; 68 } 69 inline matrix power(matrix x,ll y) 70 { 71 matrix t=I; 72 for(;y;y>>=1,x=x*x) 73 if(y&1)t=t*x; 74 return t; 75 } 76 inline int gen() 77 { 78 int t=p-1,x=sqrt(t);pr[0]=0; 79 for2(i,2,x)if(t%i==0) 80 { 81 pr[++pr[0]]=i; 82 while(t%i==0)t/=i; 83 } 84 if(t>1)pr[++pr[0]]=t; 85 for2(i,2,inf) 86 { 87 bool flag=0; 88 for1(j,pr[0])if(power(i,(p-1)/pr[j])==1){flag=1;break;} 89 if(!flag)return i; 90 } 91 } 92 int main() 93 { 94 freopen("input.txt","r",stdin); 95 freopen("output.txt","w",stdout); 96 I.d[1][1]=I.d[2][2]=1; 97 T.d[1][1]=T.d[1][2]=T.d[2][1]=1; 98 int cs=read(); 99 while(cs--) 100 { 101 n=read();k=read();p=read();g=gen();w=power(g,(p-1)/k);ans=A; 102 for3(i,0,-k+1) 103 { 104 int x=power(w,k+i); 105 ans=ans+power(I*x+T,n)*power(x,((-n)%k+k)%k); 106 } 107 cout<<(ll)ans.d[1][1]*power(k,p-2)%p<<endl; 108 } 109 return 0; 110 }
