Aizu 2560 Point Distance FFT

68 篇文章 0 订阅
37 篇文章 1 订阅

一般的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--;
		}
}
  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值