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

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