题解

题解

对于每个 $i$，令 $x_i=a_i+d_i,\ y_i=b_i+e_i,\ z_i=c_i+f_i$，其中 $d_i,e_i,f_i\ge0$ 且 $d_i+e_i+f_i\le r_i$。三维模和的计数可以用卷积完成，但模长分别是 $X,Y,Z$，直接 DP 需要 $\mathcal{O}(nXYZ^2)$，我们改为做三维离散傅里叶变换：卷积在频域变成逐点乘。

记 $\\omega_X=e^{2\\pi i/X}$（模下用原根构造），对频率 $(p,q,r)$ 有

$$\\widehat{F_i}(p,q,r)=\\omega_X^{pa_i}\\omega_Y^{qb_i}\\omega_Z^{rc_i}\\sum_{d+e+f\\le r_i}(\\omega_X^p)^d(\\omega_Y^q)^e(\\omega_Z^r)^f . $$

核心是快速求截断几何级数

$$S_r(\\alpha,\\beta,\\gamma)=\\sum_{d+e+f\\le r}\\alpha^d\\beta^e\\gamma^f. $$

• 当 $\\alpha,\\beta,\\gamma$ 互不相等时，部分分式可得

$$S_r=A\\cdot\\frac{1-\\alpha^{r+1}}{1-\\alpha}+B\\cdot\\frac{1-\\beta^{r+1}}{1-\\beta}+C\\cdot\\frac{1-\\gamma^{r+1}}{1-\\gamma}, $$

其中 $A=\\frac1{(1-\\beta/\\alpha)(1-\\gamma/\\alpha)}$，其他类似。

• 只有频率为 $1$ 时会出现根相同的退化情形；此时用线性递推补足：\(T_k=\sum_{d+e+f=k}\alpha^d\beta^e\gamma^f\) 满足三阶线性递推，增广一维前缀和即可用 $4\\times4$ 矩阵快速幂在 $\\mathcal{O}(\\log r)$ 求出 $S_r=\\sum_{k=0}^rT_k$。
• 全为 $1$ 的频率只需直接用组合数：$\\binom{r+3}{3}$。

频域乘积：$\\widehat F(p,q,r)=\\prod_i \\widehat{F_i}(p,q,r)$。由于 $X\\cdot Y\\cdot Z\\le2000$，三维逆变换直接用定义式 $\\mathcal{O}((XYZ)^2)$ 完全可行，最后乘以 $(XYZ)^{-1}$ 得到 $F(u,v,w)$。

答案按题意异或累加：

$$\\bigoplus_{u,v,w}\\bigl((uYZ+vZ+w)\\cdot (F(u,v,w)\\bmod \\text{MOD})\\bigr). $$

复杂度：频域部分 $\\mathcal{O}(nXYZ)$，逆变换 $\\mathcal{O}((XYZ)^2)$（上界约四百万次运算），可轻松通过约束。内存 $\\mathcal{O}(XYZ)$。

#include <bits/stdc++.h>
using namespace std;

constexpr int MOD = 466560001; // prime, 2^10 * 3^6 * 5^4 + 1
using i64 = long long;

int addmod(int a, int b) { int s = a + b; if (s >= MOD) s -= MOD; return s; }
int submod(int a, int b) { int s = a - b; if (s < 0) s += MOD; return s; }
int mulmod(i64 a, i64 b) { return int(a * b % MOD); }
int mod_pow(int a, i64 e) { int r = 1; while (e) { if (e & 1) r = mulmod(r, a); a = mulmod(a, a); e >>= 1; } return r; }

int geosum(int base, i64 r) {
    if (base == 1) return int((r + 1) % MOD);
    int num = submod(mod_pow(base, r + 1), 1);
    int den = submod(base, 1);
    return mulmod(num, mod_pow(den, MOD - 2));
}

struct Matrix4 {
    int a[4][4];
    Matrix4(bool id = false) { for (int i = 0; i < 4; ++i) for (int j = 0; j < 4; ++j) a[i][j] = (id && i == j) ? 1 : 0; }
};
Matrix4 mul(const Matrix4 &x, const Matrix4 &y) {
    Matrix4 r;
    for (int i = 0; i < 4; ++i) for (int k = 0; k < 4; ++k) if (x.a[i][k])
        for (int j = 0; j < 4; ++j) r.a[i][j] = addmod(r.a[i][j], mulmod(x.a[i][k], y.a[k][j]));
    return r;
}
array<int, 4> mat_pow_apply(Matrix4 base, i64 e, array<int, 4> v) {
    Matrix4 res(true);
    while (e) { if (e & 1) res = mul(res, base); base = mul(base, base); e >>= 1; }
    array<int, 4> r{};
    for (int i = 0; i < 4; ++i) for (int j = 0; j < 4; ++j) r[i] = addmod(r[i], mulmod(res.a[i][j], v[j]));
    return r;
}

int sum_leq_r_recurrence(int alpha, int beta, int gamma, i64 r) {
    int c1 = addmod(addmod(alpha, beta), gamma);
    int c2 = addmod(mulmod(alpha, beta), addmod(mulmod(alpha, gamma), mulmod(beta, gamma)));
    int c3 = mulmod(mulmod(alpha, beta), gamma);
    Matrix4 trans(false);
    trans.a[0][0] = c1; trans.a[0][1] = submod(0, c2); trans.a[0][2] = c3;
    trans.a[1][0] = 1;  trans.a[2][1] = 1;
    trans.a[3][0] = c1; trans.a[3][1] = submod(0, c2); trans.a[3][2] = c3; trans.a[3][3] = 1;
    array<int, 4> v0{1, 0, 0, 1}; // [T0, T-1, T-2, S0]
    return mat_pow_apply(trans, r, v0)[3];
}

int sum_leq_r(int alpha, int beta, int gamma, i64 r) {
    if (alpha == beta || alpha == gamma || beta == gamma)
        return sum_leq_r_recurrence(alpha, beta, gamma, r);
    int inv_alpha = mod_pow(alpha, MOD - 2), inv_beta = mod_pow(beta, MOD - 2), inv_gamma = mod_pow(gamma, MOD - 2);
    int A = mod_pow(mulmod(submod(1, mulmod(beta, inv_alpha)), submod(1, mulmod(gamma, inv_alpha))), MOD - 2);
    int B = mod_pow(mulmod(submod(1, mulmod(alpha, inv_beta)), submod(1, mulmod(gamma, inv_beta))), MOD - 2);
    int C = mod_pow(mulmod(submod(1, mulmod(alpha, inv_gamma)), submod(1, mulmod(beta, inv_gamma))), MOD - 2);
    int res = 0;
    res = addmod(res, mulmod(A, geosum(alpha, r)));
    res = addmod(res, mulmod(B, geosum(beta, r)));
    res = addmod(res, mulmod(C, geosum(gamma, r)));
    return res;
}

int main() {
    ios::sync_with_stdio(false); cin.tie(nullptr);
    int n, X, Y, Z; if (!(cin >> n >> X >> Y >> Z)) return 0;
    vector<int> a(n), b(n), c(n); vector<i64> r(n);
    for (int i = 0; i < n; ++i) { i64 ai, bi, ci, ri; cin >> ai >> bi >> ci >> ri; a[i] = ai % X; b[i] = bi % Y; c[i] = ci % Z; r[i] = ri; }
    int root = 13;
    int wX = mod_pow(root, (MOD - 1) / X), wY = mod_pow(root, (MOD - 1) / Y), wZ = mod_pow(root, (MOD - 1) / Z);
    vector<int> powX(X, 1), powY(Y, 1), powZ(Z, 1);
    for (int i = 1; i < X; ++i) powX[i] = mulmod(powX[i - 1], wX);
    for (int i = 1; i < Y; ++i) powY[i] = mulmod(powY[i - 1], wY);
    for (int i = 1; i < Z; ++i) powZ[i] = mulmod(powZ[i - 1], wZ);
    vector<vector<int>> powXtab(X, vector<int>(X, 1)), powYtab(Y, vector<int>(Y, 1)), powZtab(Z, vector<int>(Z, 1));
    for (int p = 0; p < X; ++p) for (int k = 1; k < X; ++k) powXtab[p][k] = mulmod(powXtab[p][k - 1], powX[p]);
    for (int q = 0; q < Y; ++q) for (int k = 1; k < Y; ++k) powYtab[q][k] = mulmod(powYtab[q][k - 1], powY[q]);
    for (int s = 0; s < Z; ++s) for (int k = 1; k < Z; ++k) powZtab[s][k] = mulmod(powZtab[s][k - 1], powZ[s]);

    auto idx = [Y, Z](int p, int q, int s) { return (p * Y + q) * Z + s; };
    int totalSize = X * Y * Z;
    vector<int> G(totalSize, 1);
    for (int p = 0; p < X; ++p) for (int q = 0; q < Y; ++q) for (int s = 0; s < Z; ++s) {
        int alpha = powX[p], beta = powY[q], gamma = powZ[s];
        int prod = 1;
        for (int i = 0; i < n; ++i) {
            int part = mulmod(powXtab[p][a[i]], mulmod(powYtab[q][b[i]], powZtab[s][c[i]]));
            int ways = (alpha == 1 && beta == 1 && gamma == 1)
                       ? mulmod(int(((r[i] + 1) % MOD) * ((r[i] + 2) % MOD) % MOD * ((r[i] + 3) % MOD) % MOD), mod_pow(6, MOD - 2))
                       : sum_leq_r(alpha, beta, gamma, r[i]);
            prod = mulmod(prod, mulmod(part, ways));
        }
        G[idx(p, q, s)] = prod;
    }

    vector<vector<int>> invPowX(X, vector<int>(X)), invPowY(Y, vector<int>(Y)), invPowZ(Z, vector<int>(Z));
    for (int u = 0; u < X; ++u) for (int p = 0; p < X; ++p) {
        int exp = (X - (i64)p * u % X) % X; invPowX[u][p] = powXtab[p][exp];
    }
    for (int v = 0; v < Y; ++v) for (int q = 0; q < Y; ++q) {
        int exp = (Y - (i64)q * v % Y) % Y; invPowY[v][q] = powYtab[q][exp];
    }
    for (int w = 0; w < Z; ++w) for (int s = 0; s < Z; ++s) {
        int exp = (Z - (i64)s * w % Z) % Z; invPowZ[w][s] = powZtab[s][exp];
    }

    int invXYZ = mod_pow(totalSize % MOD, MOD - 2);
    vector<int> fvals(totalSize, 0);
    for (int u = 0; u < X; ++u) for (int v = 0; v < Y; ++v) for (int w = 0; w < Z; ++w) {
        int sum = 0;
        for (int p = 0; p < X; ++p) for (int q = 0; q < Y; ++q) for (int s = 0; s < Z; ++s) {
            int term = mulmod(G[idx(p, q, s)], mulmod(invPowX[u][p], mulmod(invPowY[v][q], invPowZ[w][s])));
            sum = addmod(sum, term);
        }
        fvals[idx(u, v, w)] = mulmod(sum, invXYZ);
    }

    long long ans = 0;
    for (int u = 0; u < X; ++u) for (int v = 0; v < Y; ++v) for (int w = 0; w < Z; ++w) {
        long long coef = 1LL * (u * Y * Z + v * Z + w);
        ans ^= (coef * fvals[idx(u, v, w)]) % MOD;
    }
    cout << ans << '\\n';
    return 0;
}