一般的FFT重点是找不变量,这里发现如果AB两点,我们把他们标号为 2 n x + y 2nx+y 2nx+y那么通过两个标号相减的值可以反推出横坐标差和纵坐标差,那么直接把所有点按标号放入多项式卷起来即可。
AC Code:
#include<bits/stdc++.h>
#define maxn (1<<23)
#define LL long long
using namespace std;
const double Pi = 3.1415926535897932384626433832795;
int n;
struct cplx
{
double r,i;
cplx(double r=0,double i=0):r(r),i(i){}
cplx operator +(const cplx &B)const{ return cplx(r+B.r,i+B.i); }
cplx operator -(const cplx &B)const{ return cplx(r-B.r,i-B.i); }
cplx operator *(const cplx &B)const{ return cplx(r*B.r-i*B.i,r*B.i+i*B.r); }
cplx conj()const{ return cplx(r,-i); }
}w[maxn]={1},A[maxn],B[maxn];
int r[maxn];
LL arr[maxn];
inline void FFT(cplx A[maxn],int lgn,int tp)
{
int n = 1<<lgn;
if(tp==1) for(int i=1;i<n;i++) r[i] = (r[i>>1]>>1) + ((i&1)<<(lgn-1));
for(int i=1;i<n;i++) if(i < r[i]) swap(A[i],A[r[i]]);
for(int L=2;L<=n;L<<=1)
{
int l = L>>1;w[1] = cplx(cos(Pi/l),sin(Pi/l)*tp);
for(int i=2;i<l;i++) w[i] = w[i-1] * w[1];
for(int st=0;st<=n;st+=L)
for(int k=0;k<l;k++)
{
cplx tmp = w[k] * A[st + k + l];
A[st + k + l] = A[st + k] - tmp , A[st + k] = tmp + A[st + k];
}
}
if(tp==-1)
for(int i=0;i<n;i++) A[i].r/=n;
}
int main()
{
scanf("%d",&n);
int lgn = 0;
for(;4*n*n>(1<<lgn);lgn++);
int len = 1<<lgn;
LL sumcnt = 0;
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
{
int tmp;scanf("%d",&tmp);
arr[0] += tmp * (tmp - 1) / 2;
sumcnt += tmp;
A[2*n*i+j].r = tmp,A[len/2-2*n*i-j].i=tmp;
}
FFT(A,lgn,1);
for(int i=0;i<len;i++)
{
cplx a = A[i] , b = A[(len-i) & (len-1)].conj();
B[i] = (a+b)*(a-b)*cplx(0,-0.25);
}
FFT(B,lgn,-1);
double ans = 0 , sum = 0;
for(int i=len/2+1;i<len;i++)
{
int dx = (i-len/2) % (2 * n) > n ? (i-len/2) % (2 * n) - 2 * n :(i-len/2) % (2 * n),
dy = (i-len/2-dx) / (2 * n);
LL val = (LL)(B[i].r+0.5);
ans += sqrt(dx*dx+dy*dy) * val;
arr[dx*dx+dy*dy] += val;
}
printf("%.10lf\n",ans/(1ll*(sumcnt)*(sumcnt-1)/2 ));
for(int i=0,k=10000;i<=n*n*2 && k;i++)
if(arr[i] > 0)
{
printf("%d %lld\n",i,arr[i]);
k--;
}
}