src/2geom/poly.h

   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  */
  34
  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>
  43
  44 namespace Geom {
  45
  46 class Poly : public std::vector<double>{
  47 public:
  48     // coeff; // sum x^i*coeff[i]
  49
  50     //unsigned size() const { return coeff.size();}
  51     unsigned degree() const { return size()-1;}
  52
  53     //double operator[](const int i) const { return (*this)[i];}
  54     //double& operator[](const int i) { return (*this)[i];}
  55
  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);
  61
  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);
  77
  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);
  92
  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);
 104
 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());
 114
 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);
 124
 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);
 136
 137         result.resize(terms, 0.0);
 138         result.insert(result.end(), this->begin(), this->end());
 139
 140         assert(result.size() == out_size);
 141         return result;
 142     }
 143     Poly operator*(const Poly& p) const;
 144
 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     }
 153
 154     template <typename T>
 155     T operator()(T t) const { return (T)eval(t);}
 156
 157     void normalize();
 158
 159     void monicify();
 160     Poly() {}
 161     Poly(const Poly& p) : std::vector<double>(p) {}
 162     Poly(const double a) {push_back(a);}
 163
 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     }
 184
 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 };
 192
 193 inline Poly operator*(double a, Poly const & b) { return b * a;}
 194
 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);
 201
 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);
 207
 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);
 214
 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];
 228
 229         }
 230     }
 231     return out_file;
 232 }
 233
 234 } // namespace Geom
 235
 236 #endif //LIB2GEOM_SEEN_POLY_H
 237
 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 :