From Wikipedia, the free encyclopedia
#include <iostream>
#include <vector>
#include <math.h>
#include <complex>
#include <algorithm>
#include <iterator>
using namespace std;
typedef complex<double> C;
C fast_exp(C x, int y)
{
if (y == 0) return 1.0;
if (y == 1) return x;
int h = y / 2;
return fast_exp(x, h) * fast_exp(x, y-h);
}
typedef vector<double> Poly;
C eval_poly(const Poly& p, C x)
{
C y = 0.0;
for (int i = 0; i < p.size(); ++i)
y += p[i] * fast_exp(x, i);
return y;
}
int main()
{
int n;
cin >> n;
Poly p(n, 0.0);
for (int i = 0; i < n; ++i)
cin >> p[i];
vector<C> q(n-1, C(0.4, 0.9));
for (int i = 0; i < n-1; ++i)
q[i] = fast_exp(C(0.4, 0.9), i);
for (int k = 0; k < 10; ++k)
{
for (int i = 0; i < n-1; ++i)
{
C denom(1.0, 0.0);
for (int j = 0; j < n-1; ++j)
if (i != j)
denom *= q[i] - q[j];
q[i] -= eval_poly(p, q[i]) / denom;
}
}
copy(q.begin(), q.end(), ostream_iterator<C>(cout));
cout << endl;
}