Problem Statement

Write a program to find the roots of a quadratic equation, i.e., solve the equation

a

x

2

b x + c

0

{\displaystyle ax^{2}+bx+c=0}

. Your program must correctly handle non-real roots, but it need not check that

a ≠ 0

{\displaystyle a\neq 0}

. The problem of solving a quadratic equation is a good example of how dangerous it can be to ignore the peculiarities of floating-point arithmetic. The obvious way to implement the quadratic formula suffers catastrophic loss of accuracy when one of the roots to be found is much closer to 0 than the other. In their classic textbook on numeric methods Computer Methods for Mathematical Computations, George Forsythe, Michael Malcolm, and Cleve Moler suggest trying the naive algorithm with

a

1

{\displaystyle a=1}

,

b

−

10

5

{\displaystyle b=-10^{5}}

, and

c

1

{\displaystyle c=1}

. (For double-precision floats, set

b

−

10

9

{\displaystyle b=-10^{9}}

.) Consider the following implementation in Ada: As we can see, the second root has lost all significant figures. The right answer is that X2 is about

10

− 6

{\displaystyle 10^{-6}}

. The naive method is numerically unstable. Suggested by Middlebrook (D-OA), a better numerical method: to define two parameters

q

a c

/

b

{\displaystyle q={\sqrt {ac}}/b}

and

f

1

/

2 +

1 − 4

q

2

/

2

{\displaystyle f=1/2+{\sqrt {1-4q^{2}}}/2}

and the two roots of the quardratic are:

− b

a

f

{\displaystyle {\frac {-b}{a}}f}

and

− c

b f

{\displaystyle {\frac {-c}{bf}}}

Task: do it better. This means that given

a

1

{\displaystyle a=1}

,

b

−

10

9

{\displaystyle b=-10^{9}}

, and

c

1

{\displaystyle c=1}

, both of the roots your program returns should be greater than

10

− 11

{\displaystyle 10^{-11}}

. Or, if your language can't do floating-point arithmetic any more precisely than single precision, your program should be able to handle

b

−

10

6

{\displaystyle b=-10^{6}}

. Either way, show what your program gives as the roots of the quadratic in question. See page 9 of "What Every Scientist Should Know About Floating-Point Arithmetic" for a possible algorithm.

Let's start with the solution:

import "/complex" for Complex

var quadratic = Fn.new { |a, b, c|
    var d = b*b - 4*a*c
    if (d == 0) {
        // single root
        return [[-b/(2*a)], null]
    }
    if (d > 0) {
        // two real roots
        var sr = d.sqrt
        d = (b < 0) ? sr - b : -sr - b
        return [[d/(2*a), 2*c/d], null]
    }
    // two complex roots
    var den = 1 / (2*a)
    var t1 = Complex.new(-b*den, 0)
    var t2 = Complex.new(0, (-d).sqrt * den)
    return [[], [t1+t2, t1-t2]]
}

var test = Fn.new { |a, b, c|
    System.write("coefficients: %(a), %(b), %(c) -> ")
    var roots = quadratic.call(a, b, c)
    var r = roots[0]
    if (r.count == 1) {
        System.print("one real root: %(r[0])")
    } else if (r.count == 2) {
        System.print("two real roots: %(r[0]) and %(r[1])")
    } else {
        var i = roots[1]
        System.print("two complex roots: %(i[0]) and %(i[1])")
    }
}

var coeffs = [
    [1, -2, 1],
    [1,  0, 1],
    [1, -10, 1],
    [1, -1000, 1],
    [1, -1e9, 1]
]

for (c in coeffs) test.call(c[0], c[1], c[2])