您现在的位置：程式碼之家 >> 編程語言 > >> C語言 >> C++ >> 關於C++

BZOJ 3527 [ZJOI2014] 力：FFT 卷積求場強

編輯：關於C++

題目大意

給定 n 個電荷 q[i] 位於座標 i (1≤i≤n)，每個電荷產生電場 E[j] = Σ_{i≠j} q[i] / (i-j)²。求每個位置 j 的電場強度 E[j]。

n ≤ 250,000

數學建模：卷積形式

定義兩數列：

a[i] = q[i]
b[i] = 1/i² (i≥1), b[0] = 0

左側貢獻：E_left[j] = Σ_{i<j} a[i] × b[j-i] → 標準卷積

右側貢獻：E_right[j] = Σ_{i>j} a[i] × b[i-j]

令 a' 為 a 反序，則 E_right[j] = Σ_{k} a'[k] × b[j-k] → 同樣是卷積

最終 E[j] = E_left[j] - E_right[j]

FFT 實作細節

長度擴展至 2 的冪 ≥ 2n
雙精度 complex 足夠（答案保留 5 位小數）
兩次 FFT：正序卷積 + 反序卷積

完整程式碼 (C++)

#include <bits/stdc++.h>
using namespace std;

const double PI = acos(-1.0);

struct Complex {
    double x, y;
    Complex(double _x=0, double _y=0): x(_x), y(_y) {}
    Complex operator+(const Complex& o) const { return Complex(x+o.x, y+o.y); }
    Complex operator-(const Complex& o) const { return Complex(x-o.x, y-o.y); }
    Complex operator*(const Complex& o) const { return Complex(x*o.x - y*o.y, x*o.y + y*o.x); }
};

void fft(vector<Complex>& a, bool invert) {
    int n = a.size();
    // 位元反序排列
    for (int i=1, j=0; i<n; ++i) {
        int bit = n >> 1;
        for (; j & bit; bit >>= 1) j ^= bit;
        j ^= bit;
        if (i < j) swap(a[i], a[j]);
    }
    // 蝴蝶運算
    for (int len=2; len<=n; len<<=1) {
        double ang = 2*PI/len * (invert ? -1 : 1);
        Complex wlen(cos(ang), sin(ang));
        for (int i=0; i<n; i+=len) {
            Complex w(1, 0);
            for (int j=0; j<len/2; ++j) {
                Complex u = a[i+j];
                Complex v = a[i+j+len/2] * w;
                a[i+j] = u + v;
                a[i+j+len/2] = u - v;
                w = w * wlen;
            }
        }
    }
    if (invert) for (auto& x : a) x.x /= n, x.y /= n;
}

vector<double> convolution(const vector<double>& a, const vector<double>& b) {
    int n = 1;
    while (n < (int)a.size() + (int)b.size()) n <<= 1;
    vector<Complex> fa(a.begin(), a.end()), fb(b.begin(), b.end());
    fa.resize(n); fb.resize(n);
    fft(fa, false); fft(fb, false);
    for (int i=0; i<n; ++i) fa[i] = fa[i] * fb[i];
    fft(fa, true);
    vector<double> res(a.size() + b.size() - 1);
    for (size_t i=0; i<res.size(); ++i) res[i] = fa[i].x;
    return res;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    
    int n; cin >> n;
    vector<double> q(n);
    for (int i=0; i<n; ++i) cin >> q[i];
    
    // 建立 b[i] = 1/i^2
    vector<double> b(n, 0);
    for (int i=1; i<n; ++i) b[i] = 1.0 / (i * i);
    
    // 左側卷積：q * b
    auto left = convolution(q, b);
    
    // 右側卷積：反序 q * b
    vector<double> q_rev = q;
    reverse(q_rev.begin(), q_rev.end());
    auto right = convolution(q_rev, b);
    
    // 輸出結果
    cout << fixed << setprecision(5);
    for (int i=0; i<n; ++i) {
        double ans = left[i] - right[n-1-i];
        cout << ans << '\n';
    }
    return 0;
}