1 /**
2 * \file
3 * \brief \todo brief description
4 *
5 * Authors:
6 * ? <?@?.?>
7 *
8 * Copyright ?-? authors
9 *
10 * This library is free software; you can redistribute it and/or
11 * modify it either under the terms of the GNU Lesser General Public
12 * License version 2.1 as published by the Free Software Foundation
13 * (the "LGPL") or, at your option, under the terms of the Mozilla
14 * Public License Version 1.1 (the "MPL"). If you do not alter this
15 * notice, a recipient may use your version of this file under either
16 * the MPL or the LGPL.
17 *
18 * You should have received a copy of the LGPL along with this library
19 * in the file COPYING-LGPL-2.1; if not, write to the Free Software
20 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21 * You should have received a copy of the MPL along with this library
22 * in the file COPYING-MPL-1.1
23 *
24 * The contents of this file are subject to the Mozilla Public License
25 * Version 1.1 (the "License"); you may not use this file except in
26 * compliance with the License. You may obtain a copy of the License at
27 * http://www.mozilla.org/MPL/
28 *
29 * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
30 * OF ANY KIND, either express or implied. See the LGPL or the MPL for
31 * the specific language governing rights and limitations.
32 *
33 */
35 #ifndef LIB2GEOM_SEEN_POLY_H
36 #define LIB2GEOM_SEEN_POLY_H
37 #include <assert.h>
38 #include <vector>
39 #include <iostream>
40 #include <algorithm>
41 #include <complex>
42 #include <2geom/utils.h>
44 namespace Geom {
46 class Poly : public std::vector<double>{
47 public:
48 // coeff; // sum x^i*coeff[i]
50 //unsigned size() const { return coeff.size();}
51 unsigned degree() const { return size()-1;}
53 //double operator[](const int i) const { return (*this)[i];}
54 //double& operator[](const int i) { return (*this)[i];}
56 Poly operator+(const Poly& p) const {
57 Poly result;
58 const unsigned out_size = std::max(size(), p.size());
59 const unsigned min_size = std::min(size(), p.size());
60 //result.reserve(out_size);
62 for(unsigned i = 0; i < min_size; i++) {
63 result.push_back((*this)[i] + p[i]);
64 }
65 for(unsigned i = min_size; i < size(); i++)
66 result.push_back((*this)[i]);
67 for(unsigned i = min_size; i < p.size(); i++)
68 result.push_back(p[i]);
69 assert(result.size() == out_size);
70 return result;
71 }
72 Poly operator-(const Poly& p) const {
73 Poly result;
74 const unsigned out_size = std::max(size(), p.size());
75 const unsigned min_size = std::min(size(), p.size());
76 result.reserve(out_size);
78 for(unsigned i = 0; i < min_size; i++) {
79 result.push_back((*this)[i] - p[i]);
80 }
81 for(unsigned i = min_size; i < size(); i++)
82 result.push_back((*this)[i]);
83 for(unsigned i = min_size; i < p.size(); i++)
84 result.push_back(-p[i]);
85 assert(result.size() == out_size);
86 return result;
87 }
88 Poly operator-=(const Poly& p) {
89 const unsigned out_size = std::max(size(), p.size());
90 const unsigned min_size = std::min(size(), p.size());
91 resize(out_size);
93 for(unsigned i = 0; i < min_size; i++) {
94 (*this)[i] -= p[i];
95 }
96 for(unsigned i = min_size; i < out_size; i++)
97 (*this)[i] = -p[i];
98 return *this;
99 }
100 Poly operator-(const double k) const {
101 Poly result;
102 const unsigned out_size = size();
103 result.reserve(out_size);
105 for(unsigned i = 0; i < out_size; i++) {
106 result.push_back((*this)[i]);
107 }
108 result[0] -= k;
109 return result;
110 }
111 Poly operator-() const {
112 Poly result;
113 result.resize(size());
115 for(unsigned i = 0; i < size(); i++) {
116 result[i] = -(*this)[i];
117 }
118 return result;
119 }
120 Poly operator*(const double p) const {
121 Poly result;
122 const unsigned out_size = size();
123 result.reserve(out_size);
125 for(unsigned i = 0; i < out_size; i++) {
126 result.push_back((*this)[i]*p);
127 }
128 assert(result.size() == out_size);
129 return result;
130 }
131 // equivalent to multiply by x^terms, negative terms are disallowed
132 Poly shifted(unsigned const terms) const {
133 Poly result;
134 size_type const out_size = size() + terms;
135 result.reserve(out_size);
137 result.resize(terms, 0.0);
138 result.insert(result.end(), this->begin(), this->end());
140 assert(result.size() == out_size);
141 return result;
142 }
143 Poly operator*(const Poly& p) const;
145 template <typename T>
146 T eval(T x) const {
147 T r = 0;
148 for(int k = size()-1; k >= 0; k--) {
149 r = r*x + T((*this)[k]);
150 }
151 return r;
152 }
154 template <typename T>
155 T operator()(T t) const { return (T)eval(t);}
157 void normalize();
159 void monicify();
160 Poly() {}
161 Poly(const Poly& p) : std::vector<double>(p) {}
162 Poly(const double a) {push_back(a);}
164 public:
165 template <class T, class U>
166 void val_and_deriv(T x, U &pd) const {
167 pd[0] = back();
168 int nc = size() - 1;
169 int nd = pd.size() - 1;
170 for(unsigned j = 1; j < pd.size(); j++)
171 pd[j] = 0.0;
172 for(int i = nc -1; i >= 0; i--) {
173 int nnd = std::min(nd, nc-i);
174 for(int j = nnd; j >= 1; j--)
175 pd[j] = pd[j]*x + operator[](i);
176 pd[0] = pd[0]*x + operator[](i);
177 }
178 double cnst = 1;
179 for(int i = 2; i <= nd; i++) {
180 cnst *= i;
181 pd[i] *= cnst;
182 }
183 }
185 static Poly linear(double ax, double b) {
186 Poly p;
187 p.push_back(b);
188 p.push_back(ax);
189 return p;
190 }
191 };
193 inline Poly operator*(double a, Poly const & b) { return b * a;}
195 Poly integral(Poly const & p);
196 Poly derivative(Poly const & p);
197 Poly divide_out_root(Poly const & p, double x);
198 Poly compose(Poly const & a, Poly const & b);
199 Poly divide(Poly const &a, Poly const &b, Poly &r);
200 Poly gcd(Poly const &a, Poly const &b, const double tol=1e-10);
202 /*** solve(Poly p)
203 * find all p.degree() roots of p.
204 * 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?
205 */
206 std::vector<std::complex<double> > solve(const Poly & p);
208 /*** solve_reals(Poly p)
209 * find all real solutions to Poly p.
210 * currently we just use solve and pick out the suitably real looking values, there may be a better algorithm.
211 */
212 std::vector<double> solve_reals(const Poly & p);
213 double polish_root(Poly const & p, double guess, double tol);
215 inline std::ostream &operator<< (std::ostream &out_file, const Poly &in_poly) {
216 if(in_poly.size() == 0)
217 out_file << "0";
218 else {
219 for(int i = (int)in_poly.size()-1; i >= 0; --i) {
220 if(i == 1) {
221 out_file << "" << in_poly[i] << "*x";
222 out_file << " + ";
223 } else if(i) {
224 out_file << "" << in_poly[i] << "*x^" << i;
225 out_file << " + ";
226 } else
227 out_file << in_poly[i];
229 }
230 }
231 return out_file;
232 }
234 } // namespace Geom
236 #endif //LIB2GEOM_SEEN_POLY_H
238 /*
239 Local Variables:
240 mode:c++
241 c-file-style:"stroustrup"
242 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
243 indent-tabs-mode:nil
244 fill-column:99
245 End:
246 */
247 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :