A Practical Guide to Periodic Boundary Conditions

This article is aimed at programmers, physicists, mathematicians, and anyone else who needs to work with periodic spaces – whether those are angles, particle systems, colors, musical notes or anything else with a periodic structure.

dx = x2 - x1

if dx < -w/2:
    dx += w
else if dx >= w/2:
    dx -= w

Initially, working with periodic spaces seems easy. Just use the modulo operator, right? (By the way: I will advise against using the modulo operator – there are better alternatives.) But as soon as the question of "shortest connections" across boundaries arises, things get confusing: What's the smallest distance between two values? Should I go across the left or the right border? What about multiple dimensions?

This article will introduce you to a helpful framework of concepts that will allow you to answer those questions with ease. In this process, I will show you different algorithms for implementing periodic boundaries and shortest connections in your code.

The first two parts Periodic Boundaries and Connections will look at one-dimensional spaces, while the third part Multidimensional will show you how to generalize these concepts to higher dimensions.

If you spot an error in this article, please let me know! Cite this article if you plan to use any of this in your own publications.

Cite this article
BibTeX:
@online{mohr2024,
    author = {Tom Mohr},
    year = {2024},
    title = {A Practical Guide to Periodic Boundary Conditions},
    url = {https://tommohr.dev/pbc}
}
Table of Contents

Part I

Periodic Boundaries

This section introduces basic terms needed for talking about periodic boundaries.

Interval $[a, b)$

A periodic space is defined by a half-open interval $[a, b)$ with lower bound $a$ and upper bound $b$. It includes all real numbers $x$ with $a \leq x < b$. The upper bound $b$ is not part of the interval.

For angles, the interval could be $[0, 2\pi)$, $[-\pi, \pi)$ or $[0, 360^{\circ})$. For musical notes, the interval could be $[0, 12)$.

Interval Size \(w\)

The value $w = b - a$ is called the size of the interval $[a, b)$.

For angles, the interval size would be $w = 2 \pi$ or $w=360^{\circ}$. For a particle system, this could be the width or height of your canvas. For musical notes, this could be $w=12$ (an octave).

Instead of defining an interval via $a$ and $b$, we can also define it via $a$ and $w$, writing the interval as $[a, a + w)$. This is just a matter of taste.

Note that if we shift such an interval by its size, we get a new interval $$[a, b) + w = [a + w, b + w)$$ which is disjoint from the original one: $$ [a, b) + w \hspace{0.5em} \cap \hspace{0.5em} [a, b) = \emptyset $$

This is generally true for any non-zero multiple of $w$: $$ [a, b) + kw \hspace{0.5em} \cap \hspace{0.5em} [a, b) = \emptyset $$ (for any $k \in \mathbb{Z}$ with $k \neq 0$).

Proof

Proof by contradiction.
Let's assume that there is an intersection: $$[a, b) + kw \cap [a, b) \neq \emptyset.$$

Then there must be a number $x$ that is in both intervals: $$\begin{align*} &x \in [a, b) \\ &x \in [a, b) + kw \end{align*}$$ Rewritten in terms of $a$ and $w$: $$\begin{align*} &x \in [a, a + w) \\ &x \in [a + kw, a + kw + w) \end{align*}$$ Rewritten as a collection of inequalities: $$\begin{align*} a &\leq x \\ x &< a + w \\ a + kw &\leq x \\ x &< a + kw + w \end{align*}$$ Subtract $a$ from each inequality: $$\begin{align*} 0 &\leq x - a \\ x - a &< w \\ kw &\leq x - a \\ x - a &< kw + w \end{align*}$$ Combining 1 and 4, as well as 2 and 3 yields: $$\begin{align*} -w &< kw \\ kw &< w \end{align*}$$ Combine again: $$-w < kw < w$$ Divide by $w$: $$-1 < k < 1$$ This is only possible if $k = 0$, which contradicts $k \neq 0$.

Winding Number \(n\)

The winding number \(n \in \mathbb{Z}\) of a number \(x \in \mathbb{R}\) is the number of times that the interval size \(w\) has to be subtracted from \(x\) in order to bring it into the range \([a, b)\).

The value of \(n\) can, for example, be computed via \(n = \left\lfloor \frac{x - a}{w} \right\rfloor\) (normalizing, then flooring).

Therefore, \( x - nw \in [a, b)\).

If [a, b) = [0, 1), then \(n \leq x < n+1\). So here, \(n = \lfloor x \rfloor\) (floor).

Absolute Winding Number \(|n|\)

The absolute winding number \(|n| \in \{0, 1, 2, \ldots\}\) is the absolute value of the winding number \(n\). It counts how many times the interval size \(w\) has to be added to or subtracted from \(x\) in order to bring it into the range \([a, b)\).

Periodic Boundary Operator \(p\)

The function $p_{a,b}$ is called the periodic boundary operator on $[a, b)$.

It is defined as follows:

$$\begin{align*} p_{a,b}&: \mathbb{R} \to [a, b) \\ p_{a,b}(x) &= x - \left\lfloor \frac{x - a}{b - a} \right\rfloor (b - a) \end{align*}$$

This is just a complicated way of writing $$p(x) = x - nw$$ where $n = \left\lfloor \frac{x - a}{w} \right\rfloor$ is the winding number of $x$ and $w = b - a$ is the interval size.

We often omit the subscripts $a, b$ and just write $p$ if the interval is clear from the context.

That last form leads us to a visual intuition of the periodic boundary operator: It unwinds the number $x$ back into the interval $[a, b)$. We could also call it the unwinding operator.

We will now discuss some often needed properties of the periodic boundary operator:

Identity on $[a, b)$: $p \mid_{[a, b)} \hspace{0.25em} = \mathtt{id}$
Proof

In other words: $p(x) = x$ for all $x \in [a, b)$.
Meaning, on the interval $[a, b)$, the periodic boundary operator $p$ is the identity.

If $x \in [a, b)$, we have a winding number of $n = 0$.
Therefore, $p(x) = x - nw = x - 0w = x$.

Idempotency: $p^2 = p$
Proof

In other words: $p(p(x)) = p(x)$. Applying the boundary operator twice is the same as applying it once.

First, remember that $p(x) \in [a, b)$.
Then, since $p \mid_{[a, b)} \hspace{0.25em} = \text{id},$ we get $p(p(x)) = \text{id}(p(x)) = p(x).$

Period $w$: $p(x + w) = p(x)$
Proof

The periodic boundary operator has a period of $w$. In other words: $p$ is immune to addition / subtraction of the interval size.

First, a quick reminder of the definition of $p$: $$p(x) = x - \left\lfloor \frac{x - a}{w} \right\rfloor w$$

We need to be aware of the fact that the floor function is compatible with addition and subtraction of integers. In other words: $$\lfloor x \rfloor + 1 = \lfloor x + 1 \rfloor$$ Thus, $\lfloor x \rfloor + m = \lfloor x + m \rfloor$ for any $m \in \mathbb{Z}$.

$$\begin{align*} p(x) &= x - \left\lfloor \frac{x - a}{w} \right\rfloor w \\ &= x + w - \left\lfloor \frac{x - a}{w} \right\rfloor w - w \\ &= x + w - (\left\lfloor \frac{x - a}{w} \right\rfloor + 1)w \\ &= x + w - \left\lfloor \frac{x - a}{w} + 1\right\rfloor w \\ &= x + w - \left\lfloor \frac{x + w - a}{w}\right\rfloor w \\ &= p(x + w) \end{align*}$$

Thus, $p(x + mw) = p(x) \hspace{1em}$ for any $m \in \mathbb{Z}$.

Boundary Operator Implementation

There are many possible algorithms that compute \(p(x)\). We can generally differentiate between three implementations:

Cost Restrictions
Single Wind \(O(1)\) \(|n| \leq 1\).
Linear \(O(|n|)\)
Constant \(O(1)\) but probably slow.

I will now discuss these implementations in detail, because each one has its own advantages and disadvantages.

1. Single Wind:

If you can guarantee that the absolute winding number \(|n|\) is never greater than 1, you can use this implementation. For example if

With time complexity \(O(1)\) and exclusive use of cheap operations, it is the fastest implementation:

if x < a:
    x += w
else if x >= b:
    x -= w

However, it can only handle absolute winding numbers \(|n| \leq 1\). This is therefore only a partial implementation of \(p\) on the interval \([a-w, b+w)\).

2. Linear

This implementation has linear time complexity \(O(|n|)\):

while x < a:
    x += w
while x >= b:
    x -= w

This implementation depends on the absolute winding number \(|n|\) being low. In real-world applications, this is often the case: Imagine a particle system, where the particles usually only move a few pixels per frame – this will almost never result in a high absolute winding number.

3. Constant

This implementation has constant time complexity \(O(1)\). This is good if you also want to handle high absolute winding numbers, in which case the Single Wind and Linear implementations would fail (produce wrong results or take very long).

For low absolute winding numbers, this implementation is usually slower than the Single Wind and Linear implementations. This is because it involves expensive operations such as division, floor or modulo.

However, even if you can guarantee that the absolute winding number is low, the Constant implementation might still be the best choice on some platforms. For example, highly parallelized targets like GPUs usually don't like programs that branch unpredictably – making the Constant implementation perform the best.

There are two ways to implement this. I recommend using a floor-based implementation $x - \lfloor \frac{x-a}{w} \rfloor w$.

x - floor((x - a) / w) * w

You could also use a modulo-based implementation $a + ((x-a) \mod w)$, although I advise against it.

a + ((x - a) % w)

In the appendix, I explain why I advise against using the modulo operator for this.

So, which one should you choose?

Nested Boundary Operators

The application of a boundary operator can be omitted if the result is part of a sum that is passed to another boundary operator of equal size.

In other words:
Let \(p_1\) and \(p_2\) be boundary operators in two (possibly different) periodic spaces with identical size $w$. Then, \[p_2(p_1(x) + y) = p_2(x + y).\]

Proof

First, a quick reminder of the definition of \(p\): $$p(x) = x - nw,$$ where $n \in \mathbb{Z}$ is the winding number of $x$.

When we introduced the periodic boundary operator, we showed that the addition / subtraction of $w$ inside $p$ can be omitted: $$p(x + mw) = p(x)$$ for any $m \in \mathbb{Z}$.

Now, the actual proof is simple: $$ \begin{align*} p_2(p_1(x) + y) &= p_2(x - n_1 w + y) \\ &= p_2(x + y - n_1 w) \\ &= p_2(x + y). \end{align*} $$

In the following, we will use the special case \(p_1=p_2=p\) a lot: \[p(p(x) + y) = p(x + y)\]

We can also derive the idempotency of the boundary operator from this by setting \(y = 0\): \[p(p(x)) = p(x).\]

Part II

Connections

In this part, we will further explore the mathematical properties of periodic space to derive a definition of a connection in periodic space.

Periodic Addition $\oplus$

We define the addition \(\oplus\) on \([a, b)\) as \[x \oplus y = p(x + y).\]

This yields a commutative group:

Associativity: \(x \oplus y \oplus z\)
Proof
Uses the redundancy of nested boundary operators.
\[ \begin{align*} x \oplus (y \oplus z) &= x (x \oplus y) \oplus z \\ \Leftrightarrow p(x + p(y + z)) &= p(p(x + y) + z) \\ \Leftrightarrow p(x + y + z) &= p(x + y + z) \end{align*} \]
Commutativity: \(x \oplus y= y \oplus x\)
Proof
Trivial because real addition is commutative.
\(x \oplus y = p(x + y) = p(y + x) = y \oplus x\).
Neutral Element: \(p(0)\)
Proof
Uses the redundancy of nested boundary operators and the fact that \(p\) on \([a, b)\) is the identity.
\(x \oplus p(0) = p(x + p(0)) = p(x + 0) = p(x) = x\).
Inverse: \(p(-x)\)
Proof
Uses the redundancy of nested boundary operators.
\(x \oplus p(-x) = p(x + p(-x)) = p(x - x) = p(0)\).
If \(0 \in [a, b)\), then the neutral element is \(0\). On \([0, 100)\):
\(60 \oplus 53 = p(60 + 53) = p(113) = 13\).
The result of the addition remains in the interval.

Periodic Subtraction $\ominus$

What's the meaning of subtraction in periodic space? Generally, subtraction is defined as the addition with the inverse. So let's define subtraction as \(x \ominus y = x \oplus p(-y)\). We can show that this is equivalent to \[x \ominus y = p(x - y).\]

Proof
Uses the redundancy of nested boundary operators. \[ \begin{align*} &x \oplus p(-y) && \\ =& p(x + p(-y)) && \text{(nested boundary operators)} \\ =& p(x + (-y)). && \\ =& p(x - y). && \end{align*} \]

Differences in an Interval

Before we continue, let's make a simple observation:

The difference of two numbers in an interval is always smaller than the size of that interval.

In other words: For any numbers \(x, y \in [a, b)\): \[|x - y| < w\]

In other words: \[x - y \in (-w, w)\]

It's good to have this written down for the following sections.

Symmetrical Periodic Spaces $[-\frac{w}{2}, \frac{w}{2})$

Among all periodic spaces $[a, b)$, there is a nice special case: Symmetrical periodic spaces. A symmetrical periodic space is a periodic space $[a, b)$ with $-a = b$. It can always be written in the form $[-\frac{w}{2}, \frac{w}{2})$, where $w$ is the size of the space.

Proof

$[a, b)$ being symmetrical simply means that $[a, b) = [-b, b)$, so $-a = b$.
Thus, $b = \frac{2b}{2} = \frac{b+b}{2} = \frac{b-a}{2} = \frac{w}{2}$ .

Now, we can rewrite $[a, b) = [-b, b) = [-\frac{w}{2}, \frac{w}{2})$.

We can express a symmetrical boundary operator more concisely using the $\mathtt{round}()$ operation: $$p(x) = x - \mathtt{round}(\frac{x}{w}) w$$

Proof
Because generally $\mathtt{round}(x) = \lfloor x + 0.5 \rfloor$, we can rewrite the symmetrical boundary operator on $[a, b) = [-\frac{w}{2}, \frac{w}{2})$ as follows: $$ \begin{align*} p(x) &= x - \left\lfloor \frac{x - a}{w} \right\rfloor w \\ &= x - \left\lfloor \frac{x + \frac{w}{2}}{w} \right\rfloor w \\ &= x - \left\lfloor \frac{x}{w} + \frac{1}{2} \right\rfloor w \\ &= x - \mathtt{round}(\frac{x}{w}) w \end{align*} $$

Furthermore, symmetrical periodic spaces have several nice mathematical properties:

  1. The neutral element is 0.

    Proof
    As already established in the section about periodic addition, the neutral element is always \(p(0)\).
    Because the space is symmetrical, here \(0 \in [a, b)\).
    Since the periodic boundary operator acts as the identity on \([a, b)\), the neutral element is \(p(0) = 0\).
  2. Inside \((a, b)\), the inverse of a number \(x\) is its real negative \(-x\).

    Proof

    It's generally true that \(-x\) is an inverse of \(x\) because \(x \oplus -x = p(x + (-x)) = p(x - x) = p(0)\). That is if we allow the inverse to be outside the interval \([a, b)\).
    However, we might be interested in inverses only from the interval \([a, b)\). We already established in the section about periodic addition that such an inverse of \(x\) is equal to \(p(-x)\). (Also note that this is the only inverse of \(x\) which lies in the interval \([a, b)\).)

    Now consider that in the symmetrical case, if \(x \in (a, b)\), then also \(-x \in (a, b)\).
    This is because the symmetrical space is centered around 0 (i.e. \(-a=b\)):

    1. If \(x > a\), then \(-x < -a = b\).
    2. If \(x < b\), then \(-x > -b = a\).

    So if \(x \in (a, b)\), then both conditions are true and therefore \(-x \in (a, b)\).

    Also, the boundary operator generally (not just in the symmetrical case) acts as identity on \((a, b)\): \[p\mid_{(a, b)}(x) = x.\] Therefore, \(p(-x) = -x\) for any \(x \in (a, b)\), because \(-x \in (a, b)\).

These properties make symmetrical spaces mathematically very similar to the real numbers $\mathbb{R}$. After the following two sections about connections, we will see that there are even more advantages to symmetrical spaces.

Connections $C(x, y)$

A connection between two values \(x, y\) is a value \(c\) such that \[x \oplus c = y.\] The meaning of "connection" therefore completely depends on the definition of \(\oplus\) in the space that you are working in.

Depending on the definition of \(\oplus\), there can be multiple connections between two values: \[C(x, y) = \{c \mid x \oplus c = y\}\]

In infinite space, there is only one connection between two values: \[C(x, y) = \{y - x\}\]

In periodic space, there are infinitely many connections between two values: $$C(x, y) = \{(y + nw) - x \mid n \in \mathbb{Z} \}$$

Sometimes, $y + nw$ is called an image of $y$.
Thus, the set of connections $C(x, y)$ can be seen as the set of connections between $x$ and all images of $y$.

Shortest Connection $\overline{xy}$

The shortest connection \(\overline{xy}\) between two values \(x, y\) is defined as the connection \(c \in C(x, y)\) which has the smallest length $|c|$ among all connections $C(x, y)$: \[\overline{xy} = \underset{c \in C(x, y)}{\mathtt{arg min}} \hspace{1em} |c| \]

Recall that in infinite space, there is only one connection between two values: $$C(x, y) = \{y - x\}$$ The only connection is therefore also the shortest connection. We can thus write $$\overline{xy} = y - x.$$

In periodic space however, there are infinitely many connections between two values: $$C(x, y) = \{(y + nw) - x \mid n \in \mathbb{Z} \}$$

We therefore have to find the shortest one among all of these. Instead of brute-forcing all connections, we can use the following formula to directly compute the shortest connection: $$\overline{xy} = \overline{p}(y - x)$$ Here, \(\overline{p}\) is the periodic boundary operator on \([-\frac{w}{2}, \frac{w}{2})\), which is possibly different from the original boundary operator $p$ operating on $[a, b)$.

We are using the same interval size $w = b - a$ for $p$ and $\overline{p}$.
Proof

We can rewrite the set of connections as $$C(x, y) = \{y - x - mw \mid m \in \mathbb{Z} \}.$$

We are using the letter $m$ instead of $n$ to avoid confusion with the winding number $n$, which we will use later in this proof.

We are stating that among those elements, the one with the smallest absolute value is $$\overline{p}(y - x).$$

To prove this, we have to show two things:

Existence: $\overline{p}(y - x) \in C(x, y)$
Minimality: $|\overline{p}(y - x)| \leq |c| \hspace{1em} \forall c \in C(x, y)$

To prove existence, we simply recall the definition of the periodic boundary operator: $$\overline{p}(y - x) = y - x - nw$$ Since the winding number $n \in \mathbb{Z}$ is an integer, we can choose $m = n$ and get $$\overline{p}(y - x) = y - x - mw \in C(x, y).$$

To prove minimality, we will first represent all connections $c \in C(x, y)$ in terms of $\overline{p}(y - x)$.
For a given $m$ (that is, a given $c$), we define the integer $k = n - m$, where $n$ is the winding number of $y - x$.
Since then, $m = n - k$, we can rewrite $c$ as $$\begin{align*} c &= y - x - mw \\ &= y - x -(n-k)w \\ &= y - x - nw + kw \\ &= \overline{p}(y - x) + kw .\end{align*}$$

The rest of the proof consists of differentiating based on whether of $k = 0$ and showing for both cases that $|\overline{p}(y - x)| \leq |c|$.

Let's first discuss the trivial case $k=0$:
$$c = \overline{p}(y - x) + 0w = \overline{p}(y - x)$$ Thus, $$|\overline{p}(y - x)| = |c| \leq |c|.$$

For the other case $k \neq 0$, we first note that generally, $$|\overline{p}(y - x)| \leq \frac{w}{2}.$$ That is because $\overline{p}(y - x) \in [-\frac{w}{2}, \frac{w}{2})$.

This means that if $|c|$ was to be smaller than $|\overline{p}(y - x)|$, i.e. $|c| < \frac{w}{2}$, then $c$ would have to lie in the interval $(-\frac{w}{2}, \frac{w}{2})$.

However, since $k \neq 0$, the connection $c = \overline{p}(y - x) + kw$ lies in the shifted interval $[-\frac{w}{2}, \frac{w}{2}) + kw.$

In the section about the interval size $w$ we showed that this interval is disjoint to the interval $[-\frac{w}{2}, \frac{w}{2})$ because it is shifted by a non-zero multiple of its size $w$.
In other words, $c$ cannot be in $[-\frac{w}{2}, \frac{w}{2})$, therefore not in $(-\frac{w}{2}, \frac{w}{2})$, and therefore $|c|$ cannot be smaller than $\frac{w}{2}$.

Thus, $|c| \geq \frac{w}{2}$ and therefore $|\overline{p}(y - x)| \leq |c|$.

Instead of always saying the shortest connection, we could simply say the connection (singular), but mean the shortest connection $\overline{xy}$ among all possible connections $C(x, y)$. This way of speaking is sometimes called the Nearest Image Convention.

The Advantage of Symmetrical Spaces

In a symmetrical space $[-\frac{w}{2}, \frac{w}{2})$, the following is true: $$\overline{p} = p$$ In other words: The periodic boundary operator $\overline{p}$ for connections is identical to $p$, boundary operator of the original space!

Proof

$p$ uses the boundaries $[a, b) = [-\frac{w}{2}, \frac{w}{2})$, while $\overline{p}$ uses $[-\frac{w}{2}, \frac{w}{2})$ with the same interval size $w = b - a$ as $p$.
$\Rightarrow p$ uses the same boundaries as $\overline{p}$.
$\Rightarrow \overline{p} = p$.

This fact has two advantages.
The first advantage is conceptional simplicity. Look at the following identity: $$\overline{xy} = y \ominus x$$

Proof
Because \(\overline{p} = p\), $$\begin{align*} \overline{xy} &= \overline{p}(y - x) \\ &= p(y - x) \\ &= y \ominus x. \end{align*}$$

This means that in a symmetrical space, the periodic subtraction $\ominus$ suddenly has a very intuitive meaning: It computes the shortest connection between two values – just like subtraction in infinite space $\mathbb{R}$.

Previously, in asymmetrical spaces, the shortest connections could be from a different space.
In symmetrical spaces, we can completely forget about $\mathbb{R}$ and work exclusively with the group $([-\frac{w}{2}, \frac{w}{2}), \oplus)$.

The second advantage is simplicity regarding the implementation: You only need to implement a single boundary operator which can be used for values as well as for connections.

However, you might still want to use two different implementations for $\overline{p}$ and $p$. Imagine the following scenario:

Part III

Multidimensional

It's very easy to generalize all of this to multi-dimensional space. We can simply apply all operations element-wise, assuming we are using a $\mathtt{p}$-norm such as the Euclidean norm.

Periodic Space (N-Dimensional)

A multi-dimensional periodic space is defined by a list of \(N\) intervals:

\(\vdots\)

The elements of this space are \(N\)-dimensional vectors \(\mathbf{x} = \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_N \end{pmatrix}\) with \(x_i \in [a_i, b_i)\).

Boundary Operator (N-Dimensional)

The boundary operator \(\mathbf{p}\) in \(N\)-dimensional space maps real vectors \(\mathbf{x} \in \mathbb{R}^N\) to the periodic space: \[\mathbf{p}(\mathbf{x}) = \begin{pmatrix} p_1(x_1) \\ p_2(x_2) \\ \vdots \\ p_N(x_N) \end{pmatrix}\]

This is done by using 1-dimensional periodic boundary operators element-wise, where \(p_i\) operates on the interval \([a_i, b_i)\).

Shortest Connections (N-Dimensional)

In the multi-dimensional space, connections are vectors $\mathbf{c} \in \mathbb{R}^N$.

When talking about shortest connections, we need some notion of length of a connection \(\mathbf{c}\). In 1D, we used the absolute value $|c|$ for this.
In multi-dimensional space, we need a vector norm for that. We can use the Euclidean norm or other $\mathtt{p}$-norms such as the maximum norm or the sum norm (aka. Manhattan norm).

The $\mathtt{p}$-norm is defined as: $$\lVert \mathbf{c} \rVert_\mathtt{p} = \left( |c_1|^\mathtt{p} + \ldots + |c_N|^\mathtt{p} \right)^{1/\mathtt{p}}$$

The variable $\mathtt{p} \geq 1$ in the norm is not the same as the periodic boundary operator $p$.

Let's quickly look at two special cases:

  1. For $\mathtt{p} = 2$, we obtain the Euclidean norm: $$\lVert \mathbf{c} \rVert_2 = \sqrt{c_1^2 + \ldots + c_N^2}$$
  2. For $N=1$, we obtain the absolute value, regardless of the value of $\mathtt{p}$: $$\lVert \mathbf{c} \rVert_\mathtt{p} = |c_1|$$

$\mathtt{p}$-norms have a nice property: They are strictly monotone in the absolute value of each coordinate, meaning that an increase in the absolute value of a coordinate always leads to an increase in the overall norm of that vector. In other words: The following statements are equivalent:

Therefore, we can simply compute the shortest connection element-wise: $$\overline{\mathbf{xy}} = \begin{pmatrix} \overline{x_1 y_1} \\ \overline{x_2 y_2} \\ \vdots \\ \overline{x_N y_N} \end{pmatrix}$$

Because those 1D-shortest connections $\overline{x_i y_i}$ are (per definition) guaranteed to be minimal in their absolute value $|\overline{x_i y_i}|$, the $\mathtt{p}$-norm of this vector $$\lVert \overline{\mathbf{xy}} \rVert_\mathtt{p} = \left( |\overline{x_1 y_1}|^\mathtt{p} + \ldots + |\overline{x_N y_N}|^\mathtt{p} \right)^{1/\mathtt{p}}$$ is also minimal.

Summary: To obtain the shortest connection in multi-dimensional space, we compute the shortest connection in each coordinate.

Test Your Knowledge

Suppose we have values \(x, y\) from a periodic interval \([a, b)\).

  1. In which interval do the differences \(x-y\) between those values lie?
    Answer

    \((-w, w)\).
    Note: \(w = b - a\) is the size of the interval.

  2. How do I compute the shortest connection \(\overline{xy}\)?
    Answer
    \(\overline{xy} = \overline{p}(y - x)\).
    \(\overline{p}\) is the periodic boundary operator on \([-\frac{w}{2}, \frac{w}{2})\).
  3. What is the simplest implementation of \(\overline{xy}\)?
    Answer
    We can use the Single Wind implementation:
    function shortest_connection(x, y):
        diff = y - x
        if diff < -w/2:
            return diff + w
        else if diff >= w/2:
            return diff - w
        else:
            return diff
    Explanation:
    We want to implement \(\overline{p}(y - x)\).
    Because \(y - x\) are differences in an interval, our implementation of the boundary operator \(\overline{p}\) only needs to handle values in the range \((-w, w)\). Because \(\overline{p}\) works on the interval \([-\frac{w}{2}, \frac{w}{2})\), this means that we only need to handle values with an absolute winding number \(|n| \leq 1\).
    However, in some scenarios, you want to avoid branching statements. Fortunately, you can always use the Constant implementation:
    function shortest_connection(x, y):
        diff = y - x
        return diff - floor((diff + w/2) / w) * w
    Attention: Since we are implementing \(\overline{p}\) on the interval \([-\frac{w}{2}, \frac{w}{2})\), we need to use \(-\frac{w}{2}\) as value for \(a\) in the implementation, not the \(a\) from our value-interval \([a, b)\)!

Now suppose that we have angles \(\alpha, \beta\) with arbitrary values from \(\mathbb{R}\).

  1. How do I map an angle \(\alpha\) to the periodic interval \([-\pi, \pi)\)?
    Answer

    Use the periodic boundary operator \(p(\alpha)\) on the interval \([-\pi, \pi)\).

    Using the default formula for the periodic boundary operator yields: $$p(x) = x - \left\lfloor \frac{x + \pi}{2\pi} \right\rfloor 2\pi $$

    As discussed in the section about symmetrical spaces, we can express this more concisely using the $\mathtt{round}()$ operation: $$p(x) = x - \mathtt{round}(\frac{x}{2\pi}) 2\pi$$

  2. What's the simplest implementation of this operation?
    Answer
    We can always use the Constant implementation:
    (Here with round() instead of floor() since the interval is symmetrical.)
    function map_to_interval(alpha):
        return alpha - round(alpha / (2*pi)) * 2*pi
    If we can guarantee that the absolute winding number \(|n|\) is always rather low (i.e. \(\alpha\) is close to 0), we could also try the Linear implementation:
    function map_to_interval(alpha):
        while alpha < -pi:
            alpha += 2*pi
        while alpha >= pi:
            alpha -= 2*pi
        return alpha

    The choice between the Constant and Linear implementation depends on the expected size of the angles and the expected performance of the implementation (and therefore the platform you are targeting).

    In any case, we can't use the Single Wind implementation, because \(\alpha \in \mathbb{R}\) can be arbitrarily large, and we can therefore not generally assume that \(|n| \leq 1\).

  3. How do I compute the shortest connection \(\overline{\alpha\beta}\)?
    Answer
    \(\overline{\alpha\beta} = \overline{p}(\beta - \alpha)\).
    This is because the interval \([-\pi, \pi)\) is symmetrical and therefore $\overline{\alpha\beta} = \alpha \ominus \beta$, as discussed in the section about the advantage of symmetrical spaces.
  4. What is the simplest implementation of this operation?
    Answer

    Since we are in a symmetrical space, $\overline{p} = p$. Therefore, we can use the same implementation for the shortest connection as for the mapping to the interval.

    However, we could discuss whether we actually want to use the same implementation for both cases:

    Maybe we can make different assumptions about the differences $\beta - \alpha$ than about the values \(\alpha, \beta\).
    Certainly, we can't use the Single Wind implementation here either, because the difference $\beta - \alpha$ can be arbitrarily large ($\alpha, \beta$ are not from an interval).
    However, it could still be beneficial to use a different implementation for $p(\beta - \alpha)$ than for $p(\alpha)$. For example, assume that angles could be arbitrarily large, but, for some reason, we always only have small differences $\beta - \alpha$ in our application. In that case, we would use the Constant implementation for $p(\alpha)$ and the Linear implementation for $p(\beta - \alpha)$.

Appendix

Modulo-Based Implementations

The problem with modulo-based Constant implementations of the boundary operator is that the modulo operator is often behaving unexpectedly for negative numbers.

For example, in the following languages, the modulo-operator is giving wrong results for negative numbers: There are a few good ones though: You can check a modulo operator for correctness with this test:
-0.1 % 1 == 0.9
You can quick-fix an incorrect modulo operator like this:
def correct_mod(x, w):
    return ((x % w) + w) % w
Alternatively, you can search through this list on Wikipedia to find out if your language's modulo operator is behaving "correctly" (look for floored or Euclidean). The Euclidean modulo-definition only differs from the floor definition for negative divisors. We don't care about that difference because we only use \(w\) as a divisor, which is always positive.