AIZU 2560 [想法题]

表示敲完多项式乘法&高精度乘法两道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\n",ans/((long long)sum*(sum-1)));
    int flag=0;
    for(int i=0;i<m*m*2;++i)
        if(cnt[i])
        {
            printf("%d %lld\n",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("\n");
    }*/
    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\n",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\n",(*it).first,(*it).second);
    }
}