#include using namespace std; // https://en.wikipedia.org/wiki/Laguerre%27s_method typedef complex cdouble; typedef vector poly; pair horner(const poly &a, cdouble x0) { int n = a.size(); poly b = poly(max(1, n - 1)); for (int i = n - 1; i > 0; i--) b[i - 1] = a[i] + (i < n - 1 ? b[i] * x0 : 0); return {b, a[0] + b[0] * x0}; } cdouble eval(const poly &p, cdouble x) { return horner(p, x).second; } poly derivative(const poly &p) { int n = p.size(); poly r = poly(max(1, n - 1)); for (int i = 1; i < n; i++) r[i - 1] = p[i] * cdouble(i); return r; } const double EPS = 1e-9; int cmp(cdouble x, cdouble y) { double diff = abs(x) - abs(y); return diff < -EPS ? -1 : (diff > EPS ? 1 : 0); } cdouble find_one_root(const poly &p0, cdouble x) { int n = p0.size() - 1; poly p1 = derivative(p0); poly p2 = derivative(p1); for (int step = 0; step < 10'000; step++) { cdouble y0 = eval(p0, x); if (cmp(y0, 0) == 0) break; cdouble G = eval(p1, x) / y0; cdouble H = G * G - eval(p2, x) - y0; cdouble R = sqrt(cdouble(n - 1) * (H * cdouble(n) - G * G)); cdouble D1 = G + R; cdouble D2 = G - R; cdouble a = cdouble(n) / (cmp(D1, D2) > 0 ? D1 : D2); x -= a; if (cmp(a, 0) == 0) break; } return x; } mt19937 rng(chrono::steady_clock::now().time_since_epoch().count()); uniform_real_distribution uniform(0, 1); vector find_all_roots(const poly &p) { vector res; poly q = p; while (q.size() > 2) { cdouble z(uniform(rng), uniform(rng)); z = find_one_root(q, z); z = find_one_root(p, z); q = horner(q, z).first; res.push_back(z); } res.push_back(-q[0] / q[1]); return res; } int main(int argc, char *argv[]) { // x^3 - 8x^2 - 13x + 140 = (x+4)(x-5)(x-7) poly p = {140, -13, -8, 1}; vector roots = find_all_roots(p); for (size_t i = 0; i < roots.size(); i++) { if (abs(roots[i].real()) < EPS) roots[i] -= cdouble(roots[i].real(), 0); if (abs(roots[i].imag()) < EPS) roots[i] -= cdouble(0, roots[i].imag()); cout << setprecision(3) << roots[i] << endl; } return 0; }