学步园

-----------------------------------------------------------------------------------------------------------------------------

    typedef double Number

    int CesiumMath::sign(Number value){

        if (value > 0) {
              return 1;
          }
          if (value < 0) {
              return -1;
          }

          return 0;
    }

------------------------------------------------------------------------------------------------------------------------------

class QuadraticRealPolynomial{
public:
    static Number computeDiscriminant(Number a, Number b, Number c);
    static std::vector<Number>* computeRealRoots(Number a, Number b, Number c, std::vector<Number>& roots);
private:
    static Number addWithCancellationCheck(Number left, Number right, Number tolerance);
public:
    Number _root1;
    Number _root2;
------------------------------------------------------------------------------------------------------------------------------

    Number QuadraticRealPolynomial::computeDiscriminant(Number a, Number b, Number c)

    {
        Number discriminant = b * b - 4.0 * a * c;
        return discriminant;
    }

    std::vector<Number>* QuadraticRealPolynomial::computeRealRoots(Number a, Number b, Number c, std::vector<Number>& roots)
    {
         Number root1, root2;

         roots.clear();

         Number ratio;
         if (a == 0.0)
         {
             if (b == 0.0)
             {
                 // Constant function: c = 0.
                 return &roots;
             }

             // Linear function: b * x + c = 0.
             root1 = -c/b;
             roots.push_back(root1);
             return &roots;
         }
         else if (b == 0.0)
         {
             if (c == 0.0)
             {
                 // 2nd order monomial: a * x^2 = 0.
                 root1 = 0.0;
                 root2 = 0.0;

                 roots.push_back(root1);
                 roots.push_back(root2);
                 return &roots;
             }

             Number cMagnitude = std::abs(c);
             Number aMagnitude = std::abs(a);

             if ((cMagnitude < aMagnitude) && (cMagnitude / aMagnitude < CesiumMath::_EPSILON14))
             {
                 // c ~= 0.0.
                 // 2nd order monomial: a * x^2 = 0.
                 root1 = 0.0;
                 root2 = 0.0;

                 roots.push_back(root1);
                 roots.push_back(root2);
                 return &roots;

             } else if ((cMagnitude > aMagnitude) && (aMagnitude / cMagnitude < CesiumMath::_EPSILON14))
             {
                 // a ~= 0.0.
                 // Constant function: c = 0.
                 return &roots;
             }

             // a * x^2 + c = 0
             ratio = -c / a;

             if (ratio < 0.0)
             {
                 // Both roots are complex.
                 return &roots;
             }

             // Both roots are real.
             Number root = std::sqrt(ratio);
             root1 = -root;
             root2 = root;

             roots.push_back(root1);
             roots.push_back(root2);
             return &roots;
         }
         else if (c == 0.0)
         {
             // a * x^2 + b * x = 0
             ratio = -b / a;
             if (ratio < 0.0)
             {
                 root1 = ratio;
                 root2 = 0.0;

                 roots.push_back(root1);
                 roots.push_back(root2);
                 return &roots;
             }

             root1 = 0.0;
             root2 = ratio;

             roots.push_back(root1);
             roots.push_back(root2);
             return &roots;
     }

         // a * x^2 + b * x + c = 0
         Number b2 = b * b;
         Number four_ac = 4.0 * a * c;
         Number radicand = addWithCancellationCheck(b2, -four_ac, CesiumMath::_EPSILON14);

         if (radicand < 0.0)
         {
             // Both roots are complex.
             return &roots;
         }

         Number q = -0.5 * addWithCancellationCheck(b, CesiumMath::sign(b) * std::sqrt(radicand), CesiumMath::_EPSILON14);
         if (b > 0.0)
         {
             root1 = q / a;
             root2 = c / q;

             roots.push_back(root1);
             roots.push_back(root2);
             return &roots;
         }

         root1 = c / q;
         root2 = q / a;

         roots.push_back(root1);
         roots.push_back(root2);
         return &roots;
    }

    Number QuadraticRealPolynomial::addWithCancellationCheck(Number left, Number right, Number tolerance)
    {
        Number difference = left + right;
        if ((CesiumMath::sign(left) != CesiumMath::sign(right)) &&
             std::abs(difference / std::max(std::abs(left), std::abs(right))) < tolerance)
        {
            return 0.0;
        }

        return difference;
    }