Dual Numbers & Automatic Differentiation

In the last post, I talked about imaginary numbers, complex numbers, and how to use them to rotate vectors in 2d.

In this post, I want to share another interesting type of number called a “Dual Number” that uses the symbol ε (epsilon) and has a neat trick of automatically calculating the derivative of a function while you calculate the value of the function at the same time.

Dual numbers are pretty similar to imaginary numbers but there is one important difference. With imaginary numbers, i^2 = -1, but with dual numbers, ε^2 = 0 (and ε is not 0!). That may seem like a small difference, but oddly, that opens up a whole interesting world of mathematical usefulness.

Before we dig into automatic differentiation, I want to go over the mathematical basics for how dual numbers behave.

Basic Dual Number Math

Adding dual numbers is the same as adding complex numbers; you just add the real and dual parts separately:

(3 + 4ε) + (1 + 2ε) = 4 + 6ε

Subtraction works the same way as well:

(3 + 4ε) – (1 + 2ε) = 2 + 2ε

To multiply dual numbers, you use F.O.I.L. just like you do with complex numbers:

(3 + 4ε) * (1 + 2ε) =
3 + 6ε + 4ε + 8ε^2 =
3 + 10ε + 8ε^2

However, since ε^2 is zero, the last term 8ε^2 disappears:
3 + 10ε

It’s interesting to note that with complex numbers, the i^2 became -1, so the last term changed from imaginary to real, meaning that the imaginary numbers fed back into the real numbers during multiplication. With dual numbers, that isn’t the case, the dual numbers don’t feed back into the real numbers during multiplication.

In both complex and dual numbers the real terms do affect the non real terms during multiplication.

The division operator relates to the conjugate. I have source code for it below, and some of the links at the end of the post go into the details of that and other operations.

Quick Review: Derivatives (Slope)

If you know the line formula y=mx+b, but you don’t know what a derivative is you are in luck. Remember how “m” is the slope of the line, specifying how steep it is? That is what the derivative is too, it’s just the slope.

Below is a graph of y=2x+1. At every point on that line, the derivative (or slope) is 2. That means that for every step we make on the x axis to the right (positive direction), we make 2 steps up on the y axis (positive direction).

Now, check out this graph of y=x^2-0.2

The derivative (or slope) at every point on this graph is 2x. That means that the slope changes depending on where the x coordinate is!

So, when x=0, the slope is 0. You can see that in the graph where x=0, that it is horizontal, meaning that a step on the x axis becomes no steps on the y axis (only at that point where x is 0, and only if you take an infinitely small step).

When x is 1, the slope is 2, when x is 2, the slope is 4, when x is 3, the slope is 6. Since the numbers increase as we increase x from 0, that tells us that the graph gets steeper as we go to the right, which you can see in the graph.

Alternately, when x is -1, the slope is -2, when x is -2, the slope is -4, and when x is -3, the slope is -6. This shows us that as we decrease x from 0, the graph gets steeper in the opposite direction, which you can see in the graph as well.

What is Automatic Differentiation?

Let’s say you have a function (possibly a curve) describing the path of a rocket, and you want to make the rocket point down the path that it’s traveling.

One way you might do this is to evaluate your function f(T) to get the current location of your rocket (where T is how long the rocket has been flying), and then calculate the derivative f'(T) to find the slope of the graph at that point so that you can orient the rocket in that direction.

You could calculate the value and slope of the function at time T independently easily enough if you know how to get the derivative of a function (a calculus topic), or use wolframalpha.com.

However, if you have a complex equation, or maybe if the equation is controlled by user input, or game data, it might not be so easy to figure out what the derivative is at run time.

For instance… imagine having a function that rolled random numbers to figure out what mathematical operation it should preform on a number next (if we roll a 0, add 3, if we roll a 1 multiply by 2, if we roll a 2, square the number… etc). It isn’t going to be simple to take the derivative of the same mathematical function.

Here enters automatic differentiation (or AD). AD lets you calculate the derivative WHILE you are calculating the value of the function.

That way, you can do whatever math operations you want on your number, and in the end you will have both the value of f(T) as well as the derivative f'(T).

Using ε for Automatic Differentiation

You can use dual number operations on numbers to calculate the value of f(x) while also calculating f'(x) at the same time. I’ll show you how with a simple example using addition and multiplication like we went over above.

We’ll start with the function f(x)=3x+2, and calculate f(4) and f'(4).

the first thing we do is convert our 4 into a dual number, using 1 for the dual component, since we are plugging it in for the value of x, which has a derivative of 1.

4+1ε

Next, we want to multiply that by the constant 3, using 0 for the dual component since it is just a constant (and the derivative of a constant is 0)

(4+1ε) * (3 + 0ε) =
12 + 0ε + 3ε + 0ε^2 =
12 + 3e

Lastly, we need to add the constant 2, using 0 again for the dual component since it’s just a constant.
(12 + 3ε) + (2 + 0ε) =
14 + 3ε

In our result, the real number component (14) is the value of f(4) and the dual component (3) is the derivative f'(4), which is correct if you work it out!

Let’s try f(5). First we convert 5 to a dual number, with the dual component being 1.

5 + 1ε

Next we need to multiply it by the constant 3 (which has a dual component of 0)

(5 + 1ε) * (3 + 0e) =
15 + 0ε + 3ε + 0ε^2 =
15 + 3ε

Now, we add the constant 2 (which has a dual component of 0 again since it’s just a constant)
(15 + 3ε) + (2 + 0ε) =
17 + 3ε

So, our answer says that f(5) = 17, and f'(5) = 3, which again you can verify is true!

Quadratic Example

The example above worked well but it was a linear function. What if we want to do a function like f(x) = 5x^2 + 4x + 1?

Let’s calculate f(2). We are going to first calculate the 5x^2 term, so we need to start by making a dual number for the function parameter x:
(2 + 1ε)

Next, we need to multiply it by itself to make x^2:
(2 + 1ε) * (2 + 1ε) =
4 + 2ε + 2ε + 1ε^2 =
4 + 4ε

(remember that ε^2 is 0, so the last term disappears)

next, we multiply that by the constant 5 to finish making the 5x^2 term:
(4 + 4ε) * (5 + 0ε) =
20 + 0ε + 20ε + 0ε^2 =
20 + 20ε

Now, putting that number aside for a second we need to calculate the “4x” term by multiplying the value we plugged in for x by the constant 4
(2 + 1ε) * (4 + 0ε) =
8 + 0ε + 4ε + 0ε^2 =
8 + 4ε

Next, we need to add the last 2 values together (the 5x^2 term and the 4x term):
(20 + 20ε) + (8 + 4ε) =
28 + 24ε

Lastly, we need to add in the last term, the constant 1
(28 + 24ε) + (1 + 0ε) =
29 + 24e

There is our answer! For the equation y = 5x^2 + 4x + 1, f(2) = 29 and f'(2) = 24. Check it, it’s correct (:

As one last example let’s calculate f(10) and f'(10) with the same function above y = 5x^2 + 4x + 1.

First, to start calculating the 5x^2 term, we need to make 10 into a dual number and multiply it by itself to make x^2:
(10 + 1ε) * (10 + 1ε) =
100 + 10ε + 10ε + 1ε^2 =
100 + 20ε

Next, we multiply by the constant 5 to finish making the 5x^2 term:
(100 + 20ε) * (5 + 0ε) =
500 + 0ε + 100ε + 0ε^2 =
500 + 100ε

Putting that aside, let’s calculate the 4x term by multiplying our x value by the constant 4:
(10 + 1ε) * (4 + 0ε) =
40 + 0ε + 4ε + 0ε^2 =
40 + 4ε

Lastly, let’s add our terms: 5x^2, 4x and the constant 1
(500 + 100ε) + (40 + 4ε) + (1 + 0ε) =
541 + 104ε

The answer tells us that for the equation y = 5x^2 + 4x + 1, f(10) = 541 and f'(10) = 104.

Sample Code

There are lots of other mathematical operations that you can do with dual numbers. I’ve collected as many as I was able to find and made up some sample code that uses them. The sample code is below, as well as the program output.

#include 
#include 

#define PI 3.14159265359f

// In production code, this class should probably take a template parameter for
// it's scalar type instead of hard coding to float
class CDualNumber
{
public:
	CDualNumber (float real = 0.0f, float dual = 0.0f)
		: m_real(real)
		, m_dual(dual)
	{
	}

	float Real () const { return m_real; }
	float Dual () const { return m_dual; }

private:
	float m_real;
	float m_dual;
};

//----------------------------------------------------------------------
// Math Operations
//----------------------------------------------------------------------
inline CDualNumber operator + (const CDualNumber &a, const CDualNumber &b)
{
	return CDualNumber(a.Real() + b.Real(), a.Dual() + b.Dual());
}

inline CDualNumber operator - (const CDualNumber &a, const CDualNumber &b)
{
	return CDualNumber(a.Real() - b.Real(), a.Dual() - b.Dual());
}

inline CDualNumber operator * (const CDualNumber &a, const CDualNumber &b)
{
	return CDualNumber(
		a.Real() * b.Real(),
		a.Real() * b.Dual() + a.Dual() * b.Real()
	);
}

inline CDualNumber operator / (const CDualNumber &a, const CDualNumber &b)
{
	return CDualNumber(
		a.Real() / b.Real(),
		(a.Dual() * b.Real() - a.Real() * b.Dual()) / (b.Real() * b.Real())
	);
}

inline CDualNumber sqrt (const CDualNumber &a)
{
	float sqrtReal = ::sqrt(a.Real());
	return CDualNumber(
		sqrtReal,
		0.5f * a.Dual() / sqrtReal
	);
}

inline CDualNumber pow (const CDualNumber &a, float y)
{
	return CDualNumber(
		::pow(a.Real(), y),
		y * a.Dual() * ::pow(a.Real(), y - 1.0f)
	);
}

inline CDualNumber sin (const CDualNumber &a)
{
	return CDualNumber(
		::sin(a.Real()),
		a.Dual() * ::cos(a.Real())
	);
}

inline CDualNumber cos (const CDualNumber &a)
{
	return CDualNumber(
		::cos(a.Real()),
		-a.Dual() * ::sin(a.Real())
	);
}

inline CDualNumber tan (const CDualNumber &a)
{
	return CDualNumber(
		::tan(a.Real()),
		a.Dual() / (::cos(a.Real()) * ::cos(a.Real()))
	);
}

inline CDualNumber atan (const CDualNumber &a)
{
	return CDualNumber(
		::atan(a.Real()),
		a.Dual() / (1.0f + a.Real() * a.Real())
	);
}

inline CDualNumber SmoothStep (CDualNumber x)
{
	// f(x) = 3x^2 - 2x^3
	// f'(x) = 6x - 6x^2
	return x * x * (CDualNumber(3) - CDualNumber(2) * x);
}

//----------------------------------------------------------------------
// Test Functions
//----------------------------------------------------------------------

void TestSmoothStep (float x)
{
	CDualNumber y = SmoothStep(CDualNumber(x, 1.0f));
	printf("smoothstep 3x^2-2x^3(%0.4f) = %0.4fn", x, y.Real());
	printf("smoothstep 3x^2-2x^3'(%0.4f) = %0.4fnn", x, y.Dual());
}

void TestTrig (float x)
{
	CDualNumber y = sin(CDualNumber(x, 1.0f));
	printf("sin(%0.4f) = %0.4fn", x, y.Real());
	printf("sin'(%0.4f) = %0.4fnn", x, y.Dual());

	y = cos(CDualNumber(x, 1.0f));
	printf("cos(%0.4f) = %0.4fn", x, y.Real());
	printf("cos'(%0.4f) = %0.4fnn", x, y.Dual());

	y = tan(CDualNumber(x, 1.0f));
	printf("tan(%0.4f) = %0.4fn", x, y.Real());
	printf("tan'(%0.4f) = %0.4fnn", x, y.Dual());

	y = atan(CDualNumber(x, 1.0f));
	printf("atan(%0.4f) = %0.4fn", x, y.Real());
	printf("atan'(%0.4f) = %0.4fnn", x, y.Dual());
}

void TestSimple (float x)
{
	CDualNumber y = CDualNumber(3.0f) / sqrt(CDualNumber(x, 1.0f));
	printf("3/sqrt(%0.4f) = %0.4fn", x, y.Real());
	printf("3/sqrt(%0.4f)' = %0.4fnn", x, y.Dual());

	y = pow(CDualNumber(x, 1.0f) + CDualNumber(1.0f), 1.337f);
	printf("(%0.4f+1)^1.337 = %0.4fn", x, y.Real());
	printf("(%0.4f+1)^1.337' = %0.4fnn", x, y.Dual());
}

int main (int argc, char **argv)
{
	TestSmoothStep(0.5f);
	TestSmoothStep(0.75f);
	TestTrig(PI * 0.25f);
	TestSimple(3.0f);
	return 0;
}

Here is the program output:

Closing Info

When you are thinking what number ε has to be so that ε^2 is 0 but ε is not 0, you may be tempted to think that it is an imaginary number, just like i (the square root of -1) that doesn’t actually exist. This is actually not how it is… I’ve seen ε described in two ways.

One way I’ve seen it described is that it’s an infinitesimal number. That sort of makes sense to me, but not in a concrete and tangible way.

The way that makes more sense to me is to describe it as a matrix like this:
[0, 1]
[0, 0]

If you multiply that matrix by itself, you will get zero(s) as a result.

In fact, an alternate way to implement the dual numbers is to treat them like a matrix like that.

I also wanted to mention that it’s possible to modify this technique to get the 2nd derivative of a function or the 3rd, or the Nth. It isn’t only limited to the 1st derivative. Check the links at the bottom of this post for more info, but essentially, if you want 1st and 2nd derivative, you need to make it so that ε^3 = 0 instead of ε^2 = 0. There is a way to do that with matrices.

Another neat thing is that you can also extend this into multiple dimensions. This is useful for situations like if you have some terrain described by mathematical functions, when you are walking the grid of terrain to make vertex information, you can get the slope / gradient / surface normal at the same time.

Lastly, I wanted to mention a different kind of number called a hyperbolic number.

The imaginary number i^2 = -1 and we can use it to do 2d rotations.

The dual number ε^2 is 0 (and ε is not 0) and we can use it to do automatic differentiation.

Hyperbolic numbers have j, and j^2 = 1 (and j is not 1). I’m not sure, but I bet they have some interesting usefulness to them too. It would be interesting to research that more sometime. If you know anything about them, please post a comment!

Links

This shadertoy is what got me started looking into dual numbers. It’s a mandelbrot viewer done by iq using dual numbers to estimate a distance from a point to the mandelbrot set (as far as I understand it anyhow, ha!). He uses that estimated distance to color the pixels.

Shadertoy: Dual Complex Numbers

I didn’t get very much into the reasons of why this works (has to do with taylor series terms disappearing if ε^2 is 0), or the rigorous math behind deriving the operators, but here are some great links I found researching this stuff and putting this blog post together.

Wikipedia: Dual Number
[Book] Dual-Number Methods in Kinematics, Statics and Dynamics By Ian Fischer
[GDC2012] Math for Game Programmers: Dual Numbers by Gino van den Bergen
Stackexchange: Implementing trig functions for dual numbers
Exact numeric nth derivatives
Automatic Differentiation with Dual numbers
Wikipedia: Automatic Differentiation

Using Imaginary Numbers To Rotate 2D Vectors

I’m a big fan of “exotic” math and programming techniques. There’s nothing I like more than seeing something common used in an uncommon way to do something that I didn’t know was possible.

In this post I share a technique that lets you use imaginary numbers (complex numbers more specifically) to be able to rotate vectors in 2d. This technique is so simple that you can even use it to rotate vectors by hand on paper!

Quick Review: Imaginary and Complex Numbers

The imaginary number “i” is the square root of -1. Without using i, you can’t take the square root of a negative number because if the answer is negative, multiplying a negative by a negative is positive, and if the answer is a positive, multiplying a positive by a positive is still a positive.

But, using i, you CAN take the square root of negative numbers.

The square root of 25 is 5.

The square root of -25 is 5*i, or just 5i, which means “5 times the square root of -1”.

You can combine a real and imaginary number into something called a complex number like this: 3 + 5i

Quick Review: Multiplying Complex Numbers

We’ll need to be able to multiply complex numbers together to do our rotations.

(3+5i)*(2+3i) = ???

Luckily, multiplying complex numbers together is as simple as using F.O.I.L. if you remember that from high school math class, it stands for First, Outer, Inner, Last.

Using that, we can multiply and then combine term, remembering that i^2 is -1:

(3+5i)*(2+3i) =
6 + 9i + 10i + 15i^2 =
6 + 19i + 15i^2 =
6 + 19i – 15 =
-9 + 19i

There we go, that’s all there is to multiplying complex numbers together!

Getting Down To Business: Rotating

Let’s say that we want to rotate the vector (0,1) by the angle of the vector (1,1). To do that, we just convert the vectors to complex numbers, using the x axis as the real number component, and the y axis as the imaginary number component, then multiply them, and convert back to vectors.

That’s a mouth full, so here it is, step by step.

1) Convert Vectors to Complex Numbers

(0,1) = 0 + 1i = i
(1,1) = 1 + 1i = 1 + i

2) Multiply the Complex Numbers

i * (1 + i) =
i + i^2 =
i – 1 =
-1 + i

In the above we change i – 1 to -1 + i to make the next step easier. The real component of the complex number is the X axis and the imaginary component is the Y axis.

3) Convert from Complex Number to Vector

-1 + i =
-1 + 1i =
(-1, 1)

The diagram below shows this operation:

As you can see we rotated the purple vector (0,1) which has an angle of 90 degrees and length of 1, by the blue vector (1,1) which has an angle of 45 degrees and a length of sqrt(2), and as a result we got the tan vector (-1,1) which has an angle of 135 degrees and a length of sqrt(2). So, we essentially rotated the 90 degree angle by the 45 degree angle and ended up with a 135 degree angle.

Note that order of operations doesn’t matter, you could rotate the 90 degree angle by 45 degrees, or you could rotate the 45 degree angle by 90 degrees, either way you are going to end up with 135 degrees.

A caveat with this technique though is that when you do the rotation, the resulting vector’s length will be the the product of the two source vectors used; you won’t get a normalized vector out as a result unless your source vectors are normalized, or you normalize after you are done with your operations.

As another example, let’s rotate the vector (4,3) by the vector (-12,5).

The first step is to convert them to complex numbers:
(4,3) = 4 + 3i
(-12,5) = -12 + 5i

Then we multiply them:
(4 + 3i) * (-12 + 5i) =
-48 + 20i – 36i + 15i^2 =
-48 – 16i – 15 =
-63 – 16i

Then convert back to a vector to get: (-63, -16)

In the image, you can see that we started with the blue vector (4,3) which is about 37 degrees and has a length of 5.

We rotated that vector by the purple vector (-12,5) which is about 157 degrees and has a length of 13.

That gave us the tan vector of (-63, -16) which is about 194 degrees and has a length of 65.

The resulting vector has the rotations added, and the lengths multiplied together.

Unrotating

If we multiply the complex numbers together to add rotations, does that mean that we divide the complex numbers to subtract rotations? Sort of…

If you want to subtract a rotation, you multiply the imaginary part of the vector you want to subtract by -1 (or just flip the sign) and then multiply as normal.

When you flip the sign of the imaginary part, this is actually called the “complex conjugate” but if that term scares you, feel free to ignore it 😛

As an example of unrotation, let’s take the vector (5,1) and subtract the vector (2,2) from it to see what we get.

The first step is to convert them into complex numbers.
(5,1) = 5 + i
(2,2) = 2 + 2i

Next up, we are going to flip the imaginary part of the vector we want to subtract.
2 + 2i becomes 2 – 2i

Then, we multiply as per normal:
(5 + i) * (2 – 2i) =
10 – 10i + 2i – 2i^2 =
10 – 8i + 2 =
12 – 8i

Finally, we convert back to a vector to get (12, -8).

In the image above, we started with the blue vector (5,1) which is about 11 degrees and has a length of sqrt(26).

We then unrotated by the purple vector (2,2) which is 45 degrees and has a length of sqrt(8).

That gave us the tan vector as a result (12,-8) which is about -34 degrees and has a length of sqrt(208).

Our resulting vector is the result of the second vector’s angle subtracted from the first vector’s angle, but the length of the vectors are still multiplied together, not divided.

Unrotating to Get Vector Length

As a fun aside, if you unrotate a vector by itself, the result will be that the imaginary components (the Y component) will disappear, and in the result, the real component (the x component) will be the squared length of the vector.

It’s easy enough to calculate the squared length of a vector by adding x^2 and y^2 but this is an interesting property.

In fact, in the last post I mentioned CORDIC math using imaginary numbers to rotate vectors to try and find sine, cosine, etc. This is related to how it does that work. It basically rotates or unrotates a vector by smaller and smaller angles iteratively, based on the sign of the y (imaginary) component to try and get y to zero, which leaves the answer in the x component. It also has some other magic sprinkled in to where it only has to deal with integer math.

Hopefully before too long I’ll be able to make a CORDIC blog post and talk more about that in detail.

Example Code

Ok so this theory stuff is all well and good but how about some code?

Before we get to that I want to give the naked formulas for rotating or unrotating vectors.

Given two vectors (A.x, A.y) and (B.x, B.y)…

Rotating A by B to get C:
C.x = A.x*B.x – A.y*B.y
C.y = A.x*B.y + A.y*B.x

Note: In the resulting vector C, the angles will be added, but the vector lengths will be multiplied.

Unrotating A by B to get C, we just flip the sign of any terms that contain B.y:
C.x = A.x*B.x + A.y*B.y
C.y = -A.x*B.y + A.y*B.x

Note: In the resulting vector C, the angles will be subtracted, but the vector lengths will be multiplied.

Below is some sample code, along with the output, showing our examples above working correctly.

#include 

struct SVector
{
	float x;
	float y;
};

void Rotate (const SVector &A, const SVector &B, SVector &C)
{
	C.x = A.x * B.x - A.y * B.y;
	C.y = A.x * B.y + A.y * B.x;
}

void Unrotate (const SVector &A, const SVector &B, SVector &C)
{
	C.x = A.x * B.x + A.y * B.y;
	C.y = -A.x * B.y + A.y * B.x;
}

void print (const SVector &V)
{
	printf("(%0.2f,%0.2f)", V.x, V.y);
}

int main(int argc, char **argv)
{
	{
		SVector testA = {0.0f, 1.0f};
		SVector testB = {1.0f, 1.0f};
		SVector testC = {0.0f, 0.0f};
		Rotate(testA, testB, testC);
		printf("Rotating ");
		print(testA);
		printf(" by ");
		print(testB);
		printf(" = ");
		print(testC);
		printf("n");
	}
	{
		SVector testA = {4.0f, 3.0f};
		SVector testB = {-12.0f, 5.0f};
		SVector testC = {0.0f, 0.0f};
		Rotate(testA, testB, testC);
		printf("Rotating ");
		print(testA);
		printf(" by ");
		print(testB);
		printf(" = ");
		print(testC);
		printf("n");
	}
	{
		SVector testA = {5.0f, 1.0f};
		SVector testB = {2.0f, 2.0f};
		SVector testC = {0.0f, 0.0f};
		Unrotate(testA, testB, testC);
		printf("Unrotating ");
		print(testA);
		printf(" by ");
		print(testB);
		printf(" = ");
		print(testC);
		printf("n");
	}
	return 0;
}

Program output:

More Info

Check out these links for more details:
Better Explained: A Visual, Intuitive Guide to Imaginary Numbers
Better Explained: Intuitive Arithmetic With Complex Numbers

Four Ways to Calculate Sine Without Trig

Is it possible to sin without trig? That is a question that has plagued theologians for centuries. As evil as trigonometry is, modern science shows us that yes, it is possible to sin without trig. Here are some ways that I’ve come across.

1 – Slope Iteration Method


The above image uses 1024 samples from 0 to 2pi to aproximate sine iteratively using it’s slope. Red is true sin, green is the aproximation, and yellow is where they overlap.

This method comes from Game Programming Gems 8, which you can snag a copy of from amazon below if you are interested. It’s mentioned in chapter 6.1 A Practical DSP Radio Effect (which is a really cool read btw!).
Amazon: Game Programming Gems 8

This method uses calculus but is actually pretty straightforward and intuitive – and very surprising to me that it works so well!

The derivative of sin(x) is cos(x). That means, for any x you plug into sin, the slope of the function at that point is the cosine value of that same x.

In other words, sin(0) is 0, but it has a slope of cos(0) which is 1. Since slope is rise over run (change in y / change in x) that means that at sin(0), if you take an infinitely small step forward on x, you need to take the same sized step on y. That will get you to the next value of sine.

Let’s test that out!
sin(0) = 0 so we start at (0,0) on the graph.

If we try a step size of 0.1, our approximation is:
sin(0.1) = 0.1

The actual value according to google is 0.09983341664. so our error was about 0.0002. That is actually pretty close!

How about 0.25?
sin(0.25) = 0.25

The real value is 0.24740395925, so our error is about 0.003. We have about 10 times the error that we did at 0.1.

what if we try it with 0.5?
sin(0.5) = 0.5

The actual value is 0.4794255386, which leaves us with an error of about 0.02. Our error is 100 times as much as it was at 0.1. As you can see, the error increases as our step size gets larger.

If we wanted to, we could get the slope (cosine value) at the new x value and take another step. we could continue doing this to get our sine approximation, knowing that the smaller the step that we use, the more accurate our sine approximation would be.

We haven’t actually won anything at this point though, because we are just using cosine to take approximated steps through sine values. We are paying the cost of calculating cosine, so we might as well just calculate sine directly instead.

Well luckily, cosine has a similar relationship with sine; the derivative of cos(x) is -sin(x).

Now, we can use cosine values to step through sine values, and use those same sine values to step through cosine values.

Since we know that cos(0) = 1.0 and sin(0) = 0.0, we can start at an angle of 0 with those values and we can iteratively step through the values.

Here is a sample function in C++

// theta is the angle we want the sine value of.
// higher resolution means smaller step size AKA more accuracy but higher computational cost.
// I used a resolution value of 1024 in the image at the top of this section.
float SineApproximation (float theta, float resolution)
{
    // calculate the stepDelta for our angle.
    // resolution is the number of samples we calculate from 0 to 2pi radians
    const float TwoPi = 6.28318530718f;
    const float stepDelta = (TwoPi / resolution);

    // initialize our starting values
    float angle = 0.0;
    float vcos = 1.0;
    float vsin = 0.0;

    // while we are less than our desired angle
    while(angle < theta) {

        // calculate our step size on the y axis for our step size on the x axis.
        float vcosscaled = vcos * stepDelta;
        float vsinscaled = vsin * stepDelta;

        // take a step on the x axis
        angle += stepDelta;

        // take a step on the y axis
        vsin += vcosscaled;
        vcos -= vsinscaled;
    }

    // return the value we calculated
    return vsin;
}

Note that the higher the angle you use, the more the error rate accumulates. One way to help this would be to make sure that theta was between 0 and 2pi, or you could even just calculate between 0 and pi/2 and mirror / reflect the values for the other quadrants.

This function is quite a bit of work to calculate a single value of sine but it’s real power comes in the fact that it’s iterative. If you ever find yourself in a situation where you need progressive values of sine, and have some fixed angle step size through the sine values, this algorithm just needs to do a couple multiplies and adds to get to the next value of sine.

One great use of this could be in DSP / audio synthesis, for sine wave generation. Another good use could be in efficiently evaluating trigonometry based splines (a future topic I plan to make a post about!).

You can see this in action in this shadertoy or look below at the screenshots:
Shadertoy: Sin without Trig

64 Samples – Red is true sine, green is our approximation, and yellow is where they are the same

128 Samples

256 Samples

1024 Samples

2 – Taylor Series Method

Another way to calculate sine is by using an infinite Taylor series. Thanks to my friend Yuval for showing me this method.

You can get the Taylor series for sine by typing “series sin(x)” into wolfram alpha. You can see that here: Wolfram Alpha: series sin(x).

Wolfram alpha says the series is: x-x^3/6+x^5/120-x^7/5040+x^9/362880-x^11/39916800 ….

what this means is that if you plug a value in for x, you will get an approximation of sine for that x value. It’s an infinite series, but you can do as few or as many terms as you want to be able to trade off speed for accuracy.

For instance check out these graphs.

Google: graph y = x, y = sin(x)

Google: graph y = x-x^3/6, y = sin(x)

Google: graph y = x-x^3/6+x^5/120, y = sin(x)

Google: graph y = x-x^3/6+x^5/120-x^7/5040, y = sin(x)

Google: graph y = x-x^3/6+x^5/120-x^7/5040+x^9/362880, y = sin(x)

Google: graph y = x-x^3/6+x^5/120-x^7/5040+x^9/362880-x^11/39916800, y = sin(x)
(Note that I had to zoom out a bit to show where it became inaccurate)

When looking at these graphs, you’ll probably notice that very early on, the approximation is pretty good between -Pi/2 and + Pi/2. I leveraged that by only using those values (with modulus) and mirroring them to be able to get a sine value of any angle with more accuracy.

When using just x-x^3/6, there was an obvious problem at the top and bottom of the sine waves:

When i boosted the equation with another term, bringing it to x-x^3/6+x^5/120, my sine approximation was much better:

You can see this method in action in this shadertoy:
Shadertoy: Sin without Trig II

3 – Smoothstep Method

The third method may be my favorite, due to it’s sheer elegance and simplicity. Thanks to P_Malin on shadertoy.com for sharing this one with me.

There’s a function in graphics called “smoothstep” that is used to take the hard linear edge off of things, and give it a smoother, more organic feel. You can read more about it here: Wikipedia: Smoothstep.

BTW if you haven’t read the last post, I talk about how smooth step is really just a 1d bezier curve with specific control points (One Dimensional Bezier Curves). Also, smoothstep is just this function: y = (3x^2 – 2x^3).

Anyhow, if you have a triangle wave that has values from 0 to 1 on the y axis, and put it through a smoothstep, then scale it to -1 to +1 on the y axis, you get a pretty darn good sine approximation out.

Here is a step by step recipe:

Step 1 – Make a triangle wave that has values on the y axis from 0 to 1

Step 2 – Put that triangle wave through smoothstep (3x^2 – 2x^3)

Step 3 – Scale the result to having values from -1 to +1 on the axis.

That is pretty good isn’t it?

You can see this in action in this shadertoy (thanks to shadertoy’s Dave_Hoskins for some help with improving the code):
Shadertoy: Sin without Trig III

After I made that shadertoy, IQ, the creator of shadertoy who is an amazing programmer and an impressive “demoscene” guy, said that he experimented with removing the error from that technique to try to get a better sin/cos aproximation.

You can see that here: Shadertoy: Sincos approximation

Also, I recommend checking out IQ’s website. He has a lot of interesting topics on there: IQ’s Website

4 – CORDIC Mathematics

This fourth way is a method that cheaper calculators use to calculate trig functions, and other things as well.

I haven’t taken a whole lot of time to understand the details, but it looks like it’s using imaginary numbers to rotate vectors iteratively, doing a binary search across the search space to find the desired values.

The benefit of this technique is that it can be implemented with VERY little hardware support.

The number of logic gates for the implementation of a CORDIC is roughly comparable to the number required for a multiplier as both require combinations of shifts and additions.
Wikipedia: Coordinate Rotation Digital Computer

Did I miss any?

If you know of something that I’ve left out, post a comment, I’d love to hear about it!

Feistel Networks – Do They Have to use XOR?

If you have no idea what a Feistel network is, but like cryptography and/or random number generation algorithms, read this link first:
Fast & Lightweight Random “Shuffle” Functionality – FIXED!

As a quick refresher, to encrypt data with a Feistel network, you break the plain text data into a left and a right side and do N rounds of this operation:

Left[i+1]  = Right[i];
Right[i+1] = Left[i] ^ RoundFunction(Right[i], key);

Where RoundFunction is ideally some chaotic function that returns some pseudo-random-esque number based on the inputs. For instance, RoundFunction could be MD5 so that it returned the MD5 hash of the data and the key, where the key could be considered the salt of the hash. The better the round function, the better your encryption algorithm will be.

To decrypt data with a Feistel network, you break the data into a left and right side and do the same number of rounds of this operation:

Right[i] = Left[i+1];
Left[i] = Right[i+1] ^ RoundFunction(Left[i+1], key);

Ok, onto the question….

Does it Have to use XOR?

Recently a friend of mine was using Feistel networks for something pretty amazing (so amazing, I can’t even talk about it), but in doing so, he asked me an interesting question. He asked “do you think this HAS to be XOR here, where we combine the round functions result back into the data?”. Well, it turns out, it doesn’t!

The operation has to be a reversible operation though, and you have to do the reverse operation when decrypting that you did while encrypting.

For instance, when encrypting you could add the round function result in, but then when decrypting, you would have to subtract the round function result out.

Or, you could do bitwise rotation left when encrypting, and right when decrypting perhaps.

Basically, anything that has a reverse operation can be used.

You have to be careful though because you might be lured into the trap of thinking that this includes something like multiplication and division.

If you multiply when you encrypt, you might get an integer overflow and lose data that can’t be corrected by doing a divide. For instance, if you multiply 255*2 in an unsigned 8 bit number you get 254 as a result. If you divide 254 by 2 to “undo” the multiplication, you get 127 which is obviously not 255, so we’ve lost some data. In an unsigned 8 bit number, ((255*2)/2) = 127.

If you go the other way and divide on encryption, and multiply on decryption, that doesn’t work either. For instance, when you divide 3 by 2, you get 1 with integer math, and when you multiply by 2, you get 2. So, with integers… ((3/2)*2) = 2.

Confusing note: you ARE able to do irreversible operations within the round function though. Feel free to do a divide or whatever you want in there. If that is difficult to understand how that could possibly work, you aren’t alone. Step through the code a bit by hand with a simple round function and a low number of rounds and you might be able to understand better how it does what it does.

I’m really not sure if anyone else out there does this variation on the traditional Feistel networks or not, but it is pretty interesting to combine the RoundFunction result back into the data with something other than XOR.

Source Code

Here’s some simple C++ code below to play with if you want to mess around with this stuff.

#include 
#include 
#include 

static const unsigned int c_numRounds = 4;

void PrimeRandomNumberPump ()
{
	// if you are curious about this, check out:
	// https://blog.demofox.org/2013/06/18/wtf-rand/
	srand((unsigned)time(NULL));
	for (unsigned int index = 0; index < 20; ++index)
		rand();
}

unsigned char RoundFunction (unsigned char value, unsigned char key)
{
	// Not a particularly effective round function, but the round function
	// isn't the point of this code.
	// If you want a better round function, try plugging in a hash function
	// or another chaotic function that has big changes in output for
	// small changes in input.  Also, you could change c_numRounds to a
	// higher number if you want better encryption.
	return value + key | (value * key) + 3;
}

void Encrypt (unsigned char &left, unsigned char &right, unsigned char key)
{
	for (unsigned int index = 0; index < c_numRounds; ++index)
	{
		// Feistel Network Encryption:
		//  Left[i+1]  = Right[i];
		//  Right[i+1] = Left[i] ^ RoundFunction(Right[i], key);

		// let's do addition to combine the value of the round function on 
		// encryption, instead of doing xor.  Xor is used in feistel networks
		// because xor is it's own inverse operation.
		unsigned char oldLeft = left;
		left = right;
		right = oldLeft + RoundFunction(right, key);
	}
}

void Decrypt (unsigned char &left, unsigned char &right, unsigned char key)
{
	for (unsigned int index = 0; index < c_numRounds; ++index)
	{
		// Feistel Network Decryption:
		//  Right[i] = Left[i+1];
		//  Left[i] = Right[i+1] ^ RoundFunction(Left[i+1], key);

		// let's do subtraction to combine the value of the round function on 
		// decryption, instead of doing xor.  Xor is used in feistel networks
		// because xor is it's own inverse operation.
		unsigned char oldRight = right;
		right = left;
		left = oldRight - RoundFunction(left, key);
	}
}

void DoTest (unsigned char plainText1, unsigned char plainText2, unsigned char key, int &tests, int &errors)
{
	// encrypt the plaintext
	unsigned char cipherText1 = plainText1;
	unsigned char cipherText2 = plainText2;
	Encrypt(cipherText1, cipherText2, key);

	// decrypt the cipher text
	unsigned char decryptedData1 = cipherText1;
	unsigned char decryptedData2 = cipherText2;
	Decrypt(decryptedData1, decryptedData2, key);

	// if the decrypted data doesn't match the plaintext data, count it as an error
	// and show the details
	tests++;
	if (decryptedData1 != plainText1 || decryptedData2 != plainText2)
	{
		errors++;
		printf("plaintext = 0x%02X%02Xrn", (unsigned int)plainText1, (unsigned int)plainText2);
		printf("ciphertext = 0x%02X%02Xrn", (unsigned int)cipherText1, (unsigned int)cipherText2);
		printf("decrypteddata = 0x%02X%02Xrnrn", (unsigned int)decryptedData1, (unsigned int)decryptedData2);
	}
}

void main (void)
{
	// generate a key
	PrimeRandomNumberPump();
	unsigned char key = (unsigned char)rand();

	// run tests with the key
	int errors = 0;
	int tests = 0;
	for (unsigned int y = 0; y < 256; ++y)
		for (unsigned int x = 0; x < 256; ++x)
			DoTest((unsigned char)y, (unsigned char)x, key, tests, errors);
		
	// display the test results
	printf("%i tests ran, %i errors encountered. key = 0x%02Xrn", tests, errors, key);
}

Soft Maximum vs Hard Maximum

The other day i stumbled on an interesting concept called a “Soft Maximum”.

If you think of the normal maximum, you might have something like this:

float maxValue = max(valueA, valueB);

if valueA and valueB come from functions, there’s usually going to be a sharp bend in the graph of the above where the maximum value changes from valueA to valueB or vice versa.

Sometimes, instead of a sharp bend, you would like a smooth transition between the two values – like when using this for graphics or advanced mathematics.

Here’s the formula for soft max:

double SoftMaximum(double x, double y)
{
	double maximum = max(x, y);
	double minimum = min(x, y);
	return maximum + log( 1.0 + exp(minimum - maximum) );
}

Here are 2 really interesting links on computing and using soft max:

Soft Maximum

How to Compute the Soft Maximum

Check out the images below for an example of when you might use this. This is from a shadertoy shader The Popular Shader. The first image is with using normal max, and the second image uses soft max.

softminOFF

softminON

How to Test Randomness of Numbers

At first i said the answer was to check this out: Diehard Battery of Tests of Randomness which is linked to by this page which may also be of interest: Tests for Random Number Generators.

But apparently that is the “old way” and there is a new program from NIST that you can get here: NIST test suite for random numbers, which subsequently is linked to from random.org: Random.org Statistical Analysis.

Getting that program from NIST to compile was a little bit of a chore for me on msvc 2010. The biggest hurdle i hit was that msvc 2010 doesnt have erf() and erfc() so i had to google “erf.cpp” and find an implementation. If you can’t find one, erf and erfc are part of gcc which is open sourced so you can always go that route if you need to!

After compiling, i was able to run the test on my numbers but couldn’t make much sense of the results very easily. There were a few p scores and presumably some chi squared scores somewhere, but the “summary file” was very cryptic (pun intended) so i wasn’t really sure…

Anyways, just wanted to put it here for myself and others if anyone’s looking for this in the future 😛

Thanks to my buddy James for the correction and links to the newer NIST program. Thanks man!

Bottom Line

Interestingly, the tests above use the source number data to do a bunch of different things, and then measure the statistics of the results.

For instance, it will use the numbers to shuffle a deck of cards, and then it will play poker and see if there is any bias of cards dealt, or players winning.

Or, it will use the numbers as the source of numbers for a roulette wheel and see if players win at the right rate statistically.

I guess the bottom line lesson for testing random numbers is that you should use the numbers how you intend to use them, and see if there’s any statistical anomalies.

There doesn’t seem to be a magic bullet test that works for generic randomness, but I’m guessing it’s just sort of… check for patterns in every way you can, or every way you care about, and if you don’t find any, consider it random. If you are using it for purposes where randomness really matters – like security or gambling – you then hope nobody else finds a pattern you didn’t! 😛

On that topic, check this out: Wikipedia: Michael Larson

External C++ Header Guards

A buddy at work said “I wish C++ had a two pass pre processor so that we could do external header guards”. It got me thinking about some random macro stuff i had seen before and i thought “hrm… you know, that actually might be possible to do… i’m going to give it a try”.

I ended up working something up tonight at home that’s semi-palatable. The way you use it is a little weird, but i think it satisfies the spirit of the challenge, and works as a proof of concept that you can do external header guards without having to type a bunch of stuff.

If you can think of a way to improve it, post a comment or something and let me know!!

Umm…External Header Guards? What are you Talking About??

Have you ever seen something like the below? it’s called a header guard:

#ifndef BLAH_H
#define BLAH_H

// code goes here

#endif BLAH_H

You might also have seen this variation:

#pragma once

Without the header guards, if you include the header file twice, it will complain that the classes etc have already been defined. Those make sure that doesn’t happen, by only including the contents of the file if it hasn’t already been included.

External header guards on the other hand would be guards in the place it’s included instead of in the header file itself. That is more typing (more work), but the benefit there is that the compiler doesn’t have to open the header file at all to see if it’s already been included, which could make for faster compile times in large projects.

Anyways, here’s the code:

Main.cpp

// include testoneblah.h which defines the typedef ProofIncluded__blah_h
// to prove that it was really included
#define FILESEQ (testone)(blah)
#include "Includer.h"

// try to include testoneblah.h again.  It won't get included again, and
// instead, ProofIncludeBlocked__blah_h will get typedef'd by includer.h
// to prove that the file did not get included.  Comment these lines out
// and you'll get a compiler error that ProofIncludeBlocked__blah_h is
// an undeclared identifier
#define FILESEQ (testone)(blah)
#include "Includer.h" 

int main(int argc, char **argv)
{
	ProofIncluded__blah_h a;
	ProofIncludeBlocked__blah_h b;
	return 0;
}

So it is a little weird… but to include a file, you define FILESEQ with a directory and filename (without .h on it), and then include “Includer.h”. Even though it’s weird to use, and doesn’t work for .inl files (and maybe other issues, some easily solved), it’s only one extra line of typing to do an external header guard, which is about as good as you can expect.

Ideally I wish the interface were like the below, but I haven’t been able to figure out how to make that work unfortunately.

IncludeFile_((testone)(blah))

Includer.h

//=====================================================================================================
// Rip off boost, hooray!!  boost_pp is really nice, you can just grab it from the boost bundle and
// start using it because it's just a bunch of includes.  you don't need to build or link with boost
// at all.  It's really nice.  http://www.boost.org/
//=====================================================================================================
# define BOOST_PP_EMPTY()
# define BOOST_PP_SEQ_ELEM(i, seq) BOOST_PP_SEQ_ELEM_I(i, seq)
# define BOOST_PP_SEQ_ELEM_I(i, seq) BOOST_PP_SEQ_ELEM_II((BOOST_PP_SEQ_ELEM_ ## i seq))
# define BOOST_PP_SEQ_ELEM_II(res) BOOST_PP_SEQ_ELEM_IV(BOOST_PP_SEQ_ELEM_III res)
# define BOOST_PP_SEQ_ELEM_III(x, _) x BOOST_PP_EMPTY()
# define BOOST_PP_SEQ_ELEM_IV(x) x

# define BOOST_PP_SEQ_ELEM_0(x) x, BOOST_PP_NIL
# define BOOST_PP_SEQ_ELEM_1(_) BOOST_PP_SEQ_ELEM_0
//=====================================================================================================
#define EB_COMBINETEXT(a, b) EB_COMBINETEXT_INTERNAL(a, b)
#define EB_COMBINETEXT_INTERNAL(a, b) a ## b
#define EB_TOTEXT(a) EB_TOTEXT_INTERNAL (a)
#define EB_TOTEXT_INTERNAL(a) #a
//=====================================================================================================

// extract the directory and file name
#define DIR BOOST_PP_SEQ_ELEM(0, FILESEQ)
#define FILE BOOST_PP_SEQ_ELEM(1, FILESEQ)

// create the full file name: /.h
#define THEFILENAME EB_TOTEXT(EB_COMBINETEXT(DIR, EB_COMBINETEXT(/, EB_COMBINETEXT(FILE, .h))))

#if !EB_COMBINETEXT(__, EB_COMBINETEXT(FILE, _h))
	//including file: YES
	//include the file
	#include THEFILENAME
#else
	//including file: NO
    //create a typedef called ProofIncludeBlocked___H to prove we didn't include the file
	typedef char EB_COMBINETEXT(ProofIncludeBlocked__, EB_COMBINETEXT(FILE, _h));
#endif

// clean up the things we created.  boost macros can stick around ::shrug::
#undef THEFILENAME
#undef DIR
#undef FILE

// defined by caller, but we're cleaning it up for convenience
#undef FILESEQ

I ripped some macros out of the boost preprocessor library (boost_pp) to help things out a little bit. In a nutshell what we are doing is this…

We test to see if the preprocessor value __<File>_h is not true (false, or undefined). If that is the case, we include the file <Directory>/<File>.h. Else, we define a typedef ProofIncludeBlocked<File>_h to prove that we blocked the include from happening.

Blah.h

#define __blah_h 1
class ProofIncluded__blah_h {};

Blah.h defines __blah_h as 1 (true). It’s important that it uses the same naming convention as Includer.h does (__<File>_h), otherwise this setup won’t work. If you screw it up, you’ll get compile errors about multiply defined symbols.

This file also defines a class ProofIncluded__blah_h to prove that this file actually got included, and also defines something that will complain if the file is included twice.

Issues

So, this is just a proof of concept and it has some issue including…

  • Duplicate file names – if you have the same file name in different folders, this setup has issues. It might be able to be helped by including the directory name into the header guard preprocessor symbol.
  • Referencing the same file different ways – if you reference the same file in different ways because there’s multiple ways to reach it in the include search paths, it won’t be able to tell that it’s the same file if you do the last fix. Maybe the real solution is to have another parameter defined specifying the header guard symbol? dont know…
  • Only supports .h files – It assumes a .h extension, but maybe another parameter could be the file extension to use so you could include .inl, .hpp etc.

Hopefully you find it interesting at least though (:

Is it Worth It?

Poday made some good points in the comments about it not being worth it, and my friend Doug also has this to say:

It’s not needed though all compilers already do that:

(GCC):
The GNU C preprocessor is programmed to notice when a header file uses this particular construct and handle it efficiently. If a header file is contained entirely in a `#ifndef’ conditional, then it records that fact. If a subsequent `#include’ specifies the same file, and the macro in the `#ifndef’ is already defined, then the file is entirely skipped, without even reading it.

(Clang)
The MultipleIncludeOpt class implements a really simple little state machine that is used to detect the standard “#ifndef XX / #define XX” idiom that people typically use to prevent multiple inclusion of headers. If a buffer uses this idiom and is subsequently #include‘d, the preprocessor can simply check to see whether the guarding condition is defined or not. If so, the preprocessor can completely ignore the include of the header.

(Clang still)
clang_isFileMultipleIncludeGuarded – Determine whether the given header is guarded against multiple inclusions, either with the conventional #ifndef/#define/#endif macro guards or with #pragma once.

For MSVC all I could find is Herb Sutter lead architect for MSVC and head of the C++ committee in his book ‘C++ Coding Standards: 101 Rules, Guidelines, and Best Practices’:
24. Always write internal #include guards. Never write external #include guards.
With a reason of:
Don’t try to be clever: Don’t put any code or comments before and after the guarded portion, and stick to the standard form as shown.
Today’s preprocessors can detect include guards, but they might have limited intelligence and expect the guard code to appear exactly at the beginning and end of the header.

Alloca and Realloc – Useful Tools, Not Ancient Relics

If you are a C/C++ programmer, you are likely familiar with malloc() and free(), the predecessors to C++’s new and delete operators, as well as the existence of the variations of malloc such as calloc, realloc and alloca.

If you are like me, you probably thought for a long while that malloc and it’s variations were relics of days gone by, only actually useful in a few very limited situations. Some of these guys still have use though, and don’t really have equivalents in C++ to replace them.

First the boring ones…
malloc – Allocates memory. Precursor to new operator.
calloc – Allocates memory and sets the contents to zero. C’s answer to the problem of uninitialized memory that constructors solve in C++.

Now the more interesting ones!

Alloca

Believe it or not, alloca actually allocates memory on the stack. When your function goes out of scope, the stack memory is automatically returned to the stack due to the nature of how the stack and stack pointer work. No need to free the memory allocated with alloca, and in fact if you tried, you’d probably get a crash 😛

If you are a programmer who writes high performance applications, you are probably familiar with the benefits of using the stack instead of allocating memory on the heap with new or malloc.

The benefits of using the stack include…

  • Quicker allocations – Allocating memory can be a relatively costly operation in terms of time, especially if you have multiple threads running using the same (thread safe) allocator. Allocating memory on the stack is essentially the same cost as defining a local variable. Under the hood, it’s just moving the stack pointer a little farther and gives you that empty space to use.
  • No memory leaks – when the function you’ve allocated the stack memory in exits, that memory is automatically freed. This is because the stack pointer just “moves back” to what it used to be. There is not really any memory to free.
  • Less memory fragmentation – When mixing large and small memory allocations and frees, sometimes you end up with your memory in a state where there is a lot of memory free, but just not all together in one place. For instance, your program might need to allocate 50MB, and there may be 300MB free on the heap total, but if there are small 16 byte allocations littered in the memory every 10MB, your program won’t be able to find a single 50MB region to allocate and the allocation will fail. One common cause of this problem is small allocations used for things like relatively small arrays or small string buffer allocations that exist temporarily to copy or transform some data, but are not meant to stick around very long. If you can put these on the stack instead of the heap, those small allocations don’t hit the heap, and your memory will be less fragmented in the end.
  • Increased performance (fewer cache misses) – the contents of the stack are likely already in the CPU cache, so putting your data there means less information for the CPU to have to gather from RAM which is a slow operation.

However, there are some dangers when allocating memory on the stack as well

  • If you keep a pointer to the memory, that memory could be “freed” and re-used, getting filled with other random data (local variables). That can cause crashes, memory corruption or other strange program behavior.
  • If you allocate too much on the stack you could run out of stack space. The stack isn’t really meant to hold large amounts of allocated data. You can adjust your programs stack size though if this is a route you want to pursue.

Alternatives

There are some common techniques I’ve seen people use in places that could have also used alloca instead. These include…

  • Small Pool Allocators – To get around the memory fragmentation problem, sometimes people will have different memory allocators based on the size of memory being allocated. This way, small temporary allocations for things like temporary string buffers will all be allocated from one place, while larger allocations for things like textures will be allocated elsewhere. This dramatically improves the memory fragmentation issue.
  • Object Pools – Object pools are similar to small pool allocators but they work by allocating some amount of memory for specific types of objects, and have a way to remember which objects are used and which ones are free. For instance, you may dynamically allocate an array of 100 SMyStruct objects and have a flag for each to know which ones are in use and which ones aren’t. This way, the program can ask for a new object, and it can find one currently not in use and return it to the caller without needing to hit the ACTUAL memory allocator to get the data (unless all objects are spoken for, at which point it can choose to fail, or allocate a new “page” of objects to be able to hand out). This also has an interesting side effect that cache misses can drop quite a bit since the same kinds of objects will be nearer to eachother in memory.
  • DIY Stack Allocator – When I was working at Midway, a friend (Hi Shawn!) profiled the animation code and found that a lot of time was spent in allocating temporary buffers to blend bone data together. To fix this, he rolled his own stack allocator, where there was one contiguous piece of memory on the heap that could be allocated from. There was an internal index keeping track of where the “top of the stack” was, and when memory was allocated, that stack index would just move up by however many bytes were asked for. At the end of the frame, the stack index was reset to zero, thus “freeing” the memory. This dramatically improved the animation system performance by making the temporary bone blend buffer allocations essentially free.
  • Thread Specific Memory – If you are having problems where multiple threads are trying to allocate memory at the same time, causing contention and slowdowns due to thread synchronization, another option is to give each thread it’s own chunk of memory and let it allocate from that. That way there is no contention and you won’t have the slowdown of thread synchronization due to memory allocation anymore. A problem here though can be figuring out how much memory each thread needs. One thread may need a lot of memory, and another thread may need none, and you may not have any way of knowing which in advance. In this case, you’d have to allocate “a lot” of memory for each thread in advance, and pay an extra cost in memory that you technically don’t have to. But hey, at least it’s fast, maybe the trade off is worth it in your situation!

Lastly, there’s another common trick to avoid dynamic allocations involving templates, check it out!

// define the CStaticArray class
template 
class CStaticArray
{
public:
  T m_values[N];

  // you could put functions in here to do operations on the array data to make it look more like a standard
  // data type, instead of a plain vanilla array
  unsigned int Count () { return N; }

  void SomeOtherFunction () { }
};

void MyFunc ()
{
  // make an array of 32 floats
  CStaticArray m_floatArray;

  // make an array of 128 SSomeStructs
  CStaticArray m_objectArray;

  for (unsigned int index = 0; index < m_objectArray.Count(); ++index)
  {
    m_objectArray.m_values[index].DoSomething();
  }
}

The above really shines if you have a standard API for strings or dynamic arrays in your code base. You can make a version like the above which works without dynamic allocations, but gives the same interface so it's easier for fellow programmers to use and swap in and out as needed.

Another nice benefit to the above technique is that it works for stack allocations, but you can also make them member variables of other objects. In this way, you can minimize dynamic allocations. Instead of having to dynamically allocate an object, and then dynamically allocate the array inside of it, you do a single allocation to get all the memory needed.

That is the closest thing in C++ that I've seen to alloca, but even so, alloca has the advantage that you can decide how much memory to allocate at run time. With the template method, you have to know at compile time which is fine for a lot of cases, but othertimes is a deal breaker, forcing you to have to go back to dynamic allocations (or perhaps now, alloca instead?)

Realloc

Realloc is the other interesting memory allocation function.

Like I was mentioning above, the fewer allocations you can do, the better off you are in terms of performance, and also memory fragmentation.

By writing smart containers (dynamic arrays, dynamic strings, etc) you can make it so when someone tries to make a container smaller, that instead of allocating new memory that’s smaller, copying the data over, and freeing the old memory, that instead it just remembers the new size but keeps the old, larger memory around.

Then later on, if the container was told to grow, if it was smaller than the larger size from the past, it could just use some of that old memory again.

However, if that container grows larger than it used to be, you are going to have to allocate, copy, and free (costly etc) to grow the container.

Down in the guts of your computer however, there may be memory right after the current memory that’s not being used by anything else. Wouldn’t it be great if you could just say “hey… use that memory too, i don’t want to reallocate!”.

Well, realloc does ALL of the above for you without you having to write special code.

When you realloc memory, you give the old pointer and the new size, and if it’s able to, it won’t do any allocations whatsoever, and will just return you your old pointer back to you. It may allocate the next memory block for you if the new size is larger, but would still return the old pointer value in this case. Or, if the new amount of memory is smaller, it may return you back the same memory without doing anything internally (it depends on your compiler’s specific implementation of realloc what it does when)

If realloc does have to allocate new memory though, it will copy over all the old data to the new memory that it returns to you and free the old memory. So, you don’t have to CARE whether the pointer returned is old or new, just store the return value and continue on with your life.

It’s pretty cool and can help reduce actual memory allocations, lowering memory fragmentation and increasing performance.

Is pre-increment really faster than post increment? Part 2

In the first part of this blog post (Is pre-increment really faster than post increment? Part 1) I showed that it really doesn’t seem to matter if you use post or pre increment with simple integer types.

I then promised an example of where the choice of pre or post increment DOES matter, and here it is.

The long and the short of it is this…

  • When you pre increment a variable, you are changing the value of a variable, and the new value will be used for whatever code the pre increment is part of.
  • When you post increment a variable, you are changing the value of a variable, but the OLD value is used for whatever code the post increment is part of.

To make post increment work that way, it essentially needs to make a copy of the variable before the increment and return the copy for use by the code of which the post increment is a part of.

A pre increment has no such need, it modifies the value and everything uses the same variable. There is no copy needed.

The compiler / optimizer apparently does a good job of figuring out when it does or does not need to make a copy of integral types, but it doesn’t do as well with complex objects. Here’s some sample code to demonstrate this and the output that it generates in both debug and release.

Test Code

#include 
#include 
#include 

//=============================================================
class CScopeMessage
{
public:
	CScopeMessage(const char *label)
	{
		PrintIndent();
		printf("%srn", label);
		s_indentLevel++;
	}

	CScopeMessage(const char *label, int objectId)
	{
		PrintIndent();
		printf("%s obj %irn", label, objectId);
		s_indentLevel++;
	}

	CScopeMessage(const char *label, int objectId, int copyObjectId)
	{
		PrintIndent();
		printf("%s (obj %i) to obj %irn", label, objectId, copyObjectId);
		s_indentLevel++;
	}

	~CScopeMessage()
	{
		s_indentLevel--;
	}

	static void StartNewTest()
	{
		s_indentLevel = 0;
		s_lineNumber = 0;
		printf("rn");
	}

private:
	void PrintIndent()
	{
		s_lineNumber++;
		printf("%2i:  ", s_lineNumber);
		for(int index = 0; index < s_indentLevel; ++index)
			printf("  ");
	}

private:
	static int s_indentLevel;
	static int s_lineNumber;
};

int CScopeMessage::s_indentLevel = 0;
int CScopeMessage::s_lineNumber = 0;

//=============================================================
class CTestClass
{
public:
	CTestClass()
	{
		m_objectID = s_objectID;
		s_objectID++;

		// this is just noise in the test, but feel free to
		// comment out if you want to see for yourself
		//CScopeMessage msg("Constructing", m_objectID);
		m_value = new char[4];
		strcpy(m_value, "one");
	}

	CTestClass(const CTestClass &other)
	{
		m_objectID = s_objectID;
		s_objectID++;

		CScopeMessage msg("Copy Constructing", other.m_objectID, m_objectID);
		m_value = new char[strlen(other.m_value) + 1];
		strcpy(m_value, other.m_value);
	}

	~CTestClass()
	{
		CScopeMessage msg("Destroying", m_objectID);
		delete[] m_value;
	}

    // preincrement
	CTestClass &operator++()
	{
		CScopeMessage msg("Pre Increment", m_objectID);
		DoIncrement();
		return *this;
	}
 
	// postincrement
	CTestClass operator++(int)
	{
		CScopeMessage msg("Post Increment", m_objectID);
		CTestClass result(*this);
		DoIncrement();
		return result;
	}

	void DoIncrement()
	{
		CScopeMessage msg("Doing Increment", m_objectID);
	}

private:
	char *m_value;
	int m_objectID;

	static int s_objectID;
};

int CTestClass::s_objectID = 0;

//=============================================================
int main (int argc, char **argv)
{
	CTestClass test;
	{
		CScopeMessage msg("--Post Increment--");
		test++;
	}

	CScopeMessage::StartNewTest();
	{
		CScopeMessage msg("--Post Increment Assign--");
		CTestClass testB = test++;
	}

	CScopeMessage::StartNewTest();
	{
		CScopeMessage msg("--Pre Increment--");
		++test;
	}

	CScopeMessage::StartNewTest();
	{
		CScopeMessage msg("--Pre Increment Assign--");
		CTestClass testB = ++test;
	}

	system("pause");
	return 0;
}

Debug

Here’s the debug output:

prepostdebug

You can see that in the post increment operator, it calls the copy constructor not once but twice! The first copy constructor is called to create the “result” object, and the second copy constructor is called to return it by value to the caller.

CTestClass operator++(int)
{
	CScopeMessage msg("Post Increment", m_objectID);
	CTestClass result(*this);
	DoIncrement();
	return result;
}

Note that it can’t return the copy by reference because it’s a local variable. C++11’s “std::move” and xvalue type functionality is there to help with this stuff, but if you can’t use that tech yet, it isn’t much help hehe.

Interestingly, we can see that 2 copy constructors get called whether or not we assign the value returned or not.

On the pre-increment side, you can see that it only does a copy construction call if you assign the result. This is nice and is what we want. We don’t want extra object copies or memory allocations and deallocations.

Release

prepostrelease

Things are a little bit better in release, but not by much. The optimizer seems to have figured out that it doesn’t really need to do 2 object copies, since it only ever wants at most one REAL copy of the object, so it gets away with doing one object copy in both situations instead of two.

That’s an improvement, but still not as good as the pre-increment case which hasn’t visibly changed in release (not sure about the assembly of these guys, but if you look and find something interesting, post a comment!)

Summary

As always, you should check your own assembled code, or test your compiler with printf output like this post does to ensure you really know what your code is doing.

But… it seems like you might want to use pre-increment if you ever use increment operators for heavy weight objects (such as iterators), but if you want to use post increment for integral types, it ought to be fine.

That said, a lot of people say “just use pre-increment” because at worst it’ll be the same as a post increment, but at best it will be a lot more efficient.

You do whatever you want, so long as you are aware of the implications of going one way or the other 😛

A Super Tiny Random Number Generator

When I posted the last blog post about shuffling on the GameProgrammer.com mailing list, someone responded back with a super tiny random number generator that is actually pretty damn good. It is this:

x+=(x*x) | 5;

The high bit of X is the source of your random numbers, so if you want to generate an 8 bit random number, you have to call it 8 times. Apparently it passes a lot of tests for randomness really well and is a pretty high quality PRNG. Check this out for more info: http://www.woodmann.com/forum/showthread.php?3100-super-tiny-PRNG

You can start x at whatever you want, but it might take a few iterations to “warm up” especially if you start with a small seed (ie you may want to throw away the first 5-10 random bits it generates as they may not be that random). I adapted it into an example program below, along with some example output. I use time() to set the initial value of x.

#include 
#include 
#include 
#include 

// A super tiny prng
// http://www.woodmann.com/forum/showthread.php?3100-super-tiny-PRNG
//
unsigned int seed = 0;
unsigned int GenerateRandomBit()
{
	seed += (seed * seed) | 5;
	return seed & 0x80000000;
}

template 
void GenerateRandom(T& value)
{
	memset(&value, 0, sizeof(T));
	const unsigned int numBits = sizeof(T) * 8;
	unsigned int* dataPointer = (unsigned int *)&value;
	for (unsigned int index = 0; index < numBits; ++index)
	{
		if(GenerateRandomBit()) {
			unsigned int pointerIndex = index / 32;
			unsigned int mask = 1 << index % 32;
			dataPointer[pointerIndex] |= mask;
		}
	}
}

int main(int argc, char **argv)
{
	seed = (unsigned int)time(NULL);
	printf("seed = %urn", seed);

	printf("9 random uints...rn");

	for (unsigned int index = 0; index < 9; ++index)
	{
		unsigned int random;
		GenerateRandom(random);
		printf("%2u: %10u (%x)rn", index, random, random);
	}

	printf("3 random floats...rn");
	for (unsigned int index = 0; index < 3; ++index)
	{
		float f;
		GenerateRandom(f);
		printf("%2u: %f (%x)rn", index, f, *((unsigned int*)&f));
	}

	printf("8 random characters...rn");
	char text[8];
	GenerateRandom(text);
	for (unsigned int index = 0; index < 8; ++index)
	{
		printf("%2u: %crn", index, text[index]);
	}
	system("pause");
	return 0;
}

tinyprng1

tinyprng2

tinyprng3

tinyprng4