Quadratic

// Solving a x^2 + b x + c = 0, the classic calculation whose *shape* depends on its input.
//
// With `a` not zero it is a quadratic, and the discriminant b^2 - 4ac decides the answer: two
// real roots when it is positive, one repeated root when it is zero, and a complex pair when it
// is negative. With `a` zero it is not a quadratic at all, and the formula would divide by zero:
//
//     a = 0, b != 0           linear, one root x = -c / b
//     a = 0, b = 0, c != 0    inconsistent, no x makes c equal 0
//     a = 0, b = 0, c = 0     an identity, every x is a solution
//
// Reading comes first, and every way it can go wrong is a case of one error variant: the input
// may end early, fail to read, not spell a number, or spell one that is not finite. The parser
// accepts "Inf" and "NaN", because that is how such values print, so finiteness is checked
// separately.
import Allocator::{ Allocator, SystemAllocator };
import Core::IsFinite;
import Format::ParseFloat64;
import Io::{ IoErrorKind, Print, PrintLine, ReadLine };
import Math::Sqrt;
import Text::StringBuilder;

variant EntryError {
    Ended,
    Unreadable,
    NotANumber,
    NotFinite
}

// Reads one coefficient. The builder is reused for every line, so it is cleared first, and the
// view is trimmed so that " 2 " reads as 2.
func ReadCoefficient(name: char8[..], builder: &var StringBuilder) -> float64 ! EntryError {
    Print("{} = ", name);
    builder.Clear();
    match ReadLine(builder) {
        .Success(_) => {},
        .Failure(error) if error.kind == IoErrorKind::EndOfStream => fail EntryError::Ended,
        .Failure(_) => fail EntryError::Unreadable
    }
    // A number past the largest float64 is an `Overflow`, not a malformed one: it is a number,
    // just not a finite float64, so it joins Inf and NaN.
    let value = ParseFloat64(builder.View().Trim()) catch {
        .Overflow(_) => fail EntryError::NotFinite,
        else => fail EntryError::NotANumber
    };
    if !IsFinite(value) {
        fail EntryError::NotFinite;
    }
    return value;
}

// Each report starts on a fresh line, since piped input leaves the cursor after the prompt.
func Explain(error: EntryError) {
    PrintLine();
    match error {
        .Ended => PrintLine("the input ended before all three coefficients were given"),
        .Unreadable => PrintLine("the input could not be read"),
        .NotANumber => PrintLine("that is not a number"),
        .NotFinite => PrintLine(
            "a coefficient must be finite: Inf, NaN and numbers past float64 are refused")
    }
}

func SolveLinear(b: float64, c: float64) {
    if b != 0.0 {
        // Adding zero turns -0 into 0, so c = 0 prints "x = 0.0" rather than "x = -0.0".
        PrintLine("linear, not quadratic: x = {}", -c / b + 0.0);
    } else if c != 0.0 {
        PrintLine("inconsistent: {} = 0 is never true, so there is no solution", c);
    } else {
        PrintLine("an identity: 0 = 0 holds for every x");
    }
}

func SolveQuadratic(a: float64, b: float64, c: float64) {
    // Finite coefficients can still overflow: b = 1e200 makes b^2 infinite.
    let discriminant = b * b - 4.0 * a * c;
    if !IsFinite(discriminant) {
        PrintLine("the coefficients are too large to solve in float64");
        return;
    }

    if discriminant > 0.0 {
        // The textbook formula subtracts two nearly equal numbers when 4ac is small, and loses
        // the smaller root to rounding. Adding numbers of the same sign never cancels, so find
        // the larger root that way, then the other from x1 * x2 = c / a.
        let root = Sqrt(discriminant);
        let q = b >= 0.0 ? -(b + root) / 2.0 : (root - b) / 2.0;
        let first = q / a;
        let second = c / q;
        let low = first < second ? first : second;
        let high = first < second ? second : first;
        PrintLine("two real roots: x1 = {}, x2 = {}", low, high);
    } else if discriminant == 0.0 {
        PrintLine("one repeated root: x = {}", -b / (2.0 * a) + 0.0);
    } else {
        // No real root, but the complex pair is real part +/- imaginary part. Dividing by |2a|
        // keeps the imaginary part positive, so the signs printed are the signs meant.
        let real = -b / (2.0 * a) + 0.0;
        let twiceA = a < 0.0 ? -2.0 * a : 2.0 * a;
        let imaginary = Sqrt(-discriminant) / twiceA;
        PrintLine("complex roots: x = {} + {}i and x = {} - {}i",
            real, imaginary, real, imaginary);
    }
}

func Main() -> int {
    var system = SystemAllocator();
    let allocator: Allocator = system;

    // The builder's destructor returns its memory on every path out of `Main`.
    var builder = StringBuilder(allocator);

    PrintLine("Solving a x^2 + b x + c = 0");
    let a = ReadCoefficient("a", builder) catch { error => { Explain(error); return 1; } };
    let b = ReadCoefficient("b", builder) catch { error => { Explain(error); return 1; } };
    let c = ReadCoefficient("c", builder) catch { error => { Explain(error); return 1; } };
    PrintLine();

    if a == 0.0 {
        SolveLinear(b, c);
    } else {
        SolveQuadratic(a, b, c);
    }
    return 0;
}