首页 > 其他 > 详细

[ZJOI2014]力

时间:2018-12-08 22:11:18      阅读:181      评论:0      收藏:0      [点我收藏+]

嘟嘟嘟


这应该是我的第一篇fft题解吧。


对于这道比较裸的多项式乘法题,主要思路就是推出卷积的式子,然后用fft加速。
卷积我在网上找了半天,只看到一篇博客上写了这个东西。再加上自己的脑补,觉得卷积好像就是这个式子:
\[c(i) = \sum_{j = 0} ^ {i} a_j *b_{i - j}\]
其中\(a, b\)是俩多项式,只要符合这个形式,\(a, b\)就可以用fft加速。(自己意会的,如果不对路过的大佬一定要帮我指出来!)


那么对于这道题,目标就是把题中的式子变成可以卷积的式子。
题意:
求序列\(E\),每一项\(E_i\)满足
\[E_i = \sum_{j < i} {\frac{q_j}{(i - j) ^ 2}} - \sum_{j >i} {\frac{q_j}{(j - i) ^ 2}}\]
然后就开始推导了……
\[E_i = \sum_{j = 0} ^ {i - 1} \frac{q_j}{(i - j) ^ 2} - \sum_{j = i + 1} ^ {n - 1} \frac{q_j}{(j - i) ^ 2}\]
为了和卷积的下标保持统一,这里的下标也从\(0\)开始。
\(A(i)\)表示第一个多项式,\(B(i)\)表示第二个多项式,且\(f(i) = q_i\)\(g(i) = \frac{1}{i ^ 2}\),规定\(g(0) = 0\),则
\[E_i = A(i) - B(i)\]
\[A(i) = \sum_{j = 0} ^ {i} f(j) *g(i - j)\]

\(B(i) = \sum _{j = i} ^ {n - 1} f(j) *g(j - i)\)
????\(= \sum _ {j = 0} ^ {n - i - 1} f(i +j) * g(j)\)

如果令\(f‘(i)\)表示\(f(i)\)翻转后的序列,则
\[B(i) = \sum _{j = 0} ^{n - i - 1} f(n - i - j - 1) * g(i)\]
这样\(A(i), B(i)\)都符合了卷积的形式,于是fft就可做了。

#include<cstdio>
#include<iostream>
#include<cmath>
#include<algorithm>
#include<cstring>
#include<cstdlib>
#include<cctype>
#include<vector>
#include<stack>
#include<queue>
#include<complex>
using namespace std;
#define enter puts("") 
#define space putchar(‘ ‘)
#define Mem(a, x) memset(a, x, sizeof(a))
#define rg register
typedef long long ll;
typedef double db;
typedef complex<db> cp;
const int INF = 0x3f3f3f3f;
const db eps = 1e-8;
const db PI = acos(-1);
const int maxn = 3e5 + 5;
inline ll read()
{
    ll ans = 0;
    char ch = getchar(), last = ‘ ‘;
    while(!isdigit(ch)) {last = ch; ch = getchar();}
    while(isdigit(ch)) {ans = (ans << 1) + (ans << 3) + ch - ‘0‘; ch = getchar();}
    if(last == ‘-‘) ans = -ans;
    return ans;
}
inline void write(ll x)
{
    if(x < 0) x = -x, putchar(‘-‘);
    if(x >= 10) write(x / 10);
    putchar(x % 10 + ‘0‘);
}

int n, len = 1;
cp f1[maxn], f2[maxn], g[maxn], omg[maxn], inv[maxn];

void init()
{
    for(int i = 0; i < len; ++i)
    {
        omg[i] = cp(cos(2 * PI * i / len), sin(2 * PI * i / len));
        inv[i] = conj(omg[i]);
    }
}
void fft(cp* a, cp* omg)
{
    int lim = 0;
    while((1 << lim) < len) lim++;
    for(int i = 0; i < len; ++i)
    {
        int t = 0;
        for(int j = 0; j < lim; ++j) if((i >> j) & 1) t |= (1 << (lim - j - 1));
        if(i < t) swap(a[i], a[t]);
    }
    for(int l = 2; l <= len; l <<= 1)
    {
        int q = l >> 1;
        for(cp* p = a; p != a + len; p += l)    
            for(int i = 0; i < q; ++i)
            {
                cp t = omg[len / l * i] * p[i + q];
                p[i + q] = p[i] - t, p[i] += t;
            }
    }
}

int main()
{
    n = read();
    for(int i = 0; i < n; ++i)
    {
        db q; scanf("%lf", &q);
        f1[i].real(q); f2[n - i - 1].real(q);
        if(i) g[i].real(1.0 / (db)i / (db)i);
    }
    while(len < (n << 1)) len <<= 1;
    init();
    fft(f1, omg); fft(f2, omg); fft(g, omg);
    for(int i = 0; i < len; ++i) f1[i] *= g[i], f2[i] *= g[i];
    fft(f1, inv); fft(f2, inv); //别忘了再变回系数表示法 
    for(int i = 0; i < n; ++i) 
        printf("%.8lf\n", f1[i].real() / len  - f2[n - i - 1].real() / len);
    return 0;
}

[ZJOI2014]力

原文:https://www.cnblogs.com/mrclr/p/10088928.html

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