[ZJOI2014]力 解题报告

[ZJOI2014]力

题意

给出 (n) 个数 (q_i),
定义

[F_j = sum_{i=1}^{j-1} frac{q_i imes q_j}{(i-j)^2} - sum_{i=j+1}^{n} frac{q_i imes q_j}{(i-j)^2} ]

[E_i = frac{F_i}{q_i} ]

求 (E_i) ((i in [1,n])).

思路

首先, 提公因式, 化简可以得到,

[ E_j = sum_{i=1}^{j-1} frac{q_i}{(j-i)^2} - sum_{i=j+1}^{n} frac{q_i}{(j-i)^2} $$ (把 $q_j$ 提出来后约掉) 注意到式子可以分为两个部分, 这两个部分中的因子可以分成 $i$ 和 $j-i$ 两个部分, 即这两个部分的编号和为 $j$, 那么我们就找到加法卷积了. 设 $F[i] = frac{1}{i^2}$, $G[i] = q_i$. 则 $$ E_j = sum_{i=1}^{j-1} F[j-i] imes G[i] - sum_{i=j+1}^{n} F[j-i] imes G[i] ]

这样, 我们就可以分两步愉快地 FFT 了.

需要注意的是, 后半部分的 (j-i) 为负数, 所以需要将整个数组往前移动 (n) 位, 最后统计答案的时候也要往前移动 (n) 位 (没听懂的话可以看一下代码).

代码

#include<bits/stdc++.h>
#define _USE_MATH_DEFINES
#define ll long long
#define db double
using namespace std;
const int _=1e5+7;
const int __=3e5+7;
const db Pi=M_PI;
struct cn{
  db a,b;
  cn operator + (const cn &x) const {
    return (cn){a+x.a,b+x.b};
  }
  cn operator - (const cn &x) const {
    return (cn){a-x.a,b-x.b};
  }
  cn operator * (const cn &x) const {
    return (cn){a*x.a-b*x.b,a*x.b+b*x.a};
  }
};
int N,n,num[__];
db q[_],E[_];
cn f[__],g[__];
void FFT(cn *f,int id){
  for(int i=0;i<n;i++)
    if(i<num[i]) swap(f[i],f[num[i]]);
  for(int len=2;len<=n;len<<=1){
    int gap=len>>1;
    cn w1=(cn){cos(2*Pi/len),sin(2*Pi/len)*id};
    for(int k=0;k<n;k+=len){
      cn w=(cn){1,0};
      for(int i=k;i<k+gap;i++,w=w*w1){
	cn tmp=w*f[i+gap];
	f[i+gap]=f[i]-tmp;
	f[i]=f[i]+tmp;
      }
    }
  }
}
int main(){
  //freopen("force.in","r",stdin);
  cin>>N;
  for(int i=1;i<=N;i++) scanf("%lf",&q[i]);
  n=1; while(n<=2*N) n<<=1;
  for(ll i=1;i<=N;i++){    // 这里需要开 long long
    f[i].a=(db)1ll/(i*i);
    g[i].a=q[i];
  }
  for(int i=0;i<n;i++)
    num[i]=(num[i>>1]>>1)|((i&1) ?n>>1 :0);
  FFT(f,1);
  FFT(g,1);
  for(int i=0;i<n;i++) f[i]=f[i]*g[i];
  FFT(f,-1);
  for(int i=1;i<=N;i++) E[i]=f[i].a/n;
  memset(f,0,sizeof(f));
  memset(g,0,sizeof(g));
  for(ll i=0;i<=N;i++){
    f[i].a=(db)1ll/((i-N)*(i-N)); // 往前移动 N 位
    g[i].a=q[i];
  }
  f[N].a=0;
  FFT(f,1);
  FFT(g,1);
  for(int i=0;i<n;i++) f[i]=f[i]*g[i];
  FFT(f,-1);
  for(int i=1;i<=N;i++) printf("%.3lf
",E[i]-f[i+N].a/n);  // 统计答案时也往前移动 N 位
  return 0;
}