src/2geom/poly.h

   1 #ifndef LIB2GEOM_SEEN_POLY_H
   2 #define LIB2GEOM_SEEN_POLY_H
   3 #include <assert.h>
   4 #include <vector>
   5 #include <iostream>
   6 #include <algorithm>
   7 #include <complex>
   8 #include <2geom/utils.h>
   9
  10 namespace Geom {
  11
  12 class Poly : public std::vector<double>{
  13 public:
  14     // coeff; // sum x^i*coeff[i]
  15
  16     //unsigned size() const { return coeff.size();}
  17     unsigned degree() const { return size()-1;}
  18
  19     //double operator[](const int i) const { return (*this)[i];}
  20     //double& operator[](const int i) { return (*this)[i];}
  21
  22     Poly operator+(const Poly& p) const {
  23         Poly result;
  24         const unsigned out_size = std::max(size(), p.size());
  25         const unsigned min_size = std::min(size(), p.size());
  26         //result.reserve(out_size);
  27
  28         for(unsigned i = 0; i < min_size; i++) {
  29             result.push_back((*this)[i] + p[i]);
  30         }
  31         for(unsigned i = min_size; i < size(); i++)
  32             result.push_back((*this)[i]);
  33         for(unsigned i = min_size; i < p.size(); i++)
  34             result.push_back(p[i]);
  35         assert(result.size() == out_size);
  36         return result;
  37     }
  38     Poly operator-(const Poly& p) const {
  39         Poly result;
  40         const unsigned out_size = std::max(size(), p.size());
  41         const unsigned min_size = std::min(size(), p.size());
  42         result.reserve(out_size);
  43
  44         for(unsigned i = 0; i < min_size; i++) {
  45             result.push_back((*this)[i] - p[i]);
  46         }
  47         for(unsigned i = min_size; i < size(); i++)
  48             result.push_back((*this)[i]);
  49         for(unsigned i = min_size; i < p.size(); i++)
  50             result.push_back(-p[i]);
  51         assert(result.size() == out_size);
  52         return result;
  53     }
  54     Poly operator-=(const Poly& p) {
  55         const unsigned out_size = std::max(size(), p.size());
  56         const unsigned min_size = std::min(size(), p.size());
  57         resize(out_size);
  58
  59         for(unsigned i = 0; i < min_size; i++) {
  60             (*this)[i] -= p[i];
  61         }
  62         for(unsigned i = min_size; i < out_size; i++)
  63             (*this)[i] = -p[i];
  64         return *this;
  65     }
  66     Poly operator-(const double k) const {
  67         Poly result;
  68         const unsigned out_size = size();
  69         result.reserve(out_size);
  70
  71         for(unsigned i = 0; i < out_size; i++) {
  72             result.push_back((*this)[i]);
  73         }
  74         result[0] -= k;
  75         return result;
  76     }
  77     Poly operator-() const {
  78         Poly result;
  79         result.resize(size());
  80
  81         for(unsigned i = 0; i < size(); i++) {
  82             result[i] = -(*this)[i];
  83         }
  84         return result;
  85     }
  86     Poly operator*(const double p) const {
  87         Poly result;
  88         const unsigned out_size = size();
  89         result.reserve(out_size);
  90
  91         for(unsigned i = 0; i < out_size; i++) {
  92             result.push_back((*this)[i]*p);
  93         }
  94         assert(result.size() == out_size);
  95         return result;
  96     }
  97     // equivalent to multiply by x^terms, negative terms are disallowed
  98     Poly shifted(unsigned const terms) const {
  99         Poly result;
 100         size_type const out_size = size() + terms;
 101         result.reserve(out_size);
 102
 103         result.resize(terms, 0.0);
 104         result.insert(result.end(), this->begin(), this->end());
 105
 106         assert(result.size() == out_size);
 107         return result;
 108     }
 109     Poly operator*(const Poly& p) const;
 110
 111     template <typename T>
 112     T eval(T x) const {
 113         T r = 0;
 114         for(int k = size()-1; k >= 0; k--) {
 115             r = r*x + T((*this)[k]);
 116         }
 117         return r;
 118     }
 119
 120     template <typename T>
 121     T operator()(T t) const { return (T)eval(t);}
 122
 123     void normalize();
 124
 125     void monicify();
 126     Poly() {}
 127     Poly(const Poly& p) : std::vector<double>(p) {}
 128     Poly(const double a) {push_back(a);}
 129
 130 public:
 131     template <class T, class U>
 132     void val_and_deriv(T x, U &pd) const {
 133         pd[0] = back();
 134         int nc = size() - 1;
 135         int nd = pd.size() - 1;
 136         for(unsigned j = 1; j < pd.size(); j++)
 137             pd[j] = 0.0;
 138         for(int i = nc -1; i >= 0; i--) {
 139             int nnd = std::min(nd, nc-i);
 140             for(int j = nnd; j >= 1; j--)
 141                 pd[j] = pd[j]*x + operator[](i);
 142             pd[0] = pd[0]*x + operator[](i);
 143         }
 144         double cnst = 1;
 145         for(int i = 2; i <= nd; i++) {
 146             cnst *= i;
 147             pd[i] *= cnst;
 148         }
 149     }
 150
 151     static Poly linear(double ax, double b) {
 152         Poly p;
 153         p.push_back(b);
 154         p.push_back(ax);
 155         return p;
 156     }
 157 };
 158
 159 inline Poly operator*(double a, Poly const & b) { return b * a;}
 160
 161 Poly integral(Poly const & p);
 162 Poly derivative(Poly const & p);
 163 Poly divide_out_root(Poly const & p, double x);
 164 Poly compose(Poly const & a, Poly const & b);
 165 Poly divide(Poly const &a, Poly const &b, Poly &r);
 166 Poly gcd(Poly const &a, Poly const &b, const double tol=1e-10);
 167
 168 /*** solve(Poly p)
 169  * find all p.degree() roots of p.
 170  * 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?
 171  */
 172 std::vector<std::complex<double> > solve(const Poly & p);
 173
 174 /*** solve_reals(Poly p)
 175  * find all real solutions to Poly p.
 176  * currently we just use solve and pick out the suitably real looking values, there may be a better algorithm.
 177  */
 178 std::vector<double> solve_reals(const Poly & p);
 179 double polish_root(Poly const & p, double guess, double tol);
 180
 181 inline std::ostream &operator<< (std::ostream &out_file, const Poly &in_poly) {
 182     if(in_poly.size() == 0)
 183         out_file << "0";
 184     else {
 185         for(int i = (int)in_poly.size()-1; i >= 0; --i) {
 186             if(i == 1) {
 187                 out_file << "" << in_poly[i] << "*x";
 188                 out_file << " + ";
 189             } else if(i) {
 190                 out_file << "" << in_poly[i] << "*x^" << i;
 191                 out_file << " + ";
 192             } else
 193                 out_file << in_poly[i];
 194
 195         }
 196     }
 197     return out_file;
 198 }
 199
 200 } // namespace Geom
 201
 202 #endif //LIB2GEOM_SEEN_POLY_H
 203
 204 /*
 205   Local Variables:
 206   mode:c++
 207   c-file-style:"stroustrup"
 208   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 209   indent-tabs-mode:nil
 210   fill-column:99
 211   End:
 212 */
 213 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :