首页 > 其他 > 详细

【BZOJ3527】[ZJOI2014] 力(FFT)

时间:2018-12-15 18:21:56      阅读:108      评论:0      收藏:0      [点我收藏+]

题目:

BZOJ3527

分析:

FFT应用第一题……

首先很明显能把\(F_j\)约掉,变成:

\[E_j=\sum _{i<j} \frac{q_i}{(i-j)^2}-\sum_{i>j}\frac{q_i}{(i-j)^2}\]

然后去膜拜题解,我们知道两个多项式相乘的方式如下:

\[C_j=\sum_{i=0}^j A_iB_{j-i}\]

那么,如果把\(E_j\)的表达式化成上面那个形式,就可以用FFT计算了。(不会FFT?戳我:【知识总结】快速傅里叶变换(FFT)

先看减号前面的部分。显然可以变成(为了叙述方便,读入的\(q\)的下标为\([0,n)\)):

\[C_j=\sum_{i=0}^j F_iG_{j-i}\]

其中\(F_i=q_i\)\(G_i=\frac{1}{i^2}\)。特别地,\(G_0=0\)

减号后面要处理\(j\)位置以后的,怎么办?大力把\(q\)数组翻过来,这样就相当于求\(n-j-1\)以前的了:

\[D_j=\sum_{i=0}^{j} F'_iG_{j-i}\]

其中\(F'_j=q_{n-j-1}\)

那么答案\(E_j=C_j-D_{n-j-1}\)

代码:

注意\(g\)要初始化……

#include <cstdio>
#include <algorithm>
#include <cstring>
#include <cctype>
#include <cmath>
using namespace std;
#define _ 0
 
namespace zyt
{
    template<typename T>
    inline void read(T &x)
    {
        char c;
        bool f = false;
        x = 0;
        do
            c = getchar();
        while (c != '-' && !isdigit(c));
        if (c == '-')
            f = true, c = getchar();
        do
            x = x * 10 + c - '0', c = getchar();
        while (isdigit(c));
        if (f)
            x = -x;
    }
    inline void read(double &x)
    {
        scanf("%lf", &x);
    }
    template<typename T>
    inline void write(T x)
    {
        static char buf[20];
        char *pos = buf;
        if (x < 0)
            putchar('-'), x = -x;
        do
            *pos++ = x % 10 + '0';
        while (x /= 10);
        while (pos > buf)
            putchar(*--pos);
    }
    inline void write(const double a, const int fix = 9)
    {
        printf("%.*f", fix, a);
    }
    const int B = 17, N = 1 << (B + 1) | 10;
    const double PI = 3.141592653589793238462643383279502884197169399375105820974944L;
    namespace FFT
    {
        int rev[N];
        struct cpx
        {
            double x, y;
            cpx(const double _x = 0.0, const double _y = 0.0)
                : x(_x), y(_y) {}
            cpx operator + (const cpx &b) const
            {
                return cpx(x + b.x, y + b.y);
            }
            cpx operator - (const cpx &b) const
            {
                return cpx(x - b.x, y - b.y);
            }
            cpx operator * (const cpx &b) const
            {
                return cpx(x * b.x - y * b.y, x * b.y + y * b.x);
            }
            cpx conj() const
            {
                return cpx(x, -y);
            }
        }omega[N], inv[N];
        void init(const int lg2)
        {
            for (int i = 0; i < (1 << lg2); i++)
            {
                rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (lg2 - 1));
                omega[i] = cpx(cos(2 * PI * i / (1 << lg2)), sin(2 * PI * i / (1 << lg2)));
                inv[i] = omega[i].conj();
            }
        }
        void fft(cpx *a, const int n, const cpx *w)
        {
            for (int i = 0; i < n; i++)
                if (i < rev[i])
                    swap(a[i], a[rev[i]]);
            for (int len = 1; len < n; len <<= 1)
                for (int i = 0; i < n; i += (len << 1))
                    for (int k = 0; k < len; k++)
                    {
                        cpx tmp = a[i + k] - w[k * (n / (len << 1))] * a[i + len + k];
                        a[i + k] = a[i + k] + w[k * (n / (len << 1))] * a[i + len + k];
                        a[i + len + k] = tmp;
                    }
        }
        void solve(cpx *a, cpx *b, const int n)
        {
            fft(a, n, omega), fft(b, n, omega);
            for (int i = 0; i < n; i++)
                a[i] = a[i] * b[i];
            fft(a, n, inv);
            for (int i = 0; i < n; i++)
                a[i].x /= n;
        }
    }
    using namespace FFT;
    int n;
    double q[N];
    cpx f[N], g[N], revf[N];
    int work()
    {
        read(n);
        for (int i = 0; i < n; i++)
        {
            read(q[i]);
            f[i] = revf[n - i - 1] = q[i];
        }
        int m = n << 1, lg2 = 0;
        for (n = 1; n < m; n <<= 1)
            ++lg2;
        init(lg2);
        for (int i = 0; i < (m >> 1); i++)
            g[i] = (i ? 1.0 / ((double)i * i) : 0.0);
        solve(f, g, n);
        for (int i = 0; i < (m >> 1); i++)
            g[i] = (i ? 1.0 / ((double)i * i) : 0.0);
        for (int i = (m >> 1); i < n; i++)
            g[i] = 0;
        solve(revf, g, n);
        for (int i = 0; i < (m >> 1); i++)
            write(f[i].x - revf[(m >> 1) - i - 1].x), putchar('\n');
        return ~~(0^_^0);
    }
}
 
int main()
{
    return zyt::work(); 
}

【BZOJ3527】[ZJOI2014] 力(FFT)

原文:https://www.cnblogs.com/zyt1253679098/p/10124111.html

(0)
(0)
   
举报
评论 一句话评论(0
关于我们 - 联系我们 - 留言反馈 - 联系我们:wmxa8@hotmail.com
© 2014 bubuko.com 版权所有
打开技术之扣,分享程序人生!