CVBB ACM Team

到此差别页面的链接

--- 2020-2021:teams:legal_string:jxm2001:万能欧几里得算法 [2021/08/19 22:22]
jxm2001 [算法原理]
+++ 2020-2021:teams:legal_string:jxm2001:万能欧几里得算法 [2021/08/22 15:16] (当前版本)
jxm2001 [例题二]
@@ 行 62: / 行 62: @@
 最终 $T(a,c)=O(\log c)-O(\log\gcd (a,c))\sim O(\log c)$。
+===== 算法例题 =====
+==== 例题一 ====
+[[https://www.luogu.com.cn/problem/P5170|洛谷p5170]]
+=== 题意 ===
+求
+$$
+f(n)=\sum_{i=0}^n \lfloor \frac {ai+b}c\rfloor\\
+g(n)=\sum_{i=0}^n (\lfloor \frac {ai+b}c\rfloor)^2\\
+h(n)=\sum_{i=0}^n i\lfloor \frac {ai+b}c\rfloor
+$$
+=== 题解 ===
+设 $F,G,H$ 分别表示序列 $S$ 对应 $f,g,h$ 的答案，于是有
+$$
+\begin{equation}\begin{split}
+F(S1S2)&=F(S1)+\sum_{i=1}^n (\lfloor \frac {ai+b}c\rfloor+dy)=F(S1)+F(S2)+n\times dy\\
+G(S1S2)&=G(S1)+\sum_{i=1}^n (\lfloor \frac {ai+b}c\rfloor+dy)^2 \\
+&=G(S1)+\sum_{i=1}^n (\lfloor \frac {ai+b}c\rfloor)^2+
+dy\sum_{i=1}^n \lfloor \frac {ai+b}c\rfloor+dy+\sum_{i=1}^n (dy)^2\\
+&=G(S1)+G(S2)+2dyF(S2)+n\times (dy)^2\\
+H(S1S2)&=H(S1)+\sum_{i=1}^n (i+dx)(\lfloor \frac {ai+b}c\rfloor+dy) \\
+&=H(S1)+\sum_{i=1}^n \left(i(\lfloor \frac {ai+b}c\rfloor)+i\times dy+dx(\lfloor \frac {ai+b}c\rfloor)+dxdy\right)\\
+&=H(S1)+H(S2)+\frac {n(n+1)}2dy+dxF(S1)+n\times dxdy
+\end{split}\end{equation}
+$$
+最后这题的 $i$ 是从 $0$ 开始的，万能欧几里得算法的 $i$ 是从 $1$ 开始的，注意补上 $i=0$ 的贡献。
+<hidden 查看代码>
+<code cpp>
+const int mod=998244353;
+struct Node{
+	LL cntr,cntu,f,g,h;
+	Node operator * (const Node &b)const{
+		Node c;
+		LL n=b.cntr,dx=cntr,dy=cntu;
+		c.cntr=(cntr+b.cntr)%mod;
+		c.cntu=(cntu+b.cntu)%mod;
+		c.f=(f+b.f+dy*n)%mod;
+		c.g=(g+b.g+dy*b.f*2+dy*dy%mod*n)%mod;
+		c.h=(h+b.h+n*(n+1)/2%mod*dy+dx*b.f+dx*dy%mod*n)%mod;
+		return c;
+	}
+};
+Node quick_pow(Node n,int k){
+	Node ans=Node{0,0,0,0,0};
+	while(k){
+		if(k&1)ans=ans*n;
+		n=n*n;
+		k>>=1;
+	}
+	return ans;
+}
+Node asgcd(int a,int b,int c,int n,Node su,Node sr){
+	if(a>=c)
+	return asgcd(a%c,b,c,n,su,quick_pow(su,a/c)*sr);
+	int m=(1LL*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-(1LL*c*m-b-1)/a);
+}
+Node cal(int a,int b,int c,int n){
+	Node su=Node{0,1,0,0,0},sr=Node{1,0,0,0,0};
+	return quick_pow(su,b/c)*asgcd(a,b%c,c,n,su,sr);
+}
+int main()
+{
+	int T=read_int();
+	while(T--){
+		int n=read_int(),a=read_int(),b=read_int(),c=read_int();
+		Node ans=cal(a,b,c,n);
+		ans.f=(ans.f+b/c)%mod;
+		ans.g=(ans.g+1LL*(b/c)*(b/c))%mod;
+		space(ans.f);space(ans.g);enter(ans.h);
+	}
+	return 0;
+}
+</code>
+</hidden>
+==== 例题二 ====
+[[https://loj.ac/p/138|LOJ#138]]
+=== 题意 ===
+求
+$$
+f(n)=\sum_{i=0}^n i^{k_1}(\lfloor \frac {ai+b}c\rfloor)^{k_2}
+$$
+=== 题解 ===
+设 $F(S,k_1,k_2)$ 表示 $f(n)=\sum_{i=1}^n i^{k_1}(\lfloor \frac {ai+b}c\rfloor)^{k_2}$ 关于 $S$ 操作序列的答案。
+\begin{equation}\begin{split}
+F(S1S2,k_1,k_2)&=F(S1,k_1,k_2)+\sum_{i=1}^n (i^{k_1}+dx)^{k_1}(\lfloor \frac {ai+b}c\rfloor+dy)^{k_2}\\
+&=F(S1,k_1,k_2)+\sum_{i=1}^n\sum_{t_1=0}^{k_1}\sum_{t_2=0}^{k_2}{k_1\choose t_1}{k_2\choose t_2}dx^{k_1-t_1}dy^{k_2-t_2}i^{t_1}(\lfloor \frac {ai+b}c\rfloor+dy)^{t_2}\\
+&=F(S1,k_1,k_2)+\sum_{t_1=0}^{k_1}\sum_{t_2=0}^{k_2}{k_1\choose t_1}{k_2\choose t_2}dx^{k_1-t_1}dy^{k_2-t_2}F(S2,t_1,t_2)
+\end{split}\end{equation}
+对每个 $S$，维护矩阵 $F(S,i,j)(0\le i\le k_1,0\le j\le k_2)$，于是可以实现 $O\left(k_1^2k_2^2\right)$ 信息合并。
+边界 $F(L,i,j)=0,F(R,i,j)=[j==0](0\le i\le k_1,0\le j\le k_2)$，另外注意对答案补上 $i=0$ 的贡献。时间复杂度 $O\left(k_1^2k_2^2\log c\right)$。
+<hidden 查看代码>
+<code cpp>
+const int mod=1e9+7,MAXK=11;
+int C[MAXK][MAXK],k1,k2;
+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];
+		Node c;
+		int dx=cntr,dy=cntu;
+		px[0]=py[0]=1;
+		_rep(i,1,k1)
+		px[i]=1LL*px[i-1]*dx%mod;
+		_rep(i,1,k2)
+		py[i]=1LL*py[i-1]*dy%mod;
+		c.cntr=(cntr+b.cntr)%mod;
+		c.cntu=(cntu+b.cntu)%mod;
+		_rep(i,0,k1)_rep(j,0,k2){
+			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]*C[i][i2]%mod*C[j][j2]%mod*px[i-i2]%mod*py[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;
+		n=n*n;
+		k>>=1;
+	}
+	return ans;
+}
+Node asgcd(int a,int b,int c,int n,Node su,Node sr){
+	if(a>=c)
+	return asgcd(a%c,b,c,n,su,quick_pow(su,a/c)*sr);
+	int m=(1LL*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-(1LL*c*m-b-1)/a);
+}
+Node cal(int a,int b,int c,int n){
+	Node su=Node(0,1),sr=Node(1,0);
+	_rep(i,0,k1)
+	sr.f[i][0]=1;
+	return quick_pow(su,b/c)*asgcd(a,b%c,c,n,su,sr);
+}
+int main()
+{
+	C[0][0]=1;
+	_for(i,1,MAXK){
+		C[i][0]=1;
+		_rep(j,1,i)
+		C[i][j]=(C[i-1][j-1]+C[i-1][j])%mod;
+	}
+	int T=read_int();
+	while(T--){
+		int n=read_int(),a=read_int(),b=read_int(),c=read_int();
+		k1=read_int(),k2=read_int();
+		Node ans=cal(a,b,c,n);
+		int base=1;
+		_rep(i,0,k2){
+			ans.f[0][i]=(ans.f[0][i]+base)%mod;
+			base=1LL*base*(b/c)%mod;
+		}
+		enter(ans.f[k1][k2]);
+	}
+	return 0;
+}
+</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>
+===== 参考资料 =====
+这两篇还有配图，方便理解。
+[[https://www.luogu.com.cn/blog/ILikeDuck/mo-neng-ou-ji-li-dei-suan-fa|C3H5ClO的博客]]
+[[https://www.luogu.com.cn/blog/ix-35/solution-p5170|ix35_ 的博客]]

CVBB ACM Team

用户工具

站点工具

页面工具