1. 程式人生 > >BZOJ3527 [Zjoi2014]力 【fft】

BZOJ3527 [Zjoi2014]力 【fft】

alt 固定 sum ++ n) 觀察 gpo tdi long

題目

給出n個數qi,給出Fj的定義如下:
技術分享圖片
令Ei=Fi/qi,求Ei.

輸入格式

第一行一個整數n。

接下來n行每行輸入一個數,第i行表示qi。

輸出格式

n行,第i行輸出Ei。與標準答案誤差不超過1e-2即可。

輸入樣例

5

4006373.885184

15375036.435759

1717456.469144

8514941.004912

1410681.345880

輸出樣例

-16838672.693

3439.793

7509018.566

4595686.886

10903040.872

題解

卷積什麽的感覺好優美~~

卷積

先普及一下離散卷積的定義【瞎編的】:

對於兩個序列\(x(n)\)\(y(n)\)


其卷積\((x*y)(n) = \sum_{-\infty}^{\infty}x(k)y(n - k)\)

即當一個序列所有i位置上的值c(i)等於所有位置之和為i的x(k)*y(i - k)乘積的和時,可以看做c()為x()和y()的卷積

就好比多項式a(n) b(n)相乘,對於次數i的系數\(c(i)=\sum a(k)*b(i - k)\)

而求離散卷積可以使用離散快速傅裏葉\(O(nlogn)\)高效求出

本題

觀察式子
\(Ei = \sum_{j<i}\frac{qj}{(i-j)^2} - \sum_{j>i}\frac{qj}{(i-j)^2}\)
我們將兩個求和分開來求

我們令\(b(i) = \frac{1}{i^2}\),特別的,\(b(0) = 0\)
我們令\(a(i) = qi\)
我們會發現左邊【即為\(L(i)\)\(L(i) = \sum a(j)*b(i - j)\),剛好就是卷積的形式
可以用fft求出

同樣的,對於右邊
\(R(i) = \sum a(j)*b(j - i)\)
誒?不對啊,\(j + j - i\)不是定值啊。
但是ta們的位置關系還是很固定,考慮變形

我們將\(a(i)\)翻轉,即令\(c(n-i)=a(i)\)
奇跡發生了:
\(R(i) = \sum c(n - j)*b(j - i)\)
這樣我們算出的卷積,\(R(i)\)

就與\(E_{n-i}\)對應

最後將算出的兩個結果相減
呼啦啦,搞完啦~~

#include<iostream>
#include<cmath>
#include<cstdio>
#include<cstring>
#include<complex>
#include<algorithm>
#define LL long long int
#define REP(i,n) for (int i = 1; i <= (n); i++)
#define Redge(u) for (int k = h[u],to; k; k = ed[k].nxt)
#define BUG(s,n) for (int i = 1; i <= (n); i++) cout<<s[i]<<‘ ‘; puts("");
using namespace std;
const int maxn = 400005,maxm = 100005,INF = 1000000000;
const double pi = acos(-1);
typedef complex<double> E;
E a[maxn],b[maxn],aa[maxn];
int n,m,L,R[maxn];
void fft(E* a,int f){
    for (int i = 0; i < n; i++) if (i < R[i]) swap(a[i],a[R[i]]);
    for (int i = 1; i < n; i <<= 1){
        E wn(cos(pi / i),f * sin(pi / i));
        for (int j = 0; j < n; j += (i << 1)){
            E w(1,0);
            for (int k = 0; k < i; k++,w *= wn){
                E x = a[j + k],y = w * a[j + k + i];
                a[j + k] = x + y; a[j + k + i] = x - y;
            }
        }
    }
    if (f == -1) for (int i = 0; i < n; i++) a[i] /= n;
}
int main(){
    scanf("%d",&n); --n; double q;
    for (int i = 0; i <= n; i++){
        scanf("%lf",&q);
        a[i] = q; aa[n - i] = q;
    }
    for (int i = 1; i <= n; i++) b[i] = 1.0 / i / i;
    m = n << 1; for (n = 1; n <= m; n <<= 1) L++;
    for (int i = 0; i < n; i++) R[i] = (R[i >> 1] >> 1) | ((i & 1) << (L - 1));
    fft(a,1); fft(aa,1); fft(b,1);
    for (int i = 0; i < n; i++) a[i] *= b[i];
    for (int i = 0; i < n; i++) aa[i] *= b[i];
    fft(a,-1); fft(aa,-1);
    for (int i = 0; i <= (m >> 1); i++) printf("%.6lf\n",a[i].real() - aa[(m >> 1) - i].real());
    return 0;
}

BZOJ3527 [Zjoi2014]力 【fft】