Сумма k-ых степеней

#include <bits/stdc++.h>
#define all(x) begin(x), end(x)
using namespace std;
typedef long long ll;
typedef long double ld;
#define all(x) begin(x), end(x)

struct poly {
    vector<ld> a;
    poly() {}
    poly(vector<ld> b) : a(b) {}
    int size() {
        return (int)a.size();
    }
    ld &operator[](int i) {
        assert(0 <= i && i < size());
        return a[i];
    }
    void pop_zero() {
        while (a.back() == 0) {
            a.pop_back();
        }
    }
    ld operator()(ld x) {
        ld pw = 1, ans = 0;
        for (int i = 0; i < size(); ++i) {
            ans += pw * a[i];
            pw *= x;
        }
        return ans;
    }
};

poly operator*(poly a, poly b) {
    poly c;
    c.a.assign(a.size() + b.size(), 0);
    for (int i = 0; i < a.size(); ++i) {
        for (int j = 0; j < b.size(); ++j) {
            c[i + j] += a[i] * b[j];
        }
    }
    c.pop_zero();
    return c;
}

poly operator+(poly a, poly b) {
    poly c;
    c.a.assign(max(a.size(), b.size()) + 1, 0);
    for (int i = 0; i < a.size(); ++i) {
        c[i] += a[i];
    }
    for (int j = 0; j < b.size(); ++j) {
        c[j] += b[j];
    }
    c.pop_zero();
    return c;
}

signed main() {
#ifdef LC
    assert(freopen("input.txt", "r", stdin));
#endif
    ios::sync_with_stdio(0);
    cin.tie(0);

    int n = 6;
    vector<ld> x(n), y(n);
    for (int i = 0; i < n; ++i) {
        x[i] = i;
        y[i] = powl(i, n - 2);
        if (i) {
            y[i] += y[i - 1];
        }
    }

    poly p;

    for (int i = 0; i < n; ++i) {
        poly cur({1});
        ld down = 1;
        for (int j = 0; j < n; ++j) {
            if (i == j) {
                continue;
            }
            cur = cur * poly({-x[j], 1});
            down *= x[i] - x[j];
        }
        for (ld &e : cur.a) {
            e *= y[i] / down;
        }
        p = p + cur;
    }

    for (int i = 0; i < p.size(); ++i) {
        cout << fixed << setprecision(6) << p[i] << " x^" << i << endl;
    }

    const ld EPS = 1e-9;
    mt19937 rnd;
    vector<ld> roots;

    while (clock() < CLOCKS_PER_SEC) {
        ld x = rnd() / ld(rnd() + 1);
        ld f = abs(p(x));
        for (ld st = 1; st > 1e-12; st *= 0.9) {
            while (1) {
                ld c = abs(p(x - st));
                if (f > c) {
                    f = c;
                    x -= st;
                } else {
                    break;
                }
            }
            while (1) {
                ld c = abs(p(x + st));
                if (f > c) {
                    f = c;
                    x += st;
                } else {
                    break;
                }
            }
        }
        assert(abs(p(x)) < EPS);
        bool already = false;
        for (ld e : roots) {
            if (abs(e - x) < EPS) {
                already = true;
            }
        }
        if (!already) {
            roots.push_back(x);
        }
    }
    sort(all(roots));
    for (ld e : roots) {
        cout << e << endl;
    }
    return 0;
}