• bzoj3328: PYXFIB(单位根反演+矩阵快速幂)


    题面

    传送门

    题解

    我们设(A=egin{bmatrix}1 & 1 \ 1 & 0end{bmatrix}),那么(A^n)的左上角就是(F)的第(n)

    所以我们现在转化为求

    [sum_{i=0}^n[k|i]{nchoose i}A^i ]

    ([k|i])单位根反演一下

    [egin{aligned} ans &=sum_{i=0}^n[k|i]{nchoose i}A^i\ &={1over k}sum_{i=0}^n{nchoose i}A^isum_{j=0}^{k-1}omega^{ij}_k\ &={1over k}sum_{j=0}^{k-1}sum_{i=0}^n{nchoose i}A^iomega^{ij}_k\ &={1over k}sum_{j=0}^{k-1}left(Aomega_k^j+E ight)^n\ end{aligned} ]

    其中(E)表示单位矩阵

    注意,这里能用二项式定理是因为(A,E)之间满足交换律,如果是普通的矩阵是不能这么做的

    然后求个原根就行了(话说我今天才知道该怎么求原根)

    //minamoto
    #include<bits/stdc++.h>
    #define R register
    #define ll long long
    #define fp(i,a,b) for(R int i=a,I=b+1;i<I;++i)
    #define fd(i,a,b) for(R int i=a,I=b-1;i>I;--i)
    #define go(u) for(int i=head[u],v=e[i].v;i;i=e[i].nx,v=e[i].v)
    using namespace std;
    const int N=5e5+5;
    int p[N],P,g,tot,k,w,wn;ll n;
    inline int add(R int x,R int y){return x+y>=P?x+y-P:x+y;}
    inline int mul(R int x,R int y){return 1ll*x*y-1ll*x*y/P*P;}
    int ksm(R int x,R int y){
    	R int res=1;
    	for(;y;y>>=1,x=mul(x,x))if(y&1)res=mul(res,x);
    	return res;
    }
    struct Matrix{
    	int a[2][2];
    	inline Matrix(){a[0][0]=a[0][1]=a[1][0]=a[1][1]=0;}
    	inline int* operator [](const int &x){return a[x];}
    	Matrix operator *(Matrix b){
    		Matrix res;
    		res[0][0]=(1ll*a[0][0]*b[0][0]+1ll*a[0][1]*b[1][0])%P;
    		res[0][1]=(1ll*a[0][0]*b[0][1]+1ll*a[0][1]*b[1][1])%P;
    		res[1][0]=(1ll*a[1][0]*b[0][0]+1ll*a[1][1]*b[1][0])%P;
    		res[1][1]=(1ll*a[1][0]*b[0][1]+1ll*a[1][1]*b[1][1])%P;
    		return res;
    	}
    	Matrix operator +(Matrix b){
    		Matrix res;
    		res[0][0]=add(a[0][0],b[0][0]);
    		res[0][1]=add(a[0][1],b[0][1]);
    		res[1][0]=add(a[1][0],b[1][0]);
    		res[1][1]=add(a[1][1],b[1][1]);
    		return res;
    	}
    	Matrix operator *(const int &x){
    		Matrix res;
    		res[0][0]=mul(a[0][0],x);
    		res[0][1]=mul(a[0][1],x);
    		res[1][0]=mul(a[1][0],x);
    		res[1][1]=mul(a[1][1],x);
    		return res;
    	}
    }E,A,ans;
    Matrix ksm(Matrix x,ll y){
    	Matrix res;res[0][0]=res[1][1]=1;
    	for(;y;y>>=1,x=x*x)if(y&1)res=res*x;
    	return res;
    }
    void getrt(){
    	int phi=P-1,x=phi;tot=0;
    	for(R int i=2;1ll*i*i<=x;++i)if(x%i==0){
    		p[++tot]=i;
    		while(x%i==0)x/=i;
    	}
    	if(x>1)p[++tot]=x;
    	fp(i,2,phi){
    		x=0;
    		fp(j,1,tot)if(ksm(i,phi/p[j])==1){x=1;break;}
    		if(!x)return g=i,void();
    	}
    }
    int main(){
    //	freopen("testdata.in","r",stdin);
    	E[0][0]=E[1][1]=A[0][0]=A[0][1]=A[1][0]=1;
    	int T;scanf("%d",&T);
    	while(T--){
    		scanf("%lld%d%d",&n,&k,&P);
    		ans=Matrix(),getrt(),w=ksm(g,(P-1)/k),wn=1;
    		fp(i,0,k-1)ans=ans+ksm(A*wn+E,n),wn=mul(wn,w);
    		printf("%d
    ",mul(ans[0][0],ksm(k,P-2)));
    	}
    	return 0;
    }
    
  • 相关阅读:
    使用supervisor做进程控制
    HDU 4272 LianLianKan [状态压缩DP]
    UVALive 4297 First Knight [期望+高斯消元(优化)]
    HDU 4269 Defend Jian Ge [模拟]
    POJ 2497 Widget Factory [高斯消元解同余方程组]
    HDU 2996 In case of failure [KD树]
    HDU 4268 Alice and Bob [贪心]
    HDU 4273 Rescue [三维几何]
    HDU 4267 A Simple Problem with Integers [树状数组]
    HDU 4271 Find Black Hand [最短编辑距离DP]
  • 原文地址:https://www.cnblogs.com/bztMinamoto/p/10474567.html
Copyright © 2020-2023  润新知