JAG 2013 Spring F 二维卷积

题意

传送门 JAG 2013 Spring F Point Distance

题解

将距离看作二元组 $(x,y)$，此时行列指标的差都是线性的。对于 $(0,0)$，数量为 $\sum c_{i,j}(c_{i,j}-1)/2$。

对于 $(x,0),(0,y)$，枚举行列分别做一维卷积。以 $(0,y)$ 为例，枚举行 $i$，求 $\sum c_{i,j}{c}_{i,j+y}$，行指标相同，转化为 $\sum a_i{b}_{i+k}$ 的形式；对于这类形式的和，可以将 $b$ 反转，即令 $j=n-1-(i+k)$，表达式转化为 $\sum a_i{b}_j[i+j=n-1-k]$。

对于 $(x,y)$，考虑将其映射为值后，转化为一维卷积。卷积时要保证 $i+x\ge n,j+y\ge n$ 的情况贡献为 $0$，那么映射函数为 $f(x,y)=x\times 2n+y$。卷积后得到距离为 $(x,y),x\ge 0,y\ge 0$ 的有序点对数量；将 $c$ 按行反转后，即可得到 $x,y$ 异号的有序点对数量。

总时间复杂度 $O(n^2\log n)$。

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
constexpr ll MOD = 998244353, PRT = 3;
ll qpow(ll x, ll n)
{
    ll res = 1;
    while (n > 0)
    {
        if (n & 1)
            res = res * x % MOD;
        x = x * x % MOD, n >>= 1;
    }
    return res;
}
vector<int> rev;
struct Poly : vector<ll>
{
    Poly() {}
    Poly(int n) : vector<ll>(n) {}
    Poly(const initializer_list<ll> &list) : vector<ll>(list) {}
    void fft(int n, bool inverse)
    {
        if ((int)rev.size() != n)
        {
            rev.resize(n);
            for (int i = 0; i < n; ++i)
                rev[i] = rev[i >> 1] >> 1 | (i & 1 ? n >> 1 : 0);
        }
        resize(n);
        for (int i = 0; i < n; ++i)
            if (i < rev[i])
                std::swap(at(i), at(rev[i]));

        for (int m = 1; m < n; m <<= 1)
        {
            int m2 = m << 1;
            ll _w = qpow(inverse ? qpow(PRT, MOD - 2) : PRT, (MOD - 1) / m2);
            for (int i = 0; i < n; i += m2)
                for (int w = 1, j = 0; j < m; ++j, w = w * _w % MOD)
                {
                    ll &x = at(i + j), &y = at(i + j + m), t = w * y % MOD;
                    y = x - t;
                    if (y < 0)
                        y += MOD;
                    x += t;
                    if (x >= MOD)
                        x -= MOD;
                }
        }
    }
    void dft(int n) { fft(n, 0); };
    void idft(int n)
    {
        fft(n, 1);
        for (int i = 0, inv = qpow(n, MOD - 2); i < n; ++i)
            at(i) = at(i) * inv % MOD;
    }
    Poly operator*(const Poly &p) const
    {
        auto a = *this, b = p;
        int k = 1, n = a.size() + b.size() - 1;
        while (k < n)
            k <<= 1;
        a.dft(k), b.dft(k);
        for (int i = 0; i < k; ++i)
            a[i] = a[i] * b[i] % MOD;
        a.idft(k);
        a.resize(n);
        return a;
    }
};
constexpr int MAXN = 1050;
int N, C[MAXN][MAXN];
ll cnt[2 * MAXN * MAXN];

void count()
{
    int n = 2 * N * N;
    Poly f(n);
    for (int i = 0; i < N; ++i)
        for (int j = 0; j < N; ++j)
            f[i * 2 * N + j] = C[i][j];
    Poly g = f;
    reverse(g.begin(), g.end());
    f = f * g;
    for (int i = 1; i < N; ++i)
        for (int j = 1; j < N; ++j)
            cnt[i * i + j * j] += f[n - 1 - (i * 2 * N + j)];
}

int main()
{
    ios::sync_with_stdio(0), cin.tie(0), cout.tie(0);
    cin >> N;
    for (int i = 0; i < N; ++i)
        for (int j = 0; j < N; ++j)
            cin >> C[i][j];
    for (int i = 0; i < N; ++i)
        for (int j = 0; j < N; ++j)
            cnt[0] += C[i][j] * (C[i][j] - 1) / 2;
    for (int i = 0; i < N; ++i)
    {
        Poly f(N);
        for (int j = 0; j < N; ++j)
            f[j] = C[i][j];
        Poly g = f;
        reverse(g.begin(), g.end());
        f = f * g;
        for (int j = 1; j < N; ++j)
            cnt[j * j] += f[N - 1 - j];
    }
    for (int j = 0; j < N; ++j)
    {
        Poly f(N);
        for (int i = 0; i < N; ++i)
            f[i] = C[i][j];
        Poly g = f;
        reverse(g.begin(), g.end());
        f = f * g;
        for (int i = 1; i < N; ++i)
            cnt[i * i] += f[N - 1 - i];
    }
    count();
    for (int i = 0; i < N; ++i)
        reverse(C[i], C[i] + N);
    count();
    vector<pair<int, ll>> res;
    int lim = 10000;
    long double sum = 0;
    ll num = 0;
    for (int i = 0; i < 2 * N * N; ++i)
    {
        if (cnt[i] > 0 && (int)res.size() < lim)
            res.push_back({i, cnt[i]});
        sum += sqrt(i) * cnt[i], num += cnt[i];
    }
    cout << fixed << setprecision(9) << sum / num << '\n';
    for (auto &p : res)
        cout << p.first << ' ' << p.second << '\n';
    return 0;
}