10 class Poly :
public std::vector<double>{
15 unsigned degree()
const {
return size()-1;}
20 Poly operator+(
const Poly& p)
const {
23 const unsigned min_size = std::min(size(), p.size());
26 for(
unsigned i = 0; i < min_size; i++) {
27 result.push_back((*
this)[i] + p[i]);
29 for(
unsigned i = min_size; i < size(); i++)
30 result.push_back((*
this)[i]);
31 for(
unsigned i = min_size; i < p.size(); i++)
32 result.push_back(p[i]);
36 Poly operator-(
const Poly& p)
const {
38 const unsigned out_size = std::max(size(), p.size());
39 const unsigned min_size = std::min(size(), p.size());
40 result.reserve(out_size);
42 for(
unsigned i = 0; i < min_size; i++) {
43 result.push_back((*
this)[i] - p[i]);
45 for(
unsigned i = min_size; i < size(); i++)
46 result.push_back((*
this)[i]);
47 for(
unsigned i = min_size; i < p.size(); i++)
48 result.push_back(-p[i]);
49 assert(result.size() == out_size);
53 const unsigned out_size = std::max(size(), p.size());
54 const unsigned min_size = std::min(size(), p.size());
57 for(
unsigned i = 0; i < min_size; i++) {
60 for(
unsigned i = min_size; i < out_size; i++)
64 Poly operator-(
const double k)
const {
66 const unsigned out_size = size();
67 result.reserve(out_size);
69 for(
unsigned i = 0; i < out_size; i++) {
70 result.push_back((*
this)[i]);
75 Poly operator-()
const {
77 result.resize(size());
79 for(
unsigned i = 0; i < size(); i++) {
80 result[i] = -(*this)[i];
84 Poly operator*(
const double p)
const {
86 const unsigned out_size = size();
87 result.reserve(out_size);
89 for(
unsigned i = 0; i < out_size; i++) {
90 result.push_back((*
this)[i]*p);
92 assert(result.size() == out_size);
96 Poly shifted(
unsigned terms)
const {
101 const size_type out_size = size() + terms;
102 result.reserve(out_size);
110 for(
unsigned i = 0; i < terms; i++) {
111 result.push_back(0.0);
113 for(
unsigned i = 0; i < size(); i++) {
114 result.push_back((*
this)[i]);
118 assert(result.size() == out_size);
121 Poly operator*(
const Poly& p)
const;
123 template <
typename T>
126 for(
int k = size()-1; k >= 0; k--) {
127 r = r*x + T((*
this)[k]);
132 template <
typename T>
133 T operator()(T t)
const {
return (T)eval(t);}
139 Poly(
const Poly& p) : std::vector<double>(p) {}
140 Poly(
const double a) {push_back(a);}
143 template <
class T,
class U>
144 void val_and_deriv(T x, U &pd)
const {
147 int nd = pd.size() - 1;
148 for(
unsigned j = 1; j < pd.size(); j++)
150 for(
int i = nc -1; i >= 0; i--) {
151 int nnd = std::min(nd, nc-i);
152 for(
int j = nnd; j >= 1; j--)
153 pd[j] = pd[j]*x +
operator[](i);
154 pd[0] = pd[0]*x + operator[](i);
157 for(
int i = 2; i <= nd; i++) {
163 static Poly linear(
double ax,
double b) {
171 inline Poly operator*(
double a,
Poly const & b) {
return b * a;}
175 Poly divide_out_root(
Poly const & p,
double x);
178 Poly gcd(
Poly const &a,
Poly const &b,
const double tol=1e-10);
184 std::vector<std::complex<double> > solve(
const Poly & p);
190 std::vector<double> solve_reals(
const Poly & p);
191 double polish_root(
Poly const & p,
double guess,
double tol);
193 inline std::ostream &
operator<< (std::ostream &out_file,
const Poly &in_poly) {
194 if(in_poly.size() == 0)
197 for(
int i = (
int)in_poly.size()-1; i >= 0; --i) {
199 out_file <<
"" << in_poly[i] <<
"*x";
202 out_file <<
"" << in_poly[i] <<
"*x^" << i;
205 out_file << in_poly[i];
std::ostream & operator<<(std::ostream &out_file, const Geom::Matrix &m)
Definition: matrix.h:109