#include <cmath>
#include <cassert>

enum domain {
    POS_INF,
    POS_FINITE,
    POS_ZERO,
    NEG_ZERO,
    NEG_FINITE,
    NEG_INF
};

template <typename T>
domain get_domain( T x ) {
    if ( std::isinf( x ) ) {
        return ( std::signbit( x ) ) ? NEG_INF : POS_INF;
    } else if ( x == 0.0 ) {
        return ( std::signbit( x ) ) ? NEG_ZERO : POS_ZERO;
    } else {
        return ( std::signbit( x ) ) ? NEG_FINITE : POS_FINITE;
    }
}

// Calculate the real roots of the quadratic
//        2
//    a  x  + a  x + a  = 0
//     2       1      0
//
// If the roots are complex or undefined, a std::pair (NaN, NaN) is returned,
// otherwise, a pair ('r1', 'r2') is returned where 'r1' <= 'r2'

template <typename T>
std::pair<T, T> quadratic_roots( T a2, T a1, T a0 ) {
    // It is critical that we assume that if any arguments are +/-0.0,
    // that coefficient may be either a very small positive or negative number

    std::pair<T, T> const nan_pair{std::make_pair( std::nan( "" ), std::nan( "" ) )};

    if ( std::isnan( a2 ) || std::isnan( a1 ) || std::isnan( a0 ) ) {
        return nan_pair;
    }

    domain da2{get_domain( a2 )};
    domain da1{get_domain( a1 )};
    domain da0{get_domain( a0 )};

    // At most can be +inf or -inf

    bool a2inf{ (da2 == POS_INF) || (da2 == NEG_INF) };
    bool a1inf{ (da1 == POS_INF) || (da1 == NEG_INF) };
    bool a0inf{ (da0 == POS_INF) || (da0 == NEG_INF) };

    // For now, we will assume that there are no solutions
    // if any of the coefficients are +/- infinity
    //  - example: if the coefficient of x^2 is inf, and
    //    the constant coefficient is small, then it can be
    //    guaranteed that at least one root is +/- 0.0

    if ( a0inf || a1inf || a2inf ) {
        return nan_pair;
    }
        
    // A special case
    //  - Note that if a0 is small (+/- 0.0), we may be able to
    //    find that +/-inf are the roots, or possibly 0.0
    //    depending on the magnitude of a1.

    if ( a2 == 0.0 ) {
        return nan_pair;
    }

    // A special case
    if ( a1 == 0.0 ) {
        T root{ -a0/a2 };

        if ( (root < 0.0) || ((root == 0.0) && std::signbit( root )) ) {
            return nan_pair;
        } else {
            root = std::sqrt( root );
            return std::make_pair( -root, root );
        }
    }

    T disc{ a1*a1 - 4.0*a0*a2 };

    if ( (disc < 0) || ((disc == 0.0) && std::signbit( disc )) ) {
        return nan_pair;
    }

    disc = std::sqrt( disc );

    double b{ -0.5 - std::sqrt( 0.25 - a0*a2/a1/a1 ) };

    if ( a1 == 0.0 ) {
        T root{ disc/2.0/a2 };
        return std::make_pair( -root, root );
    } else if ( a1 >= 0.0 ) {
        if ( a2 > 0.0 ) {
            return std::make_pair( (a1/a2)*b, a0/(a1*b) );
        } else {
            return std::make_pair( a0/(a1*b), (a1/a2)*b );
        }
    } else {
        assert( a1 < 0.0 );

        if ( a2 > 0.0 ) {
            return std::make_pair( a0/(a1*f), (a1/a2)*f );
            return std::make_pair( -2*a0/(a1 - disc), (-a1 + disc)/2.0/a2 );
        } else {
            return std::make_pair( (a1/a2)*f, a0/(a1*f) );
            return std::make_pair( (-a1 + disc)/2.0/a2, -2*a0/(a1 - disc) );
        }
    }
}