1 /* poly/test.c
2 *
3 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
4 *
5 * This program 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 3 of the License, or (at
8 * your option) any later version.
9 *
10 * This program 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 program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 */
19
20 #include <config.h>
21 #include <stdlib.h>
22 #include <gsl/gsl_math.h>
23 #include <gsl/gsl_test.h>
24 #include <gsl/gsl_ieee_utils.h>
25 #include <gsl/gsl_poly.h>
26 #include <gsl/gsl_heapsort.h>
27
28 /* sort by Re(z) then by Im(z) */
29 static int
cmp_cplx(const double * a,const double * b)30 cmp_cplx(const double *a, const double *b)
31 {
32 double r = a[0] - b[0];
33
34 if (r == 0.0)
35 {
36 double t = a[1] - b[1];
37 return t < 0.0 ? -1 : t > 0.0 ? 1 : 0;
38 }
39 else if (r < 0.0)
40 return -1;
41 else
42 return 1;
43 }
44
45 int
main(void)46 main (void)
47 {
48 const double eps = 100.0 * GSL_DBL_EPSILON;
49
50 gsl_ieee_env_setup ();
51
52 /* Polynomial evaluation */
53
54 {
55 double x, y;
56 double c[3] = { 1.0, 0.5, 0.3 };
57 x = 0.5;
58 y = gsl_poly_eval (c, 3, x);
59 gsl_test_rel (y, 1 + 0.5 * x + 0.3 * x * x, eps,
60 "gsl_poly_eval({1, 0.5, 0.3}, 0.5)");
61 }
62
63 {
64 double x, y;
65 double d[11] = { 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1 };
66 x = 1.0;
67 y = gsl_poly_eval (d, 11, x);
68 gsl_test_rel (y, 1.0, eps,
69 "gsl_poly_eval({1,-1, 1, -1, 1, -1, 1, -1, 1, -1, 1}, 1.0)");
70
71 }
72
73 {
74 gsl_complex x, y;
75 double c[1] = {0.3};
76 GSL_SET_REAL (&x, 0.75);
77 GSL_SET_IMAG (&x, 1.2);
78 y = gsl_poly_complex_eval (c, 1, x);
79
80 gsl_test_rel (GSL_REAL (y), 0.3, eps, "y.real, gsl_poly_complex_eval ({0.3}, 0.75 + 1.2i)");
81 gsl_test_rel (GSL_IMAG (y), 0.0, eps, "y.imag, gsl_poly_complex_eval ({0.3}, 0.75 + 1.2i)");
82 }
83
84 {
85 gsl_complex x, y;
86 double c[4] = {2.1, -1.34, 0.76, 0.45};
87 GSL_SET_REAL (&x, 0.49);
88 GSL_SET_IMAG (&x, 0.95);
89 y = gsl_poly_complex_eval (c, 4, x);
90
91 gsl_test_rel (GSL_REAL (y), 0.3959143, eps, "y.real, gsl_poly_complex_eval ({2.1, -1.34, 0.76, 0.45}, 0.49 + 0.95i)");
92 gsl_test_rel (GSL_IMAG (y), -0.6433305, eps, "y.imag, gsl_poly_complex_eval ({2.1, -1.34, 0.76, 0.45}, 0.49 + 0.95i)");
93 }
94
95 {
96 gsl_complex x, y;
97 gsl_complex c[1];
98 GSL_SET_REAL (&c[0], 0.674);
99 GSL_SET_IMAG (&c[0], -1.423);
100 GSL_SET_REAL (&x, -1.44);
101 GSL_SET_IMAG (&x, 9.55);
102 y = gsl_complex_poly_complex_eval (c, 1, x);
103
104 gsl_test_rel (GSL_REAL (y), 0.674, eps, "y.real, gsl_complex_poly_complex_eval ({0.674 - 1.423i}, -1.44 + 9.55i)");
105 gsl_test_rel (GSL_IMAG (y), -1.423, eps, "y.imag, gsl_complex_poly_complex_eval ({0.674 - 1.423i}, -1.44 + 9.55i)");
106 }
107
108 {
109 gsl_complex x, y;
110 gsl_complex c[4];
111 GSL_SET_REAL (&c[0], -2.31);
112 GSL_SET_IMAG (&c[0], 0.44);
113 GSL_SET_REAL (&c[1], 4.21);
114 GSL_SET_IMAG (&c[1], -3.19);
115 GSL_SET_REAL (&c[2], 0.93);
116 GSL_SET_IMAG (&c[2], 1.04);
117 GSL_SET_REAL (&c[3], -0.42);
118 GSL_SET_IMAG (&c[3], 0.68);
119 GSL_SET_REAL (&x, 0.49);
120 GSL_SET_IMAG (&x, 0.95);
121 y = gsl_complex_poly_complex_eval (c, 4, x);
122
123 gsl_test_rel (GSL_REAL (y), 1.82462012, eps, "y.real, gsl_complex_poly_complex_eval ({-2.31 + 0.44i, 4.21 - 3.19i, 0.93 + 1.04i, -0.42 + 0.68i}, 0.49 + 0.95i)");
124 gsl_test_rel (GSL_IMAG (y), 2.30389412, eps, "y.imag, gsl_complex_poly_complex_eval ({-2.31 + 0.44i, 4.21 - 3.19i, 0.93 + 1.04i, -0.42 + 0.68i}, 0.49 + 0.95i)");
125 }
126
127 /* Quadratic */
128
129 {
130 double x0, x1;
131
132 int n = gsl_poly_solve_quadratic (4.0, -20.0, 26.0, &x0, &x1);
133
134 gsl_test (n != 0, "gsl_poly_solve_quadratic, no roots, (2x - 5)^2 = -1");
135 }
136
137 {
138 double x0, x1;
139
140 int n = gsl_poly_solve_quadratic (4.0, -20.0, 25.0, &x0, &x1);
141
142 gsl_test (n != 2, "gsl_poly_solve_quadratic, one root, (2x - 5)^2 = 0");
143 gsl_test_rel (x0, 2.5, 1e-9, "x0, (2x - 5)^2 = 0");
144 gsl_test_rel (x1, 2.5, 1e-9, "x1, (2x - 5)^2 = 0");
145 gsl_test (x0 != x1, "x0 == x1, (2x - 5)^2 = 0");
146 }
147
148 {
149 double x0, x1;
150
151 int n = gsl_poly_solve_quadratic (4.0, -20.0, 21.0, &x0, &x1);
152
153 gsl_test (n != 2, "gsl_poly_solve_quadratic, two roots, (2x - 5)^2 = 4");
154 gsl_test_rel (x0, 1.5, 1e-9, "x0, (2x - 5)^2 = 4");
155 gsl_test_rel (x1, 3.5, 1e-9, "x1, (2x - 5)^2 = 4");
156 }
157
158 {
159 double x0, x1;
160
161 int n = gsl_poly_solve_quadratic (4.0, 7.0, 0.0, &x0, &x1);
162
163 gsl_test (n != 2, "gsl_poly_solve_quadratic, two roots, x(4x + 7) = 0");
164 gsl_test_rel (x0, -1.75, 1e-9, "x0, x(4x + 7) = 0");
165 gsl_test_rel (x1, 0.0, 1e-9, "x1, x(4x + 7) = 0");
166 }
167
168 {
169 double x0, x1;
170
171 int n = gsl_poly_solve_quadratic (5.0, 0.0, -20.0, &x0, &x1);
172
173 gsl_test (n != 2,
174 "gsl_poly_solve_quadratic, two roots b = 0, 5 x^2 = 20");
175 gsl_test_rel (x0, -2.0, 1e-9, "x0, 5 x^2 = 20");
176 gsl_test_rel (x1, 2.0, 1e-9, "x1, 5 x^2 = 20");
177 }
178
179
180 {
181 double x0, x1;
182
183 int n = gsl_poly_solve_quadratic (0.0, 3.0, -21.0, &x0, &x1);
184
185 gsl_test (n != 1,
186 "gsl_poly_solve_quadratic, one root (linear) 3 x - 21 = 0");
187 gsl_test_rel (x0, 7.0, 1e-9, "x0, 3x - 21 = 0");
188 }
189
190
191 {
192 double x0, x1;
193 int n = gsl_poly_solve_quadratic (0.0, 0.0, 1.0, &x0, &x1);
194
195 gsl_test (n != 0,
196 "gsl_poly_solve_quadratic, no roots 1 = 0");
197 }
198
199
200 /* Cubic */
201
202 {
203 double x0, x1, x2;
204
205 int n = gsl_poly_solve_cubic (0.0, 0.0, -27.0, &x0, &x1, &x2);
206
207 gsl_test (n != 1, "gsl_poly_solve_cubic, one root, x^3 = 27");
208 gsl_test_rel (x0, 3.0, 1e-9, "x0, x^3 = 27");
209 }
210
211 {
212 double x0, x1, x2;
213
214 int n = gsl_poly_solve_cubic (-51.0, 867.0, -4913.0, &x0, &x1, &x2);
215
216 gsl_test (n != 3, "gsl_poly_solve_cubic, three roots, (x-17)^3=0");
217 gsl_test_rel (x0, 17.0, 1e-9, "x0, (x-17)^3=0");
218 gsl_test_rel (x1, 17.0, 1e-9, "x1, (x-17)^3=0");
219 gsl_test_rel (x2, 17.0, 1e-9, "x2, (x-17)^3=0");
220 }
221
222 {
223 double x0, x1, x2;
224
225 int n = gsl_poly_solve_cubic (-57.0, 1071.0, -6647.0, &x0, &x1, &x2);
226
227 gsl_test (n != 3,
228 "gsl_poly_solve_cubic, three roots, (x-17)(x-17)(x-23)=0");
229 gsl_test_rel (x0, 17.0, 1e-9, "x0, (x-17)(x-17)(x-23)=0");
230 gsl_test_rel (x1, 17.0, 1e-9, "x1, (x-17)(x-17)(x-23)=0");
231 gsl_test_rel (x2, 23.0, 1e-9, "x2, (x-17)(x-17)(x-23)=0");
232 }
233
234 {
235 double x0, x1, x2;
236
237 int n = gsl_poly_solve_cubic (-11.0, -493.0, +6647.0, &x0, &x1, &x2);
238
239 gsl_test (n != 3,
240 "gsl_poly_solve_cubic, three roots, (x+23)(x-17)(x-17)=0");
241 gsl_test_rel (x0, -23.0, 1e-9, "x0, (x+23)(x-17)(x-17)=0");
242 gsl_test_rel (x1, 17.0, 1e-9, "x1, (x+23)(x-17)(x-17)=0");
243 gsl_test_rel (x2, 17.0, 1e-9, "x2, (x+23)(x-17)(x-17)=0");
244 }
245
246 {
247 double x0, x1, x2;
248
249 int n = gsl_poly_solve_cubic (-143.0, 5087.0, -50065.0, &x0, &x1, &x2);
250
251 gsl_test (n != 3,
252 "gsl_poly_solve_cubic, three roots, (x-17)(x-31)(x-95)=0");
253 gsl_test_rel (x0, 17.0, 1e-9, "x0, (x-17)(x-31)(x-95)=0");
254 gsl_test_rel (x1, 31.0, 1e-9, "x1, (x-17)(x-31)(x-95)=0");
255 gsl_test_rel (x2, 95.0, 1e-9, "x2, (x-17)(x-31)(x-95)=0");
256 }
257
258 {
259 double x0, x1, x2;
260
261 int n = gsl_poly_solve_cubic (-109.0, 803.0, 50065.0, &x0, &x1, &x2);
262
263 gsl_test (n != 3,
264 "gsl_poly_solve_cubic, three roots, (x+17)(x-31)(x-95)=0");
265 gsl_test_rel (x0, -17.0, 1e-9, "x0, (x+17)(x-31)(x-95)=0");
266 gsl_test_rel (x1, 31.0, 1e-9, "x1, (x+17)(x-31)(x-95)=0");
267 gsl_test_rel (x2, 95.0, 1e-9, "x2, (x+17)(x-31)(x-95)=0");
268 }
269
270 /* Quadratic with complex roots */
271
272 {
273 gsl_complex z0, z1;
274
275 int n = gsl_poly_complex_solve_quadratic (4.0, -20.0, 26.0, &z0, &z1);
276
277 gsl_test (n != 2,
278 "gsl_poly_complex_solve_quadratic, 2 roots (2x - 5)^2 = -1");
279 gsl_test_rel (GSL_REAL (z0), 2.5, 1e-9, "z0.real, (2x - 5)^2 = -1");
280 gsl_test_rel (GSL_IMAG (z0), -0.5, 1e-9, "z0.imag, (2x - 5)^2 = -1");
281
282 gsl_test_rel (GSL_REAL (z1), 2.5, 1e-9, "z1.real, (2x - 5)^2 = -1");
283 gsl_test_rel (GSL_IMAG (z1), 0.5, 1e-9, "z1.imag, (2x - 5)^2 = -1");
284 }
285
286 {
287 gsl_complex z0, z1;
288
289 int n = gsl_poly_complex_solve_quadratic (4.0, -20.0, 25.0, &z0, &z1);
290
291 gsl_test (n != 2,
292 "gsl_poly_complex_solve_quadratic, one root, (2x - 5)^2 = 0");
293 gsl_test_rel (GSL_REAL (z0), 2.5, 1e-9, "z0.real, (2x - 5)^2 = 0");
294 gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag (2x - 5)^2 = 0");
295 gsl_test_rel (GSL_REAL (z1), 2.5, 1e-9, "z1.real, (2x - 5)^2 = 0");
296 gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag (2x - 5)^2 = 0");
297 gsl_test (GSL_REAL (z0) != GSL_REAL (z1),
298 "z0.real == z1.real, (2x - 5)^2 = 0");
299 gsl_test (GSL_IMAG (z0) != GSL_IMAG (z1),
300 "z0.imag == z1.imag, (2x - 5)^2 = 0");
301 }
302
303 {
304 gsl_complex z0, z1;
305
306 int n = gsl_poly_complex_solve_quadratic (4.0, -20.0, 21.0, &z0, &z1);
307
308 gsl_test (n != 2,
309 "gsl_poly_complex_solve_quadratic, two roots, (2x - 5)^2 = 4");
310 gsl_test_rel (GSL_REAL (z0), 1.5, 1e-9, "z0.real, (2x - 5)^2 = 4");
311 gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, (2x - 5)^2 = 4");
312 gsl_test_rel (GSL_REAL (z1), 3.5, 1e-9, "z1.real, (2x - 5)^2 = 4");
313 gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag, (2x - 5)^2 = 4");
314 }
315
316 {
317 gsl_complex z0, z1;
318
319 int n = gsl_poly_complex_solve_quadratic (4.0, 7.0, 0.0, &z0, &z1);
320
321 gsl_test (n != 2,
322 "gsl_poly_complex_solve_quadratic, two roots, x(4x + 7) = 0");
323 gsl_test_rel (GSL_REAL (z0), -1.75, 1e-9, "z0.real, x(4x + 7) = 0");
324 gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, x(4x + 7) = 0");
325 gsl_test_rel (GSL_REAL (z1), 0.0, 1e-9, "z1.real, x(4x + 7) = 0");
326 gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag, x(4x + 7) = 0");
327 }
328
329 {
330 gsl_complex z0, z1;
331
332 int n = gsl_poly_complex_solve_quadratic (5.0, 0.0, -20.0, &z0, &z1);
333
334 gsl_test (n != 2,
335 "gsl_poly_complex_solve_quadratic, two roots b = 0, 5 x^2 = 20");
336 gsl_test_rel (GSL_REAL (z0), -2.0, 1e-9, "z0.real, 5 x^2 = 20");
337 gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, 5 x^2 = 20");
338 gsl_test_rel (GSL_REAL (z1), 2.0, 1e-9, "z1.real, 5 x^2 = 20");
339 gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag, 5 x^2 = 20");
340 }
341
342 {
343 gsl_complex z0, z1;
344
345 int n = gsl_poly_complex_solve_quadratic (5.0, 0.0, 20.0, &z0, &z1);
346
347 gsl_test (n != 2,
348 "gsl_poly_complex_solve_quadratic, two roots b = 0, 5 x^2 = -20");
349 gsl_test_rel (GSL_REAL (z0), 0.0, 1e-9, "z0.real, 5 x^2 = -20");
350 gsl_test_rel (GSL_IMAG (z0), -2.0, 1e-9, "z0.imag, 5 x^2 = -20");
351 gsl_test_rel (GSL_REAL (z1), 0.0, 1e-9, "z1.real, 5 x^2 = -20");
352 gsl_test_rel (GSL_IMAG (z1), 2.0, 1e-9, "z1.imag, 5 x^2 = -20");
353 }
354
355
356 {
357 gsl_complex z0, z1;
358
359 int n = gsl_poly_complex_solve_quadratic (0.0, 3.0, -21.0, &z0, &z1);
360
361 gsl_test (n != 1,
362 "gsl_poly_complex_solve_quadratic, one root (linear) 3 x - 21 = 0");
363
364 gsl_test_rel (GSL_REAL (z0), 7.0, 1e-9, "z0.real, 3x - 21 = 0");
365 gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, 3x - 21 = 0");
366 }
367
368
369 {
370 gsl_complex z0, z1;
371
372 int n = gsl_poly_complex_solve_quadratic (0.0, 0.0, 1.0, &z0, &z1);
373 gsl_test (n != 0,
374 "gsl_poly_complex_solve_quadratic, no roots 1 = 0");
375 }
376
377
378
379 /* Cubic with complex roots */
380
381 {
382 gsl_complex z0, z1, z2;
383
384 int n = gsl_poly_complex_solve_cubic (0.0, 0.0, -27.0, &z0, &z1, &z2);
385
386 gsl_test (n != 3, "gsl_poly_complex_solve_cubic, three root, x^3 = 27");
387 gsl_test_rel (GSL_REAL (z0), -1.5, 1e-9, "z0.real, x^3 = 27");
388 gsl_test_rel (GSL_IMAG (z0), -1.5 * sqrt (3.0), 1e-9,
389 "z0.imag, x^3 = 27");
390 gsl_test_rel (GSL_REAL (z1), -1.5, 1e-9, "z1.real, x^3 = 27");
391 gsl_test_rel (GSL_IMAG (z1), 1.5 * sqrt (3.0), 1e-9, "z1.imag, x^3 = 27");
392 gsl_test_rel (GSL_REAL (z2), 3.0, 1e-9, "z2.real, x^3 = 27");
393 gsl_test_rel (GSL_IMAG (z2), 0.0, 1e-9, "z2.imag, x^3 = 27");
394 }
395
396 {
397 gsl_complex z0, z1, z2;
398
399 int n = gsl_poly_complex_solve_cubic (-1.0, 1.0, 39.0, &z0, &z1, &z2);
400
401 gsl_test (n != 3,
402 "gsl_poly_complex_solve_cubic, three root, (x+3)(x^2-4x+13) = 0");
403 gsl_test_rel (GSL_REAL (z0), -3.0, 1e-9, "z0.real, (x+3)(x^2+1) = 0");
404 gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, (x+3)(x^2+1) = 0");
405 gsl_test_rel (GSL_REAL (z1), 2.0, 1e-9, "z1.real, (x+3)(x^2+1) = 0");
406 gsl_test_rel (GSL_IMAG (z1), -3.0, 1e-9, "z1.imag, (x+3)(x^2+1) = 0");
407 gsl_test_rel (GSL_REAL (z2), 2.0, 1e-9, "z2.real, (x+3)(x^2+1) = 0");
408 gsl_test_rel (GSL_IMAG (z2), 3.0, 1e-9, "z2.imag, (x+3)(x^2+1) = 0");
409 }
410
411 {
412 gsl_complex z0, z1, z2;
413
414 int n =
415 gsl_poly_complex_solve_cubic (-51.0, 867.0, -4913.0, &z0, &z1, &z2);
416
417 gsl_test (n != 3,
418 "gsl_poly_complex_solve_cubic, three roots, (x-17)^3=0");
419 gsl_test_rel (GSL_REAL (z0), 17.0, 1e-9, "z0.real, (x-17)^3=0");
420 gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, (x-17)^3=0");
421 gsl_test_rel (GSL_REAL (z1), 17.0, 1e-9, "z1.real, (x-17)^3=0");
422 gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag, (x-17)^3=0");
423 gsl_test_rel (GSL_REAL (z2), 17.0, 1e-9, "z2.real, (x-17)^3=0");
424 gsl_test_rel (GSL_IMAG (z2), 0.0, 1e-9, "z2.imag, (x-17)^3=0");
425 }
426
427 {
428 gsl_complex z0, z1, z2;
429
430 int n =
431 gsl_poly_complex_solve_cubic (-57.0, 1071.0, -6647.0, &z0, &z1, &z2);
432
433 gsl_test (n != 3,
434 "gsl_poly_complex_solve_cubic, three roots, (x-17)(x-17)(x-23)=0");
435 gsl_test_rel (GSL_REAL (z0), 17.0, 1e-9, "z0.real, (x-17)(x-17)(x-23)=0");
436 gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, (x-17)(x-17)(x-23)=0");
437 gsl_test_rel (GSL_REAL (z1), 17.0, 1e-9, "z1.real, (x-17)(x-17)(x-23)=0");
438 gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag, (x-17)(x-17)(x-23)=0");
439 gsl_test_rel (GSL_REAL (z2), 23.0, 1e-9, "z2.real, (x-17)(x-17)(x-23)=0");
440 gsl_test_rel (GSL_IMAG (z2), 0.0, 1e-9, "z2.imag, (x-17)(x-17)(x-23)=0");
441 }
442
443 {
444 gsl_complex z0, z1, z2;
445
446 int n =
447 gsl_poly_complex_solve_cubic (-11.0, -493.0, +6647.0, &z0, &z1, &z2);
448
449 gsl_test (n != 3,
450 "gsl_poly_complex_solve_cubic, three roots, (x+23)(x-17)(x-17)=0");
451 gsl_test_rel (GSL_REAL (z0), -23.0, 1e-9,
452 "z0.real, (x+23)(x-17)(x-17)=0");
453 gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, (x+23)(x-17)(x-17)=0");
454 gsl_test_rel (GSL_REAL (z1), 17.0, 1e-9, "z1.real, (x+23)(x-17)(x-17)=0");
455 gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag, (x+23)(x-17)(x-17)=0");
456 gsl_test_rel (GSL_REAL (z2), 17.0, 1e-9, "z2.real, (x+23)(x-17)(x-17)=0");
457 gsl_test_rel (GSL_IMAG (z2), 0.0, 1e-9, "z2.imag, (x+23)(x-17)(x-17)=0");
458 }
459
460
461 {
462 gsl_complex z0, z1, z2;
463
464 int n =
465 gsl_poly_complex_solve_cubic (-143.0, 5087.0, -50065.0, &z0, &z1, &z2);
466
467 gsl_test (n != 3,
468 "gsl_poly_complex_solve_cubic, three roots, (x-17)(x-31)(x-95)=0");
469 gsl_test_rel (GSL_REAL (z0), 17.0, 1e-9, "z0.real, (x-17)(x-31)(x-95)=0");
470 gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, (x-17)(x-31)(x-95)=0");
471 gsl_test_rel (GSL_REAL (z1), 31.0, 1e-9, "z1.real, (x-17)(x-31)(x-95)=0");
472 gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag, (x-17)(x-31)(x-95)=0");
473 gsl_test_rel (GSL_REAL (z2), 95.0, 1e-9, "z2.real, (x-17)(x-31)(x-95)=0");
474 gsl_test_rel (GSL_IMAG (z2), 0.0, 1e-9, "z2.imag, (x-17)(x-31)(x-95)=0");
475 }
476
477
478 {
479 /* Wilkinson polynomial: y = (x-1)(x-2)(x-3)(x-4)(x-5) */
480
481 double a[6] = { -120, 274, -225, 85, -15, 1 };
482 double z[6*2];
483
484 gsl_poly_complex_workspace *w = gsl_poly_complex_workspace_alloc (6);
485
486 int status = gsl_poly_complex_solve (a, 6, w, z);
487
488 gsl_poly_complex_workspace_free (w);
489
490 gsl_test (status,
491 "gsl_poly_complex_solve, 5th-order Wilkinson polynomial");
492 gsl_test_rel (z[0], 1.0, 1e-9, "z0.real, 5th-order polynomial");
493 gsl_test_rel (z[1], 0.0, 1e-9, "z0.imag, 5th-order polynomial");
494 gsl_test_rel (z[2], 2.0, 1e-9, "z1.real, 5th-order polynomial");
495 gsl_test_rel (z[3], 0.0, 1e-9, "z1.imag, 5th-order polynomial");
496 gsl_test_rel (z[4], 3.0, 1e-9, "z2.real, 5th-order polynomial");
497 gsl_test_rel (z[5], 0.0, 1e-9, "z2.imag, 5th-order polynomial");
498 gsl_test_rel (z[6], 4.0, 1e-9, "z3.real, 5th-order polynomial");
499 gsl_test_rel (z[7], 0.0, 1e-9, "z3.imag, 5th-order polynomial");
500 gsl_test_rel (z[8], 5.0, 1e-9, "z4.real, 5th-order polynomial");
501 gsl_test_rel (z[9], 0.0, 1e-9, "z4.imag, 5th-order polynomial");
502 }
503
504 {
505 /* : 8-th order polynomial y = x^8 + x^4 + 1 */
506
507 double a[9] = { 1, 0, 0, 0, 1, 0, 0, 0, 1 };
508 double z[8*2];
509
510 double C = 0.5;
511 double S = sqrt (3.0) / 2.0;
512
513 gsl_poly_complex_workspace *w = gsl_poly_complex_workspace_alloc (9);
514
515 int status = gsl_poly_complex_solve (a, 9, w, z);
516
517 gsl_poly_complex_workspace_free (w);
518
519 gsl_test (status, "gsl_poly_complex_solve, 8th-order polynomial");
520
521 gsl_test_rel (z[0], -S, 1e-9, "z0.real, 8th-order polynomial");
522 gsl_test_rel (z[1], C, 1e-9, "z0.imag, 8th-order polynomial");
523 gsl_test_rel (z[2], -S, 1e-9, "z1.real, 8th-order polynomial");
524 gsl_test_rel (z[3], -C, 1e-9, "z1.imag, 8th-order polynomial");
525 gsl_test_rel (z[4], -C, 1e-9, "z2.real, 8th-order polynomial");
526 gsl_test_rel (z[5], S, 1e-9, "z2.imag, 8th-order polynomial");
527 gsl_test_rel (z[6], -C, 1e-9, "z3.real, 8th-order polynomial");
528 gsl_test_rel (z[7], -S, 1e-9, "z3.imag, 8th-order polynomial");
529 gsl_test_rel (z[8], C, 1e-9, "z4.real, 8th-order polynomial");
530 gsl_test_rel (z[9], S, 1e-9, "z4.imag, 8th-order polynomial");
531 gsl_test_rel (z[10], C, 1e-9, "z5.real, 8th-order polynomial");
532 gsl_test_rel (z[11], -S, 1e-9, "z5.imag, 8th-order polynomial");
533 gsl_test_rel (z[12], S, 1e-9, "z6.real, 8th-order polynomial");
534 gsl_test_rel (z[13], C, 1e-9, "z6.imag, 8th-order polynomial");
535 gsl_test_rel (z[14], S, 1e-9, "z7.real, 8th-order polynomial");
536 gsl_test_rel (z[15], -C, 1e-9, "z7.imag, 8th-order polynomial");
537
538 }
539
540 {
541 /* 15-th order polynomial y = (x + 1) * (x^10 + x^9 + 2 * x^7 + 4 * x^2 +
542 4 * x + 8) * (x - 1)^2 * (x - 2)^2
543
544 Problem reported by Munagala Ramanath (bug #39055)
545 */
546
547 double a[16] = { 32, -48, -8, 28, -8, 16, -16, 12,
548 -16, 6, 10, -17, 10, 2, -4, 1 };
549 double z[16*2];
550
551 double expected[16*2] = {
552 -1.6078107423472359, 0.00000000000000000,
553 -1.3066982484920768, 0.00000000000000000,
554 -1.0000000000000000, 0.00000000000000000,
555 -0.65893856175240950, -0.83459757287426684,
556 -0.65893856175240950, 0.83459757287426684,
557 -0.070891117403341281, -1.1359249087587791,
558 -0.070891117403341281, 1.1359249087587791,
559 0.57284747839410854, -1.1987808988289705,
560 0.57284747839410854, 1.1987808988289705,
561 1.0000000000000000, 0.00000000000000000,
562 1.0000000000000000, 0.00000000000000000,
563 1.1142366961812986, -0.48083981203389980,
564 1.1142366961812986, 0.48083981203389980,
565 2.0000000000000000, 0.00000000000000000,
566 2.0000000000000000, 0.00000000000000000 };
567
568 int i;
569
570 gsl_poly_complex_workspace *w = gsl_poly_complex_workspace_alloc (16);
571
572 int status = gsl_poly_complex_solve (a, 16, w, z);
573
574 gsl_poly_complex_workspace_free (w);
575
576 gsl_test (status, "gsl_poly_complex_solve, 15th-order polynomial");
577
578 gsl_heapsort(z, 15, 2 * sizeof(double), (gsl_comparison_fn_t) &cmp_cplx);
579
580 for (i = 0; i<15; i++)
581 {
582 gsl_test_rel (z[2*i], expected[2*i], 1e-7, "z%d.real, 15th-order polynomial", i);
583 gsl_test_rel (z[2*i+1], expected[2*i+1], 1e-7, "z%d.imag, 15th-order polynomial", i);
584 }
585 }
586
587
588 {
589 int i;
590
591 double xa[7] = {0.16, 0.97, 1.94, 2.74, 3.58, 3.73, 4.70 };
592 double ya[7] = {0.73, 1.11, 1.49, 1.84, 2.30, 2.41, 3.07 };
593
594 double dd_expected[7] = { 7.30000000000000e-01,
595 4.69135802469136e-01,
596 -4.34737219941284e-02,
597 2.68681098870099e-02,
598 -3.22937056934996e-03,
599 6.12763259971375e-03,
600 -6.45402453527083e-03 };
601
602 double dd[7], coeff[7], work[7];
603
604 gsl_poly_dd_init (dd, xa, ya, 7);
605
606 for (i = 0; i < 7; i++)
607 {
608 gsl_test_rel (dd[i], dd_expected[i], 1e-10, "divided difference dd[%d]", i);
609 }
610
611 for (i = 0; i < 7; i++)
612 {
613 double y = gsl_poly_dd_eval(dd, xa, 7, xa[i]);
614 gsl_test_rel (y, ya[i], 1e-10, "divided difference y[%d]", i);
615 }
616
617 gsl_poly_dd_taylor (coeff, 1.5, dd, xa, 7, work);
618
619 for (i = 0; i < 7; i++)
620 {
621 double y = gsl_poly_eval(coeff, 7, xa[i] - 1.5);
622 gsl_test_rel (y, ya[i], 1e-10, "taylor expansion about 1.5 y[%d]", i);
623 }
624 }
625
626 {
627 size_t i;
628 double xa[3] = { 1.3, 1.6, 1.9 };
629 double ya[3] = { 0.6200860, 0.4554022, 0.2818186 };
630 double dya[3] = { -0.5220232, -0.5698959, -0.5811571 };
631
632 double dd_expected[6] = { 6.200860000000e-01,
633 -5.220232000000e-01,
634 -8.974266666667e-02,
635 6.636555555556e-02,
636 2.666666666662e-03,
637 -2.774691357989e-03 };
638
639 double dd[6], za[6], coeff[6], work[6];
640
641 gsl_poly_dd_hermite_init(dd, za, xa, ya, dya, 3);
642
643 for (i = 0; i < 6; i++)
644 {
645 gsl_test_rel (dd[i], dd_expected[i], 1e-10, "hermite divided difference dd[%d]", i);
646 }
647
648 for (i = 0; i < 3; i++)
649 {
650 double y = gsl_poly_dd_eval(dd, za, 6, xa[i]);
651 gsl_test_rel (y, ya[i], 1e-10, "hermite divided difference y[%d]", i);
652 }
653
654 for (i = 0; i < 3; i++)
655 {
656 gsl_poly_dd_taylor(coeff, xa[i], dd, za, 6, work);
657 gsl_test_rel (coeff[1], dya[i], 1e-10, "hermite divided difference dy/dx[%d]", i);
658 }
659 }
660
661 {
662 double c[6] = { +1.0, -2.0, +3.0, -4.0, +5.0, -6.0 };
663 double dc[6];
664 double x;
665 x = -0.5;
666 gsl_poly_eval_derivs(c, 6, x, dc, 6);
667
668 gsl_test_rel (dc[0], c[0] + c[1]*x + c[2]*x*x + c[3]*x*x*x + c[4]*x*x*x*x + c[5]*x*x*x*x*x , eps, "gsl_poly_eval_derivs({+1, -2, +3, -4, +5, -6}, 3.75)");
669 gsl_test_rel (dc[1], c[1] + 2.0*c[2]*x + 3.0*c[3]*x*x + 4.0*c[4]*x*x*x + 5.0*c[5]*x*x*x*x , eps, "gsl_poly_eval_derivs({+1, -2, +3, -4, +5, -6} deriv 1, -12.375)");
670 gsl_test_rel (dc[2], 2.0*c[2] + 3.0*2.0*c[3]*x + 4.0*3.0*c[4]*x*x + 5.0*4.0*c[5]*x*x*x , eps, "gsl_poly_eval_derivs({+1, -2, +3, -4, +5, -6} deriv 2, +48.0)");
671 gsl_test_rel (dc[3], 3.0*2.0*c[3] + 4.0*3.0*2.0*c[4]*x + 5.0*4.0*3.0*c[5]*x*x , eps,"gsl_poly_eval_derivs({+1, -2, +3, -4, +5, -6} deriv 3, -174.0)");
672 gsl_test_rel (dc[4], 4.0*3.0*2.0*c[4] + 5.0*4.0*3.0*2.0*c[5]*x, eps, "gsl_poly_eval_derivs({+1, -2, +3, -4, +5, -6} deriv 4, +480.0)");
673 gsl_test_rel (dc[5], 5.0*4.0*3.0*2.0*c[5] , eps, "gsl_poly_eval_derivs({+1, -2, +3, -4, +5, -6} deriv 5, -720.0)");
674 }
675
676
677 /* now summarize the results */
678
679 exit (gsl_test_summary ());
680 }
681