试图讲解万能欧几里得算法
神秘算法，抽象至极。
学习自 ix35大犇的题解，但我死活没看懂，大概是我太菜了吧，写一篇题解理清思路以及加以补充。
记号有部分变动。
问题
万能欧几里得算法解决的是形如，给出直线 \(y=\frac{kx+b}{m}\)，计算直线下方的整坐标点的贡献，贡献函数应于坐标有关，且满足结合律，即可合并。
比如模板题 【模板】类欧几里德算法 这样的：
给定 \(n,\,a,\,b,\,c\) ，分别求 \(\sum\limits_{i=0}^{n}\left\lfloor \frac{ai+b}{c} \right\rfloor\,,\ \sum\limits_{i=0}^{n}{\left\lfloor \frac{ai+b}{c} \right\rfloor}^2\,,\ \sum\limits_{i=0}^{n}i\left\lfloor \frac{ai+b}{c} \right\rfloor\) ，答案对 \(998244353\) 取模。多组数据。
万能欧几里得算法将计算复杂度降到了 \(\log\) 级别。
算法内容
从左往右扫这条直线，定义与平行于 \(x\) 轴且 \(y\) 坐标为整数的直线相交为 \(U\)，与平行于 \(y\) 轴且 \(x\) 坐标为整数的直线相交为 \(R\)，特别的，如果同时满足记为 \(UR\)。这里的记号表示对答案的贡献，分成两种是因为通常来说 \(U\) 和 \(R\) 对答案的贡献往往不同。
定义函数 \(f(k,b,m,n,U,R)\) 为当前每一恒行的贡献为 \(U\)，纵列贡献为 \(R\)，直线 \(y=\frac{kx+b}{m}\) 在定义域 \(x\in (0,n]\) 产生的 \(UR\) 序列。
那么直线将得到一个确定的 \(UR\) 序列。不难发现直线跨过整坐标存在周期性规律。
若 \(k\geq m\)
对应的几何意义是：不难发现每个 \(R\) 之前必然会有 \(\Big\lfloor\frac{k}{m}\Big\rfloor\) 个 \(U\)，我们不妨把这些整行打包起来放到 \(R\) 里，创建一个新的 \(R'\)。
其中 \(R'=U^{\Big\lfloor\frac{k}{m}\Big\rfloor}R\)。
若 \(k<m\)
转一下坐标系，翻转 \(x\)、\(y\) 的地位，反解得到 \(x=\frac{my-b}{k}\)。
但这样有一个问题，我们一开始定义的，如果同时与横线与竖线相交记为 \(UR\)，在这里会刚好反过来。所以我们向下偏移一个很小的量，完美避开这个问题。
仍然有很多细节问题。
定义的时候，我们的定义域是 \((0,n]\)，而由于转之前的 \(y\) 轴截距未必是整数，转过来之后没法完美贴合到 \(x\) 轴上，所以我们要把这一段裁掉，也就是转之前第一个 \(U\) 及其以下的部分。
旋转坐标系后平移直线，带入 \(y=1\) 得到有 \(\lfloor\frac{m-b-1}{k}\rfloor\) 个旋转前的 \(R\)，再拼上一个 \(U\)（因为左开右闭）。
新的 \(n'=\lfloor\frac{kn+b}{m}\rfloor\) 应为旋转前 \(U\) 的个数。
但这样又让坐标偏移了一位，提出来这一位偏移加入常数项：
根据定义应该以 \(R\) 结尾。但转过去之后 \(R\) 变成 \(U\) 就不符合定义了，所以末尾这一部分也要截掉，这一部分有 \(n-\lfloor\frac{mn'-b-1}{k}\rfloor\) 个。
\(n'\) 减一是刚才截掉的。
若我们发现 \(n'=0\)，那说明当前全是 \(R\)，直接返回 \(R^{n}\) 即可。
还有问题：初始直线是有截距的，算法处理的必须过原点，所以我们需要先把截距的 \(\lfloor\frac{b}{m}\rfloor\) 个 \(U\) 合并在前面。
至此我们完成了贡献的合并。
模板题
原题解讲的很清楚了，感觉只要理解上面的东西这些都不是啥问题。
#include <bits/stdc++.h>
#define int int64_t
//#define int __int128
//#define MOD (1000000007)
//#define eps (1e-6)
#define endl '\n'
#define debug_endl cout<<endl;
#define debug cout<<"debug"<<endl;
using namespace std;
template<typename T_,T_ mod>
struct ModT_ {
const static T_ MD = mod;
T_ x;
inline ModT_(T_ x = 0): x(x % MD) {}
T_ get() { return x; }
inline ModT_ operator + (const ModT_ &that) const { T_ x0 = x + that.x; return ModT_(x0 < MD? x0 : x0 - MD); }
inline ModT_ operator - (const ModT_ &that) const { T_ x0 = x - that.x; return ModT_(x0 < MD? x0 + MD : x0); }
inline ModT_ operator * (const ModT_ &that) const { return ModT_((long long)x * that.x % MD); }
inline ModT_ operator / (const ModT_ &that) const { return *this * that.inverse(); }
inline ModT_ operator += (const ModT_ &that) { x += that.x; if (x >= MD) x -= MD; return *this;}
inline ModT_ operator -= (const ModT_ &that) { x -= that.x; if (x < 0) x += MD; return *this;}
inline ModT_ operator *= (const ModT_ &that) { x = (long long)x * that.x % MD; return *this;}
inline ModT_ operator /= (const ModT_ &that) { *this = *this / that; return *this;}
inline bool operator == (const ModT_ &that) { return (x==that.x); }
inline bool operator != (const ModT_ &that) { return (x!=that.x); }
inline bool operator > (const ModT_ &that) {return x > that.x; }
inline bool operator < (const ModT_ &that) {return x < that.x; }
inline bool operator >= (const ModT_ &that) {return x >= that.x; }
inline bool operator <= (const ModT_ &that) {return x <= that.x; }
inline bool operator == (const int &that) { return (x==that); }
inline bool operator != (const int &that) { return (x!=that); }
inline bool operator > (const int &that) {return x > that; }
inline bool operator < (const int &that) {return x < that; }
inline bool operator >= (const int &that) {return x >= that; }
inline bool operator <= (const int &that) {return x <= that; }
inline ModT_ operator = (const int &that) {x=that;return *this;}
inline ModT_ inverse() const {
T_ a = x, b = MD, u = 1, v = 0;
while(b) {
T_ t = a / b;
a -= t * b; std::swap(a, b);
u -= t * v; std::swap(u, v);
}
if(u < 0) u += MD;
return u;
}
};
typedef ModT_<int,998244353> modint;
struct Contrib{
modint cntr,cntu,sums,sumi,sumss,sumis;
Contrib(){cntr=cntu=sums=sumi=sumss=sumis=0;}
Contrib(int _cntr=0,int _cntu=0,int _sums=0,int _sumi=0,int _sumss=0,int _sumis=0){
cntr=_cntr,cntu=_cntu,sums=_sums,sumi=_sumi,sumss=_sumss,sumis=_sumis;
}
void operator+=(const Contrib& other){
sumss+=other.sumss+other.cntr*cntu*cntu+other.sums*cntu*2;//后面的s整体抬升cntu，展开∑(cntu+s)^2
sumis+=other.sumis+other.cntr*cntu*cntr+cntr*other.sums+other.sumi*cntu;//展开∑(cntr+i)(cntu+s)
sums+=other.sums+other.cntr*cntu;
sumi+=other.sumi+other.cntr*cntr;
cntr+=other.cntr;
cntu+=other.cntu;
}
Contrib operator+(const Contrib& other){
Contrib res=*this;
res+=other;
return res;
}
};
Contrib U(0,1,0,0,0,0),R(1,0,0,1,0,0);
Contrib qpow(Contrib a,int b){
Contrib res(0,0,0,0,0,0);
while(b){
if(b&1) res+=a;
a+=a;
b>>=1;
}
return res;
}
Contrib solve(int k,int b,int m,int n,Contrib U,Contrib R){
if(!n) return Contrib(0,0,0,0,0,0);
if(k>=m) return solve(k%m,b,m,n,U,qpow(U,k/m)+R);
int nn=(k*n+b)/m;
if(!nn) return qpow(R,n);
return qpow(R,(m-b-1)/k)+U+solve(m,(m-b-1)%k,k,nn-1,R,U)+qpow(R,n-((m*nn-b-1)/k));
}
signed main(){
//freopen(".in","r",stdin);
//freopen(".out","w",stdout);
ios::sync_with_stdio(false);
cin.tie(0),cout.tie(0);
int T;
cin>>T;
while(T--){
int n,a,b,c;
cin>>n>>a>>b>>c;
Contrib res=qpow(U,b/c)+solve(a,b%c,c,n,U,R);
cout<<(res.sums+(b/c)).x<<' ';
cout<<(res.sumss+((b/c)*(b/c))).x<<' ';
cout<<res.sumis.x<<endl;
}
return 0;
}
浙公网安备 33010602011771号