1 /* Roots.cpp
2 *
3 * Copyright (C) 1993-2020 David Weenink
4 *
5 * This code is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or (at
8 * your option) any later version.
9 *
10 * This code is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this work. If not, see <http://www.gnu.org/licenses/>.
17 */
18
19 #include "NUMlapack.h"
20 #include "NUMmachar.h"
21 #include "Polynomial.h"
22 #include "Roots.h"
23
24 #include "oo_DESTROY.h"
25 #include "Roots_def.h"
26 #include "oo_COPY.h"
27 #include "Roots_def.h"
28 #include "oo_EQUAL.h"
29 #include "Roots_def.h"
30 #include "oo_CAN_WRITE_AS_ENCODING.h"
31 #include "Roots_def.h"
32 #include "oo_READ_TEXT.h"
33 #include "Roots_def.h"
34 #include "oo_WRITE_TEXT.h"
35 #include "Roots_def.h"
36 #include "oo_READ_BINARY.h"
37 #include "Roots_def.h"
38 #include "oo_WRITE_BINARY.h"
39 #include "Roots_def.h"
40 #include "oo_DESCRIPTION.h"
41 #include "Roots_def.h"
42
43
44 Thing_implement (Roots, Daata, 1);
45
v_info()46 void structRoots :: v_info () {
47 structDaata :: v_info ();
48 MelderInfo_writeLine (U"Number of roots: ", numberOfRoots);
49 }
50
Roots_create(integer numberOfRoots)51 autoRoots Roots_create (integer numberOfRoots) {
52 try {
53 autoRoots me = Thing_new (Roots);
54 my numberOfRoots = numberOfRoots;
55 my roots = zero_COMPVEC (numberOfRoots);
56 return me;
57 } catch (MelderError) {
58 Melder_throw (U"Roots not created.");
59 }
60 }
61
Roots_getNumberOfRoots(Roots me)62 integer Roots_getNumberOfRoots (Roots me) {
63 return my numberOfRoots;
64 }
65
Roots_getRoot(Roots me,integer index)66 dcomplex Roots_getRoot (Roots me, integer index) {
67 Melder_require (index > 0 && index <= my numberOfRoots,
68 U"Root index out of range.");
69 return my roots [index];
70 }
71
Roots_fixIntoUnitCircle(Roots me)72 void Roots_fixIntoUnitCircle (Roots me) {
73 dcomplex z10 { 1.0, 0.0 };
74 for (integer iroot = 1; iroot <= my numberOfRoots; iroot ++)
75 if (abs (my roots [iroot]) > 1.0)
76 my roots [iroot] = z10 / conj (my roots [iroot]);
77 }
78
NUMdcvector_extrema_re(COMPVEC const & v,integer lo,integer hi,double * out_min,double * out_max)79 static void NUMdcvector_extrema_re (COMPVEC const& v, integer lo, integer hi, double *out_min, double *out_max) {
80 double min = v [lo].real(), max = v [lo].real();
81 for (integer i = lo + 1; i <= hi; i ++)
82 if (v [i].real() < min)
83 min = v [i].real();
84 else if (v [i].real() > max)
85 max = v [i].real();
86 if (out_min)
87 *out_min = min;
88 if (out_max)
89 *out_max = max;
90 }
91
NUMdcvector_extrema_im(COMPVEC const & v,integer lo,integer hi,double * out_min,double * out_max)92 static void NUMdcvector_extrema_im (COMPVEC const& v, integer lo, integer hi, double *out_min, double *out_max) {
93 double min = v [lo].imag(), max = v [lo].imag();
94 for (integer i = lo + 1; i <= hi; i ++)
95 if (v [i].imag() < min)
96 min = v [i].imag();
97 else if (v [i].imag() > max)
98 max = v [i].imag();
99 if (out_min)
100 *out_min = min;
101 if (out_max)
102 *out_max = max;
103 }
104
Roots_draw(Roots me,Graphics g,double rmin,double rmax,double imin,double imax,conststring32 symbol,double fontSize,bool garnish)105 void Roots_draw (Roots me, Graphics g, double rmin, double rmax, double imin, double imax,
106 conststring32 symbol, double fontSize, bool garnish) {
107 const double oldFontSize = Graphics_inqFontSize (g);
108 const double eps = 1e-6;
109
110 if (rmax <= rmin)
111 NUMdcvector_extrema_re (my roots.get(), 1, my numberOfRoots, & rmin, & rmax);
112
113 double denominator = fabs (rmax) > fabs (rmin) ? fabs (rmax) : fabs (rmin);
114 if (denominator == 0.0)
115 denominator = 1.0;
116 if (fabs ((rmax - rmin) / denominator) < eps) {
117 rmin -= 1.0;
118 rmax += 1.0;
119 }
120 if (imax <= imin)
121 NUMdcvector_extrema_im (my roots.get(), 1, my numberOfRoots, & imin, & imax);
122 denominator = fabs (imax) > fabs (imin) ? fabs (imax) : fabs (imin);
123 if (denominator == 0.0)
124 denominator = 1.0;
125 if (fabs ((imax - imin) / denominator) < eps) {
126 imin -= 1;
127 imax += 1;
128 }
129 Graphics_setInner (g);
130 Graphics_setWindow (g, rmin, rmax, imin, imax);
131 Graphics_setFontSize (g, fontSize);
132 Graphics_setTextAlignment (g, kGraphics_horizontalAlignment::CENTRE, Graphics_HALF);
133 for (integer i = 1; i <= my numberOfRoots; i ++) {
134 const double re = my roots [i].real(), im = my roots [i].imag();
135 if (re >= rmin && re <= rmax && im >= imin && im <= imax)
136 Graphics_text (g, re, im, symbol);
137 }
138 Graphics_setFontSize (g, oldFontSize);
139 Graphics_unsetInner (g);
140 if (garnish) {
141 Graphics_drawInnerBox (g);
142 if (rmin * rmax < 0.0)
143 Graphics_markLeft (g, 0.0, true, true, true, U"0");
144 if (imin * imax < 0.0)
145 Graphics_markBottom (g, 0.0, true, true, true, U"0");
146 Graphics_marksLeft (g, 2, true, true, false);
147 Graphics_textLeft (g, true, U"Imaginary part");
148 Graphics_marksBottom (g, 2, true, true, false);
149 Graphics_textBottom (g, true, U"Real part");
150 }
151 }
152
Polynomial_to_Roots(Polynomial me)153 autoRoots Polynomial_to_Roots (Polynomial me) {
154 try {
155 Melder_assert (my numberOfCoefficients == my coefficients.size); // check invariant
156 integer np1 = my numberOfCoefficients, n = np1 - 1;
157 Melder_require (n > 0,
158 U"Cannot find roots of a constant function.");
159 /*
160 Allocate storage for a special upper Hessenberg matrix (n * n)
161 The roots of a polynomial are the eigenvalues of an
162 upper Hessenberg matrix with the coefficients of the polynomial.
163 See for example the introduction in:
164 G.S. Ammar, D. Calvetti, W.B. Gragg, L. Reichel (2001):
165 "Polynomial zero finders based on Szegö polynomials.",
166 Journal of Computational and Applied Mathematics 127: 1-–16.
167 */
168 autoVEC wr = raw_VEC (n);
169 autoVEC wi = raw_VEC (n);
170 autoMAT upperHessenberg = zero_MAT (n, n);
171 MATVU uh_CM (upperHessenberg.get());
172 uh_CM.rowStride = 1; uh_CM.colStride = n;
173
174 uh_CM [1] [n] = - (my coefficients [1] / my coefficients [np1]);
175 for (integer irow = 2; irow <= n; irow ++) {
176 uh_CM [irow] [n] = - (my coefficients [irow] / my coefficients [np1]);
177 uh_CM [irow] [irow - 1] = 1.0;
178 }
179 /*
180 Find out the working storage needed
181 */
182 double wtmp;
183 integer lwork = -1, info;
184 NUMlapack_dhseqr_ ("E", "N", n, 1, n, & upperHessenberg [1] [1], n, & wr [1], & wi [1], nullptr, n, & wtmp, lwork, & info);
185 lwork = Melder_roundUp (wtmp);
186 autoVEC work = raw_VEC (lwork);
187 /*
188 Find eigenvalues/roots.
189 */
190 NUMlapack_dhseqr_ ("E", "N", n, 1, n, & upperHessenberg [1] [1], n, & wr [1], & wi [1], nullptr, n, & work [1], lwork, & info);
191
192 integer numberOfEigenvaluesFound = n, ioffset = 0;
193 if (info > 0) {
194 /*
195 if INFO = i, NUMlapack_dhseqr_ failed to compute all of the eigenvalues. Elements i+1:n of
196 WR and WI contain those eigenvalues which have been successfully computed
197 */
198 numberOfEigenvaluesFound -= info;
199 Melder_require (numberOfEigenvaluesFound > 0,
200 U"No eigenvalues found.");
201 ioffset = info;
202 } else if (info < 0) {
203 Melder_throw (U"NUMlapack_dhseqr_ returns error ", info, U".");
204 }
205
206 autoRoots thee = Roots_create (numberOfEigenvaluesFound);
207 for (integer i = 1; i <= numberOfEigenvaluesFound; i ++) {
208 thy roots [i]. real (wr [ioffset + i]);
209 thy roots [i]. imag (wi [ioffset + i]);
210 }
211 Roots_Polynomial_polish (thee.get(), me);
212 return thee;
213 } catch (MelderError) {
214 Melder_throw (me, U": no roots can be calculated.");
215 }
216 }
217
218 /*
219 workspace.size >= n * n + 9 * n =
220 n * n ; for hessenberg matrix
221 + 2 * n ; for real and imaginary parts
222 + 6 * n ; the maximum for dhseqr_
223 */
Polynomial_into_Roots(Polynomial me,Roots r,VEC const & workspace)224 void Polynomial_into_Roots (Polynomial me, Roots r, VEC const& workspace) {
225 Melder_assert (my numberOfCoefficients == my coefficients.size); // check invariant
226 r -> roots.resize (0);
227 r -> numberOfRoots = r -> roots.size;
228 integer np1 = my numberOfCoefficients, n = np1 - 1;
229 if (n == 0)
230 return;
231 /*
232 Use the workspace reserve storage for Hessenberg matrix (n * n)
233 */
234
235 MAT upperHessenberg = MAT (& workspace [1], n, n);
236 upperHessenberg <<= 0.0;
237 MATVU uh_CM (upperHessenberg);
238 uh_CM.rowStride = 1; uh_CM.colStride = n;
239 uh_CM [1] [n] = - (my coefficients [1] / my coefficients [np1]);
240 for (integer irow = 2; irow <= n; irow ++) {
241 uh_CM [irow] [n] = - (my coefficients [irow] / my coefficients [np1]);
242 uh_CM [irow] [irow - 1] = 1.0;
243 }
244 /*
245 We don't need to find out size of the working storage needed because for the current version 3.1.1.1
246 of NUMlapack_dhseqr (20200313) its size equals maximally 6*n.
247 */
248 integer endIndex = n * n;
249 VEC wr = workspace.part (endIndex + 1, endIndex + n);
250 endIndex += n;
251 VEC wi = workspace.part (endIndex + 1, endIndex + n);
252 endIndex += n;
253 VEC work = workspace.part (endIndex + 1, workspace.size);
254 Melder_assert (work.size >= 6 * n);
255 integer lwork = work.size, info;
256 NUMlapack_dhseqr_ ("E", "N", n, 1, n, & uh_CM [1] [1], n, & wr [1], & wi [1], nullptr, n, & work [1], lwork, & info);
257 integer numberOfEigenvaluesFound = n, ioffset = 0;
258 if (info > 0) {
259 /*
260 if INFO = i, NUMlapack_dhseqr failed to compute all of the eigenvalues. Elements i+1:n of
261 WR and WI contain those eigenvalues which have been successfully computed
262 */
263 numberOfEigenvaluesFound -= info;
264 Melder_require (numberOfEigenvaluesFound > 0,
265 U"No eigenvalues found.");
266 ioffset = info;
267 } else if (info < 0) {
268 Melder_throw (U"NUMlapack_dhseqr_ returns error ", info, U".");
269 }
270
271 for (integer i = 1; i <= numberOfEigenvaluesFound; i ++) {
272 dcomplex *root = r -> roots . append();
273 (*root) . real (wr [ioffset + i]);
274 (*root) . imag (wi [ioffset + i]);
275 }
276 r -> numberOfRoots = r -> roots . size; // maintain invariant
277 Roots_Polynomial_polish (r, me);
278 }
279
Roots_sort(Roots me)280 void Roots_sort (Roots me) {
281 (void) me;
282 }
283
284 /* Get value and derivative */
Polynomial_evaluateWithDerivative_z(Polynomial me,dcomplex * in_z,dcomplex * out_p,dcomplex * out_dp)285 static void Polynomial_evaluateWithDerivative_z (Polynomial me, dcomplex *in_z, dcomplex *out_p, dcomplex *out_dp) {
286 longdouble pr = my coefficients [my numberOfCoefficients], pi = 0.0;
287 longdouble dpr = 0.0, dpi = 0.0, x = in_z->real(), y = in_z->imag();
288
289 for (integer i = my numberOfCoefficients - 1; i > 0; i --) {
290 longdouble tr = dpr;
291 dpr = dpr * x - dpi * y + pr;
292 dpi = tr * y + dpi * x + pi;
293 tr = pr;
294 pr = pr * x - pi * y + my coefficients [i];
295 pi = tr * y + pi * x;
296 }
297 if (out_p)
298 *out_p = { (double) pr, (double) pi };
299 if (out_dp)
300 *out_dp = { (double) dpr, (double) dpi };
301 }
302
Polynomial_polish_complexroot_nr(Polynomial me,dcomplex * root,integer maxit)303 static void Polynomial_polish_complexroot_nr (Polynomial me, dcomplex *root, integer maxit) {
304 if (! NUMfpp)
305 NUMmachar ();
306 dcomplex zbest = *root;
307 double ymin = 1e308;
308 for (integer iter = 1; iter <= maxit; iter ++) {
309 dcomplex y, dy;
310 Polynomial_evaluateWithDerivative_z (me, root, & y, & dy);
311 const double fabsy = abs (y);
312 if (fabsy > ymin || fabs (fabsy - ymin) < NUMfpp -> eps) {
313 /*
314 We stop, because the approximation is getting worse.
315 Return the previous (hitherto best) value for z.
316 */
317 *root = zbest;
318 return;
319 }
320 ymin = fabsy;
321 zbest = *root;
322 if (abs (dy) == 0.0)
323 return;
324 const dcomplex dz = y / dy; // Newton-Raphson
325 *root -= dz;
326 }
327 // Melder_throw (U"Maximum number of iterations exceeded.");
328 }
329
Polynomial_polish_realroot(Polynomial me,double x,integer maxit)330 static double Polynomial_polish_realroot (Polynomial me, double x, integer maxit) {
331 if (! NUMfpp)
332 NUMmachar ();
333 double xbest = x, ymin = 1e308;
334 for (integer iter = 1; iter <= maxit; iter ++) {
335 double y, dy;
336 Polynomial_evaluateWithDerivative (me, x, & y, & dy);
337 const double fabsy = fabs (y);
338 if (fabsy > ymin || fabs (fabsy - ymin) < NUMfpp -> eps) {
339 /*
340 We stop, because the approximation is getting worse or we cannot get any closer.
341 Return the previous (hitherto best) value for x.
342 */
343 x = xbest;
344 return x;
345 }
346 ymin = fabsy;
347 xbest = x;
348 if (fabs (dy) == 0.0)
349 return x;
350 const double dx = y / dy; // Newton-Raphson
351 x -= dx;
352 }
353 return x;
354 // Melder_throw (U"Maximum number of iterations exceeded.");
355 }
356
357 // Precondition: complex roots occur in pairs (a,bi), (a,-bi) with b>0
Roots_Polynomial_polish(Roots me,Polynomial thee)358 void Roots_Polynomial_polish (Roots me, Polynomial thee) {
359 const integer maxit = 80;
360 integer i = 1;
361 while (i <= my numberOfRoots) {
362 const double im = my roots [i].imag(), re = my roots [i].real();
363 if (im != 0.0) {
364 Polynomial_polish_complexroot_nr (thee, & my roots [i], maxit);
365 if (i < my numberOfRoots && im == - my roots [i + 1].imag() && re == my roots [i + 1].real()) {
366 my roots [i + 1]. real (my roots [i].real());
367 my roots [i + 1]. imag (- my roots [i].imag());
368 i ++;
369 }
370 } else {
371 my roots [i]. real (Polynomial_polish_realroot (thee, my roots [i].real(), maxit));
372 }
373 i ++;
374 }
375 }
376
Roots_to_Polynomial(Roots me,bool rootsAreReal)377 autoPolynomial Roots_to_Polynomial (Roots me, bool rootsAreReal) {
378 try {
379 (void) me;
380 autoPolynomial thee;
381 if (! rootsAreReal)
382 throw MelderError();
383 return thee;
384 } catch (MelderError) {
385 Melder_throw (U"Not implemented yet");
386 }
387 }
388
389 /* End of file Roots.cpp */
390