General introduction to implementing physics engines using geometric algebra

In game programming, using a physics engine often is useful.
But a problem with typical physics engines is, that they are only written for a fixed dimension, normally 2D or 3D.
So when you wrote a 2D game and want to extend it to 3D, you have to use a new physics engine. And when you wrote a 3D game, you should be able to reuse the engine, but it will use more memory and processing time than a 2D engine.

General introduction to physics engines

Before discussing, how to generalize a physics engine for 2D, 3D and even arbitrary other numbers of dimensions, you should know, how a physics engine generally works.
First you have to save a set of objects. Often these objects are saved in special data structures for performance reasons, but for simple cases it's enough to save them in some kind of list.
Each of these objects has to save multiple components:

• pos: the position
• vel: the velocity
• dir: the current orientation
• rot: the rotational speed
• shp: the shape for collision

An object may also save other kinds of data.
Three important aspects of physics engines will be discussed:

• The representation and updating of pos and vel.
• The representation and updating of dir and rot.
• Collision detection.

Physical interactions or user events may change vel and rot. Afterwards pos is updated depending on vel and dir is updated depending on rot.
This is repeated many times as long as the simulation is running. One iteration will be called "step".
In order to compute the physical interactions, you have to test, if the objects collide.
Therefore you can iterate over all pairs of objects in the set and call some test function on them.
If two objects collide may depend on pos, dir and shp.
The simplest shape for collision detection is a circle or sphere. But normally arbitrary convex shapes are used, which are constructed out of Vertices connected by line segments for 2D and by plane segments for 3D. These segments have an inner side and an outer side.
If every Vertex of one object lies at the outer side of one segment, it is outside of it, else they collide.

Now the representations of pos and vel and their updates, of dir and rot and their updates, and of the convex shapes and their collision detection will be looked at.
It is assumed, that you already know a bit about math for computer graphics or are at least able to look up something yourself, if you don't understand the short explanations here.
In 2D and in 3D, pos and vel are normally represented as vectors of a known size at compile time. The updates of pos will be a vector addition of pos and vel.

The rotation representations are more difficult. In 2D both, dir and rot can be represented as a single number. dir is the angle the object is rotated around it's center, rot is the angle, it rotates each step. So the update of dir is just an addition of dir to rot.
In 3D dir is represented as a normalized quaternion, which acts as an axis of rotation and an angle, it is rotated around the axis. rot is represented by the normal axis of the rotation. The lenght of this axis represents the speed of the rotation, the orientation of the axis represents the direction of rotation. In this case the update is a bit more complicated: First rot has to be converted to a quaternion. The scalar part of the quaternion is represented as the cosine of the length of the vector, the vector part is represented as the normalized vector multiplied by the sine of the length of the vector. Then both quaternions are multiplied.
You may already notice, it's difficult to generalize this for 2D and 3D or even higher dimensions.

Now let's look at collision detection.
The points are, as positions, represented as vectors with the specific number of elements.
In 2D the segments of the shape are represented as line segments, which can be represented by two points.
In 3D the segments of the shape are represented as plane segments. As a simplification triangles can be used, which can be represented by three points.
So in both cases it's possible to represend the surface of a shape as segments, which are represented by one more point than the dimension, so it's easy to unify.
For the collision detection you have to know on which side of a segment a point lies.
The best way to do this is to calculate the normal of a segment, and the distance vector between one point on the segment and the other point.
Then you calculate the dot product of the normal and the distance vector. If the result is positive, the point lies on the outer side of the segment, if it's negative on the inner side.
This is possible for 2D and 3D. But calculating the normal differs here.
In order to calculate the normal of a line in 2D, you have to take the distance vector (x,y) of two points on the line by subtracting them. The normal then is (-y,x).
Calculating the normal vector in 3D works different. Here you need three points. You calculate both distance vectors from one point to each of the other points. The normal is then calculated by using the cross product.
This also is not unifyable for 2D and 3D nor extendable to higher dimensions.

Implementation using geometric algebra

This is, where using geometric algebra can help.
Geometric algebra unifies many mathematical concepts. And it is capable of unifying this, too.
Now a general introduction to geometric algebra and how to implement the previously mentioned concepts with it, will follow.

In geometric algebra, there only elements are multivectors. They are composed out of mathematical operations of basis vectors. A basis vector is normally denoted by e_n, which are orthogonal to each other.
So the 2D vector (x,y) can be written as sum of two basis vectors x*e_1+y*e_2.
Similarly a 3D vector (x,y,z) can be written as sum of three basis vectors x*e_1+y*e_2+z*e_3.
So the pos and vel can be represented in a pretty simple way.

The most interesting part of geometric algebra is, how multiplication works. There are three kinds of products: The geometric product, the outer product and the inner product.
The geometric product is the sum of the inner product and the outer product.
The simplest way to introduce them is, to show, how these basis vectors work under multiplication.
The geometric product, which probably is the most important one, is normally denoted without an operator between the operands and is defined this way:

• e_i e_i=1
• e_i e_j=e_ij
• e_j e_i=-e_ij
i and j are not equal.

So when you multiply two parallel components, they result in a scalar value.
When you multiply orthogonal components, they result in a product of vectors, in this case bivectors.
That's why they are called multivectors.
In this case, the product not symmetric, you cannot change the order, but antisymmetric, the change of the order will result in a change of the sign.

Before explaining the outer and inner product, some explanations about multivectors will follow.
Vectors can represent a 1D direction in space. You can think of them as sized lines. Bivectors can represent a 2D direction in space. You can think of them as sized planes.
You even can multiply more multivectors together in order to represent more dimensional directions, but that's not often useful in 2D or 3D.
A general multivector in two dimensions can be represented as a vector with four elements.
This vector can be written as (x,x_1,x_2,x_12). As a mathematical expression of it's basis vectors e_1 and e_2, it looks like this: x+x_1 e_1+x_2 e_2+x_12 e_12.
x is the scalar part and has grade 0, (x_1,x_2) is the vector part and has grade 1, x_12 is the bivector part and has grade 2. It is also called pseudoscalar, since it also has just a single value.
A general multivector in three dimensions can be represented as a vector with eight elements.
This vector can be written as x,x_1,x_2,x_12,x_3,x_13,x_23,x_123). As a mathematical expression of it's basis vectors e_1, e_2 and e_3, it looks like this: x+x_1 e_1+x_2 e_2+x_12 e_12+x_3 e_3+x_13 e_13+x_23 e_23+x_123 e_123.
x is the scalar part again, (x_1,x_2,x_3) is the vector part, (x_12,x13,x_23) is the bivector part, x_123 is the pseudoscalar of grade 3 in this case.
A general multivector in n dimensions has 2^n elemens.

In fact, bivectors are already a good representation for rotational speed rot. 2D rotation in the previous approach needed exactly one value to represent rotational speed. A bivector in 2D has exactly one value. For 3D rotation a sized normal vector was used, which has exactly three values, as 3D bivectors. Since bivectors can be thought of planes with a size, the representation even makes more sence now, and is even suited for higher dimensions.
Representing orientation dir is a bit more complicated. In order to rotate something, you can multiply a normalized vector by another normalized vector. The resulting object exactly represents a rotation, so that the one vector is rotated to the other vector. In order to be able to represent rotations around exactly one pi, the rotation it represents is chosen twice as much.
A multiplication of two vectors everytime has a scalar part and a bivector part.
In order to compose such rotations, you can just multiply them. With more than three dimensions, they may even contain higher graded vector parts after composition with even grade. Therefore this type is also called even graded subalgebra. Since it can be used for rotations it's also called rotor.
Since the vectors were normalized, the rotor now also is normalized.
A general rotor in 2D looks like this: x+x_12 e_12.
This is equivalent to complex numbers, if you replace e_12 with i, because it squares to minus one. This will be shown in a few steps:

• e_12 e_12 This product should be minus one.
• e_1 e_2 e_1 e_2 Write them as product of basis vectors.
• - e_1 e_2 e_2 e_1 Switch two basis vectors. Because they are orthogonal, the sign is switched.
• - e_1 e_1 Multiplying two parallel vectors evaluates to one, so they can be removed.
• - 1 Now the two parallel vectors can be evaluated again and the result is minus one.

Since the complex number is normalized, both components can be interpreted as the sine and the cosine of the angle, we saved as angle before. It would be simple to calculate the angle, but normally in order to rotate, you would need the sine and cosine anyway, so this representation may even be more practical.
But the update now is not as simple as addition anymore. It works similar as the 3D approach with the quaternions. First you convert the rotation bivector to a rotor. The scalar part is represented as the cosine, the bivector part as the sine of the rot bivector.
For 3D, a general rotor looks like this: x+x_12 e_12+x_13 e_13+x_23 e_23.
If you change some signs, 3D rotors are equivalent to quaternions. So the approach for 3D doesn't change at all.

But when extending it to higher dimensions, there is a problem: Composition of rotations may result in additional higher dimensional even grade multivector components.
A simple rotor in 4D looks like this: x+x_12 e_12+x_13 e_13+x_23 e_23+x_14 e_14+x_24 e_24+x_34 e_34
When multiplying two of these, a new component x_1234 e_1234 will be created for general bivectors.
The reason is, that in higher dimensions, double rotations are possible.
That means, if you rotate around one plane, there is at least one other orthogonal plane to it, at which it will not be rotated. The extra component in the rotor now will allow this double rotatioon in one operation.
Another problem, that occurs in 4D, is, that addition of bivectors, which represents acceleration, doesn't result in something, that represents a plane, anymore. In 2D the addition of bivectors is just an addition of scalars, in 3D it can be thought of addition of normal vectors.
In 4D, when adding two orthogonal bivectors, it's not possible to think of them as a plane of rotation. That's also because of the double rotation.
But the addition will just work as acceleration. The only problem is the transformation from bivectors to rotors. When the bivector doesn't represent a plane, the previous approach using cosine of the length for the scalar part and the sine of the length times the normalized bivector as the bivector part. You would like to get the higher dimensional part.
There is another function, that also works for complex numbers and quaternions. It's the exponential function. Without scalar part, this funciton is equivalent to the representation with sine and cosine. So instead, the exponential function is a more general way to convert bivectors to rotors, and probably easier. And the best thing, it's even unifiable with higher dimensions.
How the exponential function works for general multivectors or at least bivectors won't be explained here, but you can look for the general definition of exponential functions, if you are interested.
As long as the engine should just support 2D and 3D, using the other representation is totally fine.

Next, the collision detection will be discussed.
There were two things, you need to do. Calculate the normal vector of a segment and then calculate the dot product.
Therefor it's important to know, how the dot and something like the cross product works with geometric algebra. This is, where the outer and inner product will be used.
The inner product is denoted by ^ and defined this way:

• e_i^e_i=0
• e_i^e_j=e_ij
• e_j^e_i=-e_ij

So it's almost the same as the geometric product, but it only has the part of higher grade.
So the outer product everytime results in a higher grade, or zero. It's antisymmetric, and never symmetric.
The inner product is denoted by * defined this way:

• e_i*e_i=1
• e_i*e_j=0
• e_j*e_i=0

That*s already equivalent to the dot product for vectors. You multiply the elements of the same direction and then sum them together. This operation results in a multivector of lower grade everytime.
As you may see, the sum of both operations equals the geometric product.

So first let's think about the simple example in 2D. As before a line is given and you want to know the normal vector of it.
The normal vector is the vector, that is orthogonal to the line in 2D. To get this vector, you can multiply it with the pseudoscalar using the geometric product.
This way the parallel components from the pseudoscalar will be removed and only the orthogonal comonents stay.
Assuming your vector representing the line is called (x,y) and you multiply it with the pseudoscalar, it looks like that:

• (x e_1+y e_2) e_12 The vector is represented as sum of basis vectors.
• x e_1 e_12+y e_2 e_12 Multiply it out.
• x e_1 e_1 e_2+y e_2 e_1 e_2 Convert to basis vectors.
• x e_1 e_1 e_2-y e_1 e_2 e_2 Rotate and change sign.
• x e_2-y e_1 Remove duplicates.

Written as vector it looks like this: (-y,x).
This is exactly the same as you would calculate it in 2D.
Since the other elements already are vectors, you just can use the inner product for calculating the dot product.
For 3D you can do it in a similar way. In this case, you don't have a single line, but two lines, that lie on a plane. When multiplying single lines with the pseudoscalar of 3D, the resulting orthogonal component will be a bivector. Therefore you have to create a bivector out of the plane. In order to represent the plane as bivector, you just have to multiply both lines using the outer product. Using the geometric product would result in a rotor.
Now you can just multiply by the pseudoscalar and get the normal vector.
The dot product can then just be calculated as before.
This example demonstrates, what the outer product has in common with the cross product in 3D.
You also see, that the approach is not that different with different dimensions.
With more dimensions, you just get more points for each segment, so also the segment has a higher dimension. You get the normal by multiplying some lines of the segment, which are not parallel, using the outer product. Then you multiply by the pseudoscalar of the matching dimension and can calculate the dot product of the result and the vector between one point of the segment and one other point.

Conclusion

These were just some examples, how easy it is to extend a physics engine to geometric algebra.
There are still some parts to cover, but other aspects should stay similar to current approaches or could be derived with the knowledge in this document.
Currently the biggest problem with implementing something like that is, that there are almost no good libraries for geometrc algebra.
It's possible to implement libraries, where for every kind of object, the data structures are created at compile time for every needed dimension. Using general multivectors would be memory waste for representing vectors, especially in higher dimensions, where the multivectors size goes up exponentially.
I may cover this topic in another post.
Tell me, when you are interested in something special.
I may want to write about implementing geometric algebra itself in the most efficient way, maybe even with using compile time generics, about more advanced examples or show some example code for a physics engine programmed this way, maybe even as some kind of tutorial, so you can reprogram it.
If you didn't understand somethng, feel also free to write comments.
I hope, you found interest in geometric algebra and maybe even consider using it for programming.