1 #include <2geom/poly-laguerre-solve.h>
2 #include <iterator>
5 namespace Geom {
7 typedef std::complex<double> cdouble;
9 cdouble laguerre_internal_complex(Poly const & p,
10 double x0,
11 double tol,
12 bool & quad_root) {
13 cdouble a = 2*tol;
14 cdouble xk = x0;
15 double n = p.degree();
16 quad_root = false;
17 const unsigned shuffle_rate = 10;
18 //static double shuffle[] = {0, 0.5, 0.25, 0.75, 0.125, 0.375, 0.625, 0.875, 1.0};
19 unsigned shuffle_counter = 0;
20 while(std::norm(a) > (tol*tol)) {
21 //std::cout << "xk = " << xk << std::endl;
22 cdouble b = p.back();
23 cdouble d = 0, f = 0;
24 double err = abs(b);
25 double abx = abs(xk);
26 for(int j = p.size()-2; j >= 0; j--) {
27 f = xk*f + d;
28 d = xk*d + b;
29 b = xk*b + p[j];
30 err = abs(b) + abx*err;
31 }
33 err *= 1e-7; // magic epsilon for convergence, should be computed from tol
35 cdouble px = b;
36 if(abs(b) < err)
37 return xk;
38 //if(std::norm(px) < tol*tol)
39 // return xk;
40 cdouble G = d / px;
41 cdouble H = G*G - f / px;
43 //std::cout << "G = " << G << "H = " << H;
44 cdouble radicand = (n - 1)*(n*H-G*G);
45 //assert(radicand.real() > 0);
46 if(radicand.real() < 0)
47 quad_root = true;
48 //std::cout << "radicand = " << radicand << std::endl;
49 if(G.real() < 0) // here we try to maximise the denominator avoiding cancellation
50 a = - sqrt(radicand);
51 else
52 a = sqrt(radicand);
53 //std::cout << "a = " << a << std::endl;
54 a = n / (a + G);
55 //std::cout << "a = " << a << std::endl;
56 //if(shuffle_counter % shuffle_rate == 0)
57 //a *= shuffle[shuffle_counter / shuffle_rate];
58 xk -= a;
59 shuffle_counter++;
60 if(shuffle_counter >= 90)
61 break;
62 }
63 //std::cout << "xk = " << xk << std::endl;
64 return xk;
65 }
67 double laguerre_internal(Poly const & p,
68 Poly const & pp,
69 Poly const & ppp,
70 double x0,
71 double tol,
72 bool & quad_root) {
73 double a = 2*tol;
74 double xk = x0;
75 double n = p.degree();
76 quad_root = false;
77 while(a*a > (tol*tol)) {
78 //std::cout << "xk = " << xk << std::endl;
79 double px = p(xk);
80 if(px*px < tol*tol)
81 return xk;
82 double G = pp(xk) / px;
83 double H = G*G - ppp(xk) / px;
85 //std::cout << "G = " << G << "H = " << H;
86 double radicand = (n - 1)*(n*H-G*G);
87 assert(radicand > 0);
88 //std::cout << "radicand = " << radicand << std::endl;
89 if(G < 0) // here we try to maximise the denominator avoiding cancellation
90 a = - sqrt(radicand);
91 else
92 a = sqrt(radicand);
93 //std::cout << "a = " << a << std::endl;
94 a = n / (a + G);
95 //std::cout << "a = " << a << std::endl;
96 xk -= a;
97 }
98 //std::cout << "xk = " << xk << std::endl;
99 return xk;
100 }
103 std::vector<cdouble >
104 laguerre(Poly p, const double tol) {
105 std::vector<cdouble > solutions;
106 //std::cout << "p = " << p << " = ";
107 while(p.size() > 1)
108 {
109 double x0 = 0;
110 bool quad_root = false;
111 cdouble sol = laguerre_internal_complex(p, x0, tol, quad_root);
112 //if(abs(sol) > 1) break;
113 Poly dvs;
114 if(quad_root) {
115 dvs.push_back((sol*conj(sol)).real());
116 dvs.push_back(-(sol + conj(sol)).real());
117 dvs.push_back(1.0);
118 //std::cout << "(" << dvs << ")";
119 //solutions.push_back(sol);
120 //solutions.push_back(conj(sol));
121 } else {
122 //std::cout << sol << std::endl;
123 dvs.push_back(-sol.real());
124 dvs.push_back(1.0);
125 solutions.push_back(sol);
126 //std::cout << "(" << dvs << ")";
127 }
128 Poly r;
129 p = divide(p, dvs, r);
130 //std::cout << r << std::endl;
131 }
132 return solutions;
133 }
135 std::vector<double>
136 laguerre_real_interval(Poly const & /*ply*/,
137 const double /*lo*/, const double /*hi*/,
138 const double /*tol*/)
139 {
140 /* not implemented*/
141 assert(false);
142 std::vector<double> result;
143 return result;
144 }
146 } // namespace Geom
148 /*
149 Local Variables:
150 mode:c++
151 c-file-style:"stroustrup"
152 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
153 indent-tabs-mode:nil
154 fill-column:99
155 End:
156 */
157 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :