src/2geom/poly.h

   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"
   9
  10 class Poly : public std::vector<double>{
  11 public:
  12     // coeff; // sum x^i*coeff[i]
  13
  14     //unsigned size() const { return coeff.size();}
  15     unsigned degree() const { return size()-1;}
  16
  17     //double operator[](const int i) const { return (*this)[i];}
  18     //double& operator[](const int i) { return (*this)[i];}
  19
  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);
  25
  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);
  41
  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);
  56
  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);
  68
  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());
  78
  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);
  88
  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);
 103
 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         }
 116
 117         assert(result.size() == out_size);
 118         return result;
 119     }
 120     Poly operator*(const Poly& p) const;
 121
 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     }
 130
 131     template <typename T>
 132     T operator()(T t) const { return (T)eval(t);}
 133
 134     void normalize();
 135
 136     void monicify();
 137     Poly() {}
 138     Poly(const Poly& p) : std::vector<double>(p) {}
 139     Poly(const double a) {push_back(a);}
 140
 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     }
 161
 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 };
 169
 170 inline Poly operator*(double a, Poly const & b) { return b * a;}
 171
 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);
 178
 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);
 184
 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);
 191
 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];
 205
 206         }
 207     }
 208     return out_file;
 209 }
 210
 211
 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