ZJOI2014 力

传送门

一道FFT的标准练习题。

把(q_i)除进去，就可以得到(E(j) = sum_{i<j}frac{q_j}{(i-j)^2} - sum_{i > j}frac{q_j}{(i-j)^2}).

把这个式子前后两部分分别拆开，设(f(i) = q_i),(g(i) = frac{1}{i^2}),那么前一部分就是(sum_{i=0}^j f(j) * g(i-j)),而后一部分就是(sum_{i = 0}^{n - j - 1}f(i+j) * g(j))

后半部分好像不大好维护……我们把整个序列调过来，用(f1)表示(f)调过来后的序列，那么后半部分就可以写成(sum_{i = 0}^{n - j - 1} f1(n - i - j - 1) * g(j))

那么前后两个式子都符合了卷积的形式，可以使用FFT求解。然后它问的每一个(E_i)相对应地就是每一个系数。

看一下代码。

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<iostream>
#include<cmath>
#include<set>
#include<vector>
#include<map>
#include<queue>
#define rep(i,a,n) for(int i = a;i <= n;i++)
#define per(i,n,a) for(int i = n;i >= a;i--)
#define enter putchar('
')
#define fr friend inline
#define y1 poj
#define mp make_pair
#define pr pair<int,int>
#define fi first
#define sc second
#define pb push_back

using namespace std;
typedef long long ll;
const int M = 200005;
const int INF = 1000000009;
const double eps = 1e-7;
const double pi = acos(-1);
const ll mod = 998244353;

ll read()
{
    ll ans = 0,op = 1;char ch = getchar();
    while(ch < '0' || ch > '9') {if(ch == '-') op = -1;ch = getchar();}
    while(ch >= '0' && ch <= '9') ans = ans * 10 + ch - '0',ch = getchar();
    return ans * op;
}

struct Comp
{
   double x,y;
   Comp(){}
   Comp(double px,double py){x = px,y = py;}
   fr Comp operator + (const Comp &a,const Comp &b) {return (Comp){a.x + b.x,a.y + b.y};}
   fr Comp operator - (const Comp &a,const Comp &b) {return (Comp){a.x - b.x,a.y - b.y};}
   fr Comp operator * (const Comp &a,const Comp &b) {return (Comp){a.x * b.x - a.y * b.y,a.y * b.x + a.x * b.y};}
}f[M<<1],g[M<<1],f1[M<<1],kx,ky;

int n,rev[M<<1],len = 1,L;
double c[M<<1];

void fft(Comp *a,int f)
{
   rep(i,0,len-1) if(i < rev[i]) swap(a[i],a[rev[i]]);
   for(int i = 1;i < len;i <<= 1)
   {
      Comp omg1(cos(pi / i),f * sin(pi / i));
      for(int j = 0;j < len;j += (i << 1))
      {
     Comp omg(1,0);
     rep(k,0,i-1)
     {
        kx = a[j+k],ky = omg * a[j+k+i];
        a[j+k] = kx + ky,a[j+k+i] = kx - ky,omg = omg * omg1;
     }
      }
   }
}

int main()
{
   n = read() - 1;
   rep(i,0,n)
   {
      scanf("%lf",&f[i].x),f1[n-i].x = f[i].x;
      if(i) g[i].x = 1.0 / double(i) / double(i);
   }
   while(len <= n << 1) len <<= 1,L++;
   rep(i,0,len-1) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (L-1));
   fft(f,1),fft(f1,1),fft(g,1);
   rep(i,0,len-1) f[i] = f[i] * g[i],f1[i] = f1[i] * g[i];
   fft(f,-1),fft(f1,-1);
   rep(i,0,n) c[i] = (f[i].x - f1[n - i].x) / len;
   rep(i,0,n) printf("%.3lf
",c[i]);
   return 0;
}