Heron’s formula computes the area of a triangle given the length of each side.

where

If you have a very thin triangle, one where two of the sides approximately equal *s* and the third side is much shorter, a direct implementation Heron’s formula may not be accurate. The cardinal rule of numerical programming is to avoid subtracting nearly equal numbers, and that’s exactly what Heron’s formula does if *s* is approximately equal to two of the sides, say *a* and *b*.

William Kahan’s formula is algebraically equivalent to Heron’s formula, but is more accurate in floating point arithmetic. His procedure is to first sort the sides in decreasing order, then compute

You can find this method, for example, in Nick Higham’s book Accuracy and Stability of Numerical Algorithms.

The algebraically redundant parentheses in the expression above are not numerically redundant. As we’ll demonstrate below, the method is less accurate without them.

Optimizing compilers respect the parentheses: the results are the same when the code below is compiled with gcc with no optimization (`-O0`

) and with aggressive optimization (`-O3`

). The same is true of Visual Studio in Debug and Release mode.

## C code demo

First, here is a straightforward implementation of Heron

#include <math.h> float heron1(float a, float b, float c) { float s = 0.5 * (a + b + c); return sqrt(s*(s - a)*(s - b)*(s - c)); }

And here’s an implementation of Kahan’s version.

void swap(float* a, float* b) { float t = *b; *b = *a; *a = t; } float heron2(float a, float b, float c) { // Sort a, b, c into descending order if (a < b) swap(&a, &b); if (a < c) swap(&a, &c); if (b < c) swap(&b, &c); float p = (a + (b + c))*(c - (a - b))*(c + (a - b))*(a + (b - c)); return 0.25*sqrt(p); }

Finally, here’s an **incorrect** implementation of Kahan’s method, with “unnecessary” parentheses removed.

float heron3(float a, float b, float c) { // Sort a, b, c into descending order if (a < b) swap(&a, &b); if (a < c) swap(&a, &c); if (b < c) swap(&b, &c); float p = (a + b + c)*(c - (a - b))*(c + a - b)*(a + b - c); return 0.25*sqrt(p); }

Now we call all three methods.

int main() { float a = 100000.1, b = 100000.2, c = 0.3; printf("%0.8g\n", heron1(a, b, c)); printf("%0.8g\n", heron2(a, b, c)); printf("%0.8g\n", heron3(a, b, c)); }

And for a gold standard, here is an implementation in `bc`

with 40 decimal place precision.

scale = 40 a = 10000000.01 b = 10000000.02 c = 0.03 s = 0.5*(a + b + c) sqrt(s*(s-a)*(s-b)*(s-c))

Here are the outputs of the various methods, in order of increasing accuracy.

Heron: 14363.129 Naive: 14059.268 Kahan: 14114.293 bc: 14142.157

Here “naive” means the incorrect implementation of Kahan’s method, `heron3`

above. The `bc`

result had many more decimals but was rounded to the same precision as the C results.

**Related post**: How to compute the area of a polygon

shouldn’t if (a < c) swap(&a, &c); be if (a < b) swap(&a, &c);

if (a

bif (b c and a > min(b,c)

// need to compare a with what ever is at b (could be previous c)

e.g a = 2, b = 1, c = 3

expr, your comment got a bit mangled, but yes, I think there’s an error in the code in the post. A working compare-and-swap order is:

1. a with b (a and b will be relatively ordered, c could be anything)

2. a with c (a will definitely be the largest after this step)

3. b with c (get these two correct and we’re done).

Although several other orders would also work.

Yes, I needed to swap the order of two of my swaps. :)

I updated the code. Thanks.

The bc implementation is different from the rest.

With Mathematica one gets the correct result: 14142.142

a = 1000001/10;

b = 1000002/10;

c = 3/10;

s = (a + b + c)/2;

NumberForm[N@Sqrt[a (s – a) (s – b) (s – c)], 9]