Aizu2560 Point Distance 【FFT】

题目描述：

有一个N×N的方阵，第x行第y列有$C_{x,y}$个点($0\le C_{x,y}\le 9$)。
任选两个不同的点，求两点欧几里德距离的均值(或期望)。
然后按距离从小到大输出该距离的平方di和对应的点对数目ci。

题目分析：

想办法把两点之间的距离用横纵坐标的线性组合表示成一个值，使得可以把两点值相减之后的值还原出对应的横纵坐标之差（类似于一个加密解密的过程），然后就可以FFT卷积了
容易想到使用每个点的编号作为加密规则，假设有两个点$(i,j),(x,y)(i,j,x,y\in [0,n-1])$，那么他们的差值为：
$$(i\ast n+j)-(x\ast n+y)=(i-x)\ast 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\ast 2n+j)-(x\ast 2n+y)=(i-x)\ast 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++;
}