Quadratic equation solver

quadratic.cpp

#include <cassert>
#include <cmath>
#include <iostream>
#include <iomanip>

template <typename T>
std::pair<T, T> naive(T a, T b, T c) {
    T dis = b * b - 4 * a * c;
    assert(dis >= 0);   // not handling complex root
    T sqrtdis = std::sqrt(dis);
    T x1 = (-b - sqrtdis) / (a + a);
    T x2 = (-b + sqrtdis) / (a + a);
    // std::cout << "dis          " << dis << std::endl;
    // std::cout << "sqrtdis      " << sqrtdis << std::endl;
    // std::cout << "-b - sqrtdis " << -b - sqrtdis << std::endl;
    // std::cout << "-b + sqrtdis " << -b + sqrtdis << std::endl;
    return std::make_pair(x1, x2);
}

template <typename T>
T sgn(T x) {
    if      (x > 0) return  1;
    else if (x < 0) return -1;
    else            return  0;
}

// https://en.wikipedia.org/wiki/Loss_of_significance#A_better_algorithm
template <typename T>
std::pair<T, T> better(T a, T b, T c) {
    if (b == 0)
        return naive(a, b, c);
    T dis = b * b - 4 * a * c;
    assert(dis >= 0);   // not handling complex root
    T sqrtdis = std::sqrt(dis);
    T x1 = (-b - sgn(b) * sqrtdis) / (a + a);
    T x2 = c / (a * x1);
    return std::make_pair(x1, x2);
}

template <typename T>
void output(const char* fname, const std::pair<T, T>& p) {
    std::cout << fname << "<" << typeid(T).name() << ">: " << p.first << ", " << p.second << std::endl;
}

int main() {
    std::cout << std::fixed << std::setprecision(0);
    output("naive",  naive (1.0f, -100000000000000000001.0f, 100000000000000000000.0f));
    output("naive",  naive (1.0 , -100000000000000000001.0 , 100000000000000000000.0 ));
    output("naive",  naive (1.0L, -100000000000000000001.0L, 100000000000000000000.0L));
    output("better", better(1.0f, -100000000000000000001.0f, 100000000000000000000.0f));
    output("better", better(1.0 , -100000000000000000001.0 , 100000000000000000000.0 ));
    output("better", better(1.0L, -100000000000000000001.0L, 100000000000000000000.0L));
}

Experiment code for answering a question.

$ g++ quadratic.cpp && ./a.out
naive<f>: -inf, inf
naive<d>: 0, 100000000000000000000
naive<e>: 4, 100000000000000000000
better<f>: inf, 0
better<d>: 100000000000000000000, 1
better<e>: 100000000000000000000, 1