类欧几里得入土总结 2
直接说 LOJ万能欧几里得这道题
题目要我们求这个式子
问题转换
考虑一件这样的事情,平面上一条直线(ax+bover c),关于这条直线,我们放到网格图上,从左到右,我们遇到一条横线执行一次(U)操作,遇到一次竖线执行一次(R)操作,经过一个整点,执行UR操作。U代表迭代一次( m curB*=B),R代表迭代一次( m curA*=A)然后累加一遍({ m sum+=}( m curA) imes ({ m curB})),这相当于在模拟这个(sum)的操作,我们称(UR)构成的序列叫做操作序列。
若此处的(<+, imes>)有两者有结合律,(+)要对$ imes $有分配律,x,y是操作序列,xy表示直接首尾相接起来。我们容易证明以下等式
也就是说那么操作序列具有结合律
记五元组((n,a,b,c,x,y)=s)表示我要生成一个操作序列(s),规则是对于直线(f(x)=y=lfloor {ax+bover c} floor,xin [1,n]) ,放在网格上,从左下往右上走,每经过一个竖线添加一段操作数列(y),每经过一个横线添加一段操作序列(x)(同时经过先进行(x)操作)。
注意到这个(s)总共有(n)个(y),第(i)个(y)之前有({ai+bover c})个(x)。
沿用【瞎讲】类欧几里得入土教程的思路,我们考虑将(a,b)都化到小于(c)的情况。
一些铺垫
重要补充
设(g(i))表示在第(i)个(y)之前有多少个(x),我们有
定义(f)的和定义和(g)等价,描述的是同一个情况,为了方便理解建议用差分的角度理解(i到i-1个y之间有多少个x)
先考虑(bge c)的情况,我们要使得$b o b%c $
第(i>1)个(y)之前有(f(i)-f(i-1))个(x),而(b o b\% c)对这个差分毫无影响。
然而(i=1)的时候有小问题,第一个(y)之前个(x)我们不是差分定义的,而我们使得(b o b\%c),使得在第一个(y)之前少了一些(y),具体的,少了(lfloor{a imes 1+bover c} floor -lfloor {a imes 1+b\% cover c} floor=lfloor{bover c} floor)个(y)。
根据以上的分析,我们得到了
再考虑(age c)的情况,我们要使得(a o a\% c)
此时(f(x)=lfloor{(a\%c )x+bover c} floor+lfloor{aover c} floor x),设(f'(x)=lfloor{(a\%c )x+bover c} floor),(g')是(f')的差分,我们发现(g(i)=lfloor {aover c} floor+g'(i)),即任何一个(y)之前都紧贴着(lfloor{aover c} floor)个(x)
因此
解决问题
继续沿用上面那个博客的思路,到了真正解决问题的时刻了
考虑现在我们只需要解决的问题是(a,b<c)
边界情况:
关于(a=0)的情况退化为一条平行于(x)的直线,就是填(n)个(y)
而(b=0)无所谓
(a,b<c)意味着一件事情,这个直线非常平缓,甚至填了若干个(y)才填一个(x)。
这也意味着另一件事情,那就是一个(x)之前有若干个(y)。
考虑 第(j>1)个(x)之前有多少个(y),即考虑(maxlimits_{m} lfloor {am+bover c} floor< j)。这即(m=f'(j),j>1)
解出来
也就是从第二个(x)开始到第(lfloor{an+bover c} floor)个(x)是一个子问题(要记得令(j=j-1)并把(c)放出来,即)
而第一个(x)之前还有(f''(0))个(y),这一部分补上。
还有一个问题是,最后一个(x)之后可能还有(y)没有统计到,我们手动算出来补在后面即可。
根据以上的分析,我们得到了
总结
先把公式全写出来
然而求很多操作序列的叠加要用分治,是一个(log)的,但是发现每次分治的长度是按除法递减的,仔细分析一波发现就是(log n)。(这里地方太小我写不下)
实现时注意每次调用递归时,要注意将第一个(y)之间的(x)补上($ m case::otherwise $ 的补(y)也可以理解为这个)。
这种写法的一大好处是,递归部分完全不需要改动,需要改动的只有操作序列的合并方式。如果操作可以对应到矩阵那就完全不要改动了。
这种写法的另一大好处是,任何良好定义的(<op1(+),op2( imes )>),且curB,curA的"迭代"有结合律和交换律,都可以套用这个做法。(泰拳警告:求多项式,超现实数,Floyd数组,SG函数,图连通性的sum)
如果(sum)后面是(k)个项相乘(但是有交换律)也能解决(我愿称之为k阶-类欧几里得)
代码实现上,由于(case1)和(case3)的统一性,可以每次在递归调用前直接补上$lfloor {bover c} floor $ ( (case2)除外),那么每次进入一个递归马上可以(b o b\%c)。
此外,注意到当(f(1)=0)的时候特殊处理下。
目前LOJ rk4 不知道为啥我跑得这么快qwq
//@winlere
#include<iostream>
#include<algorithm>
#include<cstdio>
#include<cstring>
#define mod 998244353
typedef long long ll;
ll qr(){
ll c=getchar(),ret=0,f=0;
while(!isdigit(c)) f|=c==45,c=getchar();
while( isdigit(c)) ret=ret*10+c-48,c=getchar();
return ret;
}
int MOD(const int&x){return x>=mod?x-mod:x;}
int MOD(const int&x,const int&y){ll g=1ll*x*y;return g-g/mod*mod;}
int N;
struct EU{
struct MAT{
int v[20][20];
MAT(){memset(v,0,sizeof(int)*400);}
int*operator[](int x){return v[x];}
MAT operator + (const MAT&x)const{
MAT ret;
for(int t=0;t<N;++t)
for(int i=0;i<N;++i)
ret[t][i]=MOD(v[t][i]+x.v[t][i]);
return ret;
}
MAT operator * (const MAT&x)const{
MAT ret;
for(int k=0;k<N;++k)
for(int t=0;t<N;++t)
for(int i=0;i<N;++i)
ret[t][i]=MOD(ret[t][i]+MOD(v[t][k],x.v[k][i]));
return ret;
}
}sum,A,B;
EU(const int&x=1):sum(),A(),B(){for(int t=0;t<N;++t) A[t][t]=B[t][t]=x;}
EU operator + (const EU&x)const{
EU ret;
ret.A=A*x.A;
ret.B=B*x.B;
ret.sum=sum+A*x.sum*B;
return ret;
}
EU operator * (const ll&x)const{
EU ret,base=*this;
for(ll t=x;t>0;t>>=1,base=base+base)
if(t&1) ret=ret+base;
return ret;
}
};
ll f(ll x,ll a,ll b,ll c){return ((__int128)x*a+b)/c;}
EU solve(ll n,ll a,ll b,ll c,const EU&x,const EU&y){
if(!n) return EU();
b%=c;
if(!f(n,a,b,c)) return y*n;
if(a>=c) return solve(n,a%c,b,c,x,x*(a/c)+y);
ll m=f(n,a,b,c),cnt=n-f(m-1,c,c-b-1,a);
return y*((c-b-1)/a)+x+solve(m-1,c,c-b-1,a,y,x)+y*cnt;
}
int main(){
ll a=qr(),c=qr(),b=qr(),n=qr(); N=qr();
EU U,R;
EU::MAT A,B;
for(int t=0;t<N;++t)
for(int i=0;i<N;++i)
A[t][i]=qr();
for(int t=0;t<N;++t)
for(int i=0;i<N;++i)
B[t][i]=qr();
U.B=B; R.A=A; R.sum=A;
EU::MAT ans=(U*(b/c)+solve(n,a,b,c,U,R)).sum;
for(int t=0;t<N;++t,putchar('
'))
for(int i=0;i<N;++i)
printf("%d ",ans[t][i]);
return 0;
}