Programmatically Calculating GCD and LCM

I recently came across a really interesting technique for calculating GCD (greatest common divisor) and then found out you can use that to calculate LCM (least common multiple).

Greatest Common Divisor

The greatest common divisor of two numbers is the largest number that divides evenly into those numbers.

For instance the GCD of 12 and 15 is 3, the GCD of 30 and 20 is 10, and the GCD of 7 and 11 is 1.

You could calculate this with brute force – starting with 1 and counting up to the smaller number, keeping track of the largest number that divides evenly into both numbers – but for larger numbers, this technique could take a long time.

Luckily Euclid came up a better way in 300 BC!

Euclid’s algorithm to find the GCD of numbers A and B:

  1. If A and B are the same number, that number is the GCD
  2. Otherwise, subtract the smaller number from the larger number
  3. Goto 1

Pretty simple right? It’s not immediately intuitive why that works, but as an example let’s say that there’s a number that goes into A fives times, and goes into B three times. That same number must go into (A-B) two times.

Try it out on paper, think about it a bit, and check out the links at the end of this section (:

A refinement on that algorithm is to use remainder (modulus) instead of possibly having to do repeated subtraction to get the same result. For instance if you had the numbers 1015 and 2, you are going to have to subtract 2 from 1015 quite a few times before the 2 becomes the larger number.

Here’s the refined algorithm:

  1. If A and B are the same number, that number is the GCD
  2. Otherwise, set the larger number to be the remainder of the larger number divided by the smaller number
  3. Goto 1

And here’s the C++ code:

#include <stdio.h>
#include <algorithm>

unsigned int CalculateGCD (unsigned int smaller, unsigned int larger)
	// make sure A <= B before starting
	if (larger < smaller)
		std::swap(smaller, larger);

	// loop
	while (1)
		// if the remainder of larger / smaller is 0, they are the same
		// so return smaller as the GCD
		unsigned int remainder = larger % smaller;
		if (remainder == 0)
			return smaller;

		// otherwise, the new larger number is the old smaller number, and
		// the new smaller number is the remainder
		larger = smaller;
		smaller = remainder;

void WaitForEnter ()
	printf("nPress Enter to quit");

void main (void)
	// Get A
	printf("Greatest Common Devisor Calculator, using Euclid's algorithm!n");
	printf("nFirst number? ");
	unsigned int A = 0;
	if (scanf("%u",&A) == 0 || A == 0) {
		printf("nMust input a positive integer greater than 0!");

	// Get B
	printf("Second number? ");
	unsigned int B = 0;
	if (scanf("%u",&B) == 0 || B == 0) {
		printf("nMust input a positive integer greater than 0!");
	// show the result
	printf("nGCD(%u,%u) = %un", A, B, CalculateGCD(A,B));

	// wait for user to press enter

I found this stuff in Michael Abrash’s Graphics Programming Black Book Special Edition: Patient Coding, Faster Code.

That book is filled with amazing treasures of knowledge and interesting stories to boot. I highly suggest flipping to a couple random chapters and reading it a bit. Very cool stuff in there (:

You might also find these links interesting or useful!
Wikipedia: Greatest Common Divisor
Wikipedia: Euclidean Algorithm

I’m sure there’s a way to extend this algorithm to work for N numbers at a time instead of only 2 numbers. I’ll leave that as a fun exercise for you if you want to play with that 😛

Least Common Multiple

The least common multiple of two numbers is the smallest number that is evenly divisible by those numbers.

Kind of an ear full so some examples: The LCM of 3 and 4 is 12, the LCM of 1 and 7 is 7, the LCM of 20 and 35 is 140. Note that in the first two examples, the LCM is just the two numbers multiplied together, but in the 3rd example it isn’t (also an interesting thing of note is that the first 2 examples have a GCD of 1, while the 3rd example has a GCD of 5).

Well interestingly, calculating the LCM is super easy if you already know how to calculate the GCD. You just multiply the numbers together and divide by the GCD.

LCM(A,B) = (A*B) / GCD(A,B)

Interestingly though, GCD(A,B) divides evenly into both A and B and will result in an integer result. This means we can multiply by A or B after the division happens and get the exact same answer. More importantly though, it helps protect you against integer overflow in the A*B calculation. Using that knowledge the equation becomes this:

LCM(A,B) = (A / GCD(A,B))*B

Pretty neat! Here’s some C++ code that calculates LCM.

#include <stdio.h>
#include <algorithm>

unsigned int CalculateGCD (unsigned int smaller, unsigned int larger)
	// make sure A <= B before starting
	if (larger < smaller)
		std::swap(smaller, larger);

	// loop
	while (1)
		// if the remainder of larger / smaller is 0, they are the same
		// so return smaller as the GCD
		unsigned int remainder = larger % smaller;
		if (remainder == 0)
			return smaller;

		// otherwise, the new larger number is the old smaller number, and
		// the new smaller number is the remainder
		larger = smaller;
		smaller = remainder;

unsigned int CalculateLCM (unsigned int A, unsigned int B)
	// LCM(A,B) = (A/GCD(A,B))*B
	return (A/CalculateGCD(A,B))*B;

void WaitForEnter ()
	printf("nPress Enter to quit");

void main (void)
	// Get A
	printf("Least Common Multiple Calculator, using Euclid's algorithm for GCD!n");
	printf("nFirst number? ");
	unsigned int A = 0;
	if (scanf("%u",&A) == 0 || A == 0) {
		printf("nMust input a positive integer greater than 0!");

	// Get B
	printf("Second number? ");
	unsigned int B = 0;
	if (scanf("%u",&B) == 0 || B == 0) {
		printf("nMust input a positive integer greater than 0!");
	// show the result
	printf("nLCM(%u,%u) = %un", A, B, CalculateLCM(A,B));

	// wait for user to press enter

Extending this to N numbers could be an interesting thing to try too (:

Here’s tasty link about LCM: Wikipedia: Least Common Multiple

Compile Time GCD and LCM

I’ve just heard that a compile time GCD and LCM implementation has been recomended for the STL. Check out the link below, kinda neat!

Greatest Common Divisor and Least Common Multiple


Bresenham’s Drawing Algorithms

If you were asked to name a line drawing algorithm, chances are you would say Bresenham. I recently needed to write my own software line drawing algorithm (CPU and regular ram, not GPU and VRAM) and Bresenham was the first to come to mind for me as well.

Believe it or not, Jack Bresenham actually came up with 2 famous line drawing algorithms. One is a run length algorithm, and the other is a run slice algorithm. People are most often familiar with the run length algorithm, which I mainly talk about in this post, but there’s some information about the run slice algorithm and links to other Bresenham primitive rendering algorithms below as well!

An interesting thing about both line algorithms is that they involve only integer mathematics, so there is no round off error, epsilons, numerical instability or anything like that. You use only integers and simple math (no divisions even!), and get the exact right answer out. I think that is pretty cool.

Bresenham’s Run Length Line Algorithm Summarized

To help understand the code, I want to give a brief summarization of how the algorithm works at a high level.

The first step of the Bresenham line algorithm is to see if the line is longer on the X axis or Y axis. Whichever one it is longer on is the major axis, and the shorter one is the minor axis.

Then, starting at the starting point of the line, you loop across all pixels of the major axis, doing some math to decide for each pixel whether you need to move the pixel on the minor axis yet or not. Basically you keep track of the error – or the distance the pixel is away from the true position of the line on the minor axis – and if the error is greater than or equal to 1 pixel, you move the pixel on the minor axis and subtract one pixel from the error. If the error is smaller than 1 pixel, you keep it at the same value it was for the last pixel, and keep the error value the same, so that it can accumulate some more error for the next pixel and test again.

The left line is longer on the X axis, so the major axis is the X axis. The right line is longer on the Y axis, so the major axis is the Y axis. Notice that there is one pixel for each value along the major axis of each line, but repeated pixel values along the minor axis of each line.

If you want the deeper details about the algorithm or the math behind it, check out this link:
Wikipedia: Bresenham’s line algorithm

One Octant Code

Wikipedia says that many implementations implement the algorithm for a single octant (half of a quadrant) and then use coordinate conversions (flipping the x and y’s and/or negating them) to make the same logic work for any quadrant.

In that vein, here is a simple implementation for the first octant – where the major axis is X, and both X and Y are increasing from the starting point to the ending point.

	void DrawLine (int x1, int y1, int x2, int y2, unsigned int color)
		const int dx = x2 - x1;
		const int dy = y2 - y1;

		int Error = 2 * dy - dx;
		int y = y1;
		DrawPixel(x1, y1, color);
		for (int x = x1+1; x < x2; ++x)
			if (Error > 0)
				Error += 2 * dy - 2 * dx;
				Error += 2 * dy;

With the point of Bresenham being that it is efficient, let’s sacrifice readability a bit and get rid of those DrawPixel calls and perhaps over-micro-optimize a bit and make some consts for frequently used values.

	void DrawLine (unsigned int* pixels, int width, int x1, int y1, int x2, int y2, unsigned int color)
		const int dx = x2 - x1;
		const int dy = y2 - y1;
		const int dx2 = dx * 2;
		const int dy2 = dy * 2;
		const int dy2Mindx2 = dy2 - dx2;

		int Error = dy2 - dx;

		unsigned int* pixel = &amp;pixels[y1*width+x1];
		*pixel = color;
		for (int x = x2 - x1 - 1; x > 0; --x)
			if (Error > 0)
				pixel += width + 1;
				Error += dy2Mindx2;
				Error += dy2;
			*pixel = color;

Cool, now let’s make it work for any line.

Final Code

	// Generic Line Drawing
	// X and Y are flipped for Y maxor axis lines, but the pixel writes are handled correctly due to
	// minor and major axis pixel movement
	void DrawLineMajorAxis(
		unsigned int* pixel,
		int majorAxisPixelMovement,
		int minorAxisPixelMovement,
		int dx,
		int dy,
		unsigned int color
		// calculate some constants
		const int dx2 = dx * 2;
		const int dy2 = dy * 2;
		const int dy2Mindx2 = dy2 - dx2;

		// calculate the starting error value
		int Error = dy2 - dx;

		// draw the first pixel
		*pixel = color;

		// loop across the major axis
		while (dx--)
			// move on major axis and minor axis
			if (Error > 0)
				pixel += majorAxisPixelMovement + minorAxisPixelMovement;
				Error += dy2Mindx2;
			// move on major axis only
				pixel += majorAxisPixelMovement;
				Error += dy2;

			// draw the next pixel
			*pixel = color;

	// Specialized Line Drawing optimized for horizontal or vertical lines
	// X and Y are flipped for Y maxor axis lines, but the pixel writes are handled correctly due to
	// minor and major axis pixel movement
	void DrawLineSingleAxis(unsigned int* pixel, int majorAxisPixelMovement, int dx, unsigned int color)
		// draw the first pixel
		*pixel = color;

		// loop across the major axis and draw the rest of the pixels
		while (dx--)
			pixel += majorAxisPixelMovement;
			*pixel = color;

	// Draw an arbitrary line.  Assumes start and end point are within valid range
	// pixels is a pointer to where the pixels you want to draw to start aka (0,0)
	// pixelStride is the number of unsigned ints to get from one row of pixels to the next.
	// Usually, that is the same as the width of the image you are drawing to, but sometimes is not.
	void DrawLine(unsigned int* pixels, int pixelStride, int x1, int y1, int x2, int y2, unsigned int color)
		// calculate our deltas
		int dx = x2 - x1;
		int dy = y2 - y1;

		// if the X axis is the major axis
		if (abs(dx) >= abs(dy))
			// if x2 < x1, flip the points to have fewer special cases
			if (dx < 0)
				dx *= -1;
				dy *= -1;
				swap(x1, x2);
				swap(y1, y2);

			// get the address of the pixel at (x1,y1)
			unsigned int* startPixel = &amp;pixels[y1 * pixelStride + x1];

			// determine special cases
			if (dy > 0)
				DrawLineMajorAxis(startPixel, 1, pixelStride, dx, dy, color);
			else if (dy < 0)
				DrawLineMajorAxis(startPixel, 1, -pixelStride, dx, -dy, color);
				DrawLineSingleAxis(startPixel, 1, dx, color);
		// else the Y axis is the major axis
			// if y2 < y1, flip the points to have fewer special cases
			if (dy < 0)
				dx *= -1;
				dy *= -1;
				swap(x1, x2);
				swap(y1, y2);

			// get the address of the pixel at (x1,y1)
			unsigned int* startPixel = &amp;pixels[y1 * pixelStride + x1];

			// determine special cases
			if (dx > 0)
				DrawLineMajorAxis(startPixel, pixelStride, 1, dy, dx, color);
			else if (dx < 0)
				DrawLineMajorAxis(startPixel, pixelStride, -1, dy, -dx, color);
				DrawLineSingleAxis(startPixel, pixelStride, dy, color);

Using the above functions, I wrote a little test to draw lines from the center of an image towards the edges at 100 different angles, and also a horizontal and vertical line. Here’s the results:

Run Slice Algorithm

So, like I said above, there is another Bresenham line algorithm that is a run slice algorithm.

What that means is that instead of looping across the major axis figuring out if each pixel needs to move on the minor axis as well, instead it loops across the minor axis, figuring out how many pixels it needs to write for that row or column of pixels.

One benefit of that algorithm is that it is less loop iterations since it’s looping across the minor axis pixels instead of the major axis.

Here’s an interesting read about Michael Abrash and his encounters with and the run slice algorithm:

Michael Abrash’s Graphics Programming Black Book Special Edition: The Good, the Bad, and the Run-Sliced

Other Bresenham Algorithms

This was origionally just going to be about the line algorithms, and I was going to make up another post about the circle algorithm, but I changed my mind.

While doing some research to see if there was a Bresenham bezier curve algorithm I stumbled on the page below, which shows that yes, there is! There are quite a few other algorithms as well:

Many Bresenham algorithms with very short c++ implementations

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

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.


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.


#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
	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; }

	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(
		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(
		a.Dual() * ::cos(a.Real())

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

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

inline CDualNumber atan (const CDualNumber &a)
	return CDualNumber(
		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)
	TestTrig(PI * 0.25f);
	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!


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.


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.


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 ");
		printf(" by ");
		printf(" = ");
		SVector testA = {4.0f, 3.0f};
		SVector testB = {-12.0f, 5.0f};
		SVector testC = {0.0f, 0.0f};
		Rotate(testA, testB, testC);
		printf("Rotating ");
		printf(" by ");
		printf(" = ");
		SVector testA = {5.0f, 1.0f};
		SVector testB = {2.0f, 2.0f};
		SVector testC = {0.0f, 0.0f};
		Unrotate(testA, testB, testC);
		printf("Unrotating ");
		printf(" by ");
		printf(" = ");
	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 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!

One Dimensional Bezier Curves

I was recently looking at the formula for bezier curves:

Quadratic Bezier curve:
A * (1-T)^2 + B * 2 * (1-T) * T + C * T ^2

Cubic Bezier curve:

(more info available at Bezier Curves Part 2 (and Bezier Surfaces))

And I wondered… since you can have a Bezier curve in any dimension, what would happen if you made the control points (A,B,C,D) scalar? In other words, what would happen if you made bezier curves 1 dimensional?

Well it turns out something pretty interesting happens. If you replace T with x, you end up with f(x) functions, like the below:

Quadratic Bezier curve:
y = A * (1-x)^2 + B * 2 * (1-x) * x + C * x ^2

Cubic Bezier curve:
y = A*(1-x)^3+3*B*(1-x)^2*x+3*C*(1-x)*x^2+D*x^3

What makes that so interesting is that most math operations you may want to do on a bezier curve are a lot easier using y=f(x), instead of the parameterized formula Point = F(S,T).

For instance, try to calculate where a line segment intersects with a parameterized bezier curve, and then try it again with a quadratic equation. Or, try calculating the indefinite integral of the parameterized curve so that you can find the area underneath it (like, for Analytic Fog Density). Or, try to even find the given Y location that the curve has for a specific X value. These things are pretty difficult with a parameterized function, but a lot easier with the y=f(x) style function.

This ease of use does come at a price though. Specifically, you can’t move control points freely, you can only move them up and down and cannot move them left and right. If you are ok with that trade off, these 1 dimensional curves can be pretty powerful.

Below is an image of a 1 dimensional cubic bezier curve that has control points A = 0.5 B = 0.25 C = 0.75 D = 0.5. The function of this curve is y = 0.5 * (1-x)^3 + 0.25 * 3x(1-x)^2 + 0.75 * 3x^2(1-x) + 0.5 * x^3.

You can ask google to graph that for you to see that it is in fact the same: Google: graph y = 0.5 * (1-x)^3 + 0.25 * 3x(1-x)^2 + 0.75 * 3x^2(1-x) + 0.5 * x^3


Another benefit to these one dimensional bezier curves is that you could kind of use them as a “curve fitting” method. If you have some data that you wanted to approximate with a quadratic function, you could adjust the control points of a one dimensional quadratic Bezier curve to match your data set. If you need more control points to get a better aproximation of the data, you can just increase the degree of the bezier curve (check this out for more info on how to do that: Bezier Curves Part 2 (and Bezier Surfaces)).

Smoothstep as a Cubic 1d Bezier Curve

BIG THANKS to for telling me about this, this is pretty rad.

It’s kind of out of the scope of this post to talk about why smoothstep is awesome, but to give you strong evidence that it is, it’s a GLSL built in function. You may have also seen it used in the post I made about distance fields (Distance Field Textures), because one of it’s common uses is to make the edges of things look smoother. Here’s a wikipedia page on it as well if you want more info: Wikipedia: Smoothstep

Anyhow, I had no idea, but apparently the smoothstep equation is the same as if you take a 1d cubic bezier curve and make the first two control points 0.0, and the second two control points 1.0.

The equation for smoothstep is: y = 3*x^2 – 2*x^3

The equation for the bezier curve i mentioned is: y = 0*(1-x)^3+3*0*(1-x)^2*x+3*1*(1-x)*x^2+1*x^3

otherwise known as: y = 3*1*(1-x)*x^2+1*x^3

If you work it out, those are the same equations! Wolfram alpha can verify that for us even: Wolfram Alpha: does 3*x^2 – 2*x^3 = 3*1*(1-x)*x^2+1*x^3.

Kinda neat 😛

Moving Control Points on the X Axis

There’s a way you could fake moving control points on the X axis (left and right) if you really needed to. What I’m thinking is basically that you could scale X before you plug it into the equation.

For instance, if you moved the last control point to the left so that it was at 0.9 instead of 1.0, the space between the 3rd and 4th control point is now .23 instead of .33 on the x axis. You could simulate that by having some code like this:

if (x > 0.66)
  x = 0.66 + (x - 0.66) / 0.33 * 0.23

Basically, we are squishing the X values that are between 0.66 and 1.0 into 0.66 to 0.9. This is the x value we want to use, but we’d still plug the raw, unscaled x value into the function to get the y value for that x.

As another example, let’s say you moved the 3rd control point left from 0.66 to 0.5. In this situation, you would squish the X values that were between 0.33 and 0.66 into 0.33 to 0.5. HOWEVER, you would also need to EXPAND the values that were between 0.66 and 1.0 to be from 0.5 and 1.0. Since you only moved the 3rd control point left, you’d have to squish the section to the left of the control point, and expand the section to the right to make up the difference to keep the 4th control point in the same place. The code for converting X values is a little more complex, but not too bad.

What happens if you move the first control point left or right? Well, basically you have to expand or squish the first section, but you also need to add or subtract an offset for the x as well.

I’ll leave the last 2 conversions as an exercise for whoever might want to give this a shot 😛

Another complication that comes up with the above is, what if you tried to move the 3rd control point to the left, past the 2nd control point? Here are a couple ways I can think of off hand to deal with it, but there are probably other ways too, and the right way to deal with it depends on your needs.

  1. Don’t let them do it – If a control point tries to move past another control point, just prevent it from happening
  2. Switch the control points – If you move control point 3 to the left past control point 2, secretly have them start controling control point 2 as they drag it left. As far as the user is concerned, it’s doing what they want, but as far as we are concerned, the control points never actually crossed
  3. Move both – if you move control point 3 to the left past control point 2, take control point 2 along for the ride, keeping it to the left of control point 3

When allowing this fake x axis movement, it does complicate the math a bit, which might bite you depending on what you are doing with the curve. Intersecting a line segment with it for example is going to be more complex.

An alternative to this would be letting the control points move on the X axis by letting a user control a curve that controls the X axis location of the control points – hopefully this would happen behind the scenes and they would just move points in X & Y, not directly editing the curve that controls X position of control points. This is a step towards making the math simpler again, by modifying the bezier curve function, instead of requiring if statements and loops, but as far as all the possibly functions I can think of, moving one control point on the X axis is probably going to move other control points around. Also, it will probably change the shape of the graph. It might change in a reasonable way, or it might be totally unreasonable and not be a viable alternative. I’m not really sure, but maybe someday I’ll play around with it – or if you do, please post a comment back and let us know how it went for you!


Here are some links to experiment with these guys and see them in action:
HTML5 Interactive 1D Quadratic Bezier Demo
HTML5 Interactive 1D Cubic Bezier Demo
Shadertoy: Interactive 1D Quadratic Bezier Demo
Shadertoy: Interactive 1D Cubic Bezier Demo

What if My Equation DOESN’T Equal Zero??

Take a simple equation such as y = 2x. You can transform that into the equation 2x – y = 0, and then could write it as f(x,y) = 2x – y = 0.

Now we have some function of x and y that equals zero. If we plug in numbers that are valid points in the original equation y = 2x (in other words, coordinates where y is double x), we’ll get zero out as a result:

f(x,y) = 2x - y

f(1,2) = 2 * 1 - 2 = 0
f(2,4) = 2 * 2 - 4 = 0
f(0.5, 1.0) = 2 * 0.5 - 1.0 = 0

If we plug numbers in that don’t fit that pattern, we get non zero answers as a result:

f(x,y) = 2x - y

f(2,1) = 2 * 2 - 1 = 3
f(5,5) = 2 * 5 - 5 = 5
f(-10, 20) = 2 * -10 - 20 = -40

What do these numbers mean? Are they useful for anything? As it turns out, they do have meaning and they are useful! They are part of what is needed to calculate the distance from the point (x,y) to the closest point on the curve. The other part of what we need is called the “gradient vector”.

Tangent Vector & Gradient Vector

As you probably know, there is a tangent vector at every point on a curve (let’s ignore things like asymptotes for now), and the tangent vector basically points from one point on the curve to the next point on the curve. The tangent vector can also be seen as a slope of that point on the curve.

A gradient vector is another property of a point, but it is different than a tangent vector in a couple ways. First off, it is perpendicular to the tangent vector, so sticks out away from the curve, and could be seen as a normal vector. Secondly, every point in the graph has a gradient vector, even if that point isn’t part of the curve and is just floating by itself in empty space.

The interesting part here is that a gradient vector of a point not on the curve will actually lead you straight to the closest point on the curve to the point. This is because of a handy property where the line connecting a point to the closest point on the curve will be perpendicular to the curve (perpendicular to the tangent).

So, if we can calculate a gradient vector for a point, we have the vector that points towards the closest point on the curve. Technically it might be pointing away from the closest point on the curve, instead of towards it but that gets cleared up in the math.

Calculating Gradient Vector

The gradient vector is just a partial derivative for each of the variables of the equation. If you have no idea what I’m talking about, here are 2 options:

1) You can watch some videos and learn about it here: Khan Academy: Multivariable calculus
2) You can have wolfram alpha calculate it for you, by typing “gradient” before the function. Click this link to see what I mean: Wolfram Alpha: gradient x^2-y

Here are some example equations, along with their gradient vector, to make sure that you either understand what’s going on, or are properly entering your equations into wolfram alpha. Note that a gradient vector may not always be constant, in which case, it changes value for each point! If it is constant, it’s the same gradient vector across the entire graph space.

f(x,y) = 2x – y
gradient(x,y) = (2, -1)

cubic function:
f(x,y) = x^3 – 10y
gradient(x,y) = (3x^2, -10)

f(x,y) = sin(x) – y
gradient(x,y) = (cos(x), -1)

f(x,y) = x^2 + y^2 – 5
gradient(x,y) = (2x, 2y)

Calculating Distance

Ok, so now that we have the value of our equation at some point (x,y) and we have the ability to find our gradient vector at that point, we have everything we need to calculate the distance from the point to the closest point on the curve.

All you need to do is divide the absolute value of the equation value at that point by the length of the gradient vector at that point.

Example 1 – Line

Let’s start with something super simple… y = x, aka x – y = 0 aka f(x,y) = x – y.

the gradient of that function is (1, -1) at every point on the graph. The magnitude of the gradient vector is the square root of 2.

What is the distance from the line to point P(2,0)?

Well.. the value of the equation at that point f(2,0) is 2, and the absolute value of that is still 2.

Next we divide that value by the magnitude of the gradient to get 2 / square root(2), which is just square root (2). So, the distance is square root of 2.

Since the value of the function at that point was positive before taking an absolute value, that means that the gradient points from the line to the point, so you can flip the gradient around to get (-1,1), and that is the direction that the closest point on the line is from the point P(2,0). If we travel square root of 2 units from P(2,0) down the vector (-1,1) we will get to the point (1,1) which is on the graph. That is the closest point on the graph to our point.

Note that it works out in this case, but usually the gradient vector won’t be the right magnitude since it’s not normalized. It happened that in this case it is, by dumb luck, but it usually isn’t.

Example 2 – Another Line

Let’s try a more complex line. y = 2x aka 2x – y = 0 aka f(x,y) = 2x – y.

the gradient of that function is (2, -1) with a magnitude of square root of 5 at all points on the graph.

What is the distance from that line to point P(1,5)?

Well, the value of the function at (1,5) is -3. When you take the absolute value of that, it becomes 3.

So, the distance is 3 / square root (5).

In this case, since the value of the function at the point was negative before we took the absolute value, it means that the gradient points from the point to the graph, so the line y=2x is 3/square root(5) units of length away from the point P(1,5) and the direction to travel from the point to get to the line is the same direction as the gradient vector at that point aka (2,1).

Example 3 – Quadratic Function

For the last example lets do y=x^2 aka x^2 – y = 0 aka f(x,y) = x^2 – y.

The gradient of that function is (2x,-1) and since it isn’t constant at all points, we can’t get the magnitude of the gradient yet.

What is the distance from that function to point P(2,0)?

Well, the value of f(2,0) is 4, which is still 4 if you take the absolute value.

The gradient at P(2,0) is (4,-1) with a magnitude of square root of 5.

So, the distance from point P(2,0) to the curve f(x,y) = x^2 – y is 4/square root(5) and to get from the point to the curve you should travel that distance down the same direction as the vector (4,-1).


Ok ok… the other shoe dropping here is that this distance calculation is only a distance ESTIMATE, and not always the right answer.

This is because if the gradient is a constant or linear function, it will be correct, but higher order gradients don’t always follow straight lines from the point to the curve.

The good news though is that this estimation will be less than or equal to the actual distance. Maybe not as useful as a for sure distance calculation, but an upper bound estimate still has it’s uses (like in ray marching, or drawing vector graphics).

Basically, the straight line tells you the closest it could possibly be, but the actual gradient may take a longer, curving path, from the point to the curve.

More Info

Here’s the article that got me started down this path, trying to find the answer to this question, while trying to understand the main topic:
IQ: Distance Estimation

And here’s a link that really helped me understand what a gradient was, and what it was all about:
Vector Calculus: Understanding the Gradient

Khan academy videos about these topics:
Khan Academy: Multivariable calculus

Analytic Fog Density


There are a number of ways to implement the effect of fog with modern real time rendered graphics. This blog post will explain how to render fog that has varying density, based on a function of X,Y,Z location in space, like in the picture above.

Faked Fog

One way is to “fake it” and do something like set the color of a pixel on an object to be based on it’s height. For instance you might say that pixels with a y axis value above 15 are unfogged, pixels with y axis values between 15 and 10 progressively get more fogged as they get closer to 10, and pixels with y axis values less than 10 are completely fogged. That can make some fog that looks like this:


A strange side effect of doing that though, is if you go down “into” the fog, and look out of the fog, things that should be fogged won’t. For instance, looking up at a mountain from inside the fog, the mountain won’t be fogged at all, even though it should be because you are inside of the fog.

A better way to do it, if you intend for the camera to be able to go into the fog, is to calculate a fogging amount for a pixel based on how far away it is from the view point, and how dense the fog is between the view point and the destination point.

If you are doing ray based rendering, like ray tracing or ray marching, you might find yourself trying to find how much fog is between points that don’t involve the view point – like if you are calculating the reflection portion of a ray. In this case, you are just finding out how much fog there is between the point where the reflection happened and the closest intersection. You can consider the point of reflection as the “view point” for the purpose of fogging.

Sometimes, the entire scene might not be in fog. In this case, you have to find where the fog begins and ends, instead of the total distance between the view point and the destination point.

In any case, the first thing you need to do when fogging is figure out the point where the fog begins, and the point where the fog ends. Then, you can figure out how much fog there is based on how the fog density works.

Constant Density Fog


The simplest sort of fog is fog that has the same density all throughout it.

What you do in this case is just multiply the fog density by the distance spent in the fog to come up with a final fog value.

As an example, your fog density might be “0.04” and if you are fogging a pixel 10 units away, you multiply density by distance. FogAmount = 0.04 * 10.0 = 0.4.

Doing this, you know the pixel should be 40% fogged, so you interpolate the pixel’s color 40% towards the fog color. You should make sure to clamp the fog amount to be between 0 and 1 to avoid strange visual anomolies.

The image below shows a constant fog density of 0.04.


Here’s an image of the same constant density fog as viewed from inside the fog:


A problem with constant fog density though, is that if you view it from edge on, you’ll get a very noticeable hard edge where the fog begins, like you can see in the image below:


Linear Density Fog


With linear fog density, the fog gets denser linearly, the farther you go into the fog.

With a fog plane, you can get the density of the fog for a specified point by doing a standard “distance from plane to point” calculation and multiplying that by how much the fog density grows per unit of distance. If your plane is defined by A*x+B*y+C*y+D = 0, and your point is defined as X,Y,Z, you just do a dot product between the plane and the point, giving the point a W component of one.

In other words…

FogDensity(Point, Plane) = (Plane.NormX * Point.X + Plane.NormY * Point.Y + Plane.NormZ * Point.Z + Plane.D * 1.0) * FogGrowthFactor

Here’s a picture of linear fog with a fog growth factor of 0.01:


The same fog viewed from the inside:


And lastly, the fog viewed edge on to show that the “hard line” problem of linear fog is gone (dramatic difference isn’t it?!):


Analytic Fog Density – Integrals


Taking a couple steps further, you might want to use equations to define fog density with some function FogDensity = f(x,y,z,).

How could you possibly figure out how much fog there is between two given points when the density between them varies based on some random function?

One way would be to take multiple samples along the line segment between the view point and the destination point, and either calculate the fog amount in each section, or maybe average the densities you calculate and multiply the result by the total distance. You might have to take a lot of samples to make this look correct, causing low frame rate, or accepting low visual quality as a compromise.

If you look at the graphs for the previous fog types, you might notice that we are trying to find the area under the graphs between points A and B. For constant density fog, the shape is a rectangle, so we just multiply width (time in fog) by height (the constant fog density) to get the fog amount. For linear density fog, the shape is a trapezoid, so we use the trapezoid area formula which is height (in this case, the distance in the fog) times the sum of the base lengths (the fog densities at points A and B) divided by 2.

How can we get the area under the graph between A and B for an arbitrary formula though? Well, a way exists luckily, using integrals (thanks to my buddy “Danny The Physicist” for educating me on the basics of integrals!).

There’s a way to transform a formula to get an “indefinite integral”, which itself is also a formula. I won’t go into the details of how to do that, but you can easily get the indefinite integral of a function by typing it into Wolfram Alpha.

Once you have the indefinite integral (let’s call it G(x)) of the fog density formula (let’s call it F(x)), if you calculate G(B) – G(A), that will give you the area under the graph in F(X) between A and B. Yes, seriously, that gives us the area under the graph between our points, thus giving us the amount of fog that exists between the two points for an arbitrary fog density function!

Note that when you plug a value into the indefinite integral and get a number out, that number is called the definite integral.

Analytic Fog Density – Implementation Details

Now that the theory is worked out let’s talk about implementation details.

First off, coming from an additive audio synthesis type of angle, I figured I might have some good luck adding together sine waves of various frequencies and amplitudes, so I started with this:

sin(x*F) * A

F is a frequency multiplier that controls how long the sine wave is. A is an amplitude multiplier that controls how dense the fog gets max.

Next, I knew that I needed a fog density function that never goes below zero, because that would mean if you looked through a patch of negative fog density, it would make the other fog you were looking through be less dense. That is just weird, and doesn’t exist in reality (but maybe there is some interesting visual effect hiding in there somewhere??), so the formula evolved to this, making sure the function never went below zero:

(1 + sin(x*F)) * A

Plugging that equation into wolfram alpha, it says the indefinite integral is:

(x – (cos(x*F)) / F) * A

You can check that out here:
Wolfram Alpha: (1 + sin(x*F)) * A.

It’s also kind of fun to ask google to graph these functions so you can see what they do to help understand how they work. Here are the graphs for A = 0.01 and F = 0.6:
Fog Density: graph (1 + sin(x*0.6)) * 0.01
Indefinite Integral: graph (x – (cos(x*0.6)) / 0.6) * 0.01

So, if you have point A and B where the fogging begins and ends, you might think you can do this to get the right answer:
FogAmount = G(B.x) – G(A.x)

Nope! There’s a catch. That would work if A and B had no difference on the y or z axis, but since they probably do, you need to jump through some hoops. In essence, you need to stretch your answer across the entire length of the line segment between A and B.

To do that, firstly you need to get that fog amount down to unit length. You do that by modifying the formula like so:
FogAmount = (G(B.x) – G(A.x)) / (B.x – A.x)

This also has a secondary benefit of making it so that your fog amount is always positive (so long as your fog density formula F(X) can’t ever go negative!), which saves an abs() call. Making it always positive ensures that this works when viewing fog both from the left and the right.

Now that we have the fog amount down to unit length, we need to scale it to be the length of the line segment, which makes the formula into this:
FogAmount = (G(B.x) – G(A.x)) * Length(B-A)/(B.x – A.x)

That formula will now give you the correct fog amount.

But, one axis of fog wasn’t enough to look very good, so I wanted to make sure and do one sine wave on each axis. I used 0.01 amplitude for each axis, but for the X axis i used a frequency of 0.6, for the Y axis i used a frequency of 1.2 and for the Z axis i used a frequency of 0.9.

Also, I wanted to give a little bit of baseline fog, so I added some constant density fog in as well, with a constant density of 0.1.

As a bonus, I also gave each axis a “movement factor” that made the sine waves move over time. X axis had a factor of 2.0, Y axis had a factor of 1.4 and Z axis had a factor of 2.2.

Putting all of this together, here is the final fog equation (GLSL pixel shader code) for finding the fog amount between any two points at a specific point in time:

float DefiniteIntegral (in float x, in float amplitude, in float frequency, in float motionFactor)
	// Fog density on an axis:
	// (1 + sin(x*F)) * A
	// indefinite integral:
	// (x - cos(F * x)/F) * A
	// ... plus a constant (but when subtracting, the constant disappears)
	x += iGlobalTime * motionFactor;
	return (x - cos(frequency * x)/ frequency) * amplitude;

float AreaUnderCurveUnitLength (in float a, in float b, in float amplitude, in float frequency, in float motionFactor)
	// we calculate the definite integral at a and b and get the area under the curve
	// but we are only doing it on one axis, so the "width" of our area bounding shape is
	// not correct.  So, we divide it by the length from a to b so that the area is as
	// if the length is 1 (normalized... also this has the effect of making sure it's positive
	// so it works from left OR right viewing).  The caller can then multiply the shape
	// by the actual length of the ray in the fog to "stretch" it across the ray like it
	// really is.
	return (DefiniteIntegral(a, amplitude, frequency, motionFactor) - DefiniteIntegral(b, amplitude, frequency, motionFactor)) / (a - b);

float FogAmount (in vec3 src, in vec3 dest)
	float len = length(dest - src);
	// calculate base fog amount (constant density over distance)	
	float amount = len * 0.1;
	// calculate definite integrals across axes to get moving fog adjustments
	float adjust = 0.0;
	adjust += AreaUnderCurveUnitLength(dest.x, src.x, 0.01, 0.6, 2.0);
	adjust += AreaUnderCurveUnitLength(dest.y, src.y, 0.01, 1.2, 1.4);
	adjust += AreaUnderCurveUnitLength(dest.z, src.z, 0.01, 0.9, 2.2);
	adjust *= len;
	// make sure and not go over 1 for fog amount!
	return min(amount+adjust, 1.0);

More Info

I ended up only using one sine wave per axis, but I think with more sine waves, or perhaps different functions entirely, you could get some more convincing looking fog.

At some point in the future, I’d like to play around with exponential fog density (instead of linear) where the exponential power is a parameter.

I also think that maybe squaring the sine waves could make them have sharper density changes perhaps…

One thing that bugs me in the above screenshots is the obvious “hard line” in both constant and linear fog where it seems fog crosses a threshold and gets a lot denser. I’m not really sure how to fix that yet. In traditional rasterized graphics you could put the fog amount on a curve, to give it a smoother transition, but in ray based rendering, that could make things a bit odd – like you could end up with an exponential curve butting up against the start of a different exponential curve (due to reflection or refraction or similar). The fog density would end up looking like log graph paper which would probably not look so great – although honestly I haven’t tried it to see yet!

If you have any questions, or feedback about improvements you know about or have discovered in any of the above, post a comment and let me know!

Here’s a good read on fog defined by a plane, that also gets into how to make branchless calculations for the fog amounts.
Unified Distance Formulas for Halfspace Fog

Interactive demo with GLSL pixel shader source code that you can also edit in real time with WebGL:

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.


static const unsigned int c_numRounds = 4;

void PrimeRandomNumberPump ()
	// if you are curious about this, check out:
	for (unsigned int index = 0; index < 20; ++index)

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
	if (decryptedData1 != plainText1 || decryptedData2 != plainText2)
		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
	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);

Bezier Curves Part 2 (and Bezier Surfaces)

This is a follow up post to Bezier Curves. My plan was to write a post about b-splines and nurbs next, but after looking into them deeper, I found out they aren’t going to work for my needs so I’m scratching that.

Here’s some basic info on b-splines and nurbs though before diving deeper into Bezier curves and surfaces.

B-Splines (Basis Splines)

Bezier curves are nice, but the more control points you add, the more complex the math gets because the degree of the curve function increases with each control point added. You can put multiple Bezier curves end to end to be able to have more intricate curves, but another option is to use B-Splines.

B-Splines are basically Bezier curves which let you specify more control points without raising the degree of the Bezier curve. They do this by having control points only affect part of the total curve.

This way, you could make a quadratic b-spline which had 10 control points. Only a few control points control any given point on the curve, so the curve stays quadratic (and so does the math), but you get a lot more control points. A “Knot Vector” is what controls which parts of the curve the control points control.

A Bezier curve is actually a special case of B-Spline where all control points affect the entire curve.

Nurbs (Non Uniform Rational B-Spline)

Sometimes when working with curves, you want some control points to be stronger that others. You can accomplish this in Bezier curves and B-splines by doubling up or trippling up control points in the same location to make that control point twice, or three times as strong respectively.

What if you want a control point to be 1.3 times stronger though? That gets a lot more complicated.

Nurbs solve that problem by letting you specify a weight per control point.

Just like Bezier curves are a special case of B-Splines, B-Splines are a special case of nurbs. A B-Spline could be thought of as nurbs that has the same weighting for all control points.

Back to Bezier!

My end goal is to find a curve / surface type that is flexible enough to be used to make a variety of shapes by artists, but is efficient at doing line segment tests against on the GPU. To this end, B-Splines and Nurbs add algorithmic and mathematical complexity over Bezier curves, and seem to be out of the running unless I can’t find anything more promising.

My best bet right now looks like a Bezier Triangle. Specifically, a quadratic Bezier triangle, where each side of the triangle is a quadratic Bezier curve that has 3 control points. When I get those details fully worked out, I’ll report back, but for now, here’s some interesting info I found about generalizing bezier curves both in order (linear, quadratic, cubic, quartic, etc) as well as in the number of dimensions (line, curve, triangle, tetrahedron, etc).

Bezier Generalized

I found the generalized equation on the wikipedia page for Bezier triangles and am super glad i found it, it is very cool!

I want to show you some specifics to explain the generalization by example.

Quadratic Curve:
(A * S + B * T) ^ 2

Expanding that gives you:
A^2 * S^2 + A * B * 2 * S * T + B^2 * T^2

In the above, S and T are Barycentric Coordinates in a 1 dimensional Simplex. Since we know that barycentric coordinates always add up to 1, we can replace S with (1-T) to get the below:

A^2 * (1-T)^2 + A * B * 2 * (1-T) * T + B^2 * T^2

Now, ignoring T and the constants, and only looking at A and B, we have 3 forms: A^2, AB and B^2. Those are our 3 control points! Let’s replace them with A,B and C to get the below:

A * (1-T)^2 + B * 2 * (1-T) * T + C * T ^2

And there we go, there’s the quadratic Bezier curve formula seen in the previous post.

Cubic Curve:
(A * S + B * T) ^ 3

To make a cubic curve, you just change the power from 2 to 3, that’s all! If you expand that equation, you get:

We can swap S with (1-T) to get:


Looking at A/B terms we see that there is more this time: A^3, A^2B, AB^2 and B^3. Those are our 4 control points that we can replace with A,B,C,D to get:

There is the cubic Bezier curve equation from the previous chapter.

Linear Curve:
(A * S + B * T) ^ 1

To expand that, we just throw away the exponent. After we replace S with (1-T) we get:
A * (1-T) + B * T

That is the formula for linear interpolation between 2 points – which you could think of as the 2 control points of the curve.

One more example before we can generalize.

Quadratic Bezier Triangle:
(A * S + B * T + C * U) ^ 2

If you expand that you get this:

Looking at combinations of A,B & C you have: A^2, AB, AC, B^2, BC, C^2. Once again, these are your control points, and their names tell you where they lie on the triangle. A Bezier triangle is a triangle where the 3 sides of the triangle are bezier curves. A quadratic bezier triangle has quadratic bezier curves for it’s edges which mean that each side has 3 control points. Those 3 control points are made up of the 3 corners of the triangle, and then 3 more control points, each one being between end points. A^2, B^2 and C^2 represent the 3 corners of the triangle. AB is the third control point for the bezier curve on the edge AB. BC and AC follow that pattern as well! Super easy to remember.

In a cubic Bezier triangle, you get a lot more control points, but a new class of control point too: ABC. This control point is in the middle of the triangle like the name would imply.

Anyways, in the expanded quadratic bezier triangle equation above, when you replace the control points with A,B,C for the triangle corner control points (the squares) and D,E,F for the inbetween control points, you get the bezier triangle equation below:


Note that we are dealing with a simplex in 3d now, so once again, instead of needing ALL Barycentric coordinates (S,T,U) we could pick one and replace it. For instance, we could replace U with (1-S-T) to have one less variable floating around.

All Done for Now

You can use this pattern to expand either in “surface dimension”, or in the dimension of adding more control points (and increasing the order of the equation). I love it because it’s super simple to remember that simple equation, and then just re-calculate the equation you need for whatever your specific usage case is.

If this stuff is confusing, check out the wiki page for Bezier Triangles, it has a great graphic that really shows you what I’m trying to explain:
Bezier Triangle

Next up I either want to make an HTML5 interactive app for messing around with Bezier triangles, or if I can figure out how to intersect a line segment with a quadratic Bezier triangle, i’ll probably just have some real cool looking screenshots to post along w/ the equation I ended up using (;

Special thanks to wolfram alpha for crunching some of these equations. Check it out, it’s really cool!
Wolfram Alpha – Cubic Bezier Curve Expansion

For more bezier fun check out my next Bezier post: One Dimensional Bezier Curves.