1 #ifndef LIB2GEOM_STURM_HEADER
2 #define LIB2GEOM_STURM_HEADER
4 #include <2geom/poly.h>
5 #include <2geom/utils.h>
7 namespace Geom {
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 }
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 }
53 unsigned n_roots_between(double l, double r) {
54 return count_signs(l) - count_signs(r);
55 }
56 };
58 } //namespace Geom
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 :