Game Programming: Math Libraries

Most game engines define their own custom set of math libraries for use in their games.  Unity, for example, has defined its own Mathf library for  performing common operations on floating point numbers.  You may wonder what is wrong with just using System.Math as apposed to UnityEngine.Mathf.  Some people use Mathf instead of Math because it saves them the trouble of having to do double to float conversions.  However, in most cases Math is considered faster than Mathf.  Standard math libraries are certainly convenient and will most often be more efficient than any math libraries that you roll yourself, so why go to the trouble of implementing a custom library.  The problem with standard math libraries is that they may not always have the kinds of functions that are most important to game programming.  Other times it may be better to forgo accuracy in favor of efficiency by replacing expensive operations such as cos, sin and ex with less expensive and less accurate approximations of these functions.  In this article I will discuss some of the most commonly used functions found in Unity’s Mathf library and how you may implement them yourself.

Approximately:

Due to the nature of floating point numbers it is very common for operations performed on them to result in some form of rounding error.  This can lead to problems when trying to compare the equality of floating point numbers.  For example we may end up comparing a number like 2.0 to a number like 1.999999999.  For all effects and purposes these numbers are practically the same but when we check for equality we get a false statement.  This is why we sometimes need to take into account some margin of error when comparing floats.  This is exactly what the Mathf.Approximately function does.  It can be implemented in the following way:

bool Approximately(float a, float b) {
    return Abs(a - b) < Epsilon;
}

Here we take the absolute value of the difference between the two input variables and check if it is less than some arbitrarily small value called Epsilon.  It should be noted however that although some standard math already have a predefined value for epsilon, it may actually be better for the programmer to define their own value for epsilon since the standard math library’s definition of epsilon may be too small to make any difference from straight equality comparison.  An alternative approach could also be to define an Approximately function that takes the tolerance as a parameter and just replace Epsilon with this tolerance value.

Clamp:

Clamp is such a useful tool that it is any wonder why it is not a common tool found in every math library, especially since it is so simple to write.  Perhaps it is just so simple that some programmers feel it is not necessary to include it in a standard math library.  Simply put, clamp ensures that its input value lies within a certain range of values.  If the number lies between the given range it simply returns the number itself, otherwise it returns the nearest endpoint.  An example where this might be useful is in the case where the programmer wants to ensure that a player character never goes past the edges of the screen.  All the programmer has to do is Clamp the player’s position to within the screen boundaries.  Clamp can be implemented in the following way:

float Clamp(float value, float min, float max) {
    return (value < min)? min : (value > max)? max : value;
}

This function is so useful that Unity has even included a Clamp01 function that specifically clamps a floating point value to between 0 and 1.

Linear Interpolation (Lerp):

Linear interpolation is the process by which we find a point somewhere between two input values.  We can visualize it as a line drawn between two points.  Using the equation for a line we can write the formula for the desired value y as y = a + (b - a)x where a and b are the start and end values respectively and x represents what fraction of the way we are between our two values.  Note that we only get a value between a and b when x lies within the range [0, 1]; any other value for x results in some value outside the interval [a, b].  We also see that y = a when x = 0 and y = b when x = 1.  However, to avoid rounding error caused by floating point arithmetic, it is best to rewrite the equation as y = a(1 - x) + bx (this ensures that y = b when x = 1).  Depending on how you want your Lerp function to behave for certain values of x you may wish to clamp it to between 0 and 1 so as to never overshoot the endpoints.  However, in some rare cases it may be better to leave the value unclamped.  The following code sample shows how we may implement the former:

float Lerp(float a, float b, float t) {
    t = Clamp01(t);
    return a * (1 - t) + b * t;
}

For more help learning the proper use of linear interpolation, see How to Lerp like a Pro by Robert Utter.

Move Towards:

While linear interpolation certainly has its uses, it is sometimes not entirely practical since it requires the programmer to know what point we started at and what fraction of the way between the two points we are at a given time.  But suppose instead we want to procedurally step some value towards some target using some kind of displacement value.  We want a way to find our new value while making sure that we do not overshoot our target.  To do this we simply compare our expected displacement to the difference between our current value and our target value.  If our expected displacement is greater than the actual distance to our target then we just return our target value.  Here is an example:

float MoveTowards(float current, float target, float maxDelta) {
    float delta = target - current;
    return (maxDelta >= delta) ? target : current + maxDelta;
}

Note that in the implementation above a negative value of maxDelta will cause our output to procedurally step away from  our target value.

Smooth Damp:

While MoveTowards generally serves its intended purpose it can often result in game mechanics that look somewhat jerky since it forces something to come to a sudden stop once it reaches its goal.  If we want our game to look somewhat cleaner we need to apply smoothing to our calculations.  One way that we can apply smoothing is to apply some kind of S-curve to give us that smooth ease-in/ease-out motion.  One commonly used formula is x = 6t^5 - 15t^4 + 10t^3 which gives us such an S-curve and can be used with little computing cost.

C2 S-curve

Another technique to achieve smoothing is through exponential decay.  This is a kind of ease-out method of smoothing and can be accomplished by the statement x = x + (desiredX - x) * 0.1f * timeDelta where timeDelta is the time step between updates.  The advantage of this type of smoothing is that you can smooth toward a changing target and you do not need to know how much time has elapsed since the start.  However, the problem with this type of smoothing is that the initial motion is sudden.

Exponential Decay Diagram

The approach we will be using, however, will rely on the model for a critically damped spring.  This model fixes the problem we have with the exponential decay model by maintaining the ease-in property found in the S-curve.

Crtically Damped Spring Diagram

A spring will always try to maintain a certain length by exerting a force that pushes it towards that desired length.  Without damping this system will oscillate back and forth past the desired length forever.  We say the system is critically damped when it returns to equilibrium in the least possible amount of time.  We can write the equation for the force acted on the spring as F = -k(x - x_d) - bv where k is the spring constant, x_d - x is the displacement from our desired position, b is the damping coefficient and v is the velocity of the end of the spring that works to slow the spring down.  Using Newton’s second law of motion, F = ma, and writing the velocity and acceleration as the first and second derivatives of the position, we can obtain the ordinary differential equation (ODE):

m\frac{d^2x}{dt^2} + b\frac{dx}{dt} + kx = kx_d    (1)

The next step requires us to know how to solve for an equation of x.  The above ODE will have a solution of the form x(t) = x_p(t) + x_c(t) where x_p(t) is a particular solution of the equation and x_c(t) is the complementary solution of the associated homogeneous DE.  We see that the particular solution is of the form x_p(t) = A where A is a constant since the right side of Equation 1 is constant.  By taking the first and second derivatives of x_p(t) and plugging them into Equation 1 we see that A = x_d.

That was our particular solution; now let us find our complementary solution.  The homogeneous DE gives us the characteristic equation mr^2 + br + k = 0 which has the roots:

r = \frac{-b \pm \sqrt{b^2 - 4mk}}{2m}

It can be shown that the system will be critically damped when our equation has the real repeated roots r = -\sqrt{\frac{k}{m}} which only occurs when b^2 = 4mk.  We can use this to simplify Equation 1:

\frac{d^2x}{dt^2} + 2\omega\frac{dx}{dt} + \omega^2(x_d - x) = 0 where \omega = \sqrt{\frac{k}{m}}    (2)

In the above formula \omega is the spring’s natural frequency and represents the stiffness or strength of the spring.  But now we also have our complementary solution which is x_c(t) = (c_0 + c_1t)e^{-\omega t}.  So our general solution is:

x(t) = x_d + (c_0 + c_1t)e^{-\omega t}

Solving for the constants c_0 and c_1 using the initial conditions x(0) = x_0 and x'(0) = v_0 gives us the actual solution:

x(t) = x_d + ((x_0 - x_d) + (v_0 + \omega(x_0 - x_d))t)e^{-\omega t}    (3)

This formula now tells us everything we need to calculate the new position and by taking the derivative we can find the new velocity of the object after after this time step:

v(t) = (v_0 - (v_0 + \omega(x_0 - x_d))\omega t)e^{-\omega t}    (4)

Now the only parameters we need for our smoothing function are our current position and velocity, our target position our elapsed time and finally our smoothness factor \omega.  But how do we figure out what we should use for our smoothness factor?  It would be much more intuitive to control our function based on how long we want it to take to reach equilibrium.  If we define this smooth time as being the expected time to reach the target when at maximum velocity, then we solve for it quite simply by setting the acceleration to zero and the maximum velocity to \frac{x_d - x_0}{t_{sm}} where t_{sm} is the smooth time and plugging these values into Equation 2.  We then see that \omega = \frac{2}{t_{sm}}.

One final note before we write our smoothing function is that Equation 3 and 4 contain e^x.  We talked before about how operations such as these are expensive so the following code shows a polynomial that can be used to approximate e^x.  While I do not know exactly how this polynomial is produced, I would hazard a guess and say that it is probably constructed by approximating a higher order Taylor series polynomial down to a third order polynomial.  Finally to make our function exactly like the SmoothDamp function of Unity’s Math library we also need to clamp the speed according to some defined max speed.

float SmoothDamp(float current, float target, float& currentVelocity, float smoothTime,
        float maxSpeed = Infinity, float deltaTime = Time.deltaTime) {
    smoothTime = Max(0.0001f, smoothTime);
    float omega = 2.0f / smoothTime;
    float x = omega * deltaTime;
    float exp = 1.0f / (1.0f + x + 0.48 * x * x + 0.235 * x * x * x);
    float deltaX = target - current;
    float maxDelta = maxSpeed * smoothTime;

    // ensure we do not exceed our max speed
    deltaX = Clamp(deltaX, -maxDelta, maxDelta);
    float temp = (currentVelocity + omega * deltaX) * deltaTime;
    float result = current - deltaX + (deltaX + temp) * exp;
    currentVelocity = (currentVelocity - omega * temp) * exp;

    // ensure that we do not overshoot our target
    if (target - current > 0.0f == result > target) {
        result = target;
        currentVelocity = 0.0f;
    }
    return result;
}

Resources:

Kirmse, Andrew.  (2004).  Game Programming Gems 4.  Charles River Media.  ISBN 9781-58450-295-1.

Zill, Dennis G.  (2009).  A First Course in Differential Equations with Modeling Applications.  Brooks Cole.  ISBN 978-0-495-10824-5

External Links:

Critically Damped Spring Smoothing

How to Lerp like a Pro

Unity Scripting API: Mathf

Unity Forum: Formula behind SmoothDamp

Leave a comment