A tetrahedron has four triangular faces. Suppose three of those faces come together like the corner of a cube, each perpendicular to the other. Let A1, A2, and A3 be the areas of the three triangles that meet at this corner and let A0 be the area of remaining face, the one opposite the right corner.

De Gua’s theorem, published in 1783, says

A0² = A1² + A2² + A3².

We will illustrate De Gua’s theorem using Python.

from numpy import *

def area(p1, p2, p3):
return 0.5*abs(linalg.norm(cross(p1 - p3, p2 - p3)))

p0 = array([ 0, 0,  0])
p1 = array([11, 0,  0])
p2 = array([ 0, 3,  0])
p3 = array([ 0, 0, 25])

A0 = area(p1, p2, p3)
A1 = area(p0, p2, p3)
A2 = area(p0, p1, p3)
A3 = area(p0, p1, p2)

print(A0**2 - A1**2 - A2**2 - A3**2)

This prints 0, as expected.

Higher dimensions

A natural question is whether De Gua’s theorem generalizes to higher dimensions, and indeed it does.

Suppose you have a 4-simplex where one corner fits in the corner of a hypercube. The 4-simplex has 5 vertices. If you leave out any vertex, the remaining 4 points determine a tetrahedron. Let V0 be the volume of the tetrahedron formed by leaving out the vertex at the corner of the hypercube, and let V1, V2, V3, and V3 be the volumes of the tetrahedra formed by dropping out each of the other vertices. Then

V0² = V1² + V2² + V3² + V4².

You can extend the theorem to even higher dimensions analogously.

Let’s illustrate the theorem for the 4-simplex in Python. The volume of a tetrahedron can be computed as

V = det(G)1/2/6

where G is the Gram matrix computed in the code below.

def volume(p1, p2, p3, p4):
u = p1 - p4
v = p2 - p4
w = p3 - p4

# construct the Gram matrix
G = array( [
[dot(u, u), dot(u, v), dot(u, w)],
[dot(v, u), dot(v, v), dot(v, w)],
[dot(w, u), dot(w, v), dot(w, w)] ])

return sqrt(linalg.det(G))/6

p0 = array([ 0, 0,  0, 0])
p1 = array([11, 0,  0, 0])
p2 = array([ 0, 3,  0, 0])
p3 = array([ 0, 0, 25, 0])
p4 = array([ 0, 0,  0, 7])

V0 = volume(p1, p2, p3, p4)
V1 = volume(p0, p2, p3, p4)
V2 = volume(p0, p1, p3, p4)
V3 = volume(p0, p1, p2, p4)
V4 = volume(p0, p1, p2, p3)

print(V0**2 - V1**2 - V2**2 - V3**2 - V4**2)

This prints -9.458744898438454e-11. The result is 0, modulo the limitations of floating point arithmetic.

Numerical analysis

Floating point arithmetic is generally accurate to 15 decimal places. Why is the numerical error above relatively large? The loss of precision is the usual suspect: subtracting nearly equal numbers. We have

V0 = 130978.7777777778

and

V1**2 + V2**2 + V3**2 + V4**2 = 130978.7777777779

Both results are correct to 16 decimal places, but when we subtract them we lose all precision.