表示敲完多项式乘法&高精度乘法两道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=1<<23,M=1<<10; const double pi=acos(-1); 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=0;i<n;++i) if(i<R[i]) swap(A[i],A[R[i]]); for(int i=1;i<n;i<<=1) { E wn(cos(pi/i),sin(pi/i)*f); for(int p=i<<1,j=0;j<n;j+=p) { E w(1,0); for(int k=0;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=1;n<m*2;n<<=1) ++L; for(int i=0;i<m;++i) for(int j=0;j<m;++j) { scanf("%d",&mp[i][j]); sum+=mp[i][j]; cnt[0]-=mp[i][j]; a[i*n+j]=mp[i][j]; b[(n-i-1)*n+(n-j-1)]=mp[i][j]; } n=n*n*2; L=L*2+1; for(int i=0;i<n;++i) R[i]=(R[i>>1]>>1)|((i&1)<<(L-1)); fft(a,1); fft(b,1); for(int i=0;i<n;++i) a[i]*=b[i]; fft(a,-1); int tn=sqrt(n/2)/2,x,y; x=tn; y=tn; for(int i=4*tn*tn-2*tn-1;i;--i) { int tmp1; if(y>=tn) tmp1=getdist(x-(tn*2-1),y-(tn*2-1)); else tmp1=getdist(x-1-(tn*2-1),tn*2-1-(tn*2-1-y-1)); double tmp2=a[x*tn*2+y].real()/n; ans+=sqrt(tmp1)*tmp2; cnt[tmp1]+=tmp2+0.5; x+=(y==tn*2-1); y=(y==tn*2-1?0:y+1); } printf("%.10f ",ans/((long long)sum*(sum-1))); int flag=0; for(int i=0;i<m*m*2;++i) if(cnt[i]) { printf("%d %lld ",i,cnt[i]/2); if(++flag>=10000) break; } return 0; }
还是把参考的那个只用了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(0, 1); void fft(int n, double theta, Complex a[]) { for (int m = n; m >= 2; m >>= 1) { int mh = m >> 1; for (int i = 0; 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 *= 2; } int i = 0; for (int j = 1; j < n - 1; j++) { for (int k = n >> 1; k > (i ^= k); k >>= 1); if (j < i) swap(a[i], a[j]); } } int b[1025][1025]; Complex p[1<<22]; Complex q[1<<22]; int main(){ int a; scanf("%d",&a); int cnt=0; for(int i=0;i<a;i++){ for(int j=0;j<a;j++)scanf("%d",&b[i][j]); } int n=1; while(n<a*2)n*=2; for(int i=0;i<a;i++)for(int j=0;j<a;j++){ cnt+=b[i][j]; p[i*n+j]=q[(n-1-i)*n+(n-1-j)]=Complex(b[i][j],0); } fft(n*n,2*PI/n/n,p); fft(n*n,2*PI/n/n,q); for(int i=0;i<n*n;i++)p[i]=p[i]*q[i]; fft(n*n,-2*PI/n/n,p); for(int i=0;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(" "); }*/ double ret=0; long long sz=0; map<int,long long>m; for(int i=0;i<n*n;i++){ int pi=i/n; int pj=i%n; if(pj>=n/2)continue; if(pj==0&&pi>=n/2)continue; long long v=(long long)(p[n*n-1-i].real()+0.5); if(i==0){v-=cnt;v/=2;} 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 ",ret/sz); int at=0; for(map<int,long long>::iterator it=m.begin();it!=m.end();it++){ if(at==10000)break; at++; printf("%d %lld ",(*it).first,(*it).second); } }