src/2geom/poly-dk-solve.cpp

   1 #include <2geom/poly-dk-solve.h>
   2 #include <iterator>
   3
   4 /*** implementation of the Durand-Kerner method.  seems buggy*/
   5
   6 namespace Geom {
   7
   8 std::complex<double> evalu(Poly const & p, std::complex<double> x) {
   9     std::complex<double> result = 0;
  10     std::complex<double> xx = 1;
  11
  12     for(unsigned i = 0; i < p.size(); i++) {
  13         result += p[i]*xx;
  14         xx *= x;
  15     }
  16     return result;
  17 }
  18
  19 std::vector<std::complex<double> > DK(Poly const & ply, const double tol) {
  20     std::vector<std::complex<double> > roots;
  21     const int N = ply.degree();
  22
  23     std::complex<double> b(0.4, 0.9);
  24     std::complex<double> p = 1;
  25     for(int i = 0; i < N; i++) {
  26         roots.push_back(p);
  27         p *= b;
  28     }
  29     assert(roots.size() == ply.degree());
  30
  31     double error = 0;
  32     int i;
  33     for( i = 0; i < 30; i++) {
  34         error = 0;
  35         for(int r_i = 0; r_i < N; r_i++) {
  36             std::complex<double> denom = 1;
  37             std::complex<double> R = roots[r_i];
  38             for(int d_i = 0; d_i < N; d_i++) {
  39                 if(r_i != d_i)
  40                     denom *= R-roots[d_i];
  41             }
  42             assert(norm(denom) != 0);
  43             std::complex<double> dr = evalu(ply, R)/denom;
  44             error += norm(dr);
  45             roots[r_i] = R - dr;
  46         }
  47         /*std::copy(roots.begin(), roots.end(), std::ostream_iterator<std::complex<double> >(std::cout, ",\t"));
  48           std::cout << std::endl;*/
  49         if(error < tol)
  50             break;
  51     }
  52     //std::cout << error << ", " << i<< std::endl;
  53     return roots;
  54 }
  55
  56 } // namespace Geom
  57
  58 /*
  59   Local Variables:
  60   mode:c++
  61   c-file-style:"stroustrup"
  62   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
  63   indent-tabs-mode:nil
  64   fill-column:99
  65   End:
  66 */
  67 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :