• bzoj 3527


    一道FFT的典型题

    首先给出结论:

    我们令所有下标从0开始,那么:

    构造两个多项式 A_1(x)=0*x^0+q_1*x^1+....+q_{n-1}*x^{n-1}+q_n*x^nA_2(x)=frac{-1}{n^2}x^0+frac{-1}{(n-1)^2}x^1+...+0x^n+frac{1}{1^2}x^{n+1}+...+frac{1}{(n-1)^2}x^{2n-1}+frac{1}{n^2}x^{2n}

    那么将这两个多项式卷积,然后输出(n,2*n]的系数即可

    为什么?

    我们从已知开始:已知E_i=frac{F_i}{q_i},所以E_i=sum_{j<i}frac{q_j}{(i-j)^2}-sum_{j>i}frac{q_j}{(j-i)^2}

    那么我们对Ei的表达式除掉一个qi,就变成了一个与qi无关的表达式!

    所以E_i=sum_{j<i}frac{q_j}{(i-j)^2}-sum_{j>i}frac{q_j}{(j-i)^2}

    那么,如果我们设f(i)=q_i,g(i)=frac{1}{i^2},那么上式可以重写为:E_j=sum_{i=1}^{j-1}f(i)g(j-i)-sum_{i=1}^{n-j}f(j+i)g(i)

    不难发现,前面那个东西是个卷积的形式,但是后面那个东西好像不大对劲啊...

    长得不太像卷积...

    所以这里就有一个方法:如果我们把g数组反向设为g',那么原式即可改写成:

    E_j=sum_{i=1}^{j-1}f(i)g(j-i)-sum_{i=1}^{n-j}f(j+i)g'(n-i)

    那右边不也是个卷积了?

    这样的话,左右都可以直接上FFT了

    但是为了简化,我们只想用一次FFT

    所以我们再改进一下:我们把g与g'直接相连,将符号变为值的正负,这样就可以直接卷积了

    而这就是上面给出的两个多项式了

    A_1(x)=0*x^0+q_1*x^1+....+q_{n-1}*x^{n-1}+q_n*x^n

    A_2(x)=frac{-1}{n^2}x^0+frac{-1}{(n-1)^2}x^1+...+0x^n+frac{1}{1^2}x^{n+1}+...+frac{1}{(n-1)^2}x^{2n-1}+frac{1}{n^2}x^{2n}

    这样直接做卷积即可

    但是答案是哪部分呢?

    留作思考

    我们发现,卷积后有价值的项x的指数一定大于n(这是很显然的,由上述的表达式即可决定)

    所以就是输出(n,2n]的系数即可

    然后就没了

    #include <cstdio>
    #include <cmath>
    #include <cstring>
    #include <cstdlib>
    #include <iostream>
    #include <algorithm>
    #include <queue>
    #include <stack>
    using namespace std;
    const double pi=acos(-1.0);
    struct cp
    {
        double x,y;
    };
    int n,lim=1,l;
    int to[(1<<20)+5];
    cp operator + (cp a,cp b)
    {
        return (cp){a.x+b.x,a.y+b.y};
    }
    cp operator - (cp a,cp b)
    {
        return (cp){a.x-b.x,a.y-b.y};
    }
    cp operator * (cp a,cp b)
    {
        return (cp){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};
    }
    void FFT(cp *a,int len,int k)
    {
        for(int i=0;i<len;i++)if(i<to[i])swap(a[i],a[to[i]]);
        for(int i=1;i<len;i<<=1)
        {
            cp w0=(cp){cos(pi/i),k*sin(pi/i)};
            for(int j=0;j<len;j+=(i<<1))
            {
                cp w=(cp){1,0};
                for(int o=0;o<i;o++,w=w*w0)
                {
                    cp w1=a[j+o],w2=a[j+i+o]*w;
                    a[j+o]=w1+w2,a[j+i+o]=w1-w2;
                }
            }    
        }
    }
    cp a[(1<<20)+5],b[(1<<20)+5],c[(1<<20)+5];
    int main()
    {
        scanf("%d",&n);
        for(int i=1;i<=n;i++)scanf("%lf",&a[i].x);
        for(int i=0;i<n;i++)b[i].x=-1.0/((double)(n-i)*(double)(n-i));
        for(int i=n+1;i<=2*n;i++)b[i].x=1.0/((double)(i-n)*(double)(i-n));
        int m=(n<<1);
        while(lim<=2*m)lim<<=1,l++;
        for(int i=1;i<lim;i++)to[i]=((to[i>>1]>>1)|((i&1)<<(l-1)));
        FFT(a,lim,1),FFT(b,lim,1);
        for(int i=0;i<lim;i++)c[i]=a[i]*b[i];
        FFT(c,lim,-1);
        for(int i=n+1;i<=n*2;i++)printf("%lf
    ",c[i].x/(double)lim);
        return 0;
    }
  • 相关阅读:
    进程、线程、协程
    python垃圾回收机制
    python变量存储和深浅拷贝
    Linux常用命令
    二叉树四种遍历,节点个数,深度
    装饰器
    ArrayList、Vector
    集合、Collection、迭代器、List
    卖票
    关于Thread和Runnable
  • 原文地址:https://www.cnblogs.com/zhangleo/p/10803239.html
Copyright © 2020-2023  润新知