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