万能欧几里得算法

## 前言 我本意是想学个类欧就完事了，结果同机房的hehezhou大爷告诉我有个nb的东西叫万能欧几里得。这个算法网上博客不多，讲得也大都不太清楚，平移来平移去搞的我头都大了。我花了好半天终于看明白了，希望尽量清楚地将这个算法介绍给大家。 文章中用$O$代替$\Theta$符号。 ## 引入：普通的类欧 类欧几里得算法用于求解以下这类问题： 求$f(a,b,c,n)=\sum_{i=0}^{n-1}\lfloor\frac{ai+b}{c}\rfloor$ 先分类讨论，当$a\ge c$或$b\ge c$时，令$a'=a\%c,b'=b\%c$ $$f(a,b,c,n)=\sum_{i=0}^{n-1}\lfloor\frac{a'i+b'}{c}\rfloor+\lfloor\frac{ai}{c}\rfloor+\lfloor\frac{b}{c}\rfloor$$ $$f(a,b,c,n)=\frac{n(n-1)}{2}\lfloor\frac{a}{c}\rfloor+n\lfloor\frac{b}{c}\rfloor+f(a',b',c,n)$$ 这样就转化到了$a,b<c$的情况。 考虑增加一维求和： $$f(a,b,c,n)=\sum_{i=0}^{n-1}\sum_{j=0}^{\lfloor\frac{ai+b}{c}\rfloor-1}1$$ 交换求和顺序，令$m=\lfloor\frac{a(n-1)+b}{c}\rfloor$： $$f(a,b,c,n)=\sum_{j=0}^{m-1}\sum_{i=0}^{n-1}[j<\lfloor\frac{ai+b}{c}\rfloor]$$ $$f(a,b,c,n)=\sum_{j=0}^{m-1}\sum_{i=0}^{n-1}[c(j+1)\le ai+b]$$ $$f(a,b,c,n)=\sum_{j=0}^{m-1}\sum_{i=0}^{n-1}[i>\lfloor\frac{cj+c-b-1}{a}\rfloor]$$ $$f(a,b,c,n)=\sum_{j=0}^{m-1}n-1-\lfloor\frac{cj+c-b-1}{a}\rfloor$$ $$f(a,b,c,n)=(n-1)m-f(c,c-b-1,a,m)$$ 可以递归。由于$f(a,b,c,n)$中$a$可以对$c$取模，这样反复递归，类比欧几里得算法可以得到时间复杂度为$O(\log\max(a,c))$。 从几何意义考虑，$\sum_{i=0}^{n-1}\lfloor\frac{ai+b}{c}\rfloor$对应于$[0,n)$范围内，直线$y=\frac{ax+b}{c}$下的整点数。而用类欧几里得算法求解的过程，则相当于每次将直线的斜率和截距转化到$[0,1)$再翻转坐标系。 类似地，设 $$g(a,b,c,n)=\sum_{i=0}^{n-1}i\lfloor\frac{ai+b}{c}\rfloor$$ $$h(a,b,c,n)=\sum_{i=0}^{n-1}\lfloor\frac{ai+b}{c}\rfloor^2$$ 可以推得： $a\ge c$或$b\ge c$时， $$g(a,b,c,n)=g(a',b',c,n)+\frac{n(n-1)(2n-1)}{6}\lfloor\frac{a}{c}\rfloor+\frac{n(n-1)}{2}\lfloor\frac{b}{c}\rfloor$$ $$h(a,b,c,n)=h(a',b',c,n)+2\lfloor\frac{a}{c}\rfloor g(a',b',c,n)+2\lfloor\frac{b}{c}\rfloor f(a',b',c,n)+\frac{n(n-1)(2n-1)}{6}\lfloor\frac{a}{c}\rfloor^2+n(n-1)\lfloor\frac{a}{c}\rfloor\lfloor\frac{b}{c}\rfloor+n\lfloor\frac{b}{c}\rfloor^2$$ $a<c$且$b<c$时， $$g(a,b,c,n)=\frac{1}{2}(mn(n-1)-h(c,c-b-1,a,m)-f(c,c-b-1,a,m))$$ $$h(a,b,c,n)=(n-1)m^2-2g(c,c-b-1,a,m)-f(c,c-b-1,a,m)$$ 将三个函数的值一起递归，可以保证复杂度正确。 code（[洛谷模板](https://www.luogu.com.cn/problem/P5170)） ``` #include<bits/stdc++.h> using namespace std; const int MOD=998244353,inv2=499122177,inv6=166374059; int n,a,b,c; struct node{int f,g,h;}ans; int mo(int x){return x>=MOD?x-MOD:x;} int sqr(int x){return 1ll*x*x%MOD;} node asgcd(int a,int b,int c,int n) { node res; if(!a) { res.f=1ll*b/c*n%MOD; res.g=1ll*n*(n-1)/2%MOD*(b/c)%MOD; res.h=1ll*sqr(b/c)*n%MOD; } else if(a>=c||b>=c) { node o=asgcd(a%c,b%c,c,n); res.f=(o.f+1ll*n*(n-1)/2%MOD*(a/c)+1ll*b/c*n)%MOD; res.g=(o.g+1ll*a/c*n%MOD*(n-1)%MOD*(n*2-1)%MOD*inv6+1ll*n*(n-1)/2%MOD*(b/c))%MOD; res.h=(o.h+2ll*(a/c)*o.g+2ll*(b/c)*o.f+1ll*n*(n-1)%MOD*(n*2-1)%MOD*inv6%MOD*sqr(a/c)+1ll*n*(n-1)%MOD*(a/c)%MOD*(b/c)+1ll*n*sqr(b/c))%MOD; } else { int m=(1ll*a*(n-1)+b)/c; node o=asgcd(c,c-b-1,a,m); res.f=(1ll*(n-1)*m-o.f+MOD)%MOD; res.g=((1ll*m*n%MOD*(n-1)-o.h-o.f)%MOD+MOD)*inv2%MOD; res.h=(1ll*(n-1)*m%MOD*m-mo(o.g*2)-o.f+MOD+MOD)%MOD; } return res; } void solve() { scanf("%d%d%d%d",&n,&a,&b,&c); ans=asgcd(a,b,c,n+1); printf("%d %d %d\n",ans.f,ans.h,ans.g); } int main() { int t; scanf("%d",&t); while(t--)solve(); } ``` 现在你已经会做洛谷模板了，但是别急，下面还要介绍更厉害的科技。 ## 万能欧几里得 万能欧几里得可以解决几乎全部类欧几里得算法能解决的问题。其优势在于，不管要求什么柿子的值，代码都类似，不需要重新推柿子；而且只处理两个信息的合并显然比类欧处理和式友善地多。~~方便我这样的蒟蒻背板子、推柿子~~ 一般性问题描述：现在有一个空字符串$S$，在平面直角坐标系中有一条不包含左端点的线段（不与y轴平行），在其定义域内，从左向右，直线每碰到一次直线$y=a,a\in\mathbb{Z}$就在$S$末尾写一个$U$，每碰到一次直线$x=a,a\in\mathbb{Z}$就写一个$R$，碰到整点就先写$U$再写$R$。再给定一个字符串函数$f$，求$f(S)$。 例题：求$ans=\sum_{i=1}^ni\lfloor\frac{ai+b}{c}\rfloor$ 解法： 首先，我们构造例题中的柿子对应的$f$：遍历$S$，维护一个变量$cnt$初始$=0$,若当前位置是$U$，则$cnt$++，否则$f(S)$+=$i\times cnt$，其中$i$表示当前位置是第$i$个$R$。 可以看出，要求的$ans=f(\lfloor\frac{b}{c}\rfloor个U+S)$，其中串$S$由线段$y=\frac{ax+b}{c},x\in(0,n]$生成。 先定义一个运算$\times$，满足$f(S_1+S_2)=f(S_1)\times f(S_2)$。对于例题中的问题，设$S_1$中的$R$数量为$x$，$U$数量为$y$，则$f(S_1+S_2)=f(S_1)+\sum_{i=1}^{n_2}(i+x)(\lfloor\frac{a_2i+b_2}{c_2}\rfloor+y)=f(S_1)+f(S_2)+x\sum_{i=1}^{n_2}\lfloor\frac{a_2i+b_2}{c_2}\rfloor+\frac{n_2(n_2-1)}{2}y+n_2xy$。从这个柿子我们可以看出，为了求$f(S)$，我们需要维护$S$的$x,y$和另一个函数$\sum_{i=1}^n\lfloor\frac{ai+b}{c}\rfloor$的值（类似地，这个函数的值也可以由两部分字符串的信息合并得到）。为了方便，把这些值和$f$都打包在结构体里，并对结构体重载$\times$运算符。扩展$f$的意义，将其定义为$S$对应的那个结构体。这样，我们定义了运算符$\times$。 有了$\times$运算，就可以定义$f^k=\begin{cases}f^{k-1}\times f&k>0\\f(空串)&k=0\end{cases}$ 利用快速幂算法，已知$f,k$，可以在$O(\log k)$时间中求出$f^k$。 现在我们考虑如何解决问题。在线段定义域为$[0,n)$时，$f(S)$仅与$a,b,c,f(U),f(R)$有关，不妨将其记为$f(a,b,c,n,f_U,f_R)$。考虑如何求它。由于直线平移整数格不改变$S$，也就不改变$f(S)$，因此$b$可以任意对$c$取模，可以保证$b\in[0,c)$。对$a$与$c$的关系分类讨论。 当$a\ge c$时，直线$y=\frac{ax+b}{c}$相比直线$\frac{(a\%c)x+b}{c}$，每个R之前会多$\lfloor\frac{a}{c}\rfloor$个U。因此可以得到： $$f(a,b,c,n,f_1,f_2)=f(a\%c,b,c,n,f_U,f_U^{\lfloor\frac{a}{c}\rfloor}\times f_2)$$ 当$a<c$时，情况会比较麻烦。  上图中，线段$AB$对应$\sum_{i=1}^7i\lfloor\frac{2i+4}{7}\rfloor$。上图仅仅是为了方便理解，我们要讨论的是一般情况。 令$m=\lfloor\frac{an+b}{c}\rfloor$。 若$m=0$，则代表串$S$中全是$R$，出现递归出口，直接得到$f(a,b,c,n,f_U,f_R)=f_R^n$。 否则，考虑翻转坐标系。既然原本斜率$<1$，翻转后斜率就会$\ge 1$，可以进一步递归。将$AB$关于直线$y=x$轴对称得到$A'B'$。轴对称后，整点先$U$后$R$的顺序变为先$R$后$U$，为了保持先$U$后$R$的顺序，还需要将$A'B'$向右平移$\frac{1}{c}$个单位，得到$A''B''$，所在直线为$y=\frac{cx-b-1}{a}$。那么经过一顿变换，$f(S)$的值等于交换$f(U),f(R)$后，线段$A''B''$所生成的字符串的$f$值。 接下来，将$A''B''$分成三段计算：横坐标$\le1$的一段，横坐标$>m$的一段，中间的一段，分别对应图中的$B''D,A''C,CD$。 先考虑中间段。由于$f(a,b,c,n,f_U,f_R)$的定义是区间$(0,n]$，需要将$CD$左移1单位。平移后$CD$所在直线为$y=\frac{c(x+1)-b-1}{a}$，$f$值对应于$f(c,c-b-1,a,m-1,f_R,f_U)$。 再考虑两头。由于$D,C$的坐标分别为$(1,\frac{c-b-1}{a}),(m,\frac{cm-b-1}{a})$，则$B''D$段$f$值对应于$f_R^{\lfloor\frac{c-b-1}{a}\rfloor}\times f_U$；$CA''$段$f$值对应于$f_R^{n-\lfloor\frac{cm-b-1}{a}\rfloor}$。 最后，我们得到： $$f(a,b,c,n,f_U,f_R)=f_R^{\lfloor\frac{c-b-1}{a}\rfloor}\times f_U\times f(c,c-b-1,a,m-1,f_R,f_U)\times f_R^{n-\lfloor\frac{cm-b-1}{a}\rfloor}$$ 接下来分析一下它的时间复杂度。显然，这只与$a,c$有关。我们假定运算符$\times$的时间复杂度为$O(1)$。 设$T(a,c)$是$f(a,b,c,n,f_U,f_R)$的时间复杂度，我们可以得到： $$T(a,c)=T(c\%a,a)+O(\log(\frac{c}{a})),a<c$$ 由于$\log(\frac{c}{a})=\log c-\log a$，将$T(c\%a,a)$的时间复杂度也递归，可以发现$\log c-\log a$与$\log a-\log (c\%a)$抵消了$\log c$项。一直递归到底，一顿抵消后可以得到：$T(a,c)=O(\log c)-O(\log \gcd(a,c))=O(\log c)$。 ## 例题 1.[洛谷P5170](https://www.luogu.com.cn/problem/P5170) 没啥可说的，纯板子，可以对照代码理解一下算法的实现。 ``` #include<bits/stdc++.h> using namespace std; const int MOD=998244353; int mo(int x){return x>=MOD?x-MOD:x;} int sqr(int x){return 1ll*x*x%MOD;} int n,a,b,c; struct node{int f,g,h,u,r;}ans,f0,fu,fr; node operator *(node a,node b) { node c; c.u=mo(a.u+b.u); c.r=mo(a.r+b.r); c.f=(1ll*a.u*b.r+a.f+b.f)%MOD; c.g=(1ll*a.r*b.f+1ll*b.r*(b.r-1)/2%MOD*a.u+1ll*b.r*a.r%MOD*a.u+a.g+b.g)%MOD; c.h=(2ll*a.u*b.f+1ll*a.u*a.u%MOD*b.r+a.h+b.h)%MOD; return c; } node operator ^(node a,int b) { node res=f0; for(;b;b>>=1,a=a*a) if(b&1)res=res*a; return res; } node asgcd(int a,int b,int c,int n,node fu,node fr) { b%=c; if(a>=c)return asgcd(a%c,b,c,n,fu,(fu^a/c)*fr); else { int m=(1ll*a*n+b)/c; if(!m)return fr^n; return (fr^(c-b-1)/a)*fu*asgcd(c,c-b-1,a,m-1,fr,fu)*(fr^n-(1ll*c*m-b-1)/a); } } void solve() { scanf("%d%d%d%d",&n,&a,&b,&c); ans=(fu^b/c)*fr*asgcd(a,b,c,n,fu,fr); printf("%d %d %d\n",ans.f,ans.h,ans.g); } int main() { f0={0,0,0,0,0}; fu={0,0,0,1,0}; fr={0,0,0,0,1}; int t; scanf("%d",&t); while(t--)solve(); } ``` 2.[LOJ#138](https://loj.ac/problem/138) 还是板子，只不过要加上一点套路的推柿子。。 考虑如何合并信息： 设要合并串$a,b$为串$c$，$a$中$U,R$的个数分别记为$a.u,a.r$，$a$中第$i$个$R$前$U$的数量记为$a.u(i)$， $$a.f(k_1,k_2)=\sum_{i=1}^{a.r}(i-1)^{k_1}a.u^{k_2}(i)$$ 那么，合并信息过程中： $$c.f(k_1,k_2)=a.f(k_1,k_2)+\sum_{i=1}^{b.r}(i-1+a.r)^{k_1}(b.u(i)+a.u)^{k_2}$$ $$c.f(k_1,k_2)=a.f(k_1,k_2)+\sum_{i=1}^{b.r}(\sum_{j=0}^{k_1}\binom{k_1}{j}(i-1)^{j}a.r^{k_1-j})(\sum_{j=0}^{k_2}\binom{k_2}{j}b.u^j(i)a.u^{k_2-j})$$ $$c.f(k_1,k_2)=a.f(k_1,k_2)+\sum_{i=0}^{k_1}\sum_{j=0}^{k_2}\binom{k_1}{i}\binom{k_2}{j}a.r^{k_1-i}a.u^{k_2-j}\sum_{k=1}^{b.r}(k-1)^{i}b.u^j(k)$$ $$c.f(k_1,k_2)=a.f(k_1,k_2)+\sum_{i=0}^{k_1}\sum_{j=0}^{k_2}\binom{k_1}{i}\binom{k_2}{j}a.r^{k_1-i}a.u^{k_2-j}b.f(i,j)$$ 单次合并时间复杂度为$O(k_1^2k_2^2)$，因此单次询问复杂度为$O(k_1^2k_2^2\log\max(a,c))$。 这题充分体现了万能欧几里得的优势：只要实现快速合并信息就行。相比普通的类欧，思维和代码难度都大大减小。 3.[LOJ6440](https://loj.ac/problem/6440) 对于一个字符串$S$，维护它的$\sum_{i=1}^nA^xB^{\lfloor\frac{ai+b}{c}\rfloor}$和$A^xB^{\lfloor\frac{ai+b}{c}\rfloor}$即可。 ## 参考资料 https://oi-wiki.org/math/euclidean/ https://www.cnblogs.com/AThousandMoons/p/13129093.html https://www.cnblogs.com/Mr-Spade/p/10370259.html