P5170 【模板】类欧几里得算法

Description

要求在$O(lgn)$的时间内求出：
${\sum }_{i=0}^n\lfloor \frac{ai+b}{c}\rfloor$

${\sum }_{i=0}^{n}i\lfloor \frac{ai+b}{c}\rfloor$

${\sum }_{i=0}^n\lfloor \frac{ai+b}{c}{\rfloor }^2$

Solution

其实也就是一种优秀的求和思想。

Part one

我们先考虑最基本的第一个：
${\sum }_{i=0}^n\lfloor \frac{ai+b}{c}\rfloor ={\sum }_{i=0}^n(\lfloor \frac{a^{\prime }i+b^{\prime }}{c}\rfloor +pi+q)$

其中$p=\lfloor \frac{a}{c}\rfloor ,q=\lfloor \frac{b}{c}\rfloor$

那么后面的$\sum pi+q$很容易计算，因此问题转化为求：

$f(a,b,c,n)={\sum }_{i=0}^n\lfloor \frac{ai+b}{c}\rfloor (a\in [0,c),b\in [0,c))$

Part two

然后使用一个极其诡妙的思路：
$f(a,b,c,n)$
$={\sum }_{i=0}^n\lfloor \frac{ai+b}{c}\rfloor$
$={\sum }_{i=0}^n{\sum }_{j=1}^m[j\le \lfloor \frac{ai+b}{c}\rfloor ](m=\lfloor \frac{an+b}{c}\rfloor )$
$={\sum }_{i=0}^n{\sum }_{j=0}^{m-1}[jc+c\le ai+b]$
$={\sum }_{j=0}^{m-1}{\sum }_{i=0}^n[jc+c\le ai+b]$
$={\sum }_{j=0}^{m-1}(n+1-{\sum }_{i=0}^n[jc+c>ai+b])$
$={\sum }_{j=0}^{m-1}(n-\lfloor \frac{jc+c-b-1}{a}\rfloor )$
$=nm-{\sum }_{j=0}^{m-1}\lfloor \frac{jc+c-b-1}{a}\rfloor$
$=nm-f(c,c-b-1,a,m-1)$

因此可以递归计算，直到$a=0$时，我们直接算出答案。

不难发现$a,c$相当于做了一个辗转相除的过程，因此时间复杂度为$O(lgn)$。

Part three

剩下两个也是同样的思路，先考虑让$a^{\prime }=a\mathrm{m}\mathrm{o}\mathrm{d}c,b^{\prime }=b\mathrm{m}\mathrm{o}\mathrm{d}c$，然后快速计算下取整式外的贡献，下取整式内的部分通过类似上面的方法转化，推导出若干个能用类欧计算的式子相加的形式，然后递归计算即可。这里不再赘述。

Code

#include <bits/stdc++.h>using namespace std;template<typename T> inline bool upmin(T &x, T y) { return y < x ? x = y, 1 : 0; }
template<typename T> inline bool upmax(T &x, T y) { return x < y ? x = y, 1 : 0; }#define MP(A,B) make_pair(A,B)
#define PB(A) push_back(A)
#define SIZE(A) ((int)A.size())
#define LEN(A) ((int)A.length())
#define FOR(i,a,b) for(int i=(a);i<(b);++i)
#define fi first
#define se secondtypedef long long ll;
typedef unsigned long long ull;
typedef long double lod;
typedef pair<int, int> PR;
typedef vector<int> VI; const lod eps = 1e-9;
const lod pi = acos(-1);
const int oo = 1 << 30;
const ll loo = 1ll << 60;
const int mods = 998244353;
const int inv2 = (mods + 1) >> 1;
const int inv6 = (mods + 1) / 6;
const int MAXN = 600005;
const int INF = 0x3f3f3f3f; //1061109567
/*--------------------------------------------------------------------*/namespace FastIO{constexpr int SIZE = (1 << 21) + 1;int num = 0, f;char ibuf[SIZE], obuf[SIZE], que[65], *iS, *iT, *oS = obuf, *oT = obuf + SIZE - 1, c;#define gc() (iS == iT ? (iT = ((iS = ibuf) + fread(ibuf, 1, SIZE, stdin)), (iS == iT ? EOF : *iS ++)) : *iS ++)inline void flush() {fwrite(obuf, 1, oS - obuf, stdout);oS = obuf;}inline void putc(char c) {*oS ++ = c;if (oS == oT) flush();}inline void getc(char &c) {for (c = gc(); !isalpha(c) && c != EOF; c = gc());}inline void reads(char *st) {char c;int n = 0;getc(st[++ n]);for (c = gc(); isalpha(c) ; c = gc()) st[++ n] = c;}template<class I>inline void read(I &x) {for (f = 1, c = gc(); c < '0' || c > '9' ; c = gc()) if (c == '-') f = -1;for (x = 0; c >= '0' && c <= '9' ; c = gc()) x = (x << 3) + (x << 1) + (c & 15);x *= f;}template<class I>inline void print(I x) {if (x < 0) putc('-'), x = -x;if (!x) putc('0');while (x) que[++ num] = x % 10 + 48, x /= 10;while (num) putc(que[num --]);}struct Flusher_{~Flusher_(){flush();}} io_Flusher_;
}
using FastIO :: read;
using FastIO :: putc;
using FastIO :: reads;
using FastIO :: print;struct Node{ int f, g, h; Node() { f = g = h = 0; }Node(int x, int y, int z):f(x), g(y), h(z){};
};
int S1(int x) { return 1ll * x * (x + 1) / 2 % mods; }
int S2(int x) { return 1ll * x * (x + 1) % mods * (x * 2 + 1) % mods * inv6 % mods; }
int upd(int x, int y) {return x + y >= mods ? x + y - mods : x + y; 
}
Node solve(int a, int b, int c, int n) {Node Ans, ans;int p = 0, q = 0;if (a >= c || b >= c) {p = a / c, q = b / c;Ans.f = upd(Ans.f, 1ll * p * S1(n) % mods);Ans.f = upd(Ans.f, 1ll * q * (n + 1) % mods);Ans.g = upd(Ans.g, 1ll * p * p % mods * S2(n) % mods);Ans.g = upd(Ans.g, 1ll * q * q % mods * (n + 1) % mods);Ans.g = upd(Ans.g, 2ll * p * q % mods * S1(n) % mods);Ans.h = upd(Ans.h, 1ll * p * S2(n) % mods);Ans.h = upd(Ans.h, 1ll * q * S1(n) % mods);a -= p * c, b -= q * c;}if (!a) return Ans;int m = ((ll)a * n + b) / c;Node t = solve(c, c - b - 1, a, m - 1);ans.f = upd(ans.f, upd(1ll * m * n % mods, mods - t.f));ans.h = upd(ans.h, upd(1ll * m * S1(n) % mods, mods - 1ll * upd(t.g, t.f) * inv2 % mods));ans.g = upd(ans.g, upd(2ll * ans.h * p % mods, 2ll * ans.f * q % mods));ans.g = upd(ans.g, 1ll * m * m % mods * n % mods);ans.g = upd(ans.g, mods - 2ll * t.h % mods);ans.g = upd(ans.g, mods - t.f);return Node(upd(ans.f, Ans.f), upd(ans.g, Ans.g), upd(ans.h, Ans.h));
}
signed main() {
#ifndef ONLINE_JUDGEfreopen("a.in", "r", stdin);
#endifint Case;read(Case);while (Case --) {int n, a, b, c;read(n), read(a), read(b), read(c);Node ans = solve(a, b, c, n);print(ans.f), putc(' '), print(ans.g), putc(' '), print(ans.h), putc('\n');}return 0;
}