【ZJOI2014】【BZOJ3527】力 (FFT)

最新推荐文章于 2019-08-31 09:29:00 发布

原创最新推荐文章于 2019-08-31 09:29:00 发布 · 354 阅读

1 ·

CC 4.0 BY-SA版权

BZOJ 同时被 3 个专栏收录

215 篇文章

订阅专栏

省选

90 篇文章

订阅专栏

FFT

19 篇文章

订阅专栏

本文介绍了一种利用快速傅立叶变换（FFT）解决特定序列求和问题的方法。通过将原问题转换为两个序列的卷积问题，利用FFT高效地求解。代码示例展示了如何实现这一算法。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

Description

给出序列 $\{q_i\}$ ，求

E i = \sum j < i q j ( i - j ) 2 - \sum j > i q j ( i - j ) 2

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

Solution

考虑：
令 $p_{n - i} = q_i$

E i = \sum j = 1 i - 1 q j g i - j - \sum j = i + 1 n q j g j - i = \sum j = 1 i - 1 q j g i - j - \sum j = 0 n - i - 1 p j g n - j - i

$E_i=\sum_{j=1}^{i-1}q_jg_{i-j}-\sum_{j=i+1}^{n}q_jg_{j-i}=\sum_{j=1}^{i-1}q_jg_{i-j}-\sum_{j=0}^{n - i - 1}p_jg_{n-j-i}$
然后FFT求卷积即可。

Code

#include<bits/stdc++.h>

using namespace std;

typedef long long ll;
typedef complex<double> cpx;
#define rep(i, j) for (register int i = 0, i##_end_ = (j); i < i##_end_; ++ i)
#define For(i , j , k) for (register int i = (j) , i##_end_ = (k) ; i <= i##_end_ ; ++ i)
#define Fordown(i , j , k) for (register int i = (j) , i##_end_ = (k) ; i >= i##_end_ ; -- i)
#define Set(a , b) memset(a , b , sizeof(a))
#define pb(a) push_back(a)
#define mp(a, b) make_pair(a, b)
#define ALL(a) (a).begin(), (a).end()
#define SZ(a) ((int)(a).size())
#define fir first
#define sec second
#define INF (0x3f3f3f3f)
#define INF1 (2139062143)
#define Mod (1000000007)
#ifdef hany01
#define debug(...) fprintf(stderr , __VA_ARGS__)
#else
#define debug(...)
#endif

template <typename T> inline bool chkmax(T &a , T b) { return a < b ? (a = b , 1) : 0; }
template <typename T> inline bool chkmin(T &a , T b) { return b < a ? (a = b , 1) : 0; }

inline void File()
{
#ifdef hany01
    freopen("3527.in" , "r" , stdin);
    freopen("3527.out" , "w" , stdout);
#endif
}

const double pi = acos(-1.0);
const int maxn = 1 << 21;

int n, m, cnt, rev[maxn];
cpx q[maxn], p[maxn], g[maxn];

inline void Init()
{
    scanf("%d", &n);
    For(i, 1, n) scanf("%lf", &q[i].real()), p[n - i].real() = q[i].real(), g[i].real() = 1.0 / (i * 1ll * i);
    for (m = n << 1, n = 1; n < m; n <<= 1) ++ cnt;
    For(i, 0, n) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (cnt - 1));
}

inline void FFT(cpx *a, int type)
{
    For(i, 0, n) if (rev[i] > i) swap(a[i], a[rev[i]]);
    for (register int i = 2; i <= n; i <<= 1)
    {
        register cpx wn(cos(2 * pi / i), sin(2 * pi / i) * type);
        for (register int j = 0; j < n; j += i)
        {
            register cpx w(1, 0);
            rep(k, i >> 1)
            {
                register cpx x = a[j + k], y = a[j + k + (i >> 1)] * w;
                a[j + k] = x + y; a[j + k + (i >> 1)] = x - y;
                w = w * wn;
            }
        }
    }
}

inline void Solve()
{
    FFT(q, 1); FFT(p, 1); FFT(g, 1);
    For(i, 0, n) p[i] *= g[i], q[i] *= g[i];
    FFT(q, -1); FFT(p, -1);
}

inline void Print()
{
    For(i, 1, m >> 1) printf("%.6lf\n", (q[i].real() - p[(m >> 1) - i].real()) / n);
}

int main()
{
    File();
    Init();
    Solve();
    Print();
    return 0;
}
//风吹柳花满店香，吴姬压酒唤客尝。
//金陵子弟来相送，欲行不行各尽觞。
//请君试问东流水，别意与之谁短长。
//--李白《金陵酒肆留别》