万能欧几里得算法

前言

我本意是想学个类欧就完事了，结果同机房的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

可以推得：

$$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

没啥可说的，纯板子，可以对照代码理解一下算法的实现。

#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

还是板子，只不过要加上一点套路的推柿子。。

考虑如何合并信息：

设要合并串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

对于一个字符串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