src/2geom/poly.h

	poly.h revision 8001ba81cb851b38d86650a2fef5817facffb763
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen#ifndef SEEN_POLY_H
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen#define SEEN_POLY_H
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen#include <assert.h>
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen#include <vector>
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen#include <iostream>
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen#include <algorithm>
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen#include <complex>
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen#include <2geom/utils.h>
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenclass Poly : public std::vector<double>{
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenpublic:
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    // coeff; // sum x^i*coeff[i]
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    //unsigned size() const { return coeff.size();}
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    unsigned degree() const { return size()-1;}
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    //double operator[](const int i) const { return (*this)[i];}
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    //double& operator[](const int i) { return (*this)[i];}
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    Poly operator+(const Poly& p) const {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        Poly result;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        const unsigned out_size = std::max(size(), p.size());
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        const unsigned min_size = std::min(size(), p.size());
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        //result.reserve(out_size);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        for(unsigned i = 0; i < min_size; i++) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen            result.push_back((*this)[i] + p[i]);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        }
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        for(unsigned i = min_size; i < size(); i++)
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen            result.push_back((*this)[i]);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        for(unsigned i = min_size; i < p.size(); i++)
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen            result.push_back(p[i]);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        assert(result.size() == out_size);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        return result;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    }
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    Poly operator-(const Poly& p) const {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        Poly result;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        const unsigned out_size = std::max(size(), p.size());
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        const unsigned min_size = std::min(size(), p.size());
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        result.reserve(out_size);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        for(unsigned i = 0; i < min_size; i++) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen            result.push_back((*this)[i] - p[i]);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        }
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        for(unsigned i = min_size; i < size(); i++)
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen            result.push_back((*this)[i]);
8001ba81cb851b38d86650a2fef5817facffb763johanengelen        for(unsigned i = min_size; i < p.size(); i++)
8001ba81cb851b38d86650a2fef5817facffb763johanengelen            result.push_back(-p[i]);
8001ba81cb851b38d86650a2fef5817facffb763johanengelen        assert(result.size() == out_size);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        return result;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    }
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    Poly operator-=(const Poly& p) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        const unsigned out_size = std::max(size(), p.size());
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        const unsigned min_size = std::min(size(), p.size());
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        resize(out_size);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        for(unsigned i = 0; i < min_size; i++) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen            (*this)[i] -= p[i];
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        }
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        for(unsigned i = min_size; i < out_size; i++)
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen            (*this)[i] = -p[i];
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        return *this;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    }
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    Poly operator-(const double k) const {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        Poly result;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        const unsigned out_size = size();
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        result.reserve(out_size);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        for(unsigned i = 0; i < out_size; i++) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen            result.push_back((*this)[i]);
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        }
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        result[0] -= k;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        return result;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    }
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    Poly operator-() const {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        Poly result;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        result.resize(size());
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        for(unsigned i = 0; i < size(); i++) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen            result[i] = -(*this)[i];
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        }
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        return result;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    }
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    Poly operator*(const double p) const {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        Poly result;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        const unsigned out_size = size();
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        result.reserve(out_size);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        for(unsigned i = 0; i < out_size; i++) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen            result.push_back((*this)[i]*p);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        }
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        assert(result.size() == out_size);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        return result;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    }
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen// equivalent to multiply by x^terms, discard negative terms
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    Poly shifted(unsigned terms) const {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        Poly result;
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen        // This was a no-op and breaks the build on x86_64, as it's trying
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        // to take maximum of 32-bit and 64-bit integers
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        //const unsigned out_size = std::max(unsigned(0), size()+terms);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        const size_type out_size = size() + terms;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        result.reserve(out_size);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        if(terms < 0) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen            for(unsigned i = 0; i < out_size; i++) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen                result.push_back((*this)[i-terms]);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen            }
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        } else {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen            for(unsigned i = 0; i < terms; i++) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen                result.push_back(0.0);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen            }
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen            for(unsigned i = 0; i < size(); i++) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen                result.push_back((*this)[i]);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen            }
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        }
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        assert(result.size() == out_size);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        return result;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    }
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    Poly operator*(const Poly& p) const;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    template <typename T>
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    T eval(T x) const {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        T r = 0;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        for(int k = size()-1; k >= 0; k--) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen            r = r*x + T((*this)[k]);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        }
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        return r;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    }
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    template <typename T>
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    T operator()(T t) const { return (T)eval(t);}
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    void normalize();
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    void monicify();
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    Poly() {}
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    Poly(const Poly& p) : std::vector<double>(p) {}
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    Poly(const double a) {push_back(a);}
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenpublic:
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    template <class T, class U>
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    void val_and_deriv(T x, U &pd) const {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        pd[0] = back();
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        int nc = size() - 1;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        int nd = pd.size() - 1;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        for(unsigned j = 1; j < pd.size(); j++)
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen            pd[j] = 0.0;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        for(int i = nc -1; i >= 0; i--) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen            int nnd = std::min(nd, nc-i);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen            for(int j = nnd; j >= 1; j--)
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen                pd[j] = pd[j]*x + operator[](i);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen            pd[0] = pd[0]*x + operator[](i);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        }
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        double cnst = 1;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        for(int i = 2; i <= nd; i++) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen            cnst *= i;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen            pd[i] *= cnst;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        }
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    }
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    static Poly linear(double ax, double b) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        Poly p;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        p.push_back(b);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        p.push_back(ax);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        return p;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    }
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen};
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen
981b809bc6ed10a21e89444d9447e5475801874fjohanengeleninline Poly operator*(double a, Poly const & b) { return b * a;}
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenPoly integral(Poly const & p);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenPoly derivative(Poly const & p);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenPoly divide_out_root(Poly const & p, double x);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenPoly compose(Poly const & a, Poly const & b);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenPoly divide(Poly const &a, Poly const &b, Poly &r);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenPoly gcd(Poly const &a, Poly const &b, const double tol=1e-10);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen/*** solve(Poly p)
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * find all p.degree() roots of p.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * 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?
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenstd::vector<std::complex<double> > solve(const Poly & p);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen/*** solve_reals(Poly p)
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * find all real solutions to Poly p.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * currently we just use solve and pick out the suitably real looking values, there may be a better algorithm.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenstd::vector<double> solve_reals(const Poly & p);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelendouble polish_root(Poly const & p, double guess, double tol);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen
981b809bc6ed10a21e89444d9447e5475801874fjohanengeleninline std::ostream &operator<< (std::ostream &out_file, const Poly &in_poly) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    if(in_poly.size() == 0)
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        out_file << "0";
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    else {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        for(int i = (int)in_poly.size()-1; i >= 0; --i) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen            if(i == 1) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen                out_file << "" << in_poly[i] << "*x";
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen                out_file << " + ";
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen            } else if(i) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen                out_file << "" << in_poly[i] << "*x^" << i;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen                out_file << " + ";
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen            } else
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen                out_file << in_poly[i];
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen        }
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    }
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen    return out_file;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen}
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen/*
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen  Local Variables:
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen  mode:c++
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen  c-file-style:"stroustrup"
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen  indent-tabs-mode:nil
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen  fill-column:99
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen  End:
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen*/
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen#endif
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen