首页 > 代码库 > BZOJ 3665: maths

BZOJ 3665: maths

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\) 每次递增的数值.

\[\begin{bmatrix}ans\\ b\\ 1\end{bmatrix}\]

就是这样的一个矩阵,我们想它变成

\[\begin{bmatrix}ans*m+b\\ b+1\\ 1\end{bmatrix}\]

显然啊.矩阵很好构造.

\[\begin{bmatrix}m&1&0\\0&1&1\\0&0&1\end{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;//  if(p<=1000000000) return a*b%p;//  return (a*b-(LL)(a/(long double)p*b+1e-3)*p+p)%p;//  for(;b;b>>=1,a=(a+a)%p) if(b&1) res=(res+a)%p;return res;}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<<"0"<<endl;//      cout<<m<<" "<<a0<<" "<<k<<" "<<n<<endl;         //      cout<<ans1<<endl;        ans1=Pow(a0,Pow(m,n,phi),p);//      cout<<ans1<<endl;         //      cout<<"1"<<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);         //      cout<<"2"<<endl;                 ans2=Pow(k,((A.a[0][0]+2*A.a[0][1])%phi+A.a[0][2])%phi,p);         //      cout<<"3"<<endl;//      cout<<ans1<<" "<<ans2<<endl;        Out(Mul(ans1,ans2,p)),putchar(‘\n‘);//      printf("%I64d\n",Mul(ans1,ans2,p));    }    return 0;}

  

BZOJ 3665: maths