万能欧几里得
BJpers2
2020-08-06 15:05:38
类欧几里得在OI当中相当冷门，因为其过于复杂的推理难度和繁琐的实现细节劝退效果实在太强。
但名为「**万能欧几里得**」的算法则可以帮我们绕过这些繁琐的推理步骤，简洁明了地触及问题的本质。
事实上，如果掌握了本算法，那么做一般的类欧几里得题时，你的草稿纸上甚至不必出现任何求和符号。
---
## 一、问题转化
类欧几里得题大概有如下的形式：
> 设$y=\left \lfloor \frac {px+r}q\right\rfloor$
>
求
$$
\sum_{x=1}^n y,\sum_{x=1}^n xy,\sum_{x=1}^n y^2,\sum_{x=1}^n a^xy,\sum_{x=1}^n a^xb^y,\dots
$$
这些问题又可以统一换作如下表述：
> 作出函数$y= \frac {px+r}q$的图像，之后从左往右考虑该线与网格线的交点，维护一个列向量，遇横线则乘矩阵$U$，遇竖线则乘矩阵$R$（同时相遇则优先乘$U$），求刚走过$n$条竖线后的最终向量。
可以发现上述求和问题，均可以通过合理选取维护的向量$\boldsymbol{v}$以及矩阵$U,R$来转化为矩阵乘积问题。
例如求$\sum y$，则可以取$(\boldsymbol{v},U,R)$为：
$$
\begin{bmatrix} \sum_x y \\ y \\ 1\end{bmatrix},
\begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 1 \\ 0 & 0 & 1\end{bmatrix},
\begin{bmatrix} 1 & 1 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1\end{bmatrix}
$$
可以看出这样做矩阵乘法实时地正确维护了左边的三个量。
再如求$\sum a^x y$，我们也可以设计$(\boldsymbol{v},U,R)$如下：
$$
\begin{bmatrix} \sum_x a^xy \\ a^xy \\ a^x \end{bmatrix},
\begin{bmatrix} 1 & 0 & 0\\ 0 & 1 & 1\\ 0 & 0 & 1\end{bmatrix},
\begin{bmatrix} 1 & a & 0\\ 0 & a & 0\\ 0 & 0 & a\end{bmatrix}
$$
又如求$\sum xa^xy^2b^y$，我们照样可以大力写出：
$$
\begin{bmatrix} \sum_x xy^2b^y \\ xa^xy^2b^y \\ xa^xyb^y \\ xa^xb^y \\ a^xy^2b^y \\ a^xyb^y \\ a^xb^y \end{bmatrix},
\begin{bmatrix}
1&0&0&0&0&0&0\\
0&b&2b&b&0&0&0\\
0&0&b&b&0&0&0\\
0&0&0&b&0&0&0\\
0&0&0&0&b&2b&b\\
0&0&0&0&0&b&b\\
0&0&0&0&0&0&b\\
\end{bmatrix},
\begin{bmatrix}
1&a&0&0&a&0&0\\
0&a&0&0&a&0&0\\
0&0&a&0&0&a&0\\
0&0&0&a&0&0&a\\
0&0&0&0&a&0&0\\
0&0&0&0&0&a&0\\
0&0&0&0&0&0&a\\
\end{bmatrix}
$$
因此，我们只需设计解决矩乘问题的算法，类欧问题便能得到解决。
## 二、解决思路
我们形式化地将这个问题表述如下：
> 给定一个含有$n$个$R$的$U,R$矩阵序列，序列末尾是$R$。第$x$个$R$前恰有$y=\left \lfloor \frac {px+r}q\right\rfloor$个$U$。
>
> 求这些矩阵的乘积。
我们把上述问题的答案记为$\text{sol}(p,q,r,n,U,R)$。
为了使整个问题的形式更为严整，我们规定$r<q$。
接下来可以用如下方式缩小问题的规模：
### 1. 合并
注意到，当$p\ge q$时，每个$R$之前都至少有$\left \lfloor \frac p q\right\rfloor$个连续的$U$，那么我们把这些$U$和$R$合并到一起形成新的$R$，这样每个$R$之前还剩的$U$个数为：
$$
\left\lfloor\frac{px+r}q\right\rfloor-x\left\lfloor\frac p q\right\rfloor=
\left\lfloor\frac{px+r}q\right\rfloor-x\frac {p-p \bmod q} q=
\left\lfloor\frac{px-px+(p \mod q )x+r}q\right\rfloor=
\left\lfloor\frac{(p \bmod q)x+r}q\right\rfloor
$$
于是我们得到
$$
\text{sol}(p,q,r,n,U,R)=\text{sol}(p\bmod q,q,r,n,U,U^{\left\lfloor\frac p q\right\rfloor}R)
$$
这一步操作的意义在于，$p$的规模能得以大幅减小。从几何的角度来讲，它相当于把直线按$y$轴缩小到斜率$<1$，缩放中心为直线$yq=r$。
### 2. 翻转
考察第$y$个$U$之前的$R$个数，它是关于$x$的不等式
$$
y>\left\lfloor\frac{px+r}q\right\rfloor
$$
的最大整数解。
我们可以在整数上的解集不变的情况下对不等式进行改写：
$$
y>\left\lfloor\frac{px+r}q\right\rfloor
\Leftrightarrow y\ge \left\lfloor\frac{px+r}q\right\rfloor+1
\Leftrightarrow y\ge \frac{px+r+1}q
\Leftrightarrow x\le \frac{qy-r-1}p
\Leftrightarrow x\le \left\lfloor\frac{qy-r-1}p\right\rfloor
$$
此时$U$的个数为$m=\left\lfloor\frac{pn+r}q\right\rfloor$。
但是我们并不能就这样开心地把问题转化为$\text{sol}(q,p,-r-1,m,R,U)$。我们还需解决以下两个问题：
- 新的$r$应当在$[0,p)$范围内。
- 原序列以$R$结尾，翻转后以$U$结尾，这导致转化后不满足序列末尾是$R$。
对于第一个问题，我们把翻转前第一个$U$以及其前面的所有$R$单独拿出来算。
这样一来，现在的第$y$个$U$是原先第$y+1$个$U$，其前的$R$的个数变为：
$$
\left\lfloor\frac{q(y+1)-r-1}p\right\rfloor-\left\lfloor\frac {q-r-1} p\right\rfloor=
\left\lfloor\frac{qy+q-r-1}p\right\rfloor-\frac {(q-r-1)-(q-r-1) \bmod p} p=
\left\lfloor\frac{qy+(q-r-1) \bmod p}p\right\rfloor
$$
这样新的$r$就恰好在$[0,p)$范围内了。
对于第二个问题，我们把翻转前最后一段$R$单独拿出来算即可，最后总共有$n-\left\lfloor\frac{qm-r-1}p\right\rfloor$个$R$，下取整中计算了最后一个$U$前$R$的数目。
于是我们又得到：
$$
\text{sol}(p,q,r,n,U,R)=R^{\left\lfloor\frac {q-r-1} p\right\rfloor}U\cdot\text{sol}(q,p,(q-r-1)\bmod p,m-1,R,U)\cdot R^{n-\left\lfloor\frac{qm-r-1}p\right\rfloor}
$$
当然，$m$等于$0$时并不符合递归条件，但这将意味着整个序列中不存在$U$，于是此时有：
$$
\text{sol}(p,q,r,n,U,R)=R^{n}
$$
这一步操作事实上将直线关于$y=x$进行了对称，随后左移$1$个单位，为了处理直线过格点时的特例，我们在将其下移$\frac 1 p$个单位。
### 3. 合并
注意到操作$1,2$必然是交替进行的，那么我们可以把它们合并起来，即：
$$
\text{sol}(p,q,r,n,U,R)=
\begin{cases}
R^n & m=0\\
R^{\left\lfloor\frac {q-r-1} p\right\rfloor}U\cdot\text{sol}(q \bmod p,p,(q-r-1)\bmod p,m-1,R,R^{\left\lfloor\frac {q} p\right\rfloor}U)\cdot R^{n-\left\lfloor\frac{qm-r-1}p\right\rfloor} & m>0
\end{cases}
$$
注意到，如果$p=0$，由$m=\left\lfloor\frac{pn+r}q\right\rfloor,r<q$可知$m$必为$0$，那么此时递归结束。
而$p,q$的变化与求$\gcd(p,q)$时它们的变化完全一致，因此递归的层数是$O(\log p)$级别。
## 三、实现细节
1. 经过上一节的讨论，我们已经将问题足够地标准化了，因此并不需要通过大量特判处理$p,q,r,n$的各种边界情况。
2. $U,R$实际处理时不必写成矩阵的形式，因为他们事实上处理的是$x,y$分别连续变化一段距离后对需要维护的各个量产生的影响。因此只需要给予它们与矩阵等价的定义即可，由矩阵的理论基础，新定义下的结合律仍然存在。这样能够降低一定代码复杂度。
3. 一般题目初始并不会保证$p,r\in[0,q)$，我们只需求$R^{\left\lfloor\frac r q\right\rfloor}\text{sol}(p\bmod q,q,r \bmod q,n,U,U^{\left\lfloor\frac p q\right\rfloor}R)$。
4. 可以通过适当封装降低代码复杂度。
算法的大致框架如下：
```cpp
ele operator ^ (ele x,ll y){ele z;for(;y;y>>=1,x=x*x)if(y&1) z=z*x;return z;}
ele sol(ll p,ll q,ll r,ll n,ele U,ele R){
ll m=(p*n+r)/q-1;
return ~m?r=q-r-1,(R^r/p)*U*sol(q%p,p,r%p,m,R,(R^q/p)*U)*(R^n-(q*m+r)/p):R^n;
}
int main(){
ans=(U^r/q)*sol(p%q,q,r%q,n,U,(U^p/q)*R);
}
```
只要合理的设计元素，并规定乘法。上述代码就能准确无误地运行任何目前已知的类欧几里得问题。
## 四、缺陷
尽管万能欧几里得具有代码好写，拓展性强的优势，它仍然无法解决某些特殊的问题，例如
$$
\sum \sqrt y,\sum x^y,\sum a^{xy},\sum 2^{2^y},\sum A\sin(\omega y+\phi)
$$
事实上，如果求和号内的部分关于$x,y$中任意一个的差分影响不是线性的，我们就难以设计出合理的$U,R$矩阵并套用万能欧几里得解决。
所以目前看起来可以解决的问题是：
$$
\sum f(x)a^{x}g(y)b^y
$$
其中$f,g$是分别关于$x,y$的多项式。
$a,b$并不必是一个常数。只要$a,b$在一个环上就行（例如矩阵，复数等）。
当然，如果跳出$\sum_x f(x,y)$的限制，直接去构建特定的$U,R$矩阵的话，可以产生的变化就多得多了。
---
最后附上$\text{LOJ}6440$的代码，该问题求得是$\sum A^xB^y$，$A,B$是矩阵。
```cpp
#include<bits/stdc++.h>
#define FF(i) for(int i=1;i<=T;i++)
using namespace std;
typedef long long ll;
const int P=998244353,N=21;
ll p,q,r,n,T;
void chk(int&x){x-=P;x+=x>>31&P;}
struct mat{int v[N][N];mat(){memset(v,0,sizeof v);}}A,B,I;
mat operator + (const mat&a,const mat&b){mat c;FF(i)FF(j) chk(c.v[i][j]=a.v[i][j]+b.v[i][j]);return c;}
mat operator * (const mat&a,const mat&b){mat c;FF(i)FF(j)FF(k) chk(c.v[i][j]+=1ll*a.v[i][k]*b.v[k][j]%P);return c;}
struct ele{
mat x,y,s;ele(){x=y=I;}
ele operator * (const ele&a){ele b;b.x=x*a.x,b.y=y*a.y,b.s=s+x*a.s*y;return b;}
}U,R,S;
ele operator ^ (ele x,ll y){ele z;for(;y;y>>=1,x=x*x)if(y&1) z=z*x;return z;}
ele sol(ll p,ll q,ll r,ll n,ele U,ele R){ll m=((__int128)p*n+r)/q-1;return ~m?r=q-r-1,(R^r/p)*U*sol(q%p,p,r%p,m,R,(R^q/p)*U)*(R^n-((__int128)q*m+r)/p):R^n;}
int main(){
scanf("%lld%lld%lld%lld%lld",&p,&q,&r,&n,&T);
FF(i)FF(j) scanf("%d",&A.v[i][j]),I.v[i][j]=i==j;
FF(i)FF(j) scanf("%d",&B.v[i][j]);
U=R=ele();
U.y=B;R.x=R.s=A;
S=(U^r/q)*sol(p%q,q,r%q,n,U,(U^p/q)*R);
FF(i)FF(j) cout<<S.s.v[i][j]<<" \n"[j==T];
}
```