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