diff --git a/toonz/sources/common/tcore/tmathutil.cpp b/toonz/sources/common/tcore/tmathutil.cpp index 4ce50a9..1d31cab 100644 --- a/toonz/sources/common/tcore/tmathutil.cpp +++ b/toonz/sources/common/tcore/tmathutil.cpp @@ -738,6 +738,114 @@ int rootFinding(const std::vector &in_poly, std::vector &sol) { //----------------------------------------------------------------------------- +int solveEquation2(std::complex *roots, double a, double b, double c) { + if (isAlmostZero(a)) { + if (isAlmostZero(b)) return 0; + roots[0] = -c/b; + return 1; + } + + double D = b*b - a*c*4; + double k = 0.5/a; + std::complex d = D < 0 + ? std::complex(0, sqrt(-D)) + : std::complex(sqrt(D)); + roots[0] = (-b - d)*k; + roots[1] = (-b + d)*k; + return 2; +} + +//----------------------------------------------------------------------------- + +int solveEquation3(std::complex *roots, double a, double b, double c, double d) { + if (isAlmostZero(a)) + return solveEquation2(roots, b, c, d); + + // x = y - b/(3*a) + // y*y*y + p*y + q = 0 + double p = (3*a*c - b*b)/(3*a*a); + double q = (27*a*a*d - 9*a*b*c + 2*b*b*b)/(27*a*a*a); + + double Q = p*p*p/27 + q*q/4; + std::complex Qs = Q < 0 + ? std::complex(0, sqrt(-Q)) + : std::complex(sqrt(Q)); + std::complex A = pow(-q/2 + Qs, 1.0/3); + std::complex B = pow(-q/2 - Qs, 1.0/3); + + // choose complimentary B for A (A*B must be equal -p/3) + std::complex rot(-0.5, sqrt(3.0)/2); + if (!isAlmostZero(A*B + p/3)) B *= rot; + if (!isAlmostZero(A*B + p/3)) B *= rot; + + std::complex Y = (A - B)*std::complex(0, sqrt(3.0)); + std::complex y0 = A + B; + std::complex y1 = (-y0 - Y)/2.0; + std::complex y2 = (-y0 + Y)/2.0; + + double dd = b/(3*a); + roots[0] = y0 - dd; + roots[1] = y1 - dd; + roots[2] = y2 - dd; + return 3; +} + +//----------------------------------------------------------------------------- + +int solveEquation4(std::complex *roots, double a, double b, double c, double d, double e) { + if (isAlmostZero(a)) + return solveEquation3(roots, b, c, d, e); + + // x = y - b/(4*a) + // y^4 + p*y^2 + q*y + r = 0 + double dd = b/(4*a); + double p = (8*a*c - 3*b*b)/(8*a*a); + double q = (8*a*a*d - 4*a*b*c + b*b*b)/(8*a*a*a); + double r = (256*a*a*a*e - 64*a*a*b*d + 16*a*b*b*c - 3*b*b*b*b)/(256*a*a*a*a); + + if (isAlmostZero(q)) { + // biquadratic equation + // y^4 + p*y^2 + r = 0 + // handling this special case will give us more accurate results + std::complex y[2]; + solveEquation2(y, 1, p, r); + y[0] = sqrt(y[0]); + y[1] = sqrt(y[1]); + roots[0] = -y[0] - dd; + roots[1] = y[0] - dd; + roots[2] = -y[1] - dd; + roots[3] = y[1] - dd; + return 4; + } + + // solve cubic equation + // z*z*z + (p/2)*z*z + ((p*p - 4*r)/16)*z - q*q/64 = 0 + double pp = p/2; + double qq = (p*p - 4*r)/16; + double rr = -q*q/64; + std::complex z[3]; + solveEquation3(z, 1, pp, qq, rr); + + z[0] = sqrt(z[0]); + z[1] = sqrt(z[1]); + z[2] = sqrt(z[2]); + + // we need to find signs combination where following is valid: + // (+-z0)*(+-z1)*(+-z2) = -q/8 + // magnitudes are always equals, we need to check signs only + std::complex zzz = z[0]*z[1]*z[2]; + if ((zzz.real() > 0) == (q > 0)) + z[0] = -z[0]; + + roots[0] = z[0] - z[1] - z[2] - dd; + roots[1] = -z[0] + z[1] - z[2] - dd; + roots[2] = -z[0] - z[1] + z[2] - dd; + roots[3] = z[0] + z[1] + z[2] - dd; + return 4; +} + +//----------------------------------------------------------------------------- + /* */ int numberOfRootsInInterval(int order, const double *polyH, double min, diff --git a/toonz/sources/include/tmathutil.h b/toonz/sources/include/tmathutil.h index df56563..65e2637 100644 --- a/toonz/sources/include/tmathutil.h +++ b/toonz/sources/include/tmathutil.h @@ -7,6 +7,7 @@ #include "texception.h" #include +#include #ifndef __sgi #include @@ -128,6 +129,48 @@ DVAPI int rootFinding(const std::vector &poly, //----------------------------------------------------------------------------- /*! + Find complex roots of a quadratic equation + a*x^2 + b*x + c = 0 + \ret count of roots found: + 0 if a and b are almost zero + 1 if a is almost zero + 2 if a is not almost zero + roots may have duplicates, count depends only from type of equation + */ +int solveEquation2(std::complex *roots, double a, double b, double c); + +//----------------------------------------------------------------------------- + +/*! + Find complex roots of a cubic equation: + a*x^3 + b*x^2 + c*x + d = 0 + \ret count of roots found: + 0 if a, b and c are almost zero + 1 if a and b are almost zero + 2 if a is almost zero + 3 if a is not almost zero + roots may have duplicates, count depends only from type of equation + */ +int solveEquation3(std::complex *roots, double a, double b, double c, double d); + +//----------------------------------------------------------------------------- + +/*! + Find complex roots of a power of four equation: + a*x^4 + b*x^3 + c*x^2 + d*x + e = 0 + \ret count of roots found: + 0 if a, b, c and d are almost zero + 1 if a, b and c are almost zero + 2 if a and b are almost zero + 3 if a is almost zero + 4 if a is not almost zero + roots may have duplicates, count depends only from type of equation + */ +int solveEquation4(std::complex *roots, double a, double b, double c, double d, double e); + +//----------------------------------------------------------------------------- + +/*! Check if val is nearest to epsilon. \par val value to test \par eps tolerance required @@ -140,6 +183,19 @@ inline bool isAlmostZero(double val, double eps = TConsts::epsilon) { //----------------------------------------------------------------------------- /*! + Check if complex val is nearest to epsilon. + \par val value to test + \par eps tolerance required + \ret true if value is nearest to zero + */ +inline bool isAlmostZero(std::complex val, double eps = TConsts::epsilon) { + return val.real()*val.real() + val.imag()*val.imag() < eps*eps; +} + + +//----------------------------------------------------------------------------- + +/*! Find number of roots of a poly in interval min max. \par order degree of polynomious \par coeff of poly in crescent order