高次类欧几里得算法

LibreOJ #138 类欧几里得算法

题意

$T \leq 1000$ 组数据。给定 $n, a, b, c\in [0, 10^9], k_1,k_2\in [0, 10], k_1+k_2\leq 10$，求取 $$\sum_{x=0}^{n}x^{k_1}\left\lfloor\frac{ax+b}{c}\right\rfloor^{k_2}$$

部分分——观察性质

当 $k_1=0,k_2=1$ 时，是为一般类欧几里得算法的模板。上述链接亦给出了 $k_1=0,k_2=2; k_1=1,k_2=1$ 时的解法。

当 $k_2=0$ 时，该式退化为 $\sum_{x=0}^{n}x^{k_1}$，即 $k_1$ 次的等幂求和。有一个模糊的结论：

设 $n$ 为正整数。则$k$ 次的等幂求和，$\sum_{x=0}^{n}x^k$，是一个关于 $n$ 的 $k+1$ 次多项式。

该结论来源于 zyw 学姐多项式专题所选题目 BZOJ 3453 – tyvj 1858 XLkxc。之所以说模糊，是因为该多项式的系数是已知的。不过我们仍然可以暴力计算 $k+2$ 个点以插值。

再观察 OI-Wiki 上对于 $k_1+k_2\leq 2$ 时的解法。在求取 $$h(n,a,b,c)=\sum_{x=0}^n\left\lfloor\frac{ax+b}{c}\right\rfloor^2$$时采取了如下转化：

$$\begin{aligned} x^2&=2\times \frac{x(x+1)}{2}-x\\ &=2\left(\sum_{y=1}^{x}y\right)-x\end{aligned} \quad\begin{aligned} \Longrightarrow\sum_{x=0}^n\left\lfloor\frac{ax+b}{c}\right\rfloor^2=\sum_{x=0}^n2\left(\sum_{y=1}^{\left\lfloor\frac{ax+b}{c}\right\rfloor}y\right)-\left\lfloor\frac{ax+b}{c}\right\rfloor\\ \end{aligned}$$

这样一来，由于 $y$ 是一个线性算子，在 $ax+b\geq yc$ 时都有 $y$ 向总和贡献，于是我们可以应用类似于 $k_1=0,k_2=1$ 时的方法转化贡献。这个过程对 $k_2$ 作了降次。

那么，对于更高次项，有无办法采取同样的办法转化呢？是否可以设计一些转化，使得要求取的函数 $f(\cdots,k_1,k_2)$ 能够由 $f(\cdots,k_1-?,k_2-?)$ 推出呢？

伯努利数

$$B_{0\dots 10}=\left\{1,-\frac{1}{2},\frac{1}{6},0,-\frac{1}{30},0,\frac{1}{42},0,-\frac{1}{30},0,\frac{5}{66}\right\}$$

记 $S_m(n)=\sum\limits_{x=0}^{\color{red}n-1}x^m$，则有如下性质： $$S_m(n)=\frac{1}{m+1}\sum_{k=0}^{m}\binom{m+1}{k}B_kn^{m+1-k}$$

用红笔标注是因为第二版推导的“等幂求和”定义是 $\sum_{x=1}^{n}x^m$（OI-Wiki 的转化——在下文中出现的那个——用的是这种）。最终检查时才发现不对劲。

（不过，将 $B_1$ 改为 $\frac{1}{2}$ 的情况下，应用“错误”的“等幂求和”公式得到的结果是正确的，似乎也能计算。但这会牵扯到 $x=0$ 时的计算 和 $0^0=1$ 的定义问题，有点麻烦。暂且按下不表。）

我们移项可得 $$\begin{aligned} (m+1)S_m(n)&=\binom{m+1}{0}B_0n^{m+1}+\sum_{k=1}^{m}\binom{m+1}{k}B_kn^{m+1-k}\\ n^{m+1}&=(m+1)S_m(n)-\sum_{k=1}^{m}\binom{m+1}{k}B_kn^{m+1-k}\\ n^{m}&=mS_{m-1}(n)-\sum_{k=1}^{m-1}\binom{m}{k}B_kn^{m-k} \end{aligned}$$

上式右侧是一个等幂求和以及若干 $n^k, \underline{k<m}$ 与常数项之积的和。当 $m>1$ 时，若能够快速处理 $mS_{m-1}(n)$ 这一部分的贡献，便能使用 $k_2$ 更小时的函数值计算。这便是处理 $\lfloor\frac{ax+b}{c}\rfloor^m$ 的转化的一般形式。

推导过程

类似于一般类欧的思路，我们尽可能将 $(a,b,c)$ 的情况化归为 $(a\bmod c,?,c)$ 再转化贡献。

推导 $f$

$$f(n,a,b,c,k_1,k_2)=\sum_{x=0}^{n}x^{k_1}\left\lfloor\frac{ax+b}{c}\right\rfloor^{k_2}$$

记 $\alpha=\lfloor\frac{a}{c}\rfloor,\beta=\lfloor\frac{b}{c}\rfloor,a’=a\bmod c,b’=b\bmod c$。则 $$\begin{aligned} \text{原式}&=\sum_{x=0}^{n}x^{k_1}\left\lfloor\frac{(\alpha c+a’)x+(\beta c+b’)}{c}\right\rfloor^{k_2}\\ &=\sum_{x=0}^{n}x^{k_1}\left\lfloor\alpha x+\beta+\frac{a’x+b’}{c}\right\rfloor^{k_2}\\ &=\sum_{x=0}^{n}x^{k_1}\left(\alpha x+\beta+\left\lfloor\frac{a’x+b’}{c}\right\rfloor\right)^{k_2}&\text{（整数自由移出取整号）}\\ &=\sum_{x=0}^{n}x^{k_1}\sum_{\substack{e_1+e_2+e_3=k2\\e_1,e_2,e_3\in \mathbb{N}}}\binom{k1}{e_1,e_2,e_3}(\alpha x)^{e_1}\beta^{e_2}\left\lfloor\frac{a’x+b’}{c}\right\rfloor^{e_3}\\ &=\sum_{\substack{e_1+e_2+e_3=k2\\e_1,e_2,e_3\in \mathbb{N}}}\alpha^{e_1}\beta^{e_2}\binom{k2}{e_1,e_2,e_3}\sum_{x=0}^{n}x^{k_1+e_1}\left\lfloor\frac{a’x+b’}{c}\right\rfloor^{e_3} \end{aligned}$$

$\binom{k_2}{e_1,e_2,e_3}$ 是多项式系数。有 $\binom{n}{m_1,m_2,m_3}=\frac{n!}{m_1!m_2!m_3!}$，是为 $(a+b+c)^n$ 展开式中 $a^{m_1}b^{m_2}c^{m_3}$ 的系数。

推导 $f’$

$$f'(n,a’,b’,c,k_1,k_2)=\sum_{x=0}^{n}x^{k_1}\left\lfloor\frac{a’x+b’}{c}\right\rfloor^{k_2}\quad(a’,b’\in[0,c))$$

记 $t_x=\lfloor\frac{a’x+b’}{c}\rfloor$，则 $$\begin{aligned} \text{原式}&=\sum_{x=0}^{n}x^{k_1}\left(k_2S_{k_2-1}(t_x)-\sum_{q=1}^{k_2-1}\binom{k_2}{q}B_qt_x^{k_2-q}\right)&\text{（等幂求和公式变形）}\\ &=k_2\left(\sum_{x=0}^{n}x^{k_1}S_{k_2-1}(t_x)\right)-\sum_{q=1}^{k_2-1}\binom{k_2}{q}B_q\sum_{x=0}^{n}x^{k_1}t_x^{k_2-q}\\ &=k_2\left(\sum_{x=0}^{n}x^{k_1}S_{k_2-1}(t_x)\right)-\sum_{q=1}^{k_2-1}\binom{k_2}{q}B_qf'(n,a’,b’,c,k_1,k_2-q) \end{aligned}$$

右侧可以由 $k_2′<k_2$ 的函数值求出；现在考虑如何转化左侧的贡献。

推导 $g’$

$$\begin{aligned} g'(n,a’,b’,c,k_1,k_2)&=\sum_{x=0}^{n}x^{k_1}S_{k_2}\left(\left\lfloor\frac{a’x+b’}{c}\right\rfloor\right)&(a’,b’\in[0,c))\\ &=\sum_{x=0}^{n}x^{k_1}\sum_{y=0}^{\color{red}\lfloor\frac{a’x+b’}{c}\rfloor-1}y^{k_2} \end{aligned}$$

发现 $S_m(n)$ 是前缀和的形式。则对于固定的 $y$，当 $$\begin{aligned} &\left\lfloor\frac{a’x+b’}{c}\right\rfloor-1\geq y\\ \Longrightarrow&\space a’x+b\geq (y+1)c\\ \Longrightarrow&\space x>\left\lfloor\frac{yc+c-b’-1}{a’}\right\rfloor \end{aligned}$$ 时，均有 $y^{k_2}$ 对总和作出贡献。

记 $p=\lfloor\frac{a’n+b’}{c}\rfloor$。则 $$\begin{aligned} \text{原式}&=\sum_{y=0}^{p-1}y^{k_2}\sum_{x=\lfloor\frac{yc+c-b’-1}{a’}\rfloor+1}^{n}x^{k_1}\\ &=\sum_{y=0}^{p-1}y^{k_2}\left(S_{k_1}(n+1)-S_{k_1}\left(\left\lfloor\frac{yc+c-b’-1}{a’}\right\rfloor{\color{blue}+1}\right)\right)\\ &=S_{k_2}(p)S_{k_1}(n+1)-\sum_{y=0}^{p-1}y^{k_2}S_{k_1}\left(\left\lfloor\frac{yc+c{\color{blue}+a’}-b’-1}{a’}\right\rfloor\right) \end{aligned}$$

右侧看起来……很眼熟。它像是 $g'(p-1,c,c+a’-b’-1,a’,k_2,k_1)$，但此时不可能有 $c\in [0, a)$。我们必须在这个唯一利用到下层递归函数值的地方做文章。

推导 $g$

$$\begin{aligned} g(n,a,b,c,k_1,k_2)&=\sum_{x=0}^{n}x^{k_1}S_{k_2}\left(\left\lfloor\frac{ax+b}{c}\right\rfloor\right)\\ &=\sum_{x=0}^{n}x^{k_1}\sum_{q=0}^{k_2}\frac{1}{k_2+1}\binom{k_2+1}{q}B_q\left\lfloor\frac{ax+b}{c}\right\rfloor^{k_2+1-q}&\text{（等幂求和公式）}\\ &=\sum_{q=0}^{k_2}\frac{1}{k_2+1}\binom{k_2+1}{q}B_q\sum_{x=0}^{n}x^{k_1}\left\lfloor\frac{ax+b}{c}\right\rfloor^{k_2+1-q} \end{aligned}$$

右侧终于出现了我们所熟知的 $f(n,a,b,c,k_1,k_2+1-q)$。可以发现在同一层中 $f$ 的计算仅仅依赖 $g’$ 而非 $g$，$g’$ 调用的是下层递归的 $g$，故我们已经完成了计算。

最终的递归式子

结合上文，我们可以立刻得到： $$\begin{aligned} f(n,a,b,c,k_1,k_2)&=\sum_{\substack{e_1+e_2+e_3=k_2\\e_1,e_2,e_3\in\mathbb{N}}}\alpha^{e_1}\beta^{e_2}\binom{k_2}{e_1,e_2,e_3}f'(n,a’,b’,c,k_1+e_1,k_2)\\ f'(n,a’,b’,c,k_1,k_2)&=k_2g'(n,a’,b’,c,k_1,k_2-1)-\sum_{q=1}^{k_2-1}\binom{k_2}{q}B_qf'(n,a’,b’,c,k_1,k_2-q)\\ g'(n,a’,b’,c,k_1,k_2)&=S_{k_2}(p)S_{k_1}(n+1)-g(p-1,c,c+a’-b’-1,a’,k_2,k_1)&(p=\left\lfloor\frac{a’n+b’}{c}\right\rfloor)\\ g(n,a,b,c,k_1,k_2)&=\sum_{q=0}^{k_2}\frac{1}{k_2+1}\binom{k_2+1}{q}B_qf(n,a,b,c,k_1,k_2+1-q) \end{aligned}$$

边界条件有 $n<0\lor c=0$ 时返回 $0$；同时 $f'(n,a,b,c,k_1,k_2=0)=S_{k_1}(n+1)$，不需要（也不能——因为 $S_{-1}(n)$ 是未定义的）递归计算。

慢着。细心的你可能发现了一个问题：我们调用了 $g'(\cdots,k_1,k_2)$，$g’$ 调用了下一层的 $g$，而 $g(\cdots,k_1,k_2)$ 的计算式里要调用 $f(\cdots,k_1,k_2+1-q)$，$q=0$ 的话岂不是造成了无限递归？！

稍加分析可得，在 $f'(\cdots,k_1,k_2)$ 中我们最多调用 $g'(\cdots,k_1,k_2-1)$，而其调用下一层的 $g(\cdots,k_2-1,k_1)$，此时套入上式得到其最多调用 $f(\cdots,k_2-1,k_1+1)$，在此过程中 $\text{两参数之和}=k_2+k_1$，保持不变。故而我们只需计算两参数之和 $\leq\text{（题面给定的）}k_1+k_2=k_\text{sum}$ 的所有函数值，它便是自洽的。

通过上述式子分析得到，单查询时间复杂度 $\operatorname{O}(\log\max(a,c)(k_1+k_2)^3)$。通过预先计算三项式系数、伯努利数与二项式系数之积，以及同层内的 $S_{k}(n),S_{k}(p),k\in[0,k_\text{sum}]$，可以优化常数。

提交记录 #1604976

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
#include <bits/stdc++.h>
using namespace std;

#define inl inline
/* 这次没有快读。 */
#define fst first
#define scd second
#define empb emplace_back
#define all(p) begin (p), end (p)
using ll = long long;
using ull = unsigned long long;
using pint = pair <int, int>;

constexpr int P = 1e9 + 7, G = 11;
constexpr ll B[] = { 1, 500000003, 166666668, 0, 766666672,
	0, 23809524, 0, 766666672, 0, 348484851 }, PP = (ll) P*P,
fact[] = { 1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800 },
finv[] = { 1, 1, 500000004, 166666668, 41666667, 808333339,
	301388891, 900198419, 487524805, 831947206, 283194722 },
inv[] = { 0, 1, 500000004, 333333336, 250000002, 400000003, 166666668,
	142857144, 125000001, 111111112, 700000005, 818181824 };
struct coeffi {
	ll co[G+1][G+1];
	constexpr coeffi () : co () {
		for (int x = 0; x <= G; ++x) {
			co[x][0] = 1;
			for (int y = 1; y <= x; ++y)
				co[x][y] = co[x-1][y] + co[x-1][y-1];
		}
		for (int x = 0; x <= G; ++x)
			for (int y = 0; y < x; ++y)
				(co[x][y] *= B[y]) %= P;
	}
	inl constexpr const ll* operator []
		(int x) const { return co[x]; }
} constexpr co;
struct trinom {
	ll c[G][G][G];
	constexpr trinom () : c () {
		for (int k = 0; k < G; ++k)
			for (int x = 0; x <= k; ++x)
				for (int y = 0; y + x <= k; ++y)
					c[k][x][y] = fact[k] * finv[x] % P
						* finv[y] % P * finv[k-x-y] % P; 
	}
	inl constexpr const auto operator []
		(int x) const { return c[x]; }
} constexpr tri;
int T, n, a, b, c, k1, k2, ksum;
ll f[G][G], f_[G][G], g[G][G], g_[G][G],
	s_p[G], s_n1[G], _g[G][G];

inl ll S (int m, int n) {
	static ll np, res; np = n, res = 0;
	for (int k = m; k >= 0; --k, (np *= n) %= P)
		res += np * co[m+1][k] % P;
	return res % P * inv[m + 1] % P;
}
void calc (int n, ll a, ll b, ll c) {
	// ? 边界条件。 
	if (n < 0 || !c) {
		memset (f, 0, sizeof f),
		memset (g, 0, sizeof g);
		return; }
	const ll alp = a / c, bta = b / c, _a = a % c,
		_b = b % c, p = (n*_a + _b) / c;
	static ll tmp, ap, bp;
	calc (p - 1, c, c+_a-_b-1, _a);
	for (int k = 0; k <= ksum; ++k)
		s_p[k] = S (k, p), s_n1[k] = S (k, n + 1);
	memcpy (_g, g, sizeof g);
	// * 计算 g'。
	for (int k1 = 0; k1 <= ksum; ++k1)
		for (int k2 = 0; k2 + k1 <= ksum; ++k2)
			g_[k1][k2] = (s_p[k2]*s_n1[k1] + P - _g[k2][k1]) % P;
	// * 计算 f'。
	for (int k1 = 0; k1 <= ksum; ++k1) {
		f_[k1][0] = s_n1[k1], tmp = 0; // 单独计算 k2=0 的情况。 
		for (int k2 = 1; k2 + k1 <= ksum; ++k2, tmp = 0) {
			for (int q = 1; q < k2; ++q)
				tmp += co[k2][q] * f_[k1][k2-q];
			f_[k1][k2] = (k2 * g_[k1][k2-1] + PP - tmp % P) % P;
		}
	}
	// * 计算 f。
	for (int k1 = 0; k1 <= ksum; ++k1)
		for (int k2 = 0; k2 + k1 <= ksum; ++k2) {
			ap = bp = 1, tmp = 0;
			for (int e1 = 0; e1 <= k2; ++e1, (ap *= alp) %= P, bp = 1)
				for (int e2 = 0; e2 + e1 <= k2; ++e2, (bp *= bta) %= P)
					tmp += tri[k2][e1][e2] * ap % P * bp
						% P * f_[k1+e1][k2-e1-e2] % P;
			f[k1][k2] = tmp % P;
		}
	// * 计算 g。
	for (int k1 = 0; k1 <= ksum; ++k1)
		for (int k2 = 0; k2 + k1 + 1 <= ksum; ++k2) {
			for (int q = tmp = 0; q <= k2; ++q)
				tmp += co[k2+1][q] * f[k1][k2+1-q] % P;
			g[k1][k2] = tmp % P * inv[k2+1] % P;
		}
}

int main () {
	/*  */

scanf ("%d", &T);
while (T--) {
	scanf ("%d%d%d%d%d%d", &n, &a, &b, &c, &k1, &k2);
	ksum = k1 + k2; calc (n, a, b, c);
	printf ("%lld\n", f[k1][k2]);
}

	return 0;
}