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