zoukankan      html  css  js  c++  java
  • [ZJOI2014]力

    推公式发现(这不是水题吗,这要推吗

    [E_i=Sigma^{i-1}_{j=1} frac{q_j}{(i-j)^2} - Sigma^{n}_{j=i+1} frac{q_j}{(i-j)^2} ]

    [设A[i] = q[i], B[i] = frac{1}{i^2},FFT将A,B相乘可以得到E_i的前半部分 ]

    [后半部分就把数组反过来FFT就可以了 ]

    也可以合起来FFT(懒得想

    # include <bits/stdc++.h>
    # define RG register
    # define IL inline
    # define Fill(a, b) memset(a, b, sizeof(a))
    using namespace std;
    typedef long long ll;
    const int _(2e6 + 10);
    const double Pi(acos(-1));
    
    struct Complex{
    	double real, image;
    	IL Complex(){  real = image = 0;  }
    	IL Complex(RG double a, RG double b){  real = a; image = b;  }
    	IL Complex operator +(RG Complex B){  return Complex(real + B.real, image + B.image);  }
    	IL Complex operator -(RG Complex B){  return Complex(real - B.real, image - B.image);  }
    	IL Complex operator *(RG Complex B){  return Complex(real * B.real - image * B.image, real * B.image + image * B.real);  }
    } A[_], B[_];
    int n, N, M, l, r[_];
    double q[_], E[_];
    
    IL void FFT(RG Complex *P, RG int opt){
    	for(RG int i = 0; i < N; ++i) if(i < r[i]) swap(P[i], P[r[i]]);
    	for(RG int i = 1; i < N; i <<= 1){
    		RG Complex W(cos(Pi / i), opt * sin(Pi / i));
    		for(RG int p = i << 1, j = 0; j < N; j += p){
    			RG Complex w(1, 0);
    			for(RG int k = 0; k < i; ++k, w = w * W){
    				RG Complex X = P[k + j], Y = w * P[k + j + i];
    				P[k + j] = X + Y; P[k + j + i] = X - Y;
    			}
    		}
    	}
    }
    
    int main(RG int argc, RG char *argv[]){
    	scanf("%d", &n);
    	for(RG int i = 1; i <= n; ++i) scanf("%lf", &q[i]);
    	for(RG int i = 1; i <= n; ++i) A[i].real = q[i], B[i].real = 1.0 / (1.0 * i * i);
    	for(N = 1, M = 2 * n; N <= M; N <<= 1) ++l;
    	for(RG int i = 0; i < N; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
    	FFT(A, 1); FFT(B, 1);
    	for(RG int i = 0; i < N; ++i) A[i] = A[i] * B[i];
    	FFT(A, -1);
    	for(RG int i = 1; i <= n; ++i) E[i] = A[i].real;
    	for(RG int i = 0; i < N; ++i) A[i].real = A[i].image = B[i].real = B[i].image = 0;
    	for(RG int i = n; i; --i) A[n - i + 1].real = q[i], B[i].real = 1.0 / (1.0 * i * i);
    	FFT(A, 1); FFT(B, 1);
    	for(RG int i = 0; i < N; ++i) A[i] = A[i] * B[i];
    	FFT(A, -1);
    	for(RG int i = 1; i <= n; ++i) E[i] -= A[n - i + 1].real;
    	for(RG int i = 1; i <= n; ++i) printf("%.3lf
    ", E[i] / N);
    	return 0;
    }
    
    
  • 相关阅读:
    学习&分享
    跳槽
    20121113:延期通知书
    2012.9.9 baocheng博客园正式与大家见面啦!
    数据库
    ASP.Net模板引擎
    javascript图片切换效果
    dockercompose环境下zookeeper单机搭建、集群搭建
    Linux服务器日常巡检脚本
    MMOS FFB伺服直驱方向盘主控板DIY
  • 原文地址:https://www.cnblogs.com/cjoieryl/p/8206699.html
Copyright © 2011-2022 走看看