一元二次方程求解C++实现

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

    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;

    }