Sol
矩阵乘法+快速幂+欧拉定理.
首先观察题目可以发现 (A_n) 可以表示成关于 (K) 和 (A_0) 的几次幂的形式.
(A_0) 就比较简单了 (m^n) 所以第一部分 (ans1=A_0^{m^n}) .
看看 (k) 找一下规律就可以发现, (K) 的次幂形式是 (m^{n-1}+2m^{n-2}+3m^{n-3}+...+nm^0) .
这个东西可以构造一个矩阵,来递推出来.矩阵里需要有 (3) 个变量 (ans,b,1) .
其中 (ans) 是当前答案, (b) 是循环变量, (1) 就是 (b) 每次递增的数值.
[egin{bmatrix}ans\ b\ 1end{bmatrix}]
就是这样的一个矩阵,我们想它变成
[egin{bmatrix}ans*m+b\ b+1\ 1end{bmatrix}]
显然啊.矩阵很好构造.
[egin{bmatrix}m&1&0\0&1&1\0&0&1end{bmatrix}]
这个矩阵快速幂就可以了,但是注意是 (n-1) 次幂.
然后这道题还有个坑点就是 (p) 不是质数,但是保证互质,如果是 (k) 或 (A_0) 我们只需要模关于 (p) 的逆元的次幂就可以了,这里可以直接用欧拉定理.
[a^{varphi(n) }equiv 1(mod n)]
然后就成了模板题了...满满的全是模板...
md!sbt卡常数.写个快速乘.变成 (log^2) T到死.需要优化这个常数,把一个 (log) 优化掉就可以.
快速乘可以分段来写先乘前 (10^6) 的数字再乘后 (10^6) 的数字就可以了.
之前没想快速乘的这个 (log) 感觉可以过...然后为了卡常数...什么都写了...
Code
/************************************************************** Problem: 3665 User: BeiYu Language: C++ Result: Accepted Time:12284 ms Memory:20836 kb ****************************************************************/ #include<cstdio> #include<cmath> #include<cstring> #include<algorithm> #include<vector> #include<iostream> using namespace std; typedef long long LL; //typedef vector<LL> Vec; //typedef vector<Vec> Mat; struct Mat{ LL a[3][3]; Mat(){ memset(a,0,sizeof(a)); } }; LL T,p,m,n,k,a0,phi; LL ans1,ans2; char *ps=(char *)malloc(20000000); inline LL in(LL x=0){ for(;*ps>'9'||*ps<'0';ps++); for(;*ps>='0'&&*ps<='9';ps++) x=(x<<3)+(x<<1)+*ps-'0';return x; } inline void Out(LL x){ int l=0;char ch[65]; if(!x){ putchar('0');return; } if(x<0) putchar('-'),x=-x; while(x) ch[++l]=x%10+'0',x/=10; for(int i=l;i;i--) putchar(ch[i]); } inline LL GetPhi(LL p){ LL res=p,m=sqrt(p)+0.5; for(int i=2;i<=m;i++) if(p%i==0){ res=res/i*(i-1); while(p%i==0) p/=i; } if(p>1) res=res/p*(p-1);return res; } inline LL Mul(LL a,LL b,LL p){ LL t1=b/1000000,t2=b%1000000; return ((t1*a)%p*1000000%p+t2*a%p)%p; } inline LL Pow(LL a,LL b,LL p,LL res=1){ for(;b;b>>=1,a=Mul(a,a,p)) if(b&1) res=Mul(res,a,p);return res; } Mat operator * (const Mat &A,const Mat &B){ Mat C; for(int i=0;i<3;i++) for(int j=0;j<3;j++) for(int k=0;k<3;k++) C.a[i][j]=(C.a[i][j]+Mul(A.a[i][k],B.a[k][j],phi))%phi; return C; } Mat operator ^ (Mat A,LL b){ Mat res; for(int i=0;i<3;i++) for(int j=0;j<3;j++) res.a[i][j]=(i==j)?1:0; for(;b;b>>=1,A=A*A) if(b&1) res=res*A; return res; } int main(){ // freopen("maths.in","r",stdin); // freopen("maths.out","w",stdout); // ios::sync_with_stdio(false); fread(ps,1,20000000,stdin); Mat A; for(T=in(),p=in(),phi=GetPhi(p);T--;){ m=in()%phi,a0=in()%p,k=in()%p,n=in(); // cout<<Mul(n,Pow(n,5))<<endl; // cout<<m<<" "<<a0<<" "<<k<<" "<<n<<endl; ans1=Pow(a0,Pow(m,n,phi),p); // cout<<ans1<<endl;
A.a[0][0]=m,A.a[0][1]=1,A.a[0][2]=0; A.a[1][0]=0,A.a[1][1]=1,A.a[1][2]=1; A.a[2][0]=0,A.a[2][1]=0,A.a[2][2]=1; A=A^(n-1); ans2=Pow(k,((A.a[0][0]+2*A.a[0][1])%phi+A.a[0][2])%phi,p); // cout<<ans1<<" "<<ans2<<endl; Out(Mul(ans1,ans2,p)),putchar(' '); // printf("%I64d ",Mul(ans1,ans2,p)); } return 0; }