【ZJOI2014】【BZOJ3527】 力 (FFT)

Description

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

$$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_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;
}
//风吹柳花满店香，吴姬压酒唤客尝。
//金陵子弟来相送，欲行不行各尽觞。
//请君试问东流水，别意与之谁短长。
//--李白《金陵酒肆留别》