13 #ifndef eli_mutil_poly_root_closed_form_hpp
14 #define eli_mutil_poly_root_closed_form_hpp
19 #include "eli/code_eli.hpp"
31 template<
typename data__>
37 return std::numeric_limits<int>::max();
42 template<
typename data__>
53 template<
typename data__>
59 data__ desc=a[1]*a[1]-4*a[2]*a[0];
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);
72 root[0]=-a[1]/(2*a[2]);
77 data__ term=-(a[1]+((a[1]<0) ? -1 : 1)*std::sqrt(desc))/2;
81 std::sort(root.begin(), root.end());
87 template<
typename data__>
93 data__ aa(a[2]/a[3]), bb(a[1]/a[3]), cc(a[0]/a[3]);
95 data__ q(aa*aa/9-bb/3), r, q3, r2;
97 r=aa*q/3-aa*bb/18+cc/2;
115 data__ cbrt_r(std::cbrt(r));
119 root[0]=-2*cbrt_r-aa/3;
123 std::sort(root.begin(), root.end());
130 data__ sqrt_q(std::sqrt(q));
133 theta=std::acos(r/std::sqrt(q3));
135 root[0]=-2*sqrt_q*std::cos(theta/3)-aa/3;
139 std::sort(root.begin(), root.end());
145 newa=-((r<0) ? -1 : 1)*std::cbrt(std::abs(r)+std::sqrt(r2-q3));
147 root[0]=(newa+q/newa)-aa/3;
152 template<
typename data__>
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)));
171 data__ alpha(0), beta(1/3), gamma(0);
172 delta=alpha*a*a+beta*b;
176 data__ as, bs, cs, ds, A, B, C, y;
178 std::vector<data__> ysroot;
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;
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;
196 y=ysroot.back()+delta;
207 if (std::abs(y)<=sz*sz*std::numeric_limits<data__>::epsilon())
208 y=
static_cast<data__
>(0);
213 data__ gprime, hprime, gdes, hdes, sigma;
214 bool gcomplex, hcomplex;
227 if (std::abs(gdes)<=3*sz*sz*std::numeric_limits<data__>::epsilon())
228 gdes=static_cast<data__>(0);
234 gdes=std::sqrt(-gdes)/2;
239 gdes=std::sqrt(gdes)/2;
253 if (std::abs(hdes)<=2*sz*sz*std::numeric_limits<data__>::epsilon())
254 hdes=
static_cast<data__
>(0);
260 hdes=std::sqrt(-hdes)/2;
265 hdes=std::sqrt(hdes)/2;
270 if ((gdes!=0) && (hdes!=0))
275 if (gprime*hprime*gdes*hdes>=0)
281 else if ((gprime==0) || (hprime==0))
290 data__ sig_cal=(gprime*hprime-c/2)/(gdes*hdes);
293 assert((std::abs(sig_cal)>0.9) && (std::abs(sig_cal)<1.1));
310 data__ x12des, x34des;
311 x12des=gprime*gprime+(gcomplex ? -1 : 1)*gdes*gdes-4*hprime;
313 if (gcomplex || hcomplex)
315 if (hcomplex && !gcomplex)
327 x12des+=2*gprime*gdes-4*sigma*hdes;
336 if (std::abs(x12des)<=90*sz*sz*std::numeric_limits<data__>::epsilon())
337 x12des=
static_cast<data__
>(0);
339 x34des+=-2*gprime*gdes+4*sigma*hdes;
348 if (std::abs(x34des)<=90*sz*sz*std::numeric_limits<data__>::epsilon())
349 x34des=static_cast<data__>(0);
354 x12des=std::sqrt(x12des);
355 root.push_back((-(gprime+gdes)+x12des)/2);
356 root.push_back((-(gprime+gdes)-x12des)/2);
360 x34des=std::sqrt(x34des);
361 root.push_back((-(gprime-gdes)+x34des)/2);
362 root.push_back((-(gprime-gdes)-x34des)/2);
369 std::sort(root.begin(), root.end());
370 return static_cast<int>(root.size());
373 template<
typename data__>
374 int closed_form(std::vector<data__> &root,
const data__ a[],
size_t degree)
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
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