万能欧几里得学习笔记
博客还是给自己看的，所以不保证看这篇博客能看懂。
主要参考 https://www.luogu.com.cn/blog/ix-35/solution-p5170 orz，但这篇博客略有一点笔误。
万能欧几里得之所以有这个名字，是因为它的复杂度分析和 \(\gcd\) 类似，而适用面很广。比类欧几里得算法不知道高到哪里去了
主要用于在 \(O(\log \max(P,Q))\) 的时间内求解带有 $\lfloor {Px+R\over Q} \rfloor $ 的有某些性质的式子，完全包括了类欧几里得算法的适用方面，且容易拓展，式子短，代码也短。
思想
考虑一条直线 \(y={Px+R\over Q},0\le R<Q\) ，在第一象限中会与 \(x=k\) 和 \(y=h\) 有若干个交点。每次经过 \(x=k\) 时往操作序列中加入 \(R\) ，经过 \(y=h\) 时加入 \(U\) 。特别地，如果经过了一个整点，那么先加 \(U\) ，再加 \(R\) 。
一个操作可能是对某些变量进行修改，只要能合并且满足结合律就行。
情况1：\(P\ge Q\)
容易发现此时 \(x\) 每增大 1 ，至少会碰到 \(\lfloor {P\over Q} \rfloor\) 个横线。
所以我们可以把这 \(\lfloor {P\over Q} \rfloor\) 个 \(U\) 和一个 \(R\) 合并起来，即递归 \(solve(P\%Q,Q,R,L,U,U^{\lfloor {P\over Q} \rfloor}\times R)\)。
情况2：\(P<Q\)
为了能够继续操作它，我们尝试把坐标系按 \(y=x\) 对称一下，于是直线变成了 \(y=\frac{Qx-R}{P}\) ，且我们要执行 swap(U,R) 。
然而此时规则变成了经过整点时先加 \(R\) 再加 \(U\) ，这不满足递归的条件，所以要把直线向下平移一点点，变成 \(y=\frac{Qx-R-1}{P}\) 。
但是此时还是不能递归，因为 \(R'=-R-1\) 不满足 \(0\le R'<P\) 。
但我们知道 \(y(1)\ge 0\) ，所以可以把 \(x=1\) 之前的操作序列拎出来，然后把直线平移一下。\(x=1\) 之前有 \(\lfloor{Q-R-1\over P}\rfloor\) 个 \(U\) 。
平移之后直线变成了 \(y=\frac{Qx+(Q-R-1)\% P}{P}\)，可以递归……了吗？
注意还有一个问题：原来的上界是 \(x\le L\) ，而翻转坐标系之后变成了 \(y\le L\) ，这也是不满足递归条件的。我们可以令 \(m={LP+R\over Q}\) ，那么递归的上界变为 \(m-1\) （这是因为把 \(x=1\) 之前的删掉了），然后再把最后剩下来的 \(L-\lfloor {mQ-R-1\over P}\rfloor\) 个 \(U\) 给加上去。
于是两种情况都讨论完了，容易发现除了快速幂以外复杂度分析和 \(\gcd\) 是一样的，而快速幂的复杂度分析我不会 /kk
例题
然而这东西究竟怎么解决要做的题呢？不妨结合一道例题来分析。
求
其中 \(1\le n,P,Q\le 10^9\) 。
（没错，这就是类欧几里得算法例题）
那么可以维护三元组 \((cntU,cntR,ans)\) ，表示已经经过的 \(U,R\) 分别有几个，以及当前的答案是多少。
于是走过 \(U\) 的时候变成 \((cntU+1,cntR,ans)\) ，走过 \(R\) 的时候变成 \((cntU,cntR+1,ans+cntU)\) 。
合并三元组 \((a,b,c),(d,e,f)\) 变为 \((a+d,b+e,c+f+ae)\) 。
于是就解决了这道题。
代码
程序的主体部分长这样：
hh ksm(hh x,ll y){hh ret;for (;y;y>>=1,x=x*x) if (y&1) ret=ret*x;return ret;}
hh work(ll P,ll Q,ll R,ll n,hh A,hh B)
{
if (!(((lll)n*P+R)/Q)) return ksm(B,n);
if (P>=Q) return work(P%Q,Q,R,n,A,ksm(A,P/Q)*B);
hh res;
swap(A,B);
res=ksm(A,(Q-R-1)/P)*B;
ll m=((lll)n*P+R)/Q;
res=res*work(Q,P,(Q-R-1)%P,m-1,A,B);
res=res*ksm(A,n-((lll)m*Q-R-1)/P);
return res;
}
其中 hh 表示一个操作序列。
这个主体部分是基本不需要修改的，至少做洛谷类欧模板和 loj 万欧模板时一模一样。
类欧模板（洛谷5170）：
struct hh
{
ll cntA,cntB;
ll sum,sq,pr,sx;
hh(ll x=0,ll y=0,ll z=0,ll a=0,ll b=0,ll c=0){cntA=x,cntB=y,sum=z,sq=a,pr=b,sx=c;}
const hh operator * (const hh &t) const
{
ll x=(cntA+t.cntA)%mod,y=cntB+t.cntB,z=(sum+t.sum+cntA*t.cntB)%mod;
ll a=(sq+t.sq+2*cntA*t.sum%mod+t.cntB*cntA%mod*cntA%mod)%mod;
ll b=(pr+t.pr+cntB*t.sum%mod+cntA*t.sx%mod+t.cntB*cntB%mod*cntA%mod)%mod;
ll c=(sx+t.sx+t.cntB*cntB%mod)%mod;
return hh(x,y,z,a,b,c);
}
};
void work()
{
ll n,a,b,c; read(n,a,b,c);
hh A=hh(1,0,0,0,0,0),B=hh(0,1,0,0,0,1);
hh res=work(a,c,b%c,n,A,B);
res=hh(b/c,0,0,0,0,0)*res;
printf("%lld %lld %lld\n",(res.sum+b/c)%mod,(res.sq+(b/c)*(b/c))%mod,res.pr);
}
万欧模板（loj6440）：
int n;
struct Mat
{
ll a[23][23];
Mat(){memset(a,0,sizeof(a));}
const Mat operator * (const Mat &x) const
{
Mat ret;
rep(i,1,n) rep(k,1,n) rep(j,1,n) (ret.a[i][j]+=a[i][k]*x.a[k][j]%mod)%=mod;
return ret;
}
const Mat operator + (const Mat &x) const
{
Mat ret;
rep(i,1,n) rep(j,1,n) ret.a[i][j]=(a[i][j]+x.a[i][j])%mod;
return ret;
}
}x,y,I;
Mat ksm(Mat x,ll y){Mat ret=I;for (;y;y>>=1,x=x*x) if (y&1) ret=ret*x;return ret;}
struct hh
{
Mat A,B,ans;
hh(Mat a=I,Mat b=I,Mat c=Mat()){A=a,B=b,ans=c;}
const hh operator * (const hh &x) const {return hh(A*x.A,B*x.B,ans+A*x.ans*B);}
};
int main()
{
file();
ll P,Q,R,L;
read(P,Q,R,L,n);
rep(i,1,n) rep(j,1,n) read(x.a[i][j]);
rep(i,1,n) rep(j,1,n) read(y.a[i][j]);
rep(i,1,n) I.a[i][i]=1;
hh A=hh(I,y,Mat()),B=hh(x,I,x);
hh res=hh(I,ksm(y,R/Q),Mat());
res=res*work(P,Q,R%Q,L,A,B);
rep(i,1,n) rep(j,1,n) printf("%lld%c",res.ans.a[i][j]," \n"[j==n]);
return 0;
}
