高次类欧几里得算法
题意
$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}]$，可以优化常数。
- 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; }
近期评论