CVBB ACM Team

到此差别页面的链接

--- 2020-2021:teams:legal_string:jxm2001:万能欧几里得算法 [2021/08/20 11:31]
jxm2001
+++ 2020-2021:teams:legal_string:jxm2001:万能欧几里得算法 [2021/08/22 15:16] (当前版本)
jxm2001 [例题二]
@@ 行 160: / 行 160: @@
 $$
-f(n)=\sum_{i=0}^n i^{k_1}(\lfloor \frac {ai+b}c\rfloor)^{k_2}\\
+f(n)=\sum_{i=0}^n i^{k_1}(\lfloor \frac {ai+b}c\rfloor)^{k_2}
 $$
@@ 行 255: / 行 255: @@
 </code>
 </hidden>
+<hidden 卡常版>
+<code cpp>
+int C[MAXK][MAXK];
+struct Node{
+	int cntr,cntu,f[MAXK][MAXK];
+	Node(int cntr=0,int cntu=0){
+		this->cntr=cntr;
+		this->cntu=cntu;
+		mem(f,0);
+	}
+	Node operator * (const Node &b)const{
+		static int px[MAXK],py[MAXK];
+		static int b1[MAXK][MAXK],b2[MAXK][MAXK];
+		Node c;
+		int dx=cntr,dy=cntu;
+		px[0]=py[0]=1;
+		_for(i,1,MAXK)
+		px[i]=1LL*px[i-1]*dx%mod;
+		_for(i,1,MAXK)
+		py[i]=1LL*py[i-1]*dy%mod;
+		_for(i,0,MAXK)_rep(j,0,i){
+			b1[i][j]=1LL*C[i][j]*px[i-j]%mod;
+			b2[i][j]=1LL*C[i][j]*py[i-j]%mod;
+		}
+		c.cntr=(cntr+b.cntr)%mod;
+		c.cntu=(cntu+b.cntu)%mod;
+		_for(i,0,MAXK)_for(j,0,MAXK){
+			c.f[i][j]=f[i][j];
+			_rep(i2,0,i)_rep(j2,0,j)
+			c.f[i][j]=(c.f[i][j]+1LL*b.f[i2][j2]*b1[i][i2]%mod*b2[j][j2])%mod;
+		}
+		return c;
+	}
+};
+Node quick_pow(Node n,int k){
+	Node ans=Node(0,0);
+	while(k){
+		if(k&1)ans=ans*n;
+		k>>=1;
+		if(k)n=n*n;
+	}
+	return ans;
+}
+</code>
+</hidden>
+==== 例题三 ====
+[[https://loj.ac/p/6440|LOJ#6440]]
+=== 题意 ===
+求
+$$
+\sum_{i=1}^n A^iB^{\lfloor \cfrac {ai+b}c\rfloor}
+$$
+其中，$A,B$ 为 $k\times k$ 的矩阵。
+=== 题解 ===
+$$
+\begin{equation}\begin{split}
+F(S1S2)&=F(S1)+\sum_{i=1}^n A^{i+dx}B^{\lfloor \cfrac {ai+b}c\rfloor+dy}\\
+&=F(S1)+A^{dx}\left(\sum_{i=1}^n B^{\lfloor \cfrac {ai+b}c\rfloor}\right)B^{dy}
+&=F(S1)+A^{dx}F(S2)B^{dy}
+\end{split}\end{equation}
+$$
+于是可以令 $F(S)$ 维护 $(A^{cntR},B^{cntU},\text{ans})$，设 $E$ 表示单位矩阵，$Z$ 表示全 $0$ 矩阵。
+边界条件 $F(U)=(E,B,Z),F(R)=(A,E,A),F()=(E,E,Z)$。时间复杂度 $O\left(k^3\log c\right)$。
+<hidden 查看代码>
+<code cpp>
+const int mod=998244353,MAXN=21;
+const __int128 One=1;
+int sz;
+struct Matrix{
+	int a[MAXN][MAXN];
+	Matrix(int type=0){
+		mem(a,0);
+		if(type){
+			_for(i,0,MAXN)
+			a[i][i]=1;
+		}
+	}
+	Matrix operator + (const Matrix &b)const{
+		Matrix c;
+		_for(i,0,sz)_for(j,0,sz)
+		c.a[i][j]=(a[i][j]+b.a[i][j])%mod;
+		return c;
+	}
+	Matrix operator * (const Matrix &b)const{
+		Matrix c;
+		_for(i,0,sz)_for(j,0,sz)_for(k,0,sz)
+		c.a[i][j]=(c.a[i][j]+1LL*a[i][k]*b.a[k][j])%mod;
+		return c;
+	}
+};
+struct Node{
+	Matrix A,B,f;
+	Node operator * (const Node &b)const{
+		Node c;
+		c.A=A*b.A;
+		c.B=B*b.B;
+		c.f=f+A*b.f*B;
+		return c;
+	}
+};
+Node quick_pow(Node n,LL k){
+	Node ans;
+	ans.A=Matrix(1);
+	ans.B=Matrix(1);
+	while(k){
+		if(k&1)ans=ans*n;
+		n=n*n;
+		k>>=1;
+	}
+	return ans;
+}
+Node asgcd(LL a,LL b,LL c,LL n,Node su,Node sr){
+	if(a>=c)
+	return asgcd(a%c,b,c,n,su,quick_pow(su,a/c)*sr);
+	LL m=(One*a*n+b)/c;
+	if(!m)
+	return quick_pow(sr,n);
+	else
+	return quick_pow(sr,(c-b-1)/a)*su*asgcd(c,(c-b-1)%a,a,m-1,sr,su)*quick_pow(sr,n-(One*c*m-b-1)/a);
+}
+Matrix A,B;
+Node cal(LL a,LL b,LL c,LL n){
+	Node su,sr;
+	su.A=Matrix(1);
+	su.B=B;
+	sr.A=A;
+	sr.B=Matrix(1);
+	sr.f=A;
+	return quick_pow(su,b/c)*asgcd(a,b%c,c,n,su,sr);
+}
+int main()
+{
+	LL a=read_LL(),c=read_LL(),b=read_LL(),n=read_LL();
+	sz=read_int();
+	_for(i,0,sz)_for(j,0,sz)
+	A.a[i][j]=read_int();
+	_for(i,0,sz)_for(j,0,sz)
+	B.a[i][j]=read_int();
+	Node ans=cal(a,b,c,n);
+	_for(i,0,sz){
+		_for(j,0,sz)
+		space(ans.f.a[i][j]);
+		putchar('\n');
+	}
+	return 0;
+}
+</code>
+</hidden>
 ===== 参考资料 =====

CVBB ACM Team

用户工具

站点工具

页面工具