表示敲完多项式乘法&高精度乘法两道FFT模板题后就开始来磕这题了
这题相对而言应该不算模板题 不过神犇们肯定还是一眼看穿
如果原OJ访问速度较慢的话可以考虑戳这里
http://acm.hust.edu.cn/vjudge/problem/viewProblem.action?id=153505
-------------------------------------------------------------------------------------------------------------
构造什么的 还是太巧妙了
比如对于这样的样例
2
1 2
3 4
可以构造出这样的矩阵（原矩阵在左上角 额外的矩阵在右下角 且两矩阵不相交且中心对称）
a矩阵 b矩阵
1 2 0 0 0 0 0 0
3 4 0 0 0 0 0 0
0 0 0 0 0 0 4 3
0 0 0 0 0 0 2 1
如果变为多项式的话 就是这样
多项式a 1 2 0 0 3 4 0 0 0 0 0 0 0 0 0 0
多项式b 0 0 0 0 0 0 0 0 0 0 4 3 0 0 2 1
当我们把它们乘起来后 会发现贡献到同一项点对的距离是一样的
当n不为2的整数幂的时候 可以像fft模板题一样 通过补零来进行构造
构造出这两个多项式后 我们就可以直接套fft模板 然后对结果多项式的有效的位数都访问一次即可统计出答案
-------------------------------------------------------------------------------------------------------------
不过A掉后还是有一个问题还没弄懂
设多项式a，b的长度为n 那么按照常规的多项式乘法 我们应该把他们后面都补上n个0
因为结果多项式的长度为n×2 （虽然实际有效位数比n还小）
然而我参考的某AC代码 并没有用n×2的长度 而是只用到n的长度让结果 “自然溢出”了 （不是标准说法）
而且貌似n位之后的部分溢出到前面几位了
如果有神犇知道是为什么的话希望能回复下
#include <bits/stdc++.h>
using namespace std;
const int N=<<,M=<<;
const double pi=acos(-);
typedef complex <double> E;
int mp[M][M];
int n,m,L;
int R[N];
long long cnt[N];
E a[N],b[N];
int sum;
double ans;
void fft(E *A,int f)
{
for(int i=;i<n;++i)
if(i<R[i])
swap(A[i],A[R[i]]);
for(int i=;i<n;i<<=)
{
E wn(cos(pi/i),sin(pi/i)*f);
for(int p=i<<,j=;j<n;j+=p)
{
E w(,);
for(int k=;k<i;++k,w*=wn)
{
E x=A[j+k],y=w*A[j+k+i];
A[j+k]=x+y;
A[j+k+i]=x-y;
}
}
}
}
int getdist(int x,int y)
{
return x*x+y*y;
}
int main()
{
scanf("%d",&m);
for(n=;n<m*;n<<=)
++L;
for(int i=;i<m;++i)
for(int j=;j<m;++j)
{
scanf("%d",&mp[i][j]);
sum+=mp[i][j];
cnt[]-=mp[i][j];
a[i*n+j]=mp[i][j];
b[(n-i-)*n+(n-j-)]=mp[i][j];
}
n=n*n*;
L=L*+;
for(int i=;i<n;++i)
R[i]=(R[i>>]>>)|((i&)<<(L-));
fft(a,);
fft(b,);
for(int i=;i<n;++i)
a[i]*=b[i];
fft(a,-);
int tn=sqrt(n/)/,x,y;
x=tn;
y=tn;
for(int i=*tn*tn-*tn-;i;--i)
{
int tmp1;
if(y>=tn)
tmp1=getdist(x-(tn*-),y-(tn*-));
else
tmp1=getdist(x--(tn*-),tn*--(tn*--y-));
double tmp2=a[x*tn*+y].real()/n;
ans+=sqrt(tmp1)*tmp2;
cnt[tmp1]+=tmp2+0.5;
x+=(y==tn*-);
y=(y==tn*-?:y+);
}
printf("%.10f\n",ans/((long long)sum*(sum-)));
int flag=;
for(int i=;i<m*m*;++i)
if(cnt[i])
{
printf("%d %lld\n",i,cnt[i]/);
if(++flag>=)
break;
}
return ;
}
还是把参考的那个只用了n的长度的多项式的神犇的代码贴一下（如果原作者有意见的话我会删掉的……）
#include<stdio.h>
#include<algorithm>
#include<complex>
#include<vector>
#include<math.h>
#include<map>
using namespace std;
const double PI = 4.0*atan(1.0);
typedef complex<double> Complex;
const Complex I(, );
void fft(int n, double theta, Complex a[]) {
for (int m = n; m >= ; m >>= ) {
int mh = m >> ;
for (int i = ; i < mh; i++) {
Complex w = exp(i*theta*I);
for (int j = i; j < n; j += m) {
int k = j + mh;
Complex x = a[j] - a[k];
a[j] += a[k];
a[k] = w * x;
}
}
theta *= ;
}
int i = ;
for (int j = ; j < n - ; j++) {
for (int k = n >> ; k > (i ^= k); k >>= );
if (j < i) swap(a[i], a[j]);
}
}
int b[][];
Complex p[<<];
Complex q[<<];
int main(){
int a;
scanf("%d",&a);
int cnt=;
for(int i=;i<a;i++){
for(int j=;j<a;j++)scanf("%d",&b[i][j]);
}
int n=;
while(n<a*)n*=;
for(int i=;i<a;i++)for(int j=;j<a;j++){
cnt+=b[i][j];
p[i*n+j]=q[(n--i)*n+(n--j)]=Complex(b[i][j],);
}
fft(n*n,*PI/n/n,p);
fft(n*n,*PI/n/n,q); for(int i=;i<n*n;i++)p[i]=p[i]*q[i];
fft(n*n,-*PI/n/n,p);
for(int i=;i<n*n;i++)p[i]/=n*n;
/*for(int i=0;i<n;i++){
for(int j=0;j<n;j++)printf("%.0f ",p[i*n+j].real());
printf("\n");
}*/
double ret=;
long long sz=;
map<int,long long>m;
for(int i=;i<n*n;i++){
int pi=i/n;
int pj=i%n;
if(pj>=n/)continue;
if(pj==&&pi>=n/)continue;
long long v=(long long)(p[n*n--i].real()+0.5);
if(i==){v-=cnt;v/=;}
if(!v)continue; int I=min(pi,n-pi);
int J=pj;
int t=I*I+J*J;
if(m.count(t)){
m[t]+=v;
}else{
m[t]=v;
}
sz+=v;
ret+=sqrt(t)*v;
}
printf("%.12f\n",ret/sz);
int at=;
for(map<int,long long>::iterator it=m.begin();it!=m.end();it++){
if(at==)break;
at++;
printf("%d %lld\n",(*it).first,(*it).second);
}
}