1 #include <2geom/poly-dk-solve.h>
2 #include <iterator>
4 /*** implementation of the Durand-Kerner method. seems buggy*/
6 namespace Geom {
8 std::complex<double> evalu(Poly const & p, std::complex<double> x) {
9 std::complex<double> result = 0;
10 std::complex<double> xx = 1;
12 for(unsigned i = 0; i < p.size(); i++) {
13 result += p[i]*xx;
14 xx *= x;
15 }
16 return result;
17 }
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();
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());
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 }
56 } // namespace Geom
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 :