Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #include <iostream>
- #include <vector>
- #include <cmath>
- using namespace std;
- const double PI = acos(-1.0);
- class FFT {
- private:
- struct Complex {
- double re, im;
- Complex(double r = 0.0, double i = 0.0) : re(r), im(i) {}
- Complex operator+(const Complex& other) const {
- return Complex(re + other.re, im + other.im);
- }
- Complex operator-(const Complex& other) const {
- return Complex(re - other.re, im - other.im);
- }
- Complex operator*(const Complex& other) const {
- return Complex(re * other.re - im * other.im,
- re * other.im + im * other.re);
- }
- Complex conj() const { // Сопряжение
- return Complex(re, -im);
- }
- };
- // Инверсия индексов
- void bit_reverse(vector<Complex>& a) {
- int n = a.size(), logn = 0;
- while ((1 << logn) < n) ++logn;
- for (int i = 0; i < n; ++i) {
- int j = 0;
- for (int k = 0; k < logn; ++k)
- if (i & (1 << k))
- j |= (1 << (logn - 1 - k));
- if (i < j)
- swap(a[i], a[j]);
- }
- }
- // Быстрое преобразование Фурье
- void fft(vector<Complex>& a, bool invert) {
- int n = a.size();
- bit_reverse(a);
- 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.re /= n;
- x.im /= n;
- }
- }
- }
- public:
- // Умножение двух полиномов через БПФ
- vector<double> multiply(const vector<double>& a, const vector<double>& b) {
- vector<Complex> fa(a.begin(), a.end()), fb(b.begin(), b.end());
- int n = 1;
- while (n < a.size() + b.size()) n <<= 1;
- 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> result(n);
- for (int i = 0; i < n; ++i)
- result[i] = round(fa[i].re); // Округляем
- // Убираем хвостовые нули
- while (!result.empty() && fabs(result.back()) < 1e-9)
- result.pop_back();
- return result;
- }
- };
- int main() {
- // Полиномы: (1 + 2x + 3x^2) * (4 + 5x)
- vector<double> A = {1, 2, 3}; // 1 + 2x + 3x²
- vector<double> B = {4, 5}; // 4 + 5x
- FFT fft;
- vector<double> result = fft.multiply(A, B);
- cout << "Result:\n";
- for (int i = 0; i < result.size(); ++i) {
- cout << result[i] << (i ? "x^" + to_string(i) : "") << " ";
- if (i + 1 < result.size()) cout << "+ ";
- }
- cout << '\n';
- return 0;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
