Код:
#include <iostream>
#include <cstdio>
#include <vector>
#include <stack>
#include <algorithm>
#include <cstdlib>
#include <string>
#pragma warning(disable:4996)
using namespace std;
#define pb push_back
const double eps = 1e-9;
const double inf = 1e15;
#define eq(x, y) (abs((x) - (y)) < eps)
#define lt(x, y) ((x) < (y) - eps)
#define le(x, y) ((x) < (y) + eps)
#define ge(x, y) ((x) > (y) - eps)
#define gt(x, y) ((x) > (y) + eps)
//equal, less than, less or equal, greater or equal, greater than
const double FREE = 0;
class poly
{
private:
vector<double> val;
public:
poly () {};
poly (double x);
poly (double x1, double x2);
int size ();
void resize (int n);
void norm ();
poly operator + (poly &P);
poly operator * (poly &P);
poly operator * (double k);
void operator += (poly &P);
void operator *= (poly &P);
void operator *= (double k);
double operator () (double x);
double& operator [] (int i);
};
poly :: poly (double x0)
{
//new polynomial looks like (x - x0) = 0, so x0 is it's root
resize(2);
val[0] = -x0;
val[1] = 1;
};
int poly :: size ()
{
return val.size();
};
void poly :: resize (int n)
{
val.resize(n);
};
double& poly :: operator [] (int i)
{
return (i < size() ? val[i] : FREE);
};
double poly :: operator () (double x)
{
double sum = 0;
for (int i = size() - 1; i >= 0; i--)
{
sum *= x;
sum += (*this)[i];
}
return sum;
};
void poly :: norm ()
{
while (!val.empty() && eq(val.back(), 0))
val.pop_back();
};
poly poly :: operator + (poly &P)
{
poly New;
int lim = max(size(), P.size());
New.resize(lim);
for (int i = 0; i < lim; i++)
New[i] = (*this)[i] + P[i];
New.norm();
return New;
};
poly poly :: operator * (poly &P)
{
poly New;
New.resize(size() + P.size());
for (int i = 0; i < P.size(); i++)
for (int j = 0; j < size(); j++)
New[i + j] += (*this)[j] * P[i];
New.norm();
return New;
};
poly poly :: operator * (double k)
{
poly New;
New.resize(size());
for (int i = 0; i < size(); i++)
New[i] = (*this)[i] * k;
New.norm();
return New;
};
void poly :: operator += (poly &P)
{
//Very dumb implementation of += and *= - to much unnecessary memory using
*this = *this + P;
};
void poly :: operator *= (poly &P)
{
*this = *this * P;
};
void poly :: operator *= (double k)
{
*this = *this * k;
};
ostream& operator << (ostream &out, poly &P)
{
if (P.size() == 0)
{
printf("%.3lf", 0);
return out;
}
int cur = P.size() - 1;
printf("%.3lf", P[cur]);
cur--;
for (; cur != -1; cur--)
if (lt(P[cur], 0))
printf(" - %.3lf", abs(P[cur]));
else
printf(" + %.3lf", abs(P[cur]));
return out;
};
void lagrange (vector<double> &X, vector<double> &Y, poly &P)
{
//minimum amount of pairs {x, y} == 2
int n = X.size();
poly tmp;
P.resize(0);
for (int i = 0; i < n; i++)
{
int start = 1;
if (i == 0)
{
tmp = poly(X[1]);
start = 2;
}
else
tmp = poly(X[0]);
for (int j = start; j < n; j++)
{
if (j == i) continue;
else
tmp *= poly(X[j]);
}
tmp *= (Y[i] / tmp(X[i]));
P += tmp;
}
};
int main()
{
int n;
vector<double> x, y;
poly ans;
cin >> n;
x.resize(n);
y.resize(n);
for (int i = 0; i < n; i++)
cin >> x[i] >> y[i];
lagrange(x, y, ans);
cout << ans << endl;
return 0;
};
________
Код нужно оформлять по правилам:
тегом [CODE]..[/СODE]
(это кнопочка на панели форматирования с решёточкой #)
Не забывайте об этом!
Модератор.