【BZOJ3527】力(ZJOI2014)-FFT

测试地址:
做法:本题需要用到FFT。
把题目所给式子中的 qj q j 除掉,我们发现题目要求:
Ej=i<jqi(ij)2i>jqi(ij)2 E j = ∑ i < j q i ( i − j ) 2 − ∑ i > j q i ( i − j ) 2
然而本人自己做到这里就并不会推了(即使在知道要使用FFT的情况下)。在本人看了题解之后,发现要根据一种新的思路来思考——向量卷积。
根据向量卷积的定义,两个向量 A,B A , B 的卷积 C C 为:
Ci=j+k=iAjbk
为了方便讨论,下述的下标可能是负数。
我们发现若令 A,B A , B 为:
Ai=qi,Bi=sign(i)×1i2 A i = q i , B i = s i g n ( i ) × 1 i 2
其中 sign(i) s i g n ( i ) i i 的符号,即当i<0时, sign(i)=1 s i g n ( i ) = − 1 ,否则 sign(i)=1 s i g n ( i ) = 1
那么不难看出 Ei=Ci E i = C i
为了方便计算向量卷积,我们自然要把 B B 的负数下标转化成非负下标,将下标均加上n1即可。于是用FFT计算向量卷积,然后取 Cn C n ~ C2n1 C 2 n − 1 这些项即可,这些就是我们要求的 E E <script type="math/tex" id="MathJax-Element-36">E</script>。
以下为本人代码(代码中A的下标从0开始,所以最后的向量中从第n-1项开始取):

#include <bits/stdc++.h>
using namespace std;
int n,r[600010];
struct Complex
{
    double x,y;
}a[600010],b[600010];
const double pi=acos(-1.0);
Complex operator + (Complex a,Complex b) {Complex s={a.x+b.x,a.y+b.y};return s;}
Complex operator - (Complex a,Complex b) {Complex s={a.x-b.x,a.y-b.y};return s;}
Complex operator * (Complex a,Complex b) {Complex s={a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};return s;}

void FFT(Complex *a,int type)
{
    for(int i=0;i<n;i++)
        if (i<r[i]) swap(a[i],a[r[i]]);
    for(int mid=1;mid<n;mid<<=1)
    {
        Complex W={cos(pi/mid),type*sin(pi/mid)};
        for(int l=0,r=mid<<1;l<n;l+=r)
        {
            Complex w={1.0,0.0};
            for(int k=0;k<mid;k++,w=w*W)
            {
                Complex x=a[l+k],y=w*a[l+mid+k];
                a[l+k]=x+y;
                a[l+mid+k]=x-y;
            }
        }
    }
    if (type==-1)
    {
        for(int i=0;i<n;i++)
            a[i].x/=n;
    }
}

int main()
{
    scanf("%d",&n);
    memset(a,0,sizeof(a));
    memset(b,0,sizeof(b));
    for(int i=0;i<n;i++) scanf("%lf",&a[i].x);
    for(int i=0;i<(n<<1)-1;i++)
    {
        if (i<n-1) b[i].x=-1.0/((double)(i-n+1))/((double)(i-n+1));
        else if (i>n-1) b[i].x=1.0/((double)(i-n+1))/((double)(i-n+1));
             else b[i].x=0.0;
    }

    int save,bit=0,x=1;
    while(x<=(n<<1)) x<<=1,bit++;
    r[0]=0;
    for(int i=1;i<x;i++)
        r[i]=(r[i>>1]>>1)|((i&1)<<(bit-1));
    save=n;n=x;

    FFT(a,1),FFT(b,1);
    for(int i=0;i<n;i++) a[i]=a[i]*b[i];
    FFT(a,-1);
    for(int i=save-1;i<(save<<1)-1;i++)
        printf("%lf\n",a[i].x);

    return 0;
}
  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值