Skip to content
scipython.com

P10.1.1Heron's formula

Heron's formula for the area of a triangle (as used in Example E2.3):

\[A = \sqrt{s(s-a)(s-b)(s-c)} \;\mathrm{where}\;s = \textstyle\frac{1}{2}(a+b+c).\]

is inaccurate if one side is very much smaller than the other two ("needle-shaped" triangles). Why? Demonstrate that the following reformulation [1] gives a more accurate result in this case by considering the triangle with sides \((10^{-13}, 1, 1)\) which has the area \(5 \times 10^{-14}\).

\[A = \frac{1}{4}\sqrt{(a+(b+c))(c-(a-b))(c+(a-b))(a+(b-c))},\]

where the sides have been relabelled so that \(a \ge b \ge c\).

What happens if you rewrite the factors in this equation to remove their inner parentheses? Why?

References

  1. This formula is due to William Kahan, one of the designers of the IEEE-754 floating point standard.
Show solution

Consider the following Python implementations of Heron's formula:

heronprecision
import math


def heron1(a, b, c):
    """
    Simplisitic implementation of Heron's formula for the area of a triangle

    """
    s = (a + b + c) / 2
    return math.sqrt(s * (s - a) * (s - b) * (s - c))


def heron2(a, b, c):
    """
    Numerically stable implementation of Heron's formula for the area of a
    triangle -- better for long, thin triangles.

    """
    # Order the side lengths so that a >= b >= c
    a, b, c = sorted((a, b, c), reverse=True)
    return (
        math.sqrt(
            (a + (b + c)) * (c - (a - b)) * (c + (a - b)) * (a + (b - c))
        )
        / 4
    )


if __name__ == "__main__":
    a, b, c = 1, 1, 1.0e-13
    print(
        f"Simple: {heron1(a, b, c)}, Stable: {heron2(a, b, c)}, Exact: 5.e-14"
    )

For the suggested triangle, the second of these functions clearly works better, and gives the area, \(5 \times 10^{-14}\), correct to machine precision whereas the simplistic implementation is out by about 0.08%.

In [x]: a, b, c = 1, 1, 1.e-13
In [x]: print('Simple: {}, Stable: {}, Exact: 5.e-14'
                                .format(heron1(a,b,c), heron2(a,b,c)))
Simple: 4.9960036108132044e-14, Stable: 5e-14, Exact: 5.e-14

The problem with the simplistic calculation is that the semiperimeter, \(s\), is very close in value to the longer sides' lengths leading to loss of precision in the subtractions. The second implementation avoids this by not subtracting a small number from a very much larger one. Removing the parentheses breaks the algorithm since addition and subtraction are evaluated left-to-right.