QOJ-12606 题解
题意
给定四个整数 $A,B,C,D$，求有多少个不同的四元组 $(a,b,c,d)$ 满足：
$$ 0\le a\le A,\quad 0\le b\le B,\quad 0\le c\le C,\quad 0\le d\le D $$
且
$$ a+c>b+d,\qquad a+d\ge b+c $$
答案对 $10^9+7$ 取模。
核心转化
令
$$ u=a-b,\qquad v=c-d $$
则两个条件等价于：
$$ u+v>0,\qquad u\ge v $$
如果 $u\le 0$，由 $v\le u$ 可得 $u+v\le 0$，矛盾，所以必须有：
$$ u>0 $$
于是对于每个固定的 $u>0$，$v$ 需要满足：
$ - u < v \le u $
因此答案为：
$$ \sum_{u=1}^{A} f(u)\cdot H(u) $$
其中：
- $f(u)$ 表示满足 $a-b=u$ 的 $(a,b)$ 对数；
- $H(u)$ 表示满足 $-u < c-d \le u$ 的 $(c,d)$ 对数。
计算 $f(u)$
对于 $u\ge 1$：
$$ f(u)=\#\{(a,b):0\le a\le A,\ 0\le b\le B,\ a-b=u\} $$
可得：
$$ f(u)=\min(A-u,B)+1 $$
它是一个分段函数：
- 若 $A>B$：
- 当 $1\le u\le A-B$ 时，$f(u)=B+1$，为常数；
- 当 $A-B < u \le A$ 时，$f(u)=A-u+1$，为一次函数。
- 若 $A\le B$：
- 当 $1 \le u \le A$ 时，$f(u)=A-u+1$，为一次函数。
计算 $H(u)$
设总对数：
$$ N=(C+1)(D+1) $$
用容斥计算 $H(u)$：
$$ H(u)=N-Q(C,D,u+1)-Q(D,C,u) $$
其中
$$ Q(M,N,k)=\#\{(x,y):0\le x\le M,\ 0\le y\le N,\ x-y\ge k\} $$
这里：
- $Q(C,D,u+1)$ 排除 $c-d>u$ 的情况；
- $Q(D,C,u)$ 排除 $c-d\le -u$ 的情况。
计算 $Q(M,N,k)$
当 $k>M$ 时，$Q=0$。
否则令：
$$ L=M-k+1 $$
有：
$$ Q(M,N,k)=\sum_{y=0}^{\min(N,L)}(L-y) $$
化简为：
$$ Q(M,N,k)=(\min(N,L)+1)L-\frac{\min(N,L)(\min(N,L)+1)}2 $$
这是一个关于 $k$ 的分段函数：
- 当 $L\ge N$，即 $k\le M-N+1$ 时，$\min(N,L)=N$，$Q$ 是 $k$ 的一次函数；
- 当 $LM-N+1$ 时，$\min(N,L)=L$，$Q$ 是 $k$ 的二次函数。
由于本题中 $k$ 是 $u+1$ 或 $u$，所以 $Q$ 也是关于 $u$ 的低次多项式。
多项式求和
$f(u)$ 是零次或一次多项式，$H(u)$ 是零次、一次或二次多项式，因此：
$$ f(u)\cdot H(u) $$
在每一段上是次数不超过 $3$ 的多项式。
把所有分段点（包括 $f$ 的断点、两个 $Q$ 函数的线性/二次分界点、归零点，以及 $1$ 和 $A+1$）排序，相邻两个断点之间是一段。对每一段直接使用幂和公式：
$$ \sum_{i=1}^{n}1,\quad \sum_{i=1}^{n}i,\quad \sum_{i=1}^{n}i^2,\quad \sum_{i=1}^{n}i^3 $$
即可 $O(1)$ 求出这一段的和。
所有计算均在模 $10^9+7$ 下进行。
复杂度
每组数据只需要常数个分段，因此：
$$ O(1)\ \text{per test case} $$
总复杂度为 $O(T)$，可以轻松通过 $T\le 1000,\ A,B,C,D\le 10^{18}$。
参考代码
#include <bits/stdc++.h>
using namespace std;
using int64 = long long;
using i128 = __int128_t;
const int64 MOD = 1000000007LL;
const int64 INV2 = 500000004LL;
const int64 INV6 = 166666668LL;
int64 md_norm(int64 x) {
x %= MOD;
if (x < 0) x += MOD;
return x;
}
int64 md_add(int64 a, int64 b) { return md_norm(a + b); }
int64 md_sub(int64 a, int64 b) { return md_norm(a - b); }
int64 md_mul(int64 a, int64 b) { return (int64)((i128)a * b % MOD); }
struct Poly {
int64 c[4];
Poly(int64 a = 0, int64 b = 0, int64 c2 = 0, int64 c3 = 0) {
c[0] = md_norm(a);
c[1] = md_norm(b);
c[2] = md_norm(c2);
c[3] = md_norm(c3);
}
};
Poly poly_add(const Poly &a, const Poly &b) {
Poly r;
for (int i = 0; i < 4; ++i) r.c[i] = md_add(a.c[i], b.c[i]);
return r;
}
Poly poly_sub(const Poly &a, const Poly &b) {
Poly r;
for (int i = 0; i < 4; ++i) r.c[i] = md_sub(a.c[i], b.c[i]);
return r;
}
Poly poly_mul(const Poly &a, const Poly &b) {
Poly r;
for (int i = 0; i < 4; ++i) {
if (!a.c[i]) continue;
for (int j = 0; j < 4; ++j) {
if (!b.c[j]) continue;
if (i + j < 4) r.c[i + j] = md_add(r.c[i + j], md_mul(a.c[i], b.c[j]));
}
}
return r;
}
int64 sum1(int64 n) {
if (n <= 0) return 0;
return md_mul(md_mul(md_norm(n), md_norm(n + 1)), INV2);
}
int64 sum2(int64 n) {
if (n <= 0) return 0;
return md_mul(md_mul(md_mul(md_norm(n), md_norm(n + 1)), md_norm(2 * n + 1)), INV6);
}
int64 sum3(int64 n) {
if (n <= 0) return 0;
int64 s = sum1(n);
return md_mul(s, s);
}
int64 poly_sum(const Poly &p, int64 n) {
if (n <= 0) return 0;
int64 res = 0;
res = md_add(res, md_mul(p.c[0], md_norm(n)));
res = md_add(res, md_mul(p.c[1], sum1(n)));
res = md_add(res, md_mul(p.c[2], sum2(n)));
res = md_add(res, md_mul(p.c[3], sum3(n)));
return res;
}
int64 range_sum(const Poly &p, int64 L, int64 R) {
if (L > R) return 0;
return md_sub(poly_sum(p, R), poly_sum(p, L - 1));
}
// Q(M,N,k) = # {(x,y): 0<=x<=M, 0<=y<=N, x-y >= k}, k >= 1.
// Returns it as a polynomial in u, where k = u + off.
Poly get_Q_poly(int64 M, int64 N, int64 off, int64 u) {
int64 k = u + off;
if (k <= 0 || k > M) return Poly(0);
if (M >= N && k <= M - N + 1) {
// linear part: (N+1)*(M-k+1) - N*(N+1)/2
int64 c1 = md_norm(-(N + 1));
int64 c0 = md_sub(md_mul(md_norm(N + 1), md_norm(M - off + 1)),
md_mul(md_mul(md_norm(N), md_norm(N + 1)), INV2));
return Poly(c0, c1);
} else {
// quadratic part: (M-k+1)*(M-k+2)/2, let s=M-off+1, then it is (s-u)(s-u+1)/2
int64 s = M - off + 1;
int64 c2 = INV2;
int64 c1 = md_mul(md_norm(-(2 * s + 1)), INV2);
int64 c0 = md_mul(md_mul(md_norm(s), md_norm(s + 1)), INV2);
return Poly(c0, c1, c2);
}
}
// f(u) = # {(a,b): 0<=a<=A, 0<=b<=B, a-b=u}, u>=1.
Poly get_f_poly(int64 A, int64 B, int64 u) {
if (A > B && u <= A - B) return Poly(B + 1);
return Poly(A + 1, -1); // A - u + 1
}
int64 solve_one(int64 A, int64 B, int64 C, int64 D) {
if (A <= 0) return 0;
vector<int64> pts;
pts.push_back(1);
pts.push_back(A + 1);
// Breakpoint of f(u).
if (A > B && A - B + 1 <= A) pts.push_back(A - B + 1);
auto add_Q_breakpoints = [&](int64 M, int64 N, int64 off) {
if (M <= 0) return;
int64 zero_start = M - off + 1; // first u with Q = 0
if (zero_start >= 1 && zero_start <= A + 1) pts.push_back(zero_start);
if (M >= N) {
int64 b = M - N + 2 - off; // first u of the quadratic part
if (b >= 1 && b <= A + 1) pts.push_back(b);
}
};
// Breakpoints of H(u) = total - Q(C,D,u+1) - Q(D,C,u).
add_Q_breakpoints(C, D, 1);
add_Q_breakpoints(D, C, 0);
sort(pts.begin(), pts.end());
pts.erase(unique(pts.begin(), pts.end()), pts.end());
int64 ans = 0;
for (size_t i = 0; i + 1 < pts.size(); ++i) {
int64 L = pts[i];
int64 R = pts[i + 1] - 1;
if (L < 1) L = 1;
if (R > A) R = A;
if (L > R) continue;
Poly f = get_f_poly(A, B, L);
int64 total = (int64)((i128)(C + 1) * (D + 1) % MOD);
Poly H = poly_sub(poly_sub(Poly(total), get_Q_poly(C, D, 1, L)),
get_Q_poly(D, C, 0, L));
ans = md_add(ans, range_sum(poly_mul(f, H), L, R));
}
return ans;
}
int main() {
ios::sync_with_stdio(false);
cin.tie(nullptr);
int T;
cin >> T;
while (T--) {
int64 A, B, C, D;
cin >> A >> B >> C >> D;
cout << solve_one(A, B, C, D) << '\n';
}
return 0;
}