Complex Numbers in R, Part I


R, like many scientific programming languages, has first-class support for complex numbers. And, just as in most other programming languages, this functionality is ignored by the vast majority of users.

Yet complex numbers can often offer surprisingly elegant formulations and solutions to problems. I want to convince you that familiarizing yourself with R’s excellent complex number functionality is well worth the effort and will pay off tenfold in two different ways: first by showing you how they are so amazingly useful you’ll want to go out of your way to use them, and then by showing you how they are so common and fundamental to modern analysis that you couldn’t avoid them if you wanted to.

Pythagorean Triples

Let’s start with a problem which could be solved in other ways, but is greatly simplified by the introduction of complex numbers that it almost seems magical.

A Pythagorean triple is an integer solution to the Pythagorean equation:

\[ a^2 + b^2 = c^2 \quad\quad a,b,c \in \mathbb{N}^+ \tag{1} \]

You probably learned at least one of these in school – the famous 3, 4, 5 triangle:

3-4-5 triangle

In general Diophantine equations – which require integer solutions – can be quite hard to solve, so it might surprise you to hear that it’s almost trivially easy to write down an infinite of Pythagorean triples. Well, it’s easy if we use complex numbers, anyway.

A Gaussian integer is a complex number where both the real and imaginary parts are integers:

\[ Z[i] = \{ x + iv \mid x,y \in \mathbb{Z} \} \tag{2} \]

So one way of stating the problem of finding all Pythagorean triples is to find all Gaussian integers Which are an integer distance away from the origin. The distance of a complex number from the origin is called its “norm” and denoted \(\lVert z \rVert\).

\[ T = \{ z \in \mathbb{C} \mid \lVert z \rVert \in \mathbb{Z} \} \tag{3} \]

Now, in general the norm of Gaussian integer will be the square root of an integer (the integer \(u^2 + v^2\) to be precise.) Therefore if we square a Gaussian integer, it will have an integer norm and therefore represent a Pythagorean triple!

\[ \forall z \in Z[i], z^2 \in T \tag{4} \]

So that’s a pretty good start: just a few minutes work, and we’ve already found an infinite number of Pythagorean triples, and we have a computationally trivial way of constructing new triples: we simply pick any two positive integers \(u\) and \(v\) and then square the complex number \(u + iv\).

Before address the more difficult question of whether or not we’ve found all possible Pythagorean triples using this construction, let’s switch over to R and do some code our solution so far.

Our algorithm first requires us to pick pairs of positive integers. Just to be thorough, we’ll take all such pairs up to an arbitrary threshold.

Now, if we wanted just one or two complex numbers, we could use the literal syntax:

triples <- c( 3+4i, 5+12i, 9+12i )

But since we want them along a grid, we’ll use the complex() constructor. The constructor is vectorized so by passing in two vectors of equal length we can construct a one-dimensional vector of complex numbers.

n = 400
grid <- expand.grid(u=1:(2*n), v=1:(2*n))
grid <- grid[ grid$u > grid$v, ]
gaussian_integers <- complex(real=grid$u, imaginary=grid$v)

Per the theoretical discussion above, we can generate Pythagorean triples by simply squaring these. All primitive math functions in R work just as well on complex numbers: exp, log, sin, cosand of course exponentiation:

triples <- gaussian_integers^2

# display the 10 with the smallest norm
cat( triples[order(Mod(triples))][1:10], sep="\n")
# 3+4i
# 8+6i
# 5+12i
# 15+8i
# 12+16i
# 7+24i
# 24+10i
# 21+20i
# 16+30i
# 35+12i

Did it work? We’re certainly seeing some familiar pairings there, like \(5+12i\) which maps to well-known triple \((5,12,13)\). To visualize them, we can simply pass our complex vector to R’s plot() function - it will conveniently plot them in the complex plane for us!

triples <- triples[ Re(triples) <= n & Im(triples) <= n ]

# helper function to colorize complex points by their angle.
argcolor <- function(z) hsv(Arg(z)*2/pi, s=0.9, v=0.8)

plot(
  triples, 
  col=argcolor(triples),
  pch=20,
  xlim=c(0,n),
  ylim=c(0,n),
  main=paste("Squared Gaussian Integers Up to", n)
)

Now it turns out that our algorithm does not, in fact, generate all possible triples. One thing it misses are multiples: if \((3,4,5)\) is a triple, then \((6,8,10)\) should be a triple, and \((9,12,15)\) should be a triple, and so on. So we have to expand our set to have all multiples.

multiples <- lapply(1:(floor(n/3)), function(m) triples*m)
triples <- unique(do.call(c, multiples))

It also turns out that in the special case where both integers are even we can divide by two and get a new triple that was missed by the initial net we cast. But that’s the end of the special cases – with this final rule in place, we’re now guaranteed to hit every Pythagorean triple.

halves <- triples[ Re(triples) %% 2 == 0 & Im(triples) %% 2 == 0 ] / 2
triples <- unique(c(triples, halves))

Now all we need to is clean up duplicates and duplicate along the mirror line of symmetry…

triples <- triples[ Re(triples) <= n & Im(triples) <= n]
triples <- c(triples, complex(real=Im(triples), imaginary=Re(triples)))

..and we’re finally ready to visualize the real solution.

plot(triples, col=argcolor(triples), pch=20)
title(paste("All Pythagorean Triples Up to", n))

That’s too many to really understand, although there are definitely patterns emerging. Let’s zoom in and just plot a small region,but with more detail.

small_n = 25
small_triples <- triples[ Re(triples) < small_n & Im(triples) < small_n ]
small_triples <- small_triples[ order(Mod(small_triples), decreasing=TRUE) ]

# plot points
plot(
  small_triples, 
  pch=20,
  ylim=c(0,small_n), 
  xlim=c(0,small_n),
  ylab="b", xlab="a")

# add triangles. Can't rely on automatic complex plane plotting here.
segments(
  Re(small_triples), Im(small_triples), 
  0, 0, 
  col=argcolor(small_triples))
segments(
  Re(small_triples), Im(small_triples), 
  Re(small_triples), 0, 
  col=argcolor(small_triples))
segments(
  Re(small_triples), 0, 
  0, 0, 
  col=argcolor(small_triples))

# points again, so that they're in the foreground.
points(small_triples, pch=20, col=argcolor(triples), cex=1)

# text label for the points
text(
  x=small_triples + 1i, 
  cex=0.8,
  labels=paste0(
    "(", 
    Re(small_triples), 
    ",",
    Im(small_triples),
    ",",
    Mod(small_triples),
    ")"
  )
)
title(paste("Pythagorean Triples Up to", small_n))

On the zoomed in view we can see each Pythagorean triple represented as a right triangle; we can see that the integer multiples of solutions form a series of similar triangles; we can see that there’s a strong symmetry with every triple \((a,b,c)\) having a partner \((b,a,c)\) which is its mirror reflection about the like \(y=x\).

From the zoomed out view we can see that the region close to either axis is essentially devoid of solutions; that their radii (e.g. close to 50) where triples are dense, and others where they are rare to non-existent; and it looks as if triples actually get less dense as we move away from the origin.

This last observation, about triples thinning out as we move away from the origin, can be understood and quantified by once again using the complex plane. Triples are more or less the squares of Gaussian integers; we can say the number of triples with norm less than \(r\) is roughly proportional to the number of Gaussian integers in the first quadrant and inside a circle with radius \(\sqrt{r}\), which is roughly proportional to the area of the quarter-circle of radius \(\sqrt{r}\), which is \(\frac{\pi r}{4}\) or very roughly just \(r\).

Next Time

In this first part of a planned series on complex numbers in R, we dipped our toes in the water by explicitly creating some complex numbers and showing how they worked with the most basic mathematical operators, functions, and plots.

Next time in Part II, we will discuss in more depth at some scenarios where complex numbers arise naturally from the problem.