Fixed point math in C

Lately I’ve been delving into C programming and more low level concepts. I was interested in how older games were made back when computers where basic and low spec. Why did PS1 graphics looked so clunky and how where 3D graphics generally made prior to the graphical APIs. One thing I found out was that these older engines have something called fixed point math.

What is fixed point math and why do we need it?

Fixed point math is calculating mathematical formulas using just integers and no floats. The main reason is that floating point operations are non-deterministic on different machines and on different range of CPUs. GPU power is measured in floating point operations because the GPU is specialized in dealing with floating points. But floating points have problems:

• They are non-deterministic – different CPU architectures will implement floating point operations a little bit differently and compilers might optimize code differently leading to discrepencies in calculations. This is usually a small margin error that won’t affect many people but there is benefit to having determinism. If you have determinism you can implement things like lock-step multiplayer.
• You can choose to implement your engine on a multitude of devices including micro controllers that do not have FPU
• Operations could be faster (or slower) – this is mostly based on developer skill and algorithms. The CPU will be able to handle most of the integer operations in its core though. With floating points some CPUs would have to deal with an FPU.

The main question now is how does one represent real numbers. Well it is pretty simple and also way more straightforward than what floating point uses. It has the weakness of having a certain precision locked in but unlike floats – precision doesn’t change with the higher the integer part grows.

Graphics Programming with Vulkan

You can get my course on Vulkan which covers the basics of graphics programming through the Vulkan specification. You can get the best price from here:

Enroll Now (14.99$)

Start from the notebook

To understant how the fixed point math will work in the code we first need to get back to the basics. What would happen if we set the limitation that you cannot have fractions?

I will start with the metric system. Imagine you have 1 meter. If you want to add another meter you will do 1m + 1m and poof you have 2m. But what you do if you want to add a half. With our limitation of not having half representation the only way would be to change the unit. This is why millimeters exist. We could define 1m = 1000mm and suddenly we can express that half a meter as 500mm. You could easily represent something like 1.5m without having a fraction by just saying we now have 1500mm.

$1m = 1000mm$

This is the priciple of everything that we will lay out in the current article. We are just scaling the number to a unit that can represent the fractions that we want. There are limitation with the accuracy you can get with this approach though. If we look at the meter example you will notice that even in decimal we have a problem of describing half a mm now. We have to scale even further to get bigger acurracy of everything lower than 1mm will get rounded to 0 and all in between values are now rounded to a mm value.

The next question we have to ask is how will this work on a computer. We can easily just multiply and divide by decimal point we want as accuracy. But you will find it is much easier to align the accuracy to a binary representation such that some of the operations are easily done with just left shifts and right shifts.

int32_t one = 1; // our regular number
int32_t one_scaled = one << 16; // our scaled number (similar to 1m = 1000mm)

In the example above you will notice that we left shift to scale by moving the one 16 bits to the left. This is actually equivalent to multiplying one by 65536 and producing this 0x00010000. This essentially moves our whole integer part into the upper 2 bytes of our 4 byte integer and leaves the lower 2 bytes to represent our sub scale fractions. You could look at it also as if you put the decimal point between those bytes and anything that is in the lower 2 bytes is a fractional part and upper 2 bytes is integer.

$\underbrace{0001}_\text{int}.\underbrace{0000_{16}}_\text{frac}$

This is a popular implementation of fracional integer also known as fixed 16 or 16.16 or fix16. You can represent a resolution of 0.000015 or in our example the lowest incerement with meters would be 15 micrometers. Everything that is in between 15 micrometers is going to be rounded out.

Here is a little calculator to experiment how the shifting affects the numbers (You can see a similar calculator for the floating point representation in the quick inverse square root article):

I hope that playing around with those numbers you can now get how we convert from integer to scaled integer (aka fixed point number).

Special cases

Well these examples were pretty simple but there are a few special cases we need to handle when we calculate these scaled integers. One of these is multiplication and division. Let’s get back to our metric system. I will start with division as it is easier to explain. You can easily divide 1000mm by half.

$1000mm / 2 = 500mm$

But what happens if you want to divide fixed point numbers by other fixed point numbers

$4000mm / 2000mm = 2units$

You will notice that we get back to a smaller scale. These 2 units are no longer in our scaled form. This is correct when talking for a metric system. Dividing mm by mm will produce a number but not a mm unit. This is not the case for general math though because you might want to be able to do the following:

$4000mm / 3000mm = 1unit \text{ (which is different than the expected 1.5)}$

What you need to do to get the correct math here is to actually scale the dividend in such a way that the fractional part is now part of the integral part. Imagine it as doing 4000mm = 40000000um and now dividing 4000000um by 3000mm. You will get as a result 1500 units. This is what we would expect our fractional representation of the result be in our fixed point format.

For code though scaling our 4B number by another 2B will actually lose our integer part. So to do this correctly we would have to convert the number from int32_t to int64_t first.

int32_t divident = 4 << 16;
int64_t divident_scaled = (int64_t)divident << 16;
int32_t divisor = 3 << 16;
int32_t result = divident_scaled / divisor; // 1.5 in fixed point format

Multiplication is similar but in reverse. Multiplying will produce a scaled number which we have to scale back down:

int32_t lhs = 4 << 16;
int32_t rhs = 3 << 16;
int64_t result_scaled = (int64_t)lhs * rhs;
int32_t result = result_scaled >> 16; // 12 in fixed point format

CMake Project Management

You can get my course on CMake where I cover configuring CMake projects and the CMake language in detail. You can get the best price from here:

Enroll Now (14.99$)

Implementation

Well you might’ve gotten so far how you could implement a library that allows you to work with fixed point numbers. But this article wouldn’t be complete without me showing you an implementaiton as well. The implementation for me is header only which means that all functions will be prefixed with static. I start by defining our two types and a few constants (using C23 by the way):

typedef int32_t fix16_t;
typedef int64_t fix48x16_t;

static constexpr int32_t fix16_bits = 32;
static constexpr int32_t fix16_int_bits = 16;
static constexpr int32_t fix48x16_int_bits = 48;
static constexpr int32_t fix16_frac_bits = fix16_bits - fix16_int_bits;
static constexpr uint32_t fix16_sign_mask = 0x8000'0000ul;
static constexpr uint32_t fix16_int_mask = 0xFFFF'0000ul;
static constexpr uint64_t fix48x16_int_mask = 0xFFFF'FFFF'FFFF'0000ull;
static constexpr uint32_t fix16_frac_mask = 0x0000'FFFFul;
static constexpr fix16_t fix16_one = 0x0001'0000;
static constexpr fix16_t fix16_neg_one = -fix16_one;
static constexpr fix16_t fix16_half = 0x0000'8000;
static constexpr fix48x16_t fix48x16_sign_mask = (fix48x16_t)0x8000'0000'0000'0000;
static constexpr fix16_t fix16_half = 0x0000'8000;
static constexpr fix16_t fix16_pi = 0x0003'243F;
static constexpr fix16_t fix16_half_pi = fix16_pi >> 1;
static constexpr fix16_t fix16_2_pi = fix16_pi << 1;

Now our first order of business is convertion between ints and floats and these fix16_t types:

static int16_t fix16_to_i16(fix16_t x)
{
    // fast conversion - just shift the number back losing all the fractional part
    return (int16_t)(x >> fix16_frac_bits);
}

static fix16_t fix16_from_i16(int16_t x)
{
    // just shift the number up
    return (fix16_t)((uint32_t)x << fix16_frac_bits);
}

static fix48x16_t fix48x16_from_i64(int64_t x)
{
    return (fix48x16_t)((uint64_t)x << fix16_frac_bits);
}

static float fix16_to_f32(fix16_t x)
{
    // Convert to a float and divide by our scale (fixed point 1 is 65536)
    return (float)x / fix16_one;
}

static fix16_t fix16_from_f32(float_t x)
{
    // Multiply our floating point number by our scale of 65536
    return (fix16_t)(x * fix16_one);
}

static fix16_t fix16_from_f64(double_t val)
{
    return fix16_from_f32((float)val);
}

I usually also implement a generic macro to convert more easily betwen those types instead of typing out the whole function name every time we need a fix16_t:

#define fix16(x) _Generic(x, \
  int8_t: fix16_from_i16(x), \
  uint8_t: fix16_from_i16(x), \
  int16_t: fix16_from_i16(x), \
  float: fix16_from_f32(x), \
  double: fix16_from_f64(x), \
  int64_t: fix48x16_from_i64(x), \
  int32_t: fix48x16_from_i64(x), \
  uint16_t: fix48x16_from_i64(x), \
  uint32_t: fix48x16_from_i64(x), \
  uint64_t: fix48x16_from_i64(x), \
  default: fix48x16_from_i64(x))

And then come the mathematical operations. Addition and subtraction is just that. Since the number is an int you do not need to do anything for addition and subtraction. You have to do something for multiplication and division though:

static fix48x16_t fix16_mul(fix48x16_t x, fix48x16_t y)
{
    // I have put fix48x16 on purpose so that the arguments get scaled automatically
    return ((x * y) >> fix16_frac_bits);
}

static fix48x16_t fix16_div(fix48x16_t x, fix48x16_t y)
{
    // I have put fix48x16 on purpose so that the arguments get scaled automatically
    return ((x << fix16_frac_bits) / y);
}

And the article wouldn’t be full if we don’t add a few utility functions from the math library:

static bool fix16_signbit(fix16_t x)
{
    return (x & fix16_sign_mask) > 0;
}

static fix16_t fix16_sign(fix16_t x)
{
    return fix16_signbit(x) ? fix16_neg_one : fix16_one;
}

static fix16_t fix16_abs(fix16_t x)
{
    return fix16_signbit(x) ? -x : x;
}

static fix16_t fix16_floor(fix16_t x)
{
    return x & fix16_int_mask;
}

static fix16_t fix16_ceil(fix16_t x)
{
    return fix16_floor(x) + ((x & fix16_frac_mask) > 0 ? fix16_one : 0);
}

static fix16_t fix16_round(fix16_t x)
{
    return fix16_floor(x) + (fix16_signbit(x) ? ((x & fix16_frac_mask) <= -fix16_half ? fix16_neg_one : 0) : ((x & fix16_frac_mask) >= fix16_half ? fix16_one : 0));
}

static fix16_t fix16_frac(fix16_t x)
{
    return fix16_abs(x) & fix16_frac_mask;
}

static fix16_t fix16_lerp(fix16_t x, fix16_t y, fix16_t t)
{
    return (fix16_t)fix16_mul(y - x, t) + x;
}

static fix48x16_t fix16_sq(fix48x16_t x)
{
    return fix16_mul(x, x);
}

static fix16_t fix16_min(fix16_t a, fix16_t b)
{
    return a < b ? a : b;
}

static fix16_t fix16_max(fix16_t a, fix16_t b)
{
    return a < b ? b : a;
}

static fix16_t fix16_clamp(fix16_t value, fix16_t min, fix16_t max)
{
    return fix16_max(fix16_min(value, max), min);
}

I also took some inspiration of some existing libraries for the trigonometry function implementations. You can always convert to float to use optimized implementations for those but as I mentioned with floats you would be loosing the determinism part of your math.

static fix16_t fix16_sin(fix16_t x)
{
    x = x % fix16_2_pi;
    if (x > fix16_pi) x -= fix16_2_pi;
    else if (x < -fix16_pi) x += fix16_2_pi;
    fix48x16_t x_sq = fix16_sq(x);

    fix48x16_t out = fix16_mul(-13, x_sq) + 546;
    out = fix16_mul(out, x_sq) - 10923;
    out = fix16_mul(out, x_sq) + fix16_one;
    out = fix16_mul(out, x);
    return (fix16_t)out;
}

static fix16_t fix16_cos(fix16_t x)
{
    return fix16_sin(x + fix16_half_pi);
}

static fix16_t fix16_sqrt(fix48x16_t x)
{
    // simple newton-rhapson usually produces an accurate result in a few iterations
    fix48x16_t num = ((fix48x16_sign_mask & x) > 0) ? -x : x;
    if (num == 0) return 0;
    fix48x16_t guess = fix16_mul(x + fix16_one, fix16_half); // x + 1 / 2
    guess = fix16_mul(guess + fix16_div(x, guess), fix16_half);
    guess = fix16_mul(guess + fix16_div(x, guess), fix16_half);
    guess = fix16_mul(guess + fix16_div(x, guess), fix16_half);
    guess = fix16_mul(guess + fix16_div(x, guess), fix16_half);
    guess = fix16_mul(guess + fix16_div(x, guess), fix16_half);
    guess = fix16_mul(guess + fix16_div(x, guess), fix16_half);
    guess = fix16_mul(guess + fix16_div(x, guess), fix16_half);
    return (fix16_t)guess;
}

static fix48x16_t fix16_rsqrt(fix48x16_t x)
{
    return fix16_div(fix16_one, fix16_sqrt(x));
}

static fix48x16_t fix16_hypot(fix48x16_t x, fix48x16_t y)
{
    return fix16_sqrt(fix16_sq(x) + fix16_sq(y));
}


static fix16_t fix16_atan2(fix16_t y, fix16_t x)
{
    constexpr fix16_t fix16_pi_div_four = 0x0000'C90F;
    constexpr fix16_t fix16_three_pi_div_four = 0x0002'5B2F;
    if (x == 0 && y == 0) return 0;
    fix16_t abs_y = fix16_abs(y);

    fix48x16_t r = 0;
    fix48x16_t r3 = 0;
    fix48x16_t angle = 0;

    if (x >= 0)
    {
        r = fix16_div(x - abs_y, x + abs_y);
        r3 = fix16_mul(fix16_mul(r, r), r);
        angle = fix16_mul(0x0000'3240, r3) - fix16_mul(0x0000'FB50, r) + fix16_pi_div_four;
    }
    else
    {
        r = fix16_div(x + abs_y, abs_y - x);
        r3 = fix16_mul(fix16_mul(r, r), r);
        angle = fix16_mul(0x0000'3240, r3) - fix16_mul(0x0000'FB50, r) + fix16_three_pi_div_four;
    }

    return (fix16_t)(y < 0 ? -angle : angle);
}

static fix16_t fix16_atan(fix16_t x)
{
    return fix16_atan2(x, fix16_one);
}

Note that this is not the fastest trigonometry function implementations and you can probably do better off with some cached values and look up tables. See how libfixmath does it for a faster implementation.

Usage

How would you then use this library

void main() {
    fix16_t number = fix16(16);
    fix16_t result = number / 15; // no need to convert or scale if you don't work with other fix16 on divisiont
    fix16_t result2 = fix16_div(number, fix16(15)); // if you have to divide another fix16 you have to use the specialized function
    // and so on
}

Common bugs

The most common problem with these fix16 numbers is the human factor or forgetting to convert them. You would get normal integer math and would have to debug why your fixed number is suddenly out of scale or why your fixed point loops are running from 1 to 100000.

This can easily be worked aroung by reworking the whole library to use a struct that wraps around an int32_t. But the pain then converts from doing human mistakes to being very explicit and calling a lot of very specific functions for math operations. Very much like you would do for a vector library implemented in C. A language like C++ or Rust might not have this problem since you can override the operators for your custom struct types and will get a clean convert between types.

Conclusion

I hope you found this article interesting and insightful. For me it took some time to understand the whole concept and I tried to relay it as simple as possible for other people out there.