gxx_slide之Point Distance

链接：http://acm.bnu.edu.cn/v3/problem_show.php?pid=39676
题意：N × N 的点阵,(x, y) 位置有 C x,y(<=9) 个点,考虑所有点对,把点对
按照 Euclidean 距离从小到大输出。
(N ≤ 1024)
分析：感觉ppt上讲的太简洁了一点。。。。
首先考虑将二维坐标一维化，就是i,j对应到i*n+j,两点间的距离就是两个一维坐标的差值再对映回一维,比如考虑
$012\\345\\678\\$
0号点和8号点的水平距离为2,竖直距离为2,就等于((8-0)/3,(8-0)%3),即一维坐标的差值能确定一个二维距离，但是考虑2号点和6号点6-2=4,4/3=1,4%3=1,而2和6显然不是(1,1)的关系，为什么会这样呢，因为跨过了一行，但是直接这样处理我们不知道这个差值是否跨行；因此我们可以考虑像论文里说的那样在右边补一片n*n的全零的元素，这样，当产生跨行的那种情况时，余数必定>=n,然后再做一次傅立叶，就能通过这道题了。

#include<stdio.h>
#include<algorithm>
#include<math.h>
using namespace std;
typedef long long Int;
const double pi=acos(-1.0);
const int Maxn=1024*2048*2;
struct cp
{
    double a,b;
    cp operator+(const cp &x)const{return (cp){a+x.a,b+x.b};}
    cp operator-(const cp &x)const{return (cp){a-x.a,b-x.b};}
    cp operator*(const cp &x)const{return (cp){a*x.a-b*x.b,a*x.b+b*x.a};}
}a[Maxn],b[Maxn];
Int rep[Maxn];
void rev(cp*a,int n)
{
    int k;
    for(int i=1,j=n>>1;i<n-1;i++)
    {
        if(i<j)swap(a[i],a[j]);
        for(k=n>>1;j>=k;j-=k,k>>=1);j+=k;
    }
}
void dft(cp*a,int n,int flag=1)
{
    rev(a,n);
    for(int m=2;m<=n;m<<=1)
    {
        cp wm=(cp){cos(flag*pi*2/m),sin(flag*pi*2/m)};
        for(int k=0;k<n;k+=m)
        {
            cp w=(cp){1.0,0.0};
            for(int j=k;j<k+(m>>1);w=w*wm,j++)
            {
                cp u=a[j],v=a[j+(m>>1)]*w;
                a[j]=u+v;
                a[j+(m>>1)]=u-v;
            }
        }
    }
}
void mul(cp*a,cp*b,int n)
{
    dft(a,n);dft(b,n);
    for(int i=0;i<n;i++)a[i]=a[i]*b[i];
    dft(a,n,-1);
    for(int i=0;i<n;i++)a[i].a/=n;
}
int main()
{
    int n;
    while(scanf("%d",&n)!=EOF)
    {
        int N=1;
        while(N<n*n*2)N<<=1;N<<=1;
        for(int i=0;i<N;i++)a[i]=b[i]=(cp){0.0,0.0};
        Int sum0=0,total=0;
        for(int i=0;i<n;i++)
            for(int j=0;j<n;j++)
            {
                int x;scanf("%d",&x);
                total+=x;
                sum0+=x*(x-1)/2;
                a[i*(n+n)+j].a+=x;
                b[n*n*2-1-i*(n+n)-j].a+=x;
            }
        mul(a,b,N);
        for(int i=0;i<=2*(n-1)*(n-1);i++)rep[i]=0;
        for(int i=n*n*2;i<N;i++)
        {
            Int val=(Int)(a[i].a+0.5);
            if(!val)continue;
            int t1=(i-(n*n*2-1))/(n+n),t2=(i-(n*n*2-1))%(n+n);
            if(t2>=n){t2=n+n-t2;t1++;}
            rep[t1*t1+t2*t2]+=val;
        }
        rep[0]=sum0;
        double ans=0;
        for(int i=0;i<=2*(n-1)*(n-1);i++)ans+=sqrt(i+0.0)*rep[i];
        printf("%.12f\n",ans/(total*(total-1)/2));
        for(int i=0,j=0;i<=2*(n-1)*(n-1)&&j<10000;i++)
            if(rep[i])
            {
                printf("%d %lld\n",i,rep[i]);
                j++;
            }
    }
}