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