File indexing completed on 2024-12-15 04:01:13 

0001 /*
0002  * SPDX-FileCopyrightText: 2019-2023 Mattia Basaglia <dev@dragon.best>
0003  *
0004  * SPDX-License-Identifier: GPL-3.0-or-later
0005  */
0006 
0007 #include "polynomial.hpp"
0008 #include "math.hpp"
0009 
0010 
0011 namespace math = glaxnimate::math;
0012 
0013 
0014 // Returns the real cube root of a value
0015 static double cuberoot(double v)
0016 {
0017     if ( v < 0 )
0018         return -math::pow(-v, 1./3);
0019     return math::pow(v, 1./3);
0020 }
0021 
0022 std::vector<double> glaxnimate::math::cubic_roots(double a, double b, double c, double d)
0023 {
0024 
0025     // If a is 0, it's a quadratic
0026     if ( qFuzzyIsNull(a) )
0027         return quadratic_roots(b, c, d);
0028 
0029     // Cardano's algorithm.
0030     b /= a;
0031     c /= a;
0032     d /= a;
0033 
0034     double p = (3*c - b * b) / 3;
0035     double p3 = p / 3;
0036     double q = (2 * b*b*b - 9 * b * c + 27 * d) / 27;
0037     double q2 = q / 2;
0038     double discriminant = q2 * q2 + p3 * p3 * p3;
0039 
0040     // and some variables we're going to use later on:
0041 
0042     // 3 real roots:
0043     if ( discriminant < 0)
0044     {
0045         double mp3  = -p / 3;
0046         double r = math::sqrt(mp3*mp3*mp3);
0047         double t = -q / (2*r);
0048         double cosphi = t < -1 ? -1 : t > 1 ? 1 : t;
0049         double phi  = math::acos(cosphi);
0050         double crtr = cuberoot(r);
0051         double t1   = 2 * crtr;
0052         double root1 = t1 * math::cos(phi / 3) - b / 3;
0053         double root2 = t1 * math::cos((phi + 2 * math::pi) / 3) - b / 3;
0054         double root3 = t1 * math::cos((phi + 4 * math::pi) / 3) - b / 3;
0055         return {root1, root2, root3};
0056     }
0057 
0058     // 2 real roots
0059     if ( qFuzzyIsNull(discriminant) )
0060     {
0061         double u1 = q2 < 0 ? cuberoot(-q2) : -cuberoot(q2);
0062         double root1 = 2*u1 - b / 3;
0063         double root2 = -u1 - b / 3;
0064         return {root1, root2};
0065     }
0066 
0067     // 1 real root, 2 complex roots
0068     double sd = math::sqrt(discriminant);
0069     double u1 = cuberoot(sd - q2);
0070     double v1 = cuberoot(sd + q2);
0071     return {u1 - v1 - b / 3};
0072 }
0073 
0074 std::vector<double> glaxnimate::math::quadratic_roots(double a, double b, double c)
0075 {
0076     // linear
0077     if ( qFuzzyIsNull(a) )
0078     {
0079         if ( qFuzzyIsNull(b) )
0080             return {};
0081 
0082         return {-c / b};
0083     }
0084 
0085     double s = b * b - 4 * a * c;
0086 
0087     // Complex roots
0088     if ( s < 0 )
0089         return {};
0090 
0091     double single_root = -b / (2 * a);
0092 
0093     // 1 root
0094     if ( qFuzzyIsNull(s) )
0095         return {single_root};
0096 
0097     double delta = math::sqrt(s) / (2 * a);
0098 
0099     // 2 roots
0100     return {single_root - delta, single_root + delta};
0101 }