1 #ifndef SEEN_POLY_H
2 #define SEEN_POLY_H
3 #include <assert.h>
4 #include <vector>
5 #include <iostream>
6 #include <algorithm>
7 #include <complex>
8 #include "utils.h"
10 class Poly : public std::vector<double>{
11 public:
12 // coeff; // sum x^i*coeff[i]
14 //unsigned size() const { return coeff.size();}
15 unsigned degree() const { return size()-1;}
17 //double operator[](const int i) const { return (*this)[i];}
18 //double& operator[](const int i) { return (*this)[i];}
20 Poly operator+(const Poly& p) const {
21 Poly result;
22 const unsigned out_size = std::max(size(), p.size());
23 const unsigned min_size = std::min(size(), p.size());
24 //result.reserve(out_size);
26 for(unsigned i = 0; i < min_size; i++) {
27 result.push_back((*this)[i] + p[i]);
28 }
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]);
33 assert(result.size() == out_size);
34 return result;
35 }
36 Poly operator-(const Poly& p) const {
37 Poly result;
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]);
44 }
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);
50 return result;
51 }
52 Poly operator-=(const Poly& p) {
53 const unsigned out_size = std::max(size(), p.size());
54 const unsigned min_size = std::min(size(), p.size());
55 resize(out_size);
57 for(unsigned i = 0; i < min_size; i++) {
58 (*this)[i] -= p[i];
59 }
60 for(unsigned i = min_size; i < out_size; i++)
61 (*this)[i] = -p[i];
62 return *this;
63 }
64 Poly operator-(const double k) const {
65 Poly result;
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]);
71 }
72 result[0] -= k;
73 return result;
74 }
75 Poly operator-() const {
76 Poly result;
77 result.resize(size());
79 for(unsigned i = 0; i < size(); i++) {
80 result[i] = -(*this)[i];
81 }
82 return result;
83 }
84 Poly operator*(const double p) const {
85 Poly result;
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);
91 }
92 assert(result.size() == out_size);
93 return result;
94 }
95 // equivalent to multiply by x^terms, discard negative terms
96 Poly shifted(unsigned terms) const {
97 Poly result;
98 // This was a no-op and breaks the build on x86_64, as it's trying
99 // to take maximum of 32-bit and 64-bit integers
100 //const unsigned out_size = std::max(unsigned(0), size()+terms);
101 const size_type out_size = size() + terms;
102 result.reserve(out_size);
104 if(terms < 0) {
105 for(unsigned i = 0; i < out_size; i++) {
106 result.push_back((*this)[i-terms]);
107 }
108 } else {
109 for(unsigned i = 0; i < terms; i++) {
110 result.push_back(0.0);
111 }
112 for(unsigned i = 0; i < size(); i++) {
113 result.push_back((*this)[i]);
114 }
115 }
117 assert(result.size() == out_size);
118 return result;
119 }
120 Poly operator*(const Poly& p) const;
122 template <typename T>
123 T eval(T x) const {
124 T r = 0;
125 for(int k = size()-1; k >= 0; k--) {
126 r = r*x + T((*this)[k]);
127 }
128 return r;
129 }
131 template <typename T>
132 T operator()(T t) const { return (T)eval(t);}
134 void normalize();
136 void monicify();
137 Poly() {}
138 Poly(const Poly& p) : std::vector<double>(p) {}
139 Poly(const double a) {push_back(a);}
141 public:
142 template <class T, class U>
143 void val_and_deriv(T x, U &pd) const {
144 pd[0] = back();
145 int nc = size() - 1;
146 int nd = pd.size() - 1;
147 for(unsigned j = 1; j < pd.size(); j++)
148 pd[j] = 0.0;
149 for(int i = nc -1; i >= 0; i--) {
150 int nnd = std::min(nd, nc-i);
151 for(int j = nnd; j >= 1; j--)
152 pd[j] = pd[j]*x + operator[](i);
153 pd[0] = pd[0]*x + operator[](i);
154 }
155 double cnst = 1;
156 for(int i = 2; i <= nd; i++) {
157 cnst *= i;
158 pd[i] *= cnst;
159 }
160 }
162 static Poly linear(double ax, double b) {
163 Poly p;
164 p.push_back(b);
165 p.push_back(ax);
166 return p;
167 }
168 };
170 inline Poly operator*(double a, Poly const & b) { return b * a;}
172 Poly integral(Poly const & p);
173 Poly derivative(Poly const & p);
174 Poly divide_out_root(Poly const & p, double x);
175 Poly compose(Poly const & a, Poly const & b);
176 Poly divide(Poly const &a, Poly const &b, Poly &r);
177 Poly gcd(Poly const &a, Poly const &b, const double tol=1e-10);
179 /*** solve(Poly p)
180 * find all p.degree() roots of p.
181 * This function can take a long time with suitably crafted polynomials, but in practice it should be fast. Should we provide special forms for degree() <= 4?
182 */
183 std::vector<std::complex<double> > solve(const Poly & p);
185 /*** solve_reals(Poly p)
186 * find all real solutions to Poly p.
187 * currently we just use solve and pick out the suitably real looking values, there may be a better algorithm.
188 */
189 std::vector<double> solve_reals(const Poly & p);
190 double polish_root(Poly const & p, double guess, double tol);
192 inline std::ostream &operator<< (std::ostream &out_file, const Poly &in_poly) {
193 if(in_poly.size() == 0)
194 out_file << "0";
195 else {
196 for(int i = (int)in_poly.size()-1; i >= 0; --i) {
197 if(i == 1) {
198 out_file << "" << in_poly[i] << "*x";
199 out_file << " + ";
200 } else if(i) {
201 out_file << "" << in_poly[i] << "*x^" << i;
202 out_file << " + ";
203 } else
204 out_file << in_poly[i];
206 }
207 }
208 return out_file;
209 }
212 /*
213 Local Variables:
214 mode:c++
215 c-file-style:"stroustrup"
216 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
217 indent-tabs-mode:nil
218 fill-column:99
219 End:
220 */
221 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
222 #endif