題目描述：

有一個N×N的方陣，第x行第y列有 $C_{x,y}$個點( $0$

≤ C x , y ≤ 9 0\le C_{x,y}\le9 )。
任選兩個不同的點，求兩點歐幾里德距離的均值(或期望)。
然後按距離從小到大輸出該距離的平方di和對應的點對數目ci。

題目分析：

想辦法把兩點之間的距離用橫縱座標的線性組合表示成一個值，使得可以把兩點值相減之後的值還原出對應的橫縱座標之差（類似於一個加密解密的過程），然後就可以FFT卷積了
容易想到使用每個點的編號作為加密規則，假設有兩個點 $(i,$

j ) , ( x , y )   ( i , j , x , y ∈ [ 0 , n − 1 ] ) (i,j),(x,y)~(i,j,x,y\in[0,n-1]) ，那麼他們的差值為：
$(i*n+j)-(x*n+y)=(i-x)*n+(j-y)$
這時候你會想，這個值模一下n不就得到縱座標之差了嗎？！
實際上是有問題的，因為此時的 $(i-x),(j-y)$不一定都是正數
舉個例子，當n=3時，i-x=1, j-y=1, 和 i-x=2, j-y=-2 對應的是同一個值，無法判斷
也就是說，i-x的變化被j-y的變化干擾了。
那麼就不要讓它干擾！ 把 i*n+j 改作 i*2n+j
此時再看： $(i*2n+j)-(x*2n+y)=(i-x)*2n+(j-y)$
再模上2n，j-y是正是負一目瞭然，記模出來的值為x,若 $x\in(-2n,-n]\cup(0,n-1]$則為j-y正，反之則為負
然後把對應的編號轉成次冪形式，多項式A為 $x^{id}$，多項式B為 $x^{-id}$，做FFT，由冪（編號之差）反推橫縱座標之差，乘上對應係數即可
ps:上述過程有點像雜湊似的呢。。。有解密的味道。。。

• 記得用long long! 點的數量是百萬級別，不是1024…
• 0要單獨處理

#include<cstdio>
#include<cmath>
#define LL long long
#include<algorithm>
using namespace std;
const int maxn = 1<<23;
const double Pi = acos(-1);
struct complex
{
	double r,i;
	complex(double _r=0,double _i=0):r(_r),i(_i){}
	complex operator + (const complex &t)const{return complex(r+t.r,i+t.i);}
	complex operator - (const complex &t)const{return complex(r-t.r,i-t.i);}
	complex operator * (const complex &t)const{return complex(r*t.r-i*t.i,r*t.i+i*t.r);}
}a1[maxn],a2[maxn],wn,w;
void change(complex *a,int len)
{
	for(int i=1,j=len/2,k;i<len-1;i++)
	{
		if(i<j) swap(a[i],a[j]);
		for(k=len/2;j>=k;j-=k,k>>=1);
		j+=k;
	}
}
void fft(complex *a,int len,int flg)
{
	change(a,len);
	for(int i=2;i<=len;i<<=1)
	{
		wn=complex(cos(2*Pi*flg/i),sin(2*Pi*flg/i));
		for(int j=0;j<len;j+=i)
		{
			w=complex(1,0);
			for(int k=j;k<j+i/2;k++)
			{
				complex u=a[k],v=w*a[k+i/2];
				a[k]=u+v,a[k+i/2]=u-v;
				w=w*wn;
			}
		}
	}
	if(flg==-1) for(int i=0;i<len;i++) a[i].r/=len;
}
int n,reset,x,sum;
LL cnt[2100000];
double ans;
int main()
{
	freopen("1.in","r",stdin);
	freopen("1.out","w",stdout);
	scanf("%d",&n);reset=(2*n+1)*(n-1);
	int len=1;while(len<2*reset+1) len<<=1;
	for(int i=0;i<n;i++)
		for(int j=0;j<n;j++)
		{
			scanf("%d",&x);sum+=x;
			a1[i*2*n+j].r=x;
			a2[reset-i*2*n-j].r=x;
			cnt[0]+=x*(x-1)/2;
		}
	fft(a1,len,1),fft(a2,len,1);
	for(int i=0;i<len;i++) a2[i]=a1[i]*a2[i];
	fft(a2,len,-1);
	for(int i=0,x,y;i<len;i++)
	{
		LL t=(LL)(a2[i].r+0.5);
		if(!t) continue;
		y=(i-reset)%(2*n);
		if(y<=-n) y+=2*n;
		if(y>=n) y-=2*n;
		x=(i-reset-y)/(2*n);
		x=x*x+y*y;
		if(x) cnt[x]+=t;
		ans+=t*sqrt(x);
	}
	printf("%.10lf\n",ans/(1ll*sum*(sum-1)));
	int num=0,top=(n-1)*(n-1)*2;
	if(cnt[0]) printf("0 %lld\n",cnt[0]),num++;
	for(int i=1;i<=top&&num<10000;i++)
		if(cnt[i]) printf("%d %lld\n",i,cnt[i]/2),num++;
}