Language: C++


#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cstring>
#include <vector>
#include <algorithm>
#include <set>
#include <queue>
#include <map>
#include <climits>
#include <complex>
#include <numeric>
using namespace std;

#define REP(i,n) for(int i=0; i<(int)(n); i++)
#define RREP(i,n) for(int i=(int)n-1; i>=0; i--)
#define FOR(i,c) for(__typeof((c).begin())i=(c).begin();i!=(c).end();++i)
#define RFOR(i,c) for(__typeof((c).rbegin())i=(c).rbegin();i!=(c).rend();++i)
#define ALL(c) (c).begin(), (c).end()

typedef long long int ll;
typedef pair<int, int> pii;
typedef pair<int, pair<int, int> > pipii;
typedef vector<int> vi;

const int INF = 1e9;
const int MOD = 1e9+7;

typedef complex<double> Comp;
typedef vector<Comp> vc;
void fft(vc &f, bool inv = false){
    int n = f.size();
    for(int i = 0, j = 1; j < n - 1; j++){
        for(int k = n>>1; k > (i ^= k); k>>=1);
        if(i > j) swap(f[i], f[j]);
    }
    for(int m = 2; m <= n; m<<=1){
        double deg = (inv ? -1 : 1) * 2 * M_PI / m;
        Comp zeta = Comp(cos(deg), sin(deg));
        for(int i = 0; i < n; i += m){   
            Comp pow_zeta = Comp(1.0, 0.0);
            for(int j1 = i, j2 = i + m / 2; j2 < i + m; j1++, j2++){
                Comp t1 = f[j1], t2 = pow_zeta * f[j2];
                f[j1] = t1 + t2, f[j2] = t1 - t2;
                pow_zeta *= zeta;
            }
        }
    }
    return;
}
void multiply(vector<vc> &g, vector<vc> &h){
    int n = 1;
    while(n < g.size() + h.size() + 1) n<<=1;
    g.resize(n, vc(n, Comp(0.0, 0.0)));
    h.resize(n, vc(n, Comp(0.0, 0.0)));
    REP(i, n){
        g[i].resize(n, Comp(0.0, 0.0));
        h[i].resize(n, Comp(0.0, 0.0));
    }
    REP(i, n){
        fft(g[i]); fft(h[i]);
    }
    vc x(n), y(n);
    REP(j, n){
        REP(i, n){
            x[i] = g[i][j];
            y[i] = h[i][j];
        }
        fft(x);
        fft(y);
        REP(i, n){
            g[i][j] = x[i];
            h[i][j] = y[i];
        }
    }
    REP(i, n){
        REP(j, n){
            g[i][j] *= h[i][j];
        }
    }
    //inverse
    REP(j, n){
        REP(i, n){
            x[i] = g[i][j];
        }
        fft(x, true);
        REP(i, n){
            g[i][j] = x[i];
        }
    }
    REP(i, n){
        fft(g[i], true);
    }
    REP(i, n){
        REP(j, n){
            g[i][j] /= n * n;
        }
    }
    return;
}
int n;
bool In(int x, int y, vector<vc> &g){
    if(x >= g.size() || y >= g.size()) return false;
    if(x < 0 || y < 0) return false;
    return true;
}

int main(void){
    cin >> n;
    vector<vc> g(n, vc(n)), g2(n, vc(n));
    ll num = 0;
    REP(i, n){
        REP(j, n){
            int x;
            scanf("%d", &x);
            g[i][j] = g2[n - i - 1][n - j - 1] = Comp((double)x, 0.0);
            num -= g[i][j];
        }
    }
    multiply(g, g2);
    
    //map<ll, ll> ans;
    vector<ll> ans(3 * n * n, 0);
    double dist = 0.0;
    ans[0] = num;
    REP(i, g.size()){
        REP(j, g[i].size()){
            ll dn;
            int x = abs(i - (n - 1)), y = abs(j - (n - 1));
            if((dn = real(g[i][j] + 0.01)) > 0){
                ll at = (ll)(x * x + y * y);
                num += dn;
                dist += dn * sqrt((double)at);
                ans[at] += (ll)(real(g[i][j]) + 0.01);
            }
        }
    }
    printf("%.10lf\n", dist / num);
    num = 0;
    REP(i, ans.size()){
        if(ans[i] / 2  == 0) continue;
        if(num++ == 10000) break;
        cout << i << " " << ans[i] / 2  << endl;
    }
    return 0;
}

Compiler output:
a.cc: In function 'int main()':
a.cc:119:17: error: no match for 'operator-=' (operand types are 'll' {aka 'long long int'} and '__gnu_cxx::__alloc_traits<std::allocator<std::complex<double> >, std::complex<double> >::value_type' {aka 'std::complex<double>'})
  119 |             num -= g[i][j];