src/2geom/poly.cpp

   1 #include <2geom/poly.h>
   2
   3 #define HAVE_GSL
   4 #ifdef HAVE_GSL
   5 #include <gsl/gsl_poly.h>
   6 #endif
   7
   8 namespace Geom {
   9
  10 Poly Poly::operator*(const Poly& p) const {
  11     Poly result;
  12     result.resize(degree() +  p.degree()+1);
  13
  14     for(unsigned i = 0; i < size(); i++) {
  15         for(unsigned j = 0; j < p.size(); j++) {
  16             result[i+j] += (*this)[i] * p[j];
  17         }
  18     }
  19     return result;
  20 }
  21
  22 /*double Poly::eval(double x) const {
  23     return gsl_poly_eval(&coeff[0], size(), x);
  24     }*/
  25
  26 void Poly::normalize() {
  27     while(back() == 0)
  28         pop_back();
  29 }
  30
  31 void Poly::monicify() {
  32     normalize();
  33
  34     double scale = 1./back(); // unitize
  35
  36     for(unsigned i = 0; i < size(); i++) {
  37         (*this)[i] *= scale;
  38     }
  39 }
  40
  41
  42 #ifdef HAVE_GSL
  43 std::vector<std::complex<double> > solve(Poly const & pp) {
  44     Poly p(pp);
  45     p.normalize();
  46     gsl_poly_complex_workspace * w
  47         = gsl_poly_complex_workspace_alloc (p.size());
  48
  49     gsl_complex_packed_ptr z = new double[p.degree()*2];
  50     double* a = new double[p.size()];
  51     for(unsigned int i = 0; i < p.size(); i++)
  52         a[i] = p[i];
  53     std::vector<std::complex<double> > roots;
  54     //roots.resize(p.degree());
  55
  56     gsl_poly_complex_solve (a, p.size(), w, z);
  57     delete[]a;
  58
  59     gsl_poly_complex_workspace_free (w);
  60
  61     for (unsigned int i = 0; i < p.degree(); i++) {
  62         roots.push_back(std::complex<double> (z[2*i] ,z[2*i+1]));
  63         //printf ("z%d = %+.18f %+.18f\n", i, z[2*i], z[2*i+1]);
  64     }
  65     delete[] z;
  66     return roots;
  67 }
  68
  69 std::vector<double > solve_reals(Poly const & p) {
  70     std::vector<std::complex<double> > roots = solve(p);
  71     std::vector<double> real_roots;
  72
  73     for(unsigned int i = 0; i < roots.size(); i++) {
  74         if(roots[i].imag() == 0) // should be more lenient perhaps
  75             real_roots.push_back(roots[i].real());
  76     }
  77     return real_roots;
  78 }
  79 #endif
  80
  81 double polish_root(Poly const & p, double guess, double tol) {
  82     Poly dp = derivative(p);
  83
  84     double fn = p(guess);
  85     while(fabs(fn) > tol) {
  86         guess -= fn/dp(guess);
  87         fn = p(guess);
  88     }
  89     return guess;
  90 }
  91
  92 Poly integral(Poly const & p) {
  93     Poly result;
  94
  95     result.reserve(p.size()+1);
  96     result.push_back(0); // arbitrary const
  97     for(unsigned i = 0; i < p.size(); i++) {
  98         result.push_back(p[i]/(i+1));
  99     }
 100     return result;
 101
 102 }
 103
 104 Poly derivative(Poly const & p) {
 105     Poly result;
 106
 107     if(p.size() <= 1)
 108         return Poly(0);
 109     result.reserve(p.size()-1);
 110     for(unsigned i = 1; i < p.size(); i++) {
 111         result.push_back(i*p[i]);
 112     }
 113     return result;
 114 }
 115
 116 Poly compose(Poly const & a, Poly const & b) {
 117     Poly result;
 118
 119     for(unsigned i = a.size(); i > 0; i--) {
 120         result = Poly(a[i-1]) + result * b;
 121     }
 122     return result;
 123
 124 }
 125
 126 /* This version is backwards - dividing taylor terms
 127 Poly divide(Poly const &a, Poly const &b, Poly &r) {
 128     Poly c;
 129     r = a; // remainder
 130
 131     const unsigned k = a.size();
 132     r.resize(k, 0);
 133     c.resize(k, 0);
 134
 135     for(unsigned i = 0; i < k; i++) {
 136         double ci = r[i]/b[0];
 137         c[i] += ci;
 138         Poly bb = ci*b;
 139         std::cout << ci <<"*" << b << ", r= " << r << std::endl;
 140         r -= bb.shifted(i);
 141     }
 142
 143     return c;
 144 }
 145 */
 146
 147 Poly divide(Poly const &a, Poly const &b, Poly &r) {
 148     Poly c;
 149     r = a; // remainder
 150     assert(b.size() > 0);
 151
 152     const unsigned k = a.degree();
 153     const unsigned l = b.degree();
 154     c.resize(k, 0.);
 155
 156     for(unsigned i = k; i >= l; i--) {
 157         //assert(i >= 0);
 158         double ci = r.back()/b.back();
 159         c[i-l] += ci;
 160         Poly bb = ci*b;
 161         //std::cout << ci <<"*(" << b.shifted(i-l) << ") = "
 162         //          << bb.shifted(i-l) << "     r= " << r << std::endl;
 163         r -= bb.shifted(i-l);
 164         r.pop_back();
 165     }
 166     //std::cout << "r= " << r << std::endl;
 167     r.normalize();
 168     c.normalize();
 169
 170     return c;
 171 }
 172
 173 Poly gcd(Poly const &a, Poly const &b, const double /*tol*/) {
 174     if(a.size() < b.size())
 175         return gcd(b, a);
 176     if(b.size() <= 0)
 177         return a;
 178     if(b.size() == 1)
 179         return a;
 180     Poly r;
 181     divide(a, b, r);
 182     return gcd(b, r);
 183 }
 184
 185
 186
 187 /*Poly divide_out_root(Poly const & p, double x) {
 188     assert(1);
 189     }*/
 190
 191 } //namespace Geom
 192
 193 /*
 194   Local Variables:
 195   mode:c++
 196   c-file-style:"stroustrup"
 197   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 198   indent-tabs-mode:nil
 199   fill-column:99
 200   End:
 201 */
 202 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :