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