The following branchless solver was proposed, to avoid cancellation error suffered by the conventional formula in the neighbourhoods of x:y = 0:1 (zero root) and 1:0 (infinite root) --- The homogeneous quadratic equation f(x,y) == a x^2 - 2 b x y + c y^2 = 0 has solutions (A) x : y = (c + b + d) : (a + b - d) , (c + b - d) : (a + b + d) , where d^2 = b^2 - a c . Exercise: show that the roots x : y = e : 1, 1 : e of the equation x^2 - (e + 1/e)x y + y^2 = 0 are given by (A) to working precision r when e is small. Nobody has yet pointed out the embarrassing fact that in filling one hole, a larger has inadvertantly been dug: given say f = x^2 - y^2 with two roots +1,-1, method (A) will locate only one of them, returning 0:0 (undefined) for the other! More generally, Exercise: show that the roots x : y = e + 1 : 1, e - 1 : 1 of the equation x^2 - 2 e x y + (e^2 - 1)y^2 = 0 are given by (A) to working precision r for one, and with error r/e for the other, when e is small. A closer look reveals that (A) suffers excessive cancellation when (a + b) and (c + b) large and (a + b) + (c + b) small; or (a - b) and (c - b) large and (a - b) + (c - b) small; equivalently, when (6) |2 b| ~ |a + c| << |a| ~ |c| . We might try to repair case (6) by generalising (A) to the form (B) x : y = ((b + d)g + c h) : (a g + (b - d)h), x : y = ((b - d)g + c h) : (a g + (b + d)h), where g = g(a,b,c), h = h(a,b,c) are variable weight functions. This breaks down (and will be ill-conditioned in the neighbourhood) where (c h + (b + d)g) = (a g + (b - d)h) = 0 or (c h + (b - d)g) = (a g + (b + d)h) = 0; rationalising, (c h + (b + d)g)(c h + (b - d)g) = 0 and (a g + (b - d)h)(a g + (b + d)h) = 0; that is, (a g^2 + 2 b g h + c h^2) c = 0 and (a g^2 + 2 b g h + c h^2) a = 0. Finally then, we require to avoid either (7) a = c = 0; or (8) roots (a,b,c) of f(g,h) = 0 . We can't do much about case (7), which in practice does not by itself seem to cause excessive cancellation --- I think the reason is connected with the fact that relative error becomes meaningless near roots 0, oo. However to avoid case (8), all real points of the homogeneous plane curve f(g,h) = 0 would be required to lie within the parabolic region b^2 < a c. But suppose g,h are polynomial in a,b,c, with the same homogeneous degree; then f(g,h) has odd total degree, so must have some real point. Now transform the parabola projectively into an ellipse containing that point: the odd-degree branch still extends to infinity, so it must intersect the ellipse; and its pre-image intersects the parabola. This shows that a branchless solution of the form (B) is impossible with homogeneous polynomial weights; and I suspect the result may be extended to analytic weights and beyond. Perhaps a more natural strategy is to reverse weights for second root: (C) x : y = ((b + d)g + c h) : (a g + (b - d)h), x : y = ((b - d)h + c g) : (a h + (b + d)g); Kahan's method (adapted to unscaled b) uses (C) with g = 1 + sign(b), h = 1 - sign(b); it breaks down at b = c = 0 where both roots are zero, although it remains well-conditioned in the neighbourhood. More generally (C) fails in much the same way as (B); now both f(g,h) and f(h,g) are required to lack real zeros outside the parabola. Incidentally, preliminary rescaling (1) is itself nontrivial to achieve branchlessly: the fastest strategy --- where possible --- may be to be to omit it altogether initially, instead setting exception traps for exponent overflow and underflow. Fred Lunnon On 5/18/10, Fred lunnon <fred.lunnon@gmail.com> wrote:

Several authors have analysed the pitfalls involved in the numerical evaluation of the traditional schoolroom solution of a quadratic equation in one variable x: a x^2 - 2 b x + c = 0 ; x = (b + d)/a , (b - d)/a , where d^2 = b^2 - a c .

Very briefly, the relevant considerations are: (1) Rescaling a,b,c, avoiding overflow or underflow in d^2; (2) Double-precision square-root, avoiding cancellation in d^2; (3) Complex roots when d^2 < 0; (4) Evaluating x" = c/(a x') where |x'| > |x"|, avoiding cancellation in x"; (5) Detection of cases a = 0, a = b = 0, a = b = c = 0 (solutions single infinite, double infinite, undefined).

The following detailed discussions are freely accessible online: George E. Forsythe "How Do You Solve a Quadratic Equation", Tech. Report CS40, Comp. Sci. Dept., U. of Stanford (1966); Henry G. Baker "Garbage In, Garbage Out" ACM SIGPLAN Notices vol 33, 30--38 (1998).

The topic was aired in math-fun a while back ("branchless programming", 22-3/10/2007), but for some reason I didn't at the time get around to contributing my own two-penn'orth. In a geometric context, rather than as a single variable, the solution is often interpreted as a ratio, or homogeneous coordinate (x,y) of some point in projective 1-space.

From this aspect, it's irritating that the traditional formulation ties itself in logical knots over infinite roots (or worse, ignores them altogether); and furthermore it's aesthetically distressing that the symmetry of the homogeneous quadratic form has been lost.

But suppose either inhomogeneous solution is written alternatively as x = (b + d)/a = c/(b - d) by symmetry; now adding numerators and denominators leads to the symmetric reformulation:

The homogeneous quadratic equation a x^2 - 2 b x y + c y^2 = 0 has solutions x : y = (c + b + d) : (a + b - d) , (c + b - d) : (a + b + d) , where d^2 = b^2 - a c .

In this form there is no need to detect and specially treat cases (4) and (5). Where parallel processing is available therefore, it should execute marginally faster than the conventional formula (though of course both versions are dominated by the square root).

Fred Lunnon [18/05/10]