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