【矩阵乘法&快速乘模板】【题解】【刷题比赛】

P1707 刷题比赛

Solution:
很明显是一道矩阵乘法加速递推，只是推转移矩阵比较复杂。
做这种题的思路是：

1. 如果递推为线性递推，且暴力递推次数很多，可以考虑矩阵优化
2. 先根据要计算的信息，列出一个行数为1的状态矩阵，设为\(F_n\)。
3. 再仿照\(F_n\)写出\(F_{n+1}\)。
4. 再对\(F_{n+1}\)的每一项分别考虑，看这一项由\(Fn\)的哪些项乘上相应的系数转移而来，构造出转移矩阵\(A\)
5. 写出\(F_1\)，则\(F_n=F_1*A^{n-1}\),用矩阵快速幂的模板计算即可

另外，本题在做乘法时有可能爆long long，需要用快速乘或龟速乘来避免

龟速乘

与快速幂的思路十分类似。若要计算\(x*y\)，设\(y=a_0*2^0+a_1*2^1+\cdots +a_n*2^n\)
那么由乘法分配律得，
\( \begin{aligned} x*y&=x*2^0*a_0+x*2^1*a_1+\cdots +x*2^n*a_n\\ &=x*2^0[a_0=1]+x*2^1[a_1=1]+\cdots +x*2^n[a_n=1] \end{aligned} \)
那么看代码就很容易了

#define int long long
int fmul(int x,int y,int mod)
{
	int res=0;
	while(y)
	{
		if(y&1) res+=x,res%=mod;
		y>>=1;
		x<<=1;
		x%=mod;
	}
	return res;
}

快速乘

龟速乘的复杂度是带个log的，因而出现了O(1)的快速乘。
其原理跟long double的精度有关，在此不做讨论，只给出代码。

#define int long long
int fmul(int x,int y,int mod)
{
	unsigned int res=
	(unsigned int)x*y-
	(unsigned int)((long double)x/mod*y+0.5L)*mod;
	return (res+mod)%mod;
}

最后给出本题代码

#include<bits/stdc++.h>
using namespace std;
#define int long long
inline int read()
{
	register int x=0,w=1;
	register char ch=getchar();
	while((ch<'0'||ch>'9')&&ch!='-') ch=getchar();
	if(ch=='-') {ch=getchar();w=-1;}
	while(ch>='0'&&ch<='9') {x=(x<<3)+(x<<1)+(ch^48);ch=getchar();}
	return x*w;
}
inline void write(int x)
{
	if(x<0) {putchar('-');x=~(x-1);}
	if(x>9) write(x/10);
	putchar('0'+x%10);
}
const int N=11;
int n,mod;
int fmul(int x,int y,int mod)
{
	unsigned int res=
	(unsigned int)x*y-
	(unsigned int)((long double)x/mod*y+0.5L)*mod;
	return (res+mod)%mod;
}
//int fmul(int x,int y,int mod)
//{
//	int res=0;
//	while(y)
//	{
//		if(y&1) res+=x,res%=mod;
//		y>>=1;
//		x<<=1;
//		x%=mod;
//	}
//	return res;
//}
struct Mat{
	int a[N][N];
	Mat(){
		memset(a,0,sizeof a);
	}
	Mat operator*(const Mat &x)const{
	    Mat res;
	    for(int i=0;i<N;++i)
	      for(int j=0;j<N;++j)
	        for(int k=0;k<N;++k)
	          res.a[i][j]+=fmul(a[i][k],x.a[k][j],mod),res.a[i][j]%=mod;
	    return res;
	}
}f,A;
Mat fpow(Mat x,int y)
{
	Mat res=x;y--;
	while(y)
	{
		if(y&1) res=res*x;
		y>>=1;
		x=x*x;
	}
	return res;
}
signed main()
{
    int p,q,r,t,u,v,w,x,y,z;
    n=read();mod=read();
    cin>>p>>q>>r>>t>>u>>v>>w>>x>>y>>z;
    f.a[0][0]=f.a[0][1]=f.a[0][2]=f.a[0][8]=f.a[0][9]=f.a[0][10]=1;
    f.a[0][3]=f.a[0][4]=f.a[0][5]=3;
    f.a[0][6]=w;f.a[0][7]=z;
    A.a[0][3]=q;
    A.a[1][4]=v;
    A.a[2][5]=y;
    A.a[3][0]=1;A.a[3][3]=p;A.a[3][4]=A.a[3][5]=1;
    A.a[4][1]=A.a[4][3]=A.a[4][5]=1;A.a[4][4]=u;
    A.a[5][2]=A.a[5][3]=A.a[5][4]=1;A.a[5][5]=x;
    A.a[6][4]=1;A.a[6][6]=w;
    A.a[7][5]=1;A.a[7][7]=z;
    A.a[8][3]=r;A.a[8][8]=1;
    A.a[9][3]=t;A.a[9][5]=A.a[9][9]=1;A.a[9][8]=2;
    A.a[10][3]=A.a[10][8]=A.a[10][9]=A.a[10][10]=1;A.a[10][5]=2;
    f=f*fpow(A,n-1);
	printf("nodgd %lld\nCiocio %lld\nNicole %lld", f.a[0][0], f.a[0][1], f.a[0][2]);
	return 0;
}