10
class Poly : public std::vector<double>{
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 {
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]);
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);
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);
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());
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);
95
// equivalent to multiply by x^terms, discard negative terms
96
Poly shifted(unsigned terms) const {
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);
105
for(unsigned i = 0; i < out_size; i++) {
106
result.push_back((*this)[i-terms]);
109
for(unsigned i = 0; i < terms; i++) {
110
result.push_back(0.0);
112
for(unsigned i = 0; i < size(); i++) {
113
result.push_back((*this)[i]);
117
assert(result.size() == out_size);
120
Poly operator*(const Poly& p) const;
122
template <typename T>
125
for(int k = size()-1; k >= 0; k--) {
126
r = r*x + T((*this)[k]);
131
template <typename T>
132
T operator()(T t) const { return (T)eval(t);}
138
Poly(const Poly& p) : std::vector<double>(p) {}
139
Poly(const double a) {push_back(a);}
142
template <class T, class U>
143
void val_and_deriv(T x, U &pd) const {
146
int nd = pd.size() - 1;
147
for(unsigned j = 1; j < pd.size(); j++)
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);
156
for(int i = 2; i <= nd; i++) {
162
static Poly linear(double ax, double b) {
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);
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?
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.
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)
196
for(int i = (int)in_poly.size()-1; i >= 0; --i) {
198
out_file << "" << in_poly[i] << "*x";
201
out_file << "" << in_poly[i] << "*x^" << i;
204
out_file << in_poly[i];
215
c-file-style:"stroustrup"
216
c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
221
// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
added
removed
src/2geom/poly.h
