Code-Eli  0.3.6
poly_root_test_suite.hpp
Go to the documentation of this file.
1 /*********************************************************************************
2 * Copyright (c) 2013 David D. Marshall <ddmarsha@calpoly.edu>
3 *
4 * All rights reserved. This program and the accompanying materials
5 * are made available under the terms of the Eclipse Public License v1.0
6 * which accompanies this distribution, and is available at
7 * http://www.eclipse.org/legal/epl-v10.html
8 *
9 * Contributors:
10 * David D. Marshall - initial code and implementation
11 ********************************************************************************/
12 
13 #ifndef poly_root_test_suite_hpp
14 #define poly_root_test_suite_hpp
15 
16 #include <typeinfo> // typeid
17 #include <string> // std::string
18 #include <sstream> // std::stringstream
19 #include <iomanip> // std::setw
20 #include <limits> // std::numeric_limits
21 
26 
27 template<typename data__>
28 class poly_root_test_suite : public Test::Suite
29 {
30  protected:
31  void AddTests(const float &)
32  {
41  }
42 
43  void AddTests(const double &)
44  {
53  }
54 
55  void AddTests(const long double &)
56  {
65  }
66 
67  public:
69  {
70  // add the tests
71  AddTests(data__());
72  }
74  {
75  }
76 
77  private:
79  {
80  // this is case where know the Sturm functions and the root count in several intervals
81  // the roots are between (-4, -3), (-3, -2), (-1, 0), (0, 1), (1, 2)
82  {
83  typename eli::mutil::poly::polynomial<data__>::coefficient_type coef(6), coef_out, coef_ref;
84 
85  // set coefficients
86  coef << 2, -10, -20, 0, 5, 1;
87 
89 
90  TEST_ASSERT(eli::mutil::poly::root::descartes_rule(f)==5);
91  TEST_ASSERT(eli::mutil::poly::root::descartes_rule(f, true)==2);
92  TEST_ASSERT(eli::mutil::poly::root::descartes_rule(f, false)==3);
93  }
94 
95  // this is case where know the Sturm functions and the real roots, -1 and 1. The
96  // other two roots are complex conjugates.
97  {
98  typename eli::mutil::poly::polynomial<data__>::coefficient_type coef(5), coef_out, coef_ref;
99 
100  // set coefficients
101  coef << -1, -1, 0, 1, 1;
102 
104 
105  TEST_ASSERT(eli::mutil::poly::root::descartes_rule(f)==4); // should be 2 but extra roots reported on neg. side
106  TEST_ASSERT(eli::mutil::poly::root::descartes_rule(f, true)==1);
107  TEST_ASSERT(eli::mutil::poly::root::descartes_rule(f, false)==3); // should be 1 but extra roots reported
108  }
109 
110  // this is case with a manually built polynomial with no repeated roots
111  {
112  data__ roots[7] = {-6, -4, -2, 2, 4, 6, 8};
113  eli::mutil::poly::polynomial<data__> f(roots, roots+7);
114 
115  TEST_ASSERT(eli::mutil::poly::root::descartes_rule(f)==7);
116  TEST_ASSERT(eli::mutil::poly::root::descartes_rule(f, true)==4);
117  TEST_ASSERT(eli::mutil::poly::root::descartes_rule(f, false)==3);
118  }
119 
120  // this is case with a manually built polynomial with repeated root -2 (twice) with
121  // -4, 2 and 5 as roots as well
122  {
123  data__ roots[5] = {-2, -2, -4, 2, 5};
124  eli::mutil::poly::polynomial<data__> f(roots, roots+5);
125 
126  TEST_ASSERT(eli::mutil::poly::root::descartes_rule(f)==5);
127  TEST_ASSERT(eli::mutil::poly::root::descartes_rule(f, true)==2);
128  TEST_ASSERT(eli::mutil::poly::root::descartes_rule(f, false)==3);
129  }
130 
131  // this is case with a manually built polynomial with repeated root -2 (thrice) with
132  // -4, 2 and 5 as roots as well
133  {
134  data__ roots[6] = {-2, -2, -2, -4, 2, 5};
135  eli::mutil::poly::polynomial<data__> f(roots, roots+6);
136 
137  TEST_ASSERT(eli::mutil::poly::root::descartes_rule(f)==6);
138  TEST_ASSERT(eli::mutil::poly::root::descartes_rule(f, true)==2);
139  TEST_ASSERT(eli::mutil::poly::root::descartes_rule(f, false)==4);
140  }
141 
142  // this is case with a manually built polynomial with repeated roots -2 (thrice) and
143  // 2 (twice) with -4 and 5 as roots as well
144  {
145  data__ roots[7] = {-2, -2, -2, 2, 2, -4, 5};
146  eli::mutil::poly::polynomial<data__> f(roots, roots+7);
147 
148  TEST_ASSERT(eli::mutil::poly::root::descartes_rule(f)==7);
149  TEST_ASSERT(eli::mutil::poly::root::descartes_rule(f, true)==3);
150  TEST_ASSERT(eli::mutil::poly::root::descartes_rule(f, false)==4);
151  }
152 
153  // this is case with a manually built polynomial with repeated roots -2 (twice) and
154  // 5 (thrice) with -4 and 2 as roots as well
155  {
156  data__ roots[7] = {-2, -2, 2, 2, 5, -4, 2};
157  eli::mutil::poly::polynomial<data__> f(roots, roots+7);
158 
159  TEST_ASSERT(eli::mutil::poly::root::descartes_rule(f)==7);
160  TEST_ASSERT(eli::mutil::poly::root::descartes_rule(f, true)==4);
161  TEST_ASSERT(eli::mutil::poly::root::descartes_rule(f, false)==3);
162  }
163  }
164 
165  void sturm_test()
166  {
167  // this is case where know the Sturm functions and the root count in several intervals
168  // the roots are between (-4, -3), (-3, -2), (-1, 0), (0, 1), (1, 2)
169  {
170  typename eli::mutil::poly::polynomial<data__>::coefficient_type coef(6), coef_out, coef_ref;
171 
172  // set coefficients
173  coef << 2, -10, -20, 0, 5, 1;
174 
176  std::vector< eli::mutil::poly::polynomial<data__> > p;
177 
178  // get the Sturm functions and manipulate to match reference values
180  p[3].multiply(3);
181  p[4].multiply(17);
182 
183  // test p[0]
184  coef_ref.resize(6);
185  coef_ref << 2, -10, -20, 0, 5, 1;
186  p[0].get_coefficients(coef_out);
187  TEST_ASSERT(coef_out==coef_ref);
188 
189  // test p[1]
190  coef_ref.resize(5);
191  coef_ref << -2, -8, 0, 4, 1;
192  p[1].get_coefficients(coef_out);
193  TEST_ASSERT(coef_out==coef_ref);
194 
195  // test p[2]
196  coef_ref.resize(4);
197  coef_ref << -1, 0, 3, 1;
198  p[2].get_coefficients(coef_out);
199  TEST_ASSERT(coef_out==coef_ref);
200 
201  // test p[3]
202  coef_ref.resize(3);
203  coef_ref << 1, 7, 3;
204  p[3].get_coefficients(coef_out);
205  if ((typeid(data__)==typeid(float)) || (typeid(data__)==typeid(double)))
206  {
207  TEST_ASSERT((coef_out-coef_ref).norm()<5*std::numeric_limits<data__>::epsilon());
208  }
209  else
210  {
211  TEST_ASSERT(coef_out==coef_ref);
212  }
213 
214  // test p[4]
215  coef_ref.resize(2);
216  coef_ref << 11, 17;
217  p[4].get_coefficients(coef_out);
218  TEST_ASSERT((coef_out-coef_ref).norm()<4986*std::numeric_limits<data__>::epsilon());
219 
220  // test p[5]
221  TEST_ASSERT(p[5].degree()==0);
222  TEST_ASSERT(p[5].f(0)>0);
223 
224  // calculate the root count and get back the Sturm functions
225  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, -5, 5)==5); // total real roots
226  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, -5, 0)==3); // negative real roots
227  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, 0, 5)==2); // positive real roots
228  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, -4, -3)==1); // root bracket
229  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, -3, -2)==1); // root bracket
230  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, -1, 0)==1); // root bracket
231  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, 0, 1)==1); // root bracket
232  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, 1, 2)==1); // root bracket
233  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, -2, -1)==0); // no root
234 
235  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p)==5); // total number of real roots
236  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, false)==3); // total number of negative real roots
237  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, true)==2); // total number of positive real roots
238  }
239 
240  // this is case where know the Sturm functions and the real roots, -1 and 1. The
241  // other two roots are complex conjugates.
242  {
243  typename eli::mutil::poly::polynomial<data__>::coefficient_type coef(5), coef_out, coef_ref;
244 
245  // set coefficients
246  coef << -1, -1, 0, 1, 1;
247 
249  std::vector< eli::mutil::poly::polynomial<data__> > p;
250 
251  // get the Sturm functions and manipulate to match reference values
253  p[1].multiply(4);
254 
255  // test p[0]
256  coef_ref.resize(5);
257  coef_ref << -1, -1, 0, 1, 1;
258  p[0].get_coefficients(coef_out);
259  TEST_ASSERT(coef_out==coef_ref);
260 
261  // test p[1]
262  coef_ref.resize(4);
263  coef_ref << -1, 0, 3, 4;
264  p[1].get_coefficients(coef_out);
265  TEST_ASSERT(coef_out==coef_ref);
266 
267  // test p[2]
268  coef_ref.resize(3);
269  coef_ref << 5, 4, 1;
270  p[2].get_coefficients(coef_out);
271  TEST_ASSERT(coef_out==coef_ref);
272 
273  // test p[3]
274  coef_ref.resize(2);
275  coef_ref << -2, -1;
276  p[3].get_coefficients(coef_out);
277  TEST_ASSERT(coef_out==coef_ref);
278 
279  // test p[4]
280  coef_ref.resize(1);
281  coef_ref << -1;
282  p[4].get_coefficients(coef_out);
283  TEST_ASSERT(coef_out==coef_ref);
284 
285  // calculate the root count and get back the Sturm functions
286  TEST_ASSERT(eli::mutil::poly::root::sturm_count(f, -5, 5)==2); // total real roots
287  TEST_ASSERT(eli::mutil::poly::root::sturm_count(f, -5, 0)==1); // negative real roots
288  TEST_ASSERT(eli::mutil::poly::root::sturm_count(f, 0, 5)==1); // positive real roots
289  TEST_ASSERT(eli::mutil::poly::root::sturm_count(f, -2, 0)==1); // root bracket
290  TEST_ASSERT(eli::mutil::poly::root::sturm_count(f, 0, 2)==1); // root bracket
291  TEST_ASSERT(eli::mutil::poly::root::sturm_count(f, -3, -2)==0); // no root
292 
293  TEST_ASSERT(eli::mutil::poly::root::sturm_count(f)==2); // total number of real roots
294  TEST_ASSERT(eli::mutil::poly::root::sturm_count(f, false)==1); // total number of negative real roots
295  TEST_ASSERT(eli::mutil::poly::root::sturm_count(f, true)==1); // total number of positive real roots
296  }
297 
298  // this is case with a manually built polynomial with no repeated roots
299  {
300  data__ roots[7] = {-6, -4, -2, 2, 4, 6, 8};
301  eli::mutil::poly::polynomial<data__> f(roots, roots+7);
302  std::vector< eli::mutil::poly::polynomial<data__> > p;
303 
304  // get the Sturm functions and manipulate to match reference values
306 
307  // calculate the root count and get back the Sturm functions
308  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, -9, 9)==7); // total real roots
309  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, -9, 0)==3); // negative real roots
310  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, 0, 9)==4); // positive real roots
311  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, -7, -5)==1); // root bracket
312  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, -5, -3)==1); // root bracket
313  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, -3, -1)==1); // root bracket
314  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, 1, 3)==1); // root bracket
315  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, 3, 5)==1); // root bracket
316  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, 5, 7)==1); // root bracket
317  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, 7, 9)==1); // root bracket
318  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, -8, -7)==0); // no root
319 
320  // switched to passing the polynomial for no performance reason, just testing variety
321  TEST_ASSERT(eli::mutil::poly::root::sturm_count(f)==7); // total number of real roots
322  TEST_ASSERT(eli::mutil::poly::root::sturm_count(f, false)==3); // total number of negative real roots
323  TEST_ASSERT(eli::mutil::poly::root::sturm_count(f, true)==4); // total number of positive real roots
324  }
325 
326  // this is case with a manually built polynomial with repeated root -2 (twice) with
327  // -4, 2 and 5 as roots as well
328  {
329  data__ roots[5] = {-2, -2, -4, 2, 5};
330  eli::mutil::poly::polynomial<data__> f(roots, roots+5);
331  std::vector< eli::mutil::poly::polynomial<data__> > p;
332 
333  // get the Sturm functions and manipulate to match reference values
335 
336  // calculate the root count and get back the Sturm functions
337  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, -9, 9)==4); // total real roots
338  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, -9, 0)==2); // negative real roots
339  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, 0, 9)==2); // positive real roots
340  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, -5, -3)==1); // root bracket
341  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, -3, -1)==1); // double root bracket
342  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, 1, 3)==1); // root bracket
343  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, 4, 6)==1); // root bracket
344  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, 3, 4)==0); // no root
345 
346  // switched to passing the polynomial for no performance reason, just testing variety
347  TEST_ASSERT(eli::mutil::poly::root::sturm_count(f)==4); // total number of real roots
348  TEST_ASSERT(eli::mutil::poly::root::sturm_count(f, false)==2); // total number of negative real roots
349  TEST_ASSERT(eli::mutil::poly::root::sturm_count(f, true)==2); // total number of positive real roots
350  }
351 
352  // this is case with a manually built polynomial with repeated root -2 (thrice) with
353  // -4, 2 and 5 as roots as well
354  {
355  data__ roots[6] = {-2, -2, -2, -4, 2, 5};
356  eli::mutil::poly::polynomial<data__> f(roots, roots+6);
357  std::vector< eli::mutil::poly::polynomial<data__> > p;
358 
359  // get the Sturm functions and manipulate to match reference values
361 
362  // calculate the root count and get back the Sturm functions
363  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, -9, 9)==4); // total real roots
364  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, -9, 0)==2); // negative real roots
365  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, 0, 9)==2); // positive real roots
366  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, -5, -3)==1); // root bracket
367  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, -3, -1)==1); // double root bracket
368  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, 1, 3)==1); // root bracket
369  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, 4, 6)==1); // root bracket
370  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, 3, 4)==0); // no root
371 
372  // switched to passing the polynomial for no performance reason, just testing variety
373  TEST_ASSERT(eli::mutil::poly::root::sturm_count(f)==4); // total number of real roots
374  TEST_ASSERT(eli::mutil::poly::root::sturm_count(f, false)==2); // total number of negative real roots
375  TEST_ASSERT(eli::mutil::poly::root::sturm_count(f, true)==2); // total number of positive real roots
376  }
377 
378  // this is case with a manually built polynomial with repeated roots -2 (thrice) and
379  // 2 (twice) with -4 and 5 as roots as well
380  {
381  data__ roots[7] = {-2, -2, -2, 2, 2, -4, 5};
382  eli::mutil::poly::polynomial<data__> f(roots, roots+7);
383  std::vector< eli::mutil::poly::polynomial<data__> > p;
384 
385  // get the Sturm functions and manipulate to match reference values
387 
388  // calculate the root count and get back the Sturm functions
389  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, -9, 9)==4); // total real roots
390  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, -9, 0)==2); // negative real roots
391  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, 0, 9)==2); // positive real roots
392  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, -5, -3)==1); // root bracket
393  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, -3, -1)==1); // triple root bracket
394  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, 1, 3)==1); // double root bracket
395  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, 4, 6)==1); // root bracket
396  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, 3, 4)==0); // no root
397 
398  // switched to passing the polynomial for no performance reason, just testing variety
399  TEST_ASSERT(eli::mutil::poly::root::sturm_count(f)==4); // total number of real roots
400  TEST_ASSERT(eli::mutil::poly::root::sturm_count(f, false)==2); // total number of negative real roots
401  TEST_ASSERT(eli::mutil::poly::root::sturm_count(f, true)==2); // total number of positive real roots
402  }
403 
404  // this is case with a manually built polynomial with repeated roots -2 (twice) and
405  // 5 (thrice) with -4 and 2 as roots as well
406  {
407  data__ roots[7] = {-2, -2, 2, 2, 5, -4, 2};
408  eli::mutil::poly::polynomial<data__> f(roots, roots+7);
409  std::vector< eli::mutil::poly::polynomial<data__> > p;
410 
411  // get the Sturm functions and manipulate to match reference values
413 
414  // calculate the root count and get back the Sturm functions
415  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, -9, 9)==4); // total real roots
416  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, -9, 0)==2); // negative real roots
417  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, 0, 9)==2); // positive real roots
418  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, -5, -3)==1); // root bracket
419  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, -3, -1)==1); // double root bracket
420  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, 1, 3)==1); // root bracket
421  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, 4, 6)==1); // triple root bracket
422  TEST_ASSERT(eli::mutil::poly::root::sturm_count(p, 3, 4)==0); // no root
423 
424  // switched to passing the polynomial for no performance reason, just testing variety
425  TEST_ASSERT(eli::mutil::poly::root::sturm_count(f)==4); // total number of real roots
426  TEST_ASSERT(eli::mutil::poly::root::sturm_count(f, false)==2); // total number of negative real roots
427  TEST_ASSERT(eli::mutil::poly::root::sturm_count(f, true)==2); // total number of positive real roots
428  }
429  }
430 
431  // test 0th order polynomials (non-zero and zero)
433  {
435 
437  int nroot;
438  std::vector<data__> root;
439 
440  // zero case
441  f=0;
442  f.get_coefficients(a);
443  nroot=eli::mutil::poly::root::closed_form(root, a.data(), 0);
444  TEST_ASSERT(nroot==std::numeric_limits<int>::max());
445  TEST_ASSERT(root.size()==0);
446 
447  // non-zero case
448  f=1;
449  f.get_coefficients(a);
450  nroot=eli::mutil::poly::root::closed_form(root, a.data(), 0);
451  TEST_ASSERT(nroot==0);
452  TEST_ASSERT(root.size()==0);
453  }
454 
455  // test 1st order polynomials
457  {
459 
461  data__ root_val;
462  int nroot, deg(1);
463  std::vector<data__> root;
464 
465  // degenerate polynomial
466  a.resize(deg+1);
467  a.setZero();
468  a[deg-1]=1;
469  f.set_coefficients(a);
470  f.get_coefficients(a);
471  nroot=eli::mutil::poly::root::closed_form(root, a.data(), deg);
472  TEST_ASSERT(nroot==0);
473  TEST_ASSERT(root.size()==0);
474 
475  // root at zero case
476  root_val=static_cast<data__>(0);
477  f.set_roots(root_val);
478  f.get_coefficients(a);
479  nroot=eli::mutil::poly::root::closed_form(root, a.data(), deg);
480  TEST_ASSERT(nroot==1);
481  TEST_ASSERT(root.size()==1);
482  if (root.size()==1)
483  {
484  TEST_ASSERT(root[0]==root_val);
485  }
486  root.resize(0);
487 
488  // root away from zero case
489  root_val=2;
490  f.set_roots(root_val);
491  f.get_coefficients(a);
492  nroot=eli::mutil::poly::root::closed_form(root, a.data(), deg);
493  TEST_ASSERT(nroot==1);
494  TEST_ASSERT(root.size()==1);
495  if (root.size()==1)
496  {
497  TEST_ASSERT(root[0]==root_val);
498  }
499  }
500 
501  // test 2nd order polynomials
502  // (no real roots, one real double root, two real distinct roots)
504  {
506 
508  int nroot, deg(2);
509  data__ root_val[2];
510  std::vector<data__> root;
511 
512  // degenerate polynomial
513  a.resize(deg+1);
514  a.setZero();
515  a[deg-1]=1;
516  f.set_coefficients(a);
517  f.get_coefficients(a);
518  nroot=eli::mutil::poly::root::closed_form(root, a.data(), deg);
519  TEST_ASSERT(nroot==1);
520  TEST_ASSERT(root.size()==1);
521  if (root.size()==1)
522  {
523  TEST_ASSERT(root[0]==0);
524  }
525  root.resize(0);
526 
527  // two real roots case
528  root_val[0]=-1;
529  root_val[1]=2;
530  f.set_roots(root_val, root_val+deg);
531  f.get_coefficients(a);
532  nroot=eli::mutil::poly::root::closed_form(root, a.data(), deg);
533  TEST_ASSERT(nroot==2);
534  TEST_ASSERT(root.size()==2);
535  if (root.size()==2)
536  {
537  TEST_ASSERT(root[0]==root_val[0]);
538  TEST_ASSERT(root[1]==root_val[1]);
539  }
540  root.resize(0);
541 
542  // double root case
543  root_val[0]=3;
544  root_val[1]=3;
545  f.set_roots(root_val, root_val+deg);
546  f.get_coefficients(a);
547  nroot=eli::mutil::poly::root::closed_form(root, a.data(), deg);
548  TEST_ASSERT(nroot==2);
549  TEST_ASSERT(root.size()==2);
550  if (root.size()==2)
551  {
552  TEST_ASSERT(root[0]==root_val[0]);
553  TEST_ASSERT(root[1]==root_val[1]);
554  }
555  root.resize(0);
556 
557  // no root case
558  a.resize(deg+1);
559  a << 1, 1, 1;
560  f.set_coefficients(a);
561  f.get_coefficients(a);
562  nroot=eli::mutil::poly::root::closed_form(root, a.data(), deg);
563  TEST_ASSERT(nroot==0);
564  TEST_ASSERT(root.size()==0);
565  }
566 
567  // test 3rd order polynomials
568  // (one real triple root, one unique real & one real double root, three real unique roots, one real root)
570  {
572 
574  int nroot, deg(3);
575  data__ root_val[3];
576  std::vector<data__> root;
577 
578  // degenerate polynomial
579  a.resize(deg+1);
580  a.setZero();
581  a[deg-1]=1;
582  f.set_coefficients(a);
583  f.get_coefficients(a);
584  nroot=eli::mutil::poly::root::closed_form(root, a.data(), deg);
585  TEST_ASSERT(nroot==2);
586  TEST_ASSERT(root.size()==2);
587  if (root.size()==2)
588  {
589  TEST_ASSERT(root[0]==0);
590  TEST_ASSERT(root[1]==0);
591  }
592  root.resize(0);
593 
594  // caused problems
595  a.resize(deg+1);
596  a << -54, -27, 0, 1;
597  root_val[0]=-3;
598  root_val[1]=-3;
599  root_val[2]=6;
600  f.set_coefficients(a);
601  f.get_coefficients(a);
602  nroot=eli::mutil::poly::root::closed_form(root, a.data(), deg);
603  TEST_ASSERT(nroot==3);
604  TEST_ASSERT(root.size()==3);
605  if (root.size()==3)
606  {
607  TEST_ASSERT(root[0]==root_val[0]);
608  TEST_ASSERT(root[1]==root_val[1]);
609  TEST_ASSERT(root[2]==root_val[2]);
610  }
611  root.resize(0);
612 
613  // triple root case
614  root_val[0]=2;
615  root_val[1]=2;
616  root_val[2]=2;
617  f.set_roots(root_val, root_val+deg);
618  f.get_coefficients(a);
619  nroot=eli::mutil::poly::root::closed_form(root, a.data(), deg);
620  TEST_ASSERT(nroot==3);
621  TEST_ASSERT(root.size()==3);
622  if (root.size()==3)
623  {
624  TEST_ASSERT(root[0]==root_val[0]);
625  TEST_ASSERT(root[1]==root_val[1]);
626  TEST_ASSERT(root[2]==root_val[2]);
627  }
628  root.resize(0);
629 
630  // one unique & double root case
631  root_val[0]=-2;
632  root_val[1]=1;
633  root_val[2]=1;
634  f.set_roots(root_val, root_val+deg);
635  f.get_coefficients(a);
636  nroot=eli::mutil::poly::root::closed_form(root, a.data(), deg);
637  TEST_ASSERT(nroot==3);
638  TEST_ASSERT(root.size()==3);
639  if (root.size()==3)
640  {
641  TEST_ASSERT(root[0]==root_val[0]);
642  TEST_ASSERT(root[1]==root_val[1]);
643  TEST_ASSERT(root[2]==root_val[2]);
644  }
645  root.resize(0);
646 
647  // three unique root case
648  root_val[0]=-1;
649  root_val[1]=2;
650  root_val[2]=1;
651  f.set_roots(root_val, root_val+deg);
652  f.get_coefficients(a);
653  nroot=eli::mutil::poly::root::closed_form(root, a.data(), deg);
654  TEST_ASSERT(nroot==3);
655  TEST_ASSERT(root.size()==3);
656  if (root.size()==3)
657  {
658  TEST_ASSERT_DELTA(root[0], root_val[0], 2*std::numeric_limits<data__>::epsilon());
659  TEST_ASSERT_DELTA(root[1], root_val[2], 2*std::numeric_limits<data__>::epsilon());
660  TEST_ASSERT_DELTA(root[2], root_val[1], 2*std::numeric_limits<data__>::epsilon());
661  }
662  root.resize(0);
663 
664  // one real root case
665  a.resize(deg+1);
666  a << 1, 2, 2, 1;
667  root_val[0]=-1;
668  f.set_coefficients(a);
669  f.get_coefficients(a);
670  nroot=eli::mutil::poly::root::closed_form(root, a.data(), deg);
671  TEST_ASSERT(nroot==1);
672  TEST_ASSERT(root.size()==1);
673  if (root.size()==1)
674  {
675  if (typeid(data__)==typeid(double))
676  {
677  TEST_ASSERT_DELTA(root[0], root_val[0], 2*std::numeric_limits<data__>::epsilon());
678  }
679  else if (typeid(data__)==typeid(long double))
680  {
681  TEST_ASSERT_DELTA(root[0], root_val[0], 513*std::numeric_limits<data__>::epsilon());
682  }
683  else
684  {
685  TEST_ASSERT_DELTA(root[0], root_val[0], std::numeric_limits<data__>::epsilon());
686  }
687  }
688  }
689 
690  // test 4th order polynomials
691  // (no real roots, one real quadruple root, one unique real & one triple root, two real double roots, two unique real & one double root, one real double root, two unique real roots)
693  {
695 
697  int nroot, deg(4);
698  data__ root_val[4];
699  std::vector<data__> root;
700 
701  // degenerate polynomial
702  a.resize(deg+1);
703  a.setZero();
704  a[deg-1]=1;
705  f.set_coefficients(a);
706  f.get_coefficients(a);
707  nroot=eli::mutil::poly::root::closed_form(root, a.data(), deg);
708  TEST_ASSERT(nroot==3);
709  TEST_ASSERT(root.size()==3);
710  TEST_ASSERT(root[0]==0);
711  TEST_ASSERT(root[1]==0);
712  TEST_ASSERT(root[2]==0);
713  root.resize(0);
714 
715  // 4 real root case sample calculation
716  a.resize(deg+1);
717  a << 4, 0, -5, 0, 1;
718  root_val[0]=-2;
719  root_val[1]=-1;
720  root_val[2]=1;
721  root_val[3]=2;
722  f.set_coefficients(a);
723  f.get_coefficients(a);
724  nroot=eli::mutil::poly::root::closed_form(root, a.data(), deg);
725  TEST_ASSERT(nroot==4);
726  TEST_ASSERT(root.size()==4);
727  if (root.size()==4)
728  {
729  TEST_ASSERT_DELTA(root[0], root_val[0], 3*std::numeric_limits<data__>::epsilon());
730  TEST_ASSERT_DELTA(root[1], root_val[1], 3*std::numeric_limits<data__>::epsilon());
731  TEST_ASSERT_DELTA(root[2], root_val[2], 3*std::numeric_limits<data__>::epsilon());
732  TEST_ASSERT_DELTA(root[3], root_val[3], 3*std::numeric_limits<data__>::epsilon());
733  }
734  root.resize(0);
735 
736  // 2 real root case sample calculation
737  a.resize(deg+1);
738  a << -1, 2, 0, 2, 1;
739  {
740  data__ one(1), two(2);
741  root_val[0]=-one-std::sqrt(two);
742  root_val[1]=-one+std::sqrt(two);
743  }
744  f.set_coefficients(a);
745  f.get_coefficients(a);
746  nroot=eli::mutil::poly::root::closed_form(root, a.data(), deg);
747  TEST_ASSERT(nroot==2);
748  TEST_ASSERT(root.size()==2);
749  if (root.size()==2)
750  {
751  if (typeid(data__)==typeid(long double))
752  {
753  TEST_ASSERT_DELTA(root[0], root_val[0], 17*std::numeric_limits<data__>::epsilon());
754  TEST_ASSERT_DELTA(root[1], root_val[1], 17*std::numeric_limits<data__>::epsilon());
755  }
756  else
757  {
758  TEST_ASSERT(root[0]==root_val[0]);
759  TEST_ASSERT(root[1]==root_val[1]);
760  }
761  }
762  root.resize(0);
763 
764  // no real root case sample calculation
765  a.resize(deg+1);
766  a << 4, 0, 5, 0, 1;
767  f.set_coefficients(a);
768  f.get_coefficients(a);
769  nroot=eli::mutil::poly::root::closed_form(root, a.data(), deg);
770  TEST_ASSERT(nroot==0);
771  TEST_ASSERT(root.size()==0);
772 
773  // quadruple root case
774  root_val[0]=2;
775  root_val[1]=2;
776  root_val[2]=2;
777  root_val[3]=2;
778  f.set_roots(root_val, root_val+deg);
779  f.get_coefficients(a);
780  nroot=eli::mutil::poly::root::closed_form(root, a.data(), deg);
781  TEST_ASSERT(nroot==4);
782  TEST_ASSERT(root.size()==4);
783  if (root.size()==4)
784  {
785  TEST_ASSERT(root[0]==root_val[0]);
786  TEST_ASSERT(root[1]==root_val[1]);
787  TEST_ASSERT(root[2]==root_val[2]);
788  TEST_ASSERT(root[3]==root_val[3]);
789  }
790  root.resize(0);
791 
792  // one unique & triple root case
793  root_val[0]=-2;
794  root_val[1]=1;
795  root_val[2]=1;
796  root_val[3]=1;
797  f.set_roots(root_val, root_val+deg);
798  f.get_coefficients(a);
799  nroot=eli::mutil::poly::root::closed_form(root, a.data(), deg);
800  TEST_ASSERT(nroot==4);
801  TEST_ASSERT(root.size()==4);
802  if (root.size()==4)
803  {
804  TEST_ASSERT(root[0]==root_val[0]);
805  TEST_ASSERT(root[1]==root_val[1]);
806  TEST_ASSERT(root[2]==root_val[2]);
807  TEST_ASSERT(root[3]==root_val[3]);
808  }
809  root.resize(0);
810 
811  // two real double root case
812  root_val[0]=-2;
813  root_val[1]=-2;
814  root_val[2]=1;
815  root_val[3]=1;
816  f.set_roots(root_val, root_val+deg);
817  f.get_coefficients(a);
818  nroot=eli::mutil::poly::root::closed_form(root, a.data(), deg);
819  TEST_ASSERT(nroot==4);
820  TEST_ASSERT(root.size()==4);
821  if (root.size()==4)
822  {
823  TEST_ASSERT(root[0]==root_val[0]);
824  TEST_ASSERT(root[1]==root_val[1]);
825  TEST_ASSERT(root[2]==root_val[2]);
826  TEST_ASSERT(root[3]==root_val[3]);
827  }
828  root.resize(0);
829 
830  // two unique real & real double root case
831  root_val[0]=-2;
832  root_val[1]=-2;
833  root_val[2]=1;
834  root_val[3]=4;
835  f.set_roots(root_val, root_val+deg);
836  f.get_coefficients(a);
837  nroot=eli::mutil::poly::root::closed_form(root, a.data(), deg);
838  TEST_ASSERT(nroot==4);
839  TEST_ASSERT(root.size()==4);
840  if (root.size()==4)
841  {
842  if ((typeid(data__)==typeid(float)) || (typeid(data__)==typeid(double)) || (typeid(data__)==typeid(long double)))
843  {
844  TEST_ASSERT(root[0]==root_val[0]);
845  TEST_ASSERT(root[1]==root_val[1]);
846  TEST_ASSERT_DELTA(root[2], root_val[2], std::numeric_limits<data__>::epsilon());
847  TEST_ASSERT(root[3]==root_val[3]);
848  }
849  else
850  {
851  TEST_ASSERT(root[0]==root_val[0]);
852  TEST_ASSERT(root[1]==root_val[1]);
853  TEST_ASSERT(root[2]==root_val[2]);
854  TEST_ASSERT(root[3]==root_val[3]);
855  }
856  }
857  root.resize(0);
858 
859  // four unique root case
860  root_val[0]=-1;
861  root_val[1]=2;
862  root_val[2]=1;
863  root_val[3]=4;
864  f.set_roots(root_val, root_val+deg);
865  f.get_coefficients(a);
866  nroot=eli::mutil::poly::root::closed_form(root, a.data(), deg);
867  TEST_ASSERT(nroot==4);
868  TEST_ASSERT(root.size()==4);
869  if (root.size()==4)
870  {
871  if ((typeid(data__)==typeid(float)) || (typeid(data__)==typeid(long double)))
872  {
873  TEST_ASSERT_DELTA(root[0], root_val[0], std::numeric_limits<data__>::epsilon());
874  TEST_ASSERT_DELTA(root[1], root_val[2], std::numeric_limits<data__>::epsilon());
875  TEST_ASSERT(root[2]==root_val[1]);
876  TEST_ASSERT(root[3]==root_val[3]);
877  }
878  else
879  {
880  TEST_ASSERT(root[0]==root_val[0]);
881  TEST_ASSERT(root[1]==root_val[2]);
882  TEST_ASSERT(root[2]==root_val[1]);
883  TEST_ASSERT(root[3]==root_val[3]);
884  }
885  }
886  root.resize(0);
887 
888  // one real double root
889  root_val[0]=1;
890  root_val[1]=1;
891  a.resize(deg+1);
892  a << root_val[0]*root_val[1],
893  root_val[0]*root_val[1]-root_val[0]-root_val[1],
894  root_val[0]*root_val[1]-root_val[0]-root_val[1]+1,
895  1-root_val[0]-root_val[1],
896  1;
897  f.set_coefficients(a);
898  f.get_coefficients(a);
899  nroot=eli::mutil::poly::root::closed_form(root, a.data(), deg);
900  TEST_ASSERT(nroot==2);
901  TEST_ASSERT(root.size()==2);
902  if (root.size()==2)
903  {
904  TEST_ASSERT(root[0]==root_val[0]);
905  TEST_ASSERT(root[1]==root_val[1]);
906  }
907  root.resize(0);
908 
909  // two real root
910  root_val[0]=2;
911  root_val[1]=1;
912  a.resize(deg+1);
913  a << root_val[0]*root_val[1],
914  root_val[0]*root_val[1]-root_val[0]-root_val[1],
915  root_val[0]*root_val[1]-root_val[0]-root_val[1]+1,
916  1-root_val[0]-root_val[1],
917  1;
918  f.set_coefficients(a);
919  f.get_coefficients(a);
920  nroot=eli::mutil::poly::root::closed_form(root, a.data(), deg);
921  TEST_ASSERT(nroot==2);
922  TEST_ASSERT(root.size()==2);
923  if (root.size()==2)
924  {
925  TEST_ASSERT(root[0]==root_val[1]);
926  TEST_ASSERT(root[1]==root_val[0]);
927  }
928  root.resize(0);
929 
930  // no real root
931  a.resize(deg+1);
932  a << 2, 4, 5, 3, 1;
933  f.set_coefficients(a);
934  f.get_coefficients(a);
935  nroot=eli::mutil::poly::root::closed_form(root, a.data(), deg);
936  TEST_ASSERT(nroot==0);
937  TEST_ASSERT(root.size()==0);
938  }
939 
941  {
942  // -1, 1, -1, 1, -1, 1 has one real root which is 1
943  }
944 };
945 
946 #endif
void AddTests(const double &)
Definition: poly_root_test_suite.hpp:43
~poly_root_test_suite()
Definition: poly_root_test_suite.hpp:73
poly_root_test_suite()
Definition: poly_root_test_suite.hpp:68
void closed_form_root_test_3()
Definition: poly_root_test_suite.hpp:569
Definition: poly_root_test_suite.hpp:28
void set_coefficients(const coefficient_type &ain)
Definition: polynomial.hpp:81
int descartes_rule(const polynomial< data__ > &f, bool positive)
Definition: descartes_rule.hpp:32
Definition: polynomial.hpp:31
Eigen::Matrix< data_type, Eigen::Dynamic, 1 > coefficient_type
Definition: polynomial.hpp:35
void set_roots(const data_type &root)
Definition: polynomial.hpp:118
int sturm_count(const std::vector< polynomial< data__ > > &sturm_fun, const data2__ &xmin, const data2__ &xmax)
Definition: sturm_count.hpp:82
void sturm_test()
Definition: poly_root_test_suite.hpp:165
void sturm_functions(std::vector< polynomial< data__ > > &sturm_fun, const polynomial< data__ > &f)
Definition: sturm_count.hpp:34
void get_coefficients(coefficient_type &aout) const
Definition: polynomial.hpp:86
int closed_form(std::vector< data__ > &root, const data__ a[], size_t degree)
Definition: closed_form.hpp:374
void AddTests(const long double &)
Definition: poly_root_test_suite.hpp:55
void closed_form_root_test_4()
Definition: poly_root_test_suite.hpp:692
void closed_form_root_test_0()
Definition: poly_root_test_suite.hpp:432
void descartes_test()
Definition: poly_root_test_suite.hpp:78
void AddTests(const float &)
Definition: poly_root_test_suite.hpp:31
void closed_form_root_test_2()
Definition: poly_root_test_suite.hpp:503
void closed_form_root_test_1()
Definition: poly_root_test_suite.hpp:456
void root_finder_test()
Definition: poly_root_test_suite.hpp:940