src/2geom/poly.cpp

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