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

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 second

typedef 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_JUDGE
	freopen("a.in", "r", stdin);
#endif
	int 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;
}