Plane geometry with Julia¶

Photo by Kaboompics .com from Pexels

When I first got interested in computer algebra systems, I came across a book

PLANE  GEOMETRY:   AN  ELEMENTARY  TEXTBOOK 
BY   SHALOSH  B.  EKHAD,  XIV   (CIRCA  2050)    
DOWNLOADED  FROM  THE  FUTURE  BY   DORON  ZEILBERGER

Of course, Shalosh is not really a time traveler from a future. He/she/it is the computer of Doron Zeilberger, a mathematician who has been advocating for the use of computers in mathematics for decades. He often writes articles and papers with Shalosh named as a coauthor.

The title his book is just a joke. What Zeilberger really wants to say is that in the future, kids won't need to learn to do (at least) plane geometry with pencil and paper anymore. Their homework will be writing codes so their computer 🤖️ will do the math for them.

The book was written with Maple. But in principle, this can be done in any programming languages. With SymPy.jl handling the symbolic computation and Plots.jl drawing nice pictures, Julia has become a good choice for this task. Writing the code in a modern language also makes things a bit more organized. Moreover, Julia is free but Maple costs 🤑️.

So I am going to show you how to repeat Zeilberger's proof of Napoleon's theorem using Julia. According to Wikipedia,

In geometry, Napoleon's theorem states that if equilateral triangles are constructed on the sides of any triangle, either all outward or all inward, the lines connecting the centres of those equilateral triangles themselves form an equilateral triangle.

The code in this notebook has been collected into a package PlaneGeometry.jl

Install needed packages¶

We will need SymPy.jl for the computation and Plot.jl for some nice plots. Run the following cell if you don't already have them.

In [ ]:

using Pkg; 
Pkg.activate("."); 
Pkg.add("SymPy")
Pkg.add("Plots")

Define some basic geometry objects¶

In [1]:

import Base: isequal, ==

abstract type GeoObject end
abstract type GeoShape <: GeoObject end

struct Point <: GeoObject
    x::Number
    y::Number
end
(==)(p1::Point, p2::Point) = p1.x==p2.x && p1.y==p2.y

struct Triangle <: GeoShape
    A::Point
    B::Point
    C::Point
end
Triangle(ax, ay, bx, by, cx, cy) = Triangle(Point(ax, ay), Point(bx, by), Point(cx, cy))
(==)(t1::Triangle, t2::Triangle) = vertices(t1) == vertices(t2)

function vertices(tri::Triangle) 
    [tri.A, tri.B, tri.C]
end

struct Edge <: GeoShape
    src::Point
    dst::Point
end

function edges(tri::Triangle)
    elist = Edge[]
    pts = vertices(tri)
    for i in 1:length(pts)-1
        push!(elist, Edge(pts[i], pts[i+1]))
    end
    push!(elist, Edge(pts[length(pts)], pts[1]))
    elist
end

struct Circle
    center::Point
    radius::Number
end
(==)(c1::Circle, c2::Circle) = c1.center == c2.center && c1.radius == c2.radius

ccenter(c::Circle) = c.center

Out[1]:

ccenter (generic function with 1 method)

How to draw a triangle¶

Let's consider the following triangle.

In [2]:

A = Point(0,0); B = Point(1, 3); C = Point(4,2);
tri = Triangle(A, B, C)

Out[2]:

Triangle(Point(0, 0), Point(1, 3), Point(4, 2))

To draw it, we create a Plots.Shape object.

In [3]:

using Plots

function shape(ptlist::Vector{Point})
    xlist = [pt.x for pt in ptlist]
    ylist = [pt.y for pt in ptlist]
    shape = Shape(xlist, ylist)
end;
shape(tri::Triangle) = shape(vertices(tri));

trishape=shape(tri)

Out[3]:

Shape([0.0, 1.0, 4.0], [0.0, 3.0, 2.0])

Then we can just feed the shape into plot(). The leg=false argument hides the unnecessary plot legend.

In [4]:

plot(trishape, leg=false, fill=(0, :green), aspect_ratio=:equal, fillalpha= 0.2)

Out[4]:

Let's point out where are the 3 points $A$, $B$ and $C$.

In [5]:

scatter!(trishape.x, trishape.y, color=:red, series_annotations = text.(["A", "B", "C"], :bottom))

Out[5]:

Equilateral triangles¶

There are two points that are can form a equilateral triangle with $A$ and $B$. Let's find them.

First we define a functions that computes the squared euclidean distance between two points.

In [6]:

function squaredist(A, B)
    (A.x-B.x)^2+(A.y-B.y)^2
end;

To find these two points, we need to solve a quadratic equation. For this we use SymPy.jl

In [7]:

using SymPy

function equipoints(A, B)
    x, y = @vars x y
    pt = Point(x, y)
    dist1 = squaredist(pt, A)
    dist2 = squaredist(pt, B)
    dist3 = squaredist(A, B)
    sol = solve([dist1-dist3, dist2-dist3], [x,y]) # Soving a quadratic eqaution.
    [Point(s[1], s[2]) for s in sol]
end;

These are the two points that can form equilateral triangles with $A$ and $B$.

In [8]:

ptAB = equipoints(A, B)

Out[8]:

2-element Array{Point,1}:
 Point(1/2 - 3*sqrt(3)/2, sqrt(3)/2 + 3/2)
 Point(1/2 + 3*sqrt(3)/2, 3/2 - sqrt(3)/2)

Let's consider only the outer equilateral triangle. So among the two points we only keep the one that is closer to $C$. The following function does that for us.

In [9]:

function outer_equitri(A, B, C)
    ptAB = equipoints(A, B)
    dist = map(pt->squaredist(pt, C), ptAB)
    if dist[1] >= dist[2]
        return Triangle(A, B, ptAB[1])
    else
        return Triangle(A, B, ptAB[2])
    end
end;

And we write function to do this for all the three edges of $\Delta ABC$.

In [10]:

function outer_equitriangles(tri)
    pts = vertices(tri)
    map(i->outer_equitri(circshift(pts, i)...), 0:2)
end

Out[10]:

outer_equitriangles (generic function with 1 method)

These are the three outer equilateral triangles for $\Delta ABC$.

In [11]:

outer_tri = outer_equitriangles(tri)

Out[11]:

3-element Array{Triangle,1}:
 Triangle(Point(0, 0), Point(1, 3), Point(1/2 - 3*sqrt(3)/2, sqrt(3)/2 + 3/2))
 Triangle(Point(4, 2), Point(0, 0), Point(sqrt(3) + 2, 1 - 2*sqrt(3)))
 Triangle(Point(1, 3), Point(4, 2), Point(sqrt(3)/2 + 5/2, 5/2 + 3*sqrt(3)/2))

Let's draw these three triangles to see if we get it correctly.

In [12]:

for t in outer_tri
    plot!(shape(t), leg=false, fill=(0, :pink), fillalpha=0.7)
end
current()  # Show the current plot

Out[12]:

Looks fine!

Circumcenter¶

Now we can compute the orthocenter) of the three triangles.

First a function that computes the circle that goes through several points.

In [13]:

function circumcircle(points)
    x, y = @vars x y
    c = Point(x, y)
    dist = [squaredist(pt, c) for pt in points]
    equations = [Eq(dist[i], dist[i+1]) for i in 1:length(points)-1]
    sol = solve(equations, [x, y])
    center = Point(simplify(sol[x]), simplify(sol[y]))
    Circle(center, sqrt(squaredist(center, points[1])))
end;

So the circumcircles of the three outer equilateral triangles are

In [14]:

outer_circles = circumcircle.(vertices.(outer_tri))

Out[14]:

3-element Array{Circle,1}:
 Circle(Point(1/2 - sqrt(3)/2, sqrt(3)/6 + 3/2), sqrt((1/2 - sqrt(3)/2)^2 + (sqrt(3)/6 + 3/2)^2))
 Circle(Point(sqrt(3)/3 + 2, 1 - 2*sqrt(3)/3), sqrt((-2 + sqrt(3)/3)^2 + (-2*sqrt(3)/3 - 1)^2))
 Circle(Point(sqrt(3)/6 + 5/2, sqrt(3)/2 + 5/2), sqrt((-1/2 + sqrt(3)/2)^2 + (sqrt(3)/6 + 3/2)^2))

Let's draw them to see if we computed them correctly. To draw a circle, again we need to create a Plot.Shape object to represent it.

In [15]:

function shape(c::Circle)
    θ = LinRange(0, 2*π, 300)
    c.center.x .+ c.radius*sin.(θ), c.center.y .+ c.radius*cos.(θ)
end;

In [16]:

for c in outer_circles
    plot!(shape(c), leg=false, fill=(0, :orange), aspect_ratio=:equal, fillalpha=0.2)
end
current() # Show current plot

Out[16]:

Now let's draw the triangle formed by the center of these 3 circles.

In [17]:

outer_centers = ccenter.(outer_circles);
napoleon_tri = Triangle(outer_centers...);
plot!(shape(napoleon_tri), leg=false, fill=(0, :blue), fillalpha=0.7, aspect_ratio=:equal)

Out[17]:

A surprise¶

According to the theorem, this blue triangle is actually equilateral? Is it? A function to check it --

In [18]:

function isequilateral(tri)
    dist = [squaredist(e.src, e.dst) for e in edges(tri)]
    if (dist[1] == dist[2]) && (dist[1] == dist[3])
        return true
    else
        return false
    end
end

Out[18]:

isequilateral (generic function with 1 method)

In [19]:

isequilateral(napoleon_tri)

Out[19]:

false

Wow! The theorem is wrong! 😱️ Of course not. Math tells us the theorem is correct. What is wrong is our computation. Let's check the lengths of each edge of napoleon_tri.

In [20]:

dist1 = squaredist(napoleon_tri.A, napoleon_tri.B)

Out[20]:

\begin{equation*}\left(\frac{1}{2} + \frac{5 \sqrt{3}}{6}\right)^{2} + \left(- \frac{3}{2} - \frac{5 \sqrt{3}}{6}\right)^{2}\end{equation*}

In [21]:

dist2 = squaredist(napoleon_tri.B, napoleon_tri.C)

Out[21]:

\begin{equation*}\left(- \frac{1}{2} + \frac{\sqrt{3}}{6}\right)^{2} + \left(- \frac{7 \sqrt{3}}{6} - \frac{3}{2}\right)^{2}\end{equation*}

Let's compare them numerically.

In [22]:

N(dist1) == N(dist2)

Out[22]:

true

The problem is, by default, SymPy expressions are not simplified. So to check these two distances are actually equal, we need to do some simplifications. So when we compute distances, we should use SymPy.simplify() on the result.

In [23]:

function squaredist(A, B)
    d = (A.x-B.x)^2+(A.y-B.y)^2
    if d isa Sym
        simplify(d)
    else
        d
    end
end;

No we can show that the theorem is at least correct in this example. 😀️

In [24]:

isequilateral(napoleon_tri)

Out[24]:

true

A coincidence?¶

Maybe we are just lucky and chose a triangle that the theorem holds. To be sure, we can try some more triangles.

First we put how we draw the picture above in a function.

In [25]:

function napoleon_draw(xA, yA, xB, yB, xC, yC)
    A = Point(xA, yA); 
    B = Point(xB, yB); 
    C = Point(xC, yC);
    tri = Triangle(A, B, C)

    trishape=shape(tri)

    plot(trishape, leg=false, fill=(0, :green), aspect_ratio=:equal, fillalpha= 0.2)

    scatter!(trishape.x, trishape.y, color=:red, series_annotations = text.(["A", "B", "C"], :bottom))

    outer_tri = outer_equitriangles(tri)
    for t in outer_tri
        plot!(shape(t), leg=false, fill=(0, :pink), aspect_ratio=:equal, fillalpha=0.7)
    end

    outer_circles = circumcircle.(vertices.(outer_tri))

    outer_centers = ccenter.(outer_circles);

    napoleon_tri = Triangle(outer_centers...);

    plt = plot!(shape(napoleon_tri), leg=false, fill=(0, :blue), fillalpha=0.7, aspect_ratio=:equal)

    hold = isequilateral(napoleon_tri)

    plt, hold
end

Out[25]:

napoleon_draw (generic function with 1 method)

Let's pick a random triangle.

In [26]:

function napoleon_rand()
    pts = rand(0:1//10:1,6);
    plt, hold = napoleon_draw(pts...)
    if(hold == true)
        println("Theorem holds! 😀️")
    else
        println("Theorem does not holds! 😱️")
    end
    plt
end

Out[26]:

napoleon_rand (generic function with 1 method)

In [27]:

napoleon_rand()

Theorem holds! 😀️

Out[27]:

Let's try again.

In [28]:

napoleon_rand()

Theorem holds! 😀️

Out[28]:

A symbolic proof¶

Of course, examples are not proofs. But what are proofs? Often a proof is just a computation done with symbols instead of fixed numbers.

Note that symbols, e.g., $x$, is different from Julia's variables. It does not hold any value but is just a placeholder in a computation. This means if a symbolic computation with $x$ is valid, then then the computation still holds if we replace $x$ with any number.

First, note that we can assume that $A = (0,0)$, $B$ is on the positive half of $y$-axis, and $C$ is in the right half of the plane. In other words, something like this.

In [29]:

plt, ret = napoleon_draw(0, 0, 0, 3, 3, 1)
plt

Out[29]:

We can always move a triangle so it satisfies these conditions.

Let represent their coordinates by symbols.

In [30]:

@vars by cx positive=true;
@vars cy;

A = Point(0, 0); B = Point(0, by); C = Point(cx, cy);
tri = Triangle(A, B, C)

Out[30]:

Triangle(Point(0, 0), Point(0, by), Point(cx, cy))

All the computations done before, we can just copy and paste!

In [31]:

function napoleon_check(tri)
    outer_tri = outer_equitriangles(tri);
    outer_circles = circumcircle.(vertices.(outer_tri));
    outer_centers = ccenter.(outer_circles);
    npt = Triangle(outer_centers...)
    isequilateral(npt)
end

napoleon_check(tri)

Out[31]:

false

We get a false. But again, this is a programming error. The problem is, for SymPy to know which of two symbolic expressions is larger, it needs a bit help. So when we choose among the two points which one is for the outer equilateral triangle, we chose the wrong one. We can again fix this by using simplify().

In [32]:

function outer_equitri(A, B, C)
    ptAB = equipoints(A, B)
    dist = map(pt->squaredist(pt, C), ptAB)
    d = simplify(dist[1] - dist[2])
    if d >= 0
        return Triangle(A, B, ptAB[1])
    else
        return Triangle(A, B, ptAB[2])
    end
end;

This time it works!

In [33]:

napoleon_check(tri)

Out[33]:

true

So, Napoleon is right! 😀️ But do you think future kids will actually prove this theorem like this in class? 🤔️