Code-Eli  0.3.6
closed_form.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 eli_mutil_poly_root_closed_form_hpp
14 #define eli_mutil_poly_root_closed_form_hpp
15 
16 #include <vector>
17 #include <limits>
18 
19 #include "eli/code_eli.hpp"
20 
21 #include "eli/constants/math.hpp"
22 
23 namespace eli
24 {
25  namespace mutil
26  {
27  namespace poly
28  {
29  namespace root
30  {
31  template<typename data__>
32  int closed_form_0(std::vector<data__> &root, const data__ a[])
33  {
34  root.resize(0);
35 
36  if (a[0]==0)
37  return std::numeric_limits<int>::max();
38 
39  return 0;
40  }
41 
42  template<typename data__>
43  int closed_form_1(std::vector<data__> &root, const data__ a[])
44  {
45  if (a[1]==0)
46  return closed_form_0(root, a);
47 
48  root.resize(1);
49  root[0]=-a[0]/a[1];
50  return 1;
51  }
52 
53  template<typename data__>
54  int closed_form_2(std::vector<data__> &root, const data__ a[])
55  {
56  if (a[2]==0)
57  return closed_form_1(root, a);
58 
59  data__ desc=a[1]*a[1]-4*a[2]*a[0];
60  if (desc<0)
61  {
62  root.resize(0);
63  return 0;
64  }
65  root.resize(2);
66 
67  data__ sz=std::max(std::max(std::abs(a[0]), std::abs(a[1])), std::abs(a[2]));
68  if (desc<=sz*std::numeric_limits<data__>::epsilon())
69  desc=static_cast<data__>(0);
70  if (desc==0)
71  {
72  root[0]=-a[1]/(2*a[2]);
73  root[1]=root[0];
74  return 2;
75  }
76 
77  data__ term=-(a[1]+((a[1]<0) ? -1 : 1)*std::sqrt(desc))/2;
78  root[0]=term/a[2];
79  root[1]=a[0]/term;
80 
81  std::sort(root.begin(), root.end());
82  return 2;
83  }
84 
87  template<typename data__>
88  int closed_form_3(std::vector<data__> &root, const data__ a[])
89  {
90  if (a[3]==0)
91  return closed_form_2(root, a);
92 
93  data__ aa(a[2]/a[3]), bb(a[1]/a[3]), cc(a[0]/a[3]);
94 
95  data__ q(aa*aa/9-bb/3), r, q3, r2;
96  q=aa*aa/9-bb/3;
97  r=aa*q/3-aa*bb/18+cc/2;
98  q3=q*q*q;
99  r2=r*r;
100 
101  // unique real root and double real root or triple root
102  if (r2==q3)
103  {
104  // triple real root
105  if (aa*aa==3*bb)
106  {
107  root.resize(3);
108  root[0]=-aa/3;
109  root[1]=-aa/3;
110  root[2]=-aa/3;
111 
112  return 3;
113  }
114 
115  data__ cbrt_r(std::cbrt(r));
116 
117  root.resize(3);
118 
119  root[0]=-2*cbrt_r-aa/3;
120  root[1]=cbrt_r-aa/3;
121  root[2]=root[1];
122 
123  std::sort(root.begin(), root.end());
124  return 3;
125  }
126 
127  // three real roots
128  if (r2<q3)
129  {
130  data__ sqrt_q(std::sqrt(q));
131  data__ theta;
132 
133  theta=std::acos(r/std::sqrt(q3));
134  root.resize(3);
135  root[0]=-2*sqrt_q*std::cos(theta/3)-aa/3;
136  root[1]=-2*sqrt_q*std::cos((theta+eli::constants::math<data__>::two_pi())/3)-aa/3;
137  root[2]=-2*sqrt_q*std::cos((theta-eli::constants::math<data__>::two_pi())/3)-aa/3;
138 
139  std::sort(root.begin(), root.end());
140  return 3;
141  }
142 
143  // must have only one real root
144  data__ newa;
145  newa=-((r<0) ? -1 : 1)*std::cbrt(std::abs(r)+std::sqrt(r2-q3));
146  root.resize(1);
147  root[0]=(newa+q/newa)-aa/3;
148 
149  return 1;
150  }
151 
152  template<typename data__>
153  int closed_form_4(std::vector<data__> &root, const data__ aa[])
154  {
155  // data__ zero(0), one(1), two(2), three(3), four(4), six(6), ninety(90);
156 
157  if (aa[4]==0)
158  return closed_form_3(root, aa);
159 
162  data__ a=aa[3]/aa[4], b=aa[2]/aa[4], c=aa[1]/aa[4], d=aa[0]/aa[4];
163  data__ sz=std::max(std::max(std::abs(a), std::abs(b)), std::max(std::abs(c), std::abs(d)));
164  data__ delta, DELTA;
165  // std::cout << "\ta=" << a << "\tb=" << b << "\tc=" << c << "\td=" << d << std::endl;
166 
167  // NOTE: These can actually be set to a variety of values
168  // data__ alpha(2), beta(3), gamma(0); // made up
169  // data__ alpha(0), beta(0), gamma(0); // Ferrari-Lagrange
170  // data__ alpha(-1.0/4.0), beta(1), gamma(-1.0/4.0); // Descartes-Euler-Cardano
171  data__ alpha(0), beta(1/3), gamma(0); // new scheme
172  delta=alpha*a*a+beta*b;
173  DELTA=gamma*a;
174 
175  // set the cubic coefficients and solve cubic equation for y
176  data__ as, bs, cs, ds, A, B, C, y;
177  data__ ay[4];
178  std::vector<data__> ysroot;
179  as=a+4*DELTA;
180  bs=b+(3*a+6*DELTA)*DELTA;
181  cs=c+(2*b+(3*a+4*DELTA)*DELTA)*DELTA;
182  ds=d+(c+(b+(a+DELTA)*DELTA)*DELTA)*DELTA;
183  A=3*delta-bs;
184  B=(3*delta-2*bs)*delta+as*cs-4*ds;
185  C=((delta-bs)*delta+(as*cs-4*ds))*delta+4*bs*ds-as*as*ds-cs*cs;
186  ay[0]=C; ay[1]=B; ay[2]=A; ay[3]=1;
187  closed_form_3(ysroot, ay);
188  // std::cout << "A=" << A << "\tB=" << B << "\tC=" << C << std::endl;
189  // std::cout << "A=" << 0 << "\tB=" << a*c-b*b/3-4*d << "\tC=" << a*b*c/3-a*a*d-2*b*b*b/27-c*c+8*b*d/3 << std::endl;
190  // for (size_t i=0; i<ysroot.size(); ++i)
191  // std::cout << "\tys[" << i << "]=" << ysroot[i];
192  // std::cout << std::endl;
193 
194  // get the largest root to try and minimize chance of complex roots below
195  // std::cout << "nyroots=" << ysroot.size() << std::endl;
196  y=ysroot.back()+delta;
197  // NOTE: if ys is approximately delta, then can have loss of significance and
198  // really should chose a different delta.
199  // Note: have to compensate for low precision trig function calls in getting root of cubic
200  // if (typeid(data__)==typeid(long double))
201  // {
202  // if (std::abs(y)<=sz*sz*std::numeric_limits<double>::epsilon())
203  // y=zero;
204  // }
205  // else
206  {
207  if (std::abs(y)<=sz*sz*std::numeric_limits<data__>::epsilon())
208  y=static_cast<data__>(0);
209  }
210 
211  // solve the two quadratic equations. Need to track the discriminants in case complex
212  // values for g or h
213  data__ gprime, hprime, gdes, hdes, sigma;
214  bool gcomplex, hcomplex;
215 
216  gprime=a/2;
217  gdes=a*a-4*b+4*y;
218  // std::cout << "gdes=" << gdes << "\tepsilon=" << 3*sz*sz*std::numeric_limits<data__>::epsilon() << std::endl;
219  // Note: have to compensate for low precision trig function calls in getting root of cubic
220  // if (typeid(data__)==typeid(long double))
221  // {
222  // if (std::abs(gdes)<=30*sz*sz*std::numeric_limits<double>::epsilon())
223  // gdes=zero;
224  // }
225  // else
226  {
227  if (std::abs(gdes)<=3*sz*sz*std::numeric_limits<data__>::epsilon())
228  gdes=static_cast<data__>(0);
229  }
230  if (gdes<0)
231  {
232  // std::cout << "Complex g quadratic roots!" << std::endl;
233  gcomplex=true;
234  gdes=std::sqrt(-gdes)/2;
235  }
236  else
237  {
238  gcomplex=false;
239  gdes=std::sqrt(gdes)/2;
240  }
241 
242  hprime=y/2;
243  hdes=y*y-4*d;
244  // std::cout << "hdes=" << hdes << "\tepsilon=" << 2*sz*sz*std::numeric_limits<data__>::epsilon() << std::endl;
245  // Note: have to compensate for low precision trig function calls in getting root of cubic
246  // if (typeid(data__)==typeid(long double))
247  // {
248  // if (std::abs(hdes)<=20*sz*sz*std::numeric_limits<double>::epsilon())
249  // hdes=zero;
250  // }
251  // else
252  {
253  if (std::abs(hdes)<=2*sz*sz*std::numeric_limits<data__>::epsilon())
254  hdes=static_cast<data__>(0);
255  }
256  if (hdes<0)
257  {
258  // std::cout << "Complex h quadratic roots!" << std::endl;
259  hcomplex=true;
260  hdes=std::sqrt(-hdes)/2;
261  }
262  else
263  {
264  hcomplex=false;
265  hdes=std::sqrt(hdes)/2;
266  }
267 
268  // need to find correct sign for h discriminant (sigma) that solves original equations
269  // std::cout << "gprime=" << gprime << "\tgdes=" << gdes << "\thprime=" << hprime << "\thdes=" << hdes << std::endl;
270  if ((gdes!=0) && (hdes!=0))
271  {
272  // catch special case when c is zero
273  if (c==0)
274  {
275  if (gprime*hprime*gdes*hdes>=0)
276  sigma=1;
277  else
278  sigma=-1;
279  }
280  // catch special case when gprime or hprime are zero
281  else if ((gprime==0) || (hprime==0))
282  {
283  if (gdes*hdes>0)
284  sigma=-1;
285  else
286  sigma=1;
287  }
288  else
289  {
290  data__ sig_cal=(gprime*hprime-c/2)/(gdes*hdes);
291 
292  // sig_cal should be very close to -1 or +1
293  assert((std::abs(sig_cal)>0.9) && (std::abs(sig_cal)<1.1));
294  if (sig_cal>=0)
295  sigma=1;
296  else
297  sigma=-1;
298  }
299  // std::cout << "sigma=" << sigma << std::endl;
300  // std::cout << "sigma(1)=" << (gprime*hprime-c/2)/(gdes*hdes) << std::endl;
301  // std::cout << "zero?=" << (gprime+gdes)*(hprime-hdes)+(gprime-gdes)*(hprime+hdes)-c << std::endl;
302  // std::cout << "zero?=" << (gprime+gdes)*(hprime+hdes)+(gprime-gdes)*(hprime-hdes)-c << std::endl;
303  }
304  else
305  {
306  sigma=1;
307  }
308 
309  // now calculate the 4 roots
310  data__ x12des, x34des;
311  x12des=gprime*gprime+(gcomplex ? -1 : 1)*gdes*gdes-4*hprime;
312  x34des=x12des;
313  if (gcomplex || hcomplex)
314  {
315  if (hcomplex && !gcomplex)
316  {
317  root.resize(0);
318  return 0;
319  }
320 
321  // should we ever get here?
322  assert(false);
323  }
324  else
325  {
326  // no complex math needs to be done
327  x12des+=2*gprime*gdes-4*sigma*hdes;
328  // Note: have to compensate for low precision math in long double
329  // if (typeid(data__)==typeid(long double))
330  // {
331  // if (std::abs(x12des)<=ninety*sz*sz*std::numeric_limits<double>::epsilon())
332  // x12des=zero;
333  // }
334  // else
335  {
336  if (std::abs(x12des)<=90*sz*sz*std::numeric_limits<data__>::epsilon())
337  x12des=static_cast<data__>(0);
338  }
339  x34des+=-2*gprime*gdes+4*sigma*hdes;
340  // Note: have to compensate for low precision math in long double
341  // if (typeid(data__)==typeid(long double))
342  // {
343  // if (std::abs(x34des)<=ninety*sz*sz*std::numeric_limits<double>::epsilon())
344  // x34des=zero;
345  // }
346  // else
347  {
348  if (std::abs(x34des)<=90*sz*sz*std::numeric_limits<data__>::epsilon())
349  x34des=static_cast<data__>(0);
350  }
351  // std::cout << "x12des=" << x12des << "\tx34des=" << x34des << "\tepsilon=" << 90*sz*sz*std::numeric_limits<data__>::epsilon() << std::endl;
352  if (x12des>=0)
353  {
354  x12des=std::sqrt(x12des);
355  root.push_back((-(gprime+gdes)+x12des)/2);
356  root.push_back((-(gprime+gdes)-x12des)/2);
357  }
358  if (x34des>=0)
359  {
360  x34des=std::sqrt(x34des);
361  root.push_back((-(gprime-gdes)+x34des)/2);
362  root.push_back((-(gprime-gdes)-x34des)/2);
363  }
364  }
365 
366  // for (size_t i=0; i<root.size(); ++i)
367  // std::cout << "\troot[" << i << "]=" << root[i];
368  // std::cout << std::endl;
369  std::sort(root.begin(), root.end());
370  return static_cast<int>(root.size());
371  }
372 
373  template<typename data__>
374  int closed_form(std::vector<data__> &root, const data__ a[], size_t degree)
375  {
376  switch (degree)
377  {
378  case(0):
379  {
380  return closed_form_0(root, a);
381  break;
382  }
383  case(1):
384  {
385  return closed_form_1(root, a);
386  break;
387  }
388  case(2):
389  {
390  return closed_form_2(root, a);
391  break;
392  }
393  case(3):
394  {
395  return closed_form_3(root, a);
396  break;
397  }
398  case(4):
399  {
400  return closed_form_4(root, a);
401  break;
402  }
403  default:
404  {
405  assert(false);
406  return -1;
407  break;
408  }
409  }
410  }
411  }
412  }
413  }
414 }
415 
416 #endif
Definition: math.hpp:20
int closed_form_2(std::vector< data__ > &root, const data__ a[])
Definition: closed_form.hpp:54
int closed_form_0(std::vector< data__ > &root, const data__ a[])
Definition: closed_form.hpp:32
Definition: math.hpp:25
int closed_form(std::vector< data__ > &root, const data__ a[], size_t degree)
Definition: closed_form.hpp:374
int closed_form_3(std::vector< data__ > &root, const data__ a[])
Definition: closed_form.hpp:88
int closed_form_4(std::vector< data__ > &root, const data__ aa[])
Definition: closed_form.hpp:153
int closed_form_1(std::vector< data__ > &root, const data__ a[])
Definition: closed_form.hpp:43