src/2geom/sturm.h

   1 #ifndef LIB2GEOM_STURM_HEADER
   2 #define LIB2GEOM_STURM_HEADER
   3
   4 #include <2geom/poly.h>
   5 #include <2geom/utils.h>
   6
   7 namespace Geom {
   8
   9 class sturm : public std::vector<Poly>{
  10 public:
  11     sturm(Poly const &X) {
  12         push_back(X);
  13         push_back(derivative(X));
  14         Poly Xi = back();
  15         Poly Xim1 = X;
  16         std::cout << "sturm:\n" << Xim1 << std::endl;
  17         std::cout << Xi << std::endl;
  18         while(Xi.size() > 1) {
  19             Poly r;
  20             divide(Xim1, Xi, r);
  21             std::cout << r << std::endl;
  22             assert(r.size() < Xi.size());
  23             Xim1 = Xi;
  24             Xi = -r;
  25             assert(Xim1.size() > Xi.size());
  26             push_back(Xi);
  27         }
  28     }
  29
  30     unsigned count_signs(double t) {
  31         unsigned n_signs = 0;/*  Number of sign-changes */
  32         const double big = 1e20; // a number such that practical polys would overflow on evaluation
  33         if(t >= big) {
  34             int old_sign = sgn((*this)[0].back());
  35             for (unsigned i = 1; i < size(); i++) {
  36                 int sign = sgn((*this)[i].back());
  37                 if (sign != old_sign)
  38                     n_signs++;
  39                 old_sign = sign;
  40             }
  41         } else {
  42             int old_sign = sgn((*this)[0].eval(t));
  43             for (unsigned i = 1; i < size(); i++) {
  44                 int sign = sgn((*this)[i].eval(t));
  45                 if (sign != old_sign)
  46                     n_signs++;
  47                 old_sign = sign;
  48             }
  49         }
  50         return n_signs;
  51     }
  52
  53     unsigned n_roots_between(double l, double r) {
  54         return count_signs(l) - count_signs(r);
  55     }
  56 };
  57
  58 } //namespace Geom
  59
  60 #endif
  61 /*
  62   Local Variables:
  63   mode:c++
  64   c-file-style:"stroustrup"
  65   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
  66   indent-tabs-mode:nil
  67   fill-column:99
  68   End:
  69 */
  70 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :