Table of Contents

[Back to top]

1. Foreword & Introduction

This article presents my own approach of implementing a fixed point number class. It's not the best nor most efficient implementation, but it hopefully helps in understanding fixed point numbers. And of course, no guarantee for correctness... ;-)
After summarizing shortly the theoretical basics, the fixed point number format along with the most common mathematical operations are developed and discussed step by step. This tutorial is based on 16.16 fixed point numbers, but the explanations should be general enough to easily extended them to any fixed point number format.

[Back to top]

2. Motivation

Why another article on this topic you might ask?
Well, the goal is to have a very practical tutorial with examples for each operation to quickly grasp the topic. It shall not be over-complicated, still it should be not wrong at all and should provide the basics in order to have a nice starting point to dive into the fixed point number subject.

What's the use of fixed point numbers?
Some older or low-power microcontrollers do not have a floating point number unit, thus cannot calculate decimal numbers in hardware. So the handling of decimal numbers must be performed in software. As an example the first versions of J2ME (Java Platform, Micro Edition) used on mobile devices had no native support for decimal numbers. Also for some special use cases like in DSP, fixed point numbers could turn out more practical and efficient. And by the way, it's an interesting topic.

[Back to top]

3. Basics

Fixed points numbers are represented using integer values. As the name implies, the position of the point between the integer and fractional part is fixed - meaning the number of bits representing the integer part as well as the number of bits representing the fractional part is always the same. That's a huge difference to floating point formats like IEEE 754 where the number of bits used for the integer part and fractional parts differs from number to number, gaining a much better accuracy and wider range of presentable numbers.

In the following, the values NUM_INT_BITS and NUM_FRAC_BITS define the number of bits used for representing the integer respectively the fractional part of a decimal number.
The examples and implementation uses all a 16.16 representation, meaning:

• NUM_INT_BITS = 16
• NUM_FRAC_BITS = 16
• NUM_TOTAL_BITS = NUM_INT_BITS + NUM_FRAC_BITS = 32

Thus a single fixed point number is represented by a 32 bit value as follows:

Here the weight of each bit that contributes to the actual numerical value is explicitly noted. To have a shorter representation, all subsequent binary representations of fixed point numbers are shown without indices and weights:

The bit value at position i is denoted as b_i, where i is the index in range of {15, ... 1, 0, -1, ..., -16}.

[Back to top]

3.1 Ranges and representations

The integer part can contain 2^NUM_INT_BITS values. Because we want to represent also negative values, the two-complement form is used.
This means the integer as itself can take values of range { -2^NUM_INT_BITS, 2^NUM_INT_BITS - 1 }.
More general, because we use a NUM_TOTAL_BITS = 32 bit, the range is
{ 2^{NUM_TOTAL_BITS} / 2^{NUM_FRAC_BITS} - 1, -2^{NUM_TOTAL_BITS} / 2^{NUM_FRAC_BITS} } = { 2³² / 2¹⁶ - 1, -2³² / 2¹⁶ } = { 2¹⁶ -1, -2¹⁶ } = { 32767, -32767 }.

The fractional part consists of NUM_FRAC_BITS, with the weight ranging from 2^-1 to 2^-16.
So the biggest fractional smaller than the next integer value is 2^-1 + 2^-2 + ... + 2^{-NUM_FRAC_BITS}, with NUM_FRAC_BITS = 16 this result in ≈ 0,999985. This is by the way the same as 2 - 2^-16 / (2^-16 - 1) ≈ 0,999985.
The smallest possible difference between two numbers is given by 2^-16 ≈ 0,00001526.

In general, the value of a fixed-point decimal is given by the equation:

Let's have a look at some interesting and boundary values:

• The number zero is trivial - all bits are zero.
   
• The biggest number that can be represented is as follows:

  0

  1

  1

  1

  1

  1

  1

  1

  1

  1

  1

  1

  1

  1

  1

  1

  1

  1

  1

  1

  1

  1

  1

  1

  1

  1

  1

  1

  1

  1

  1

  1

  This (0x7FFFF) has all integer bits set except the MSB bit, so it's value is 2¹⁴ + ... + 2¹ + 2⁰ + 2^-1 + ... + + 2^-16 ≈ 32767,999985.
   
• The smallest number that can be represented is as follows:

  1

  0

  0

  0

  0

  0

  0

  0

  0

  0

  0

  0

  0

  0

  0

  0

  0

  0

  0

  0

  0

  0

  0

  0

  0

  0

  0

  0

  0

  0

  0

  0

  This value (0x8000) has only the MSB bit set, so it's value is -2¹⁵ = -32768.

[Back to top]

4. Creating fixed-point numbers

To work with fixed-point numbers, the first implementation step is to create fixed-point numbers with the desired value.

[Back to top]

4.1 From integer numbers

This is the easiest case. Having an integer value, just shift it the number of fractional bits to the left to skip the fractional part. This also works for negative values as integer values are stored the same way (in two-complement from) as the integer part of our fixed-point numbers.

[Back to top]

4.2 From floating-point numbers

At first, a constant SUNFP_ONE is introduced that represents the value 1 in fixed-point format. Using the previous chapter, this is just 1 << NUM_FRAC_BITS, thus 0x00010000. This constant is used to build and create a fractional value between 0 and 1 in fixed-point number - just multiply it with SUNFP_ONE!

Summarized, creating a fixed-point number from a floating-point number can be done in two steps. First convert the integer part, then the fractional part and finally add both values.

This can be shortened to following one-liner!

[Back to top]

4.3 From strings

Basic idea is to find the decimal point: the part in front of the decimal point is the integer part and this can be directly converted from string to int. Then the integer conversion from chapter 4.1.1 can be used.
Same applies to the case if no decimal point can be found in the string - then there is no fractional part at all, only an integer part.

The handling of the fractional part is a bit more complicated. Take the fractional part as itself and convert it to an integer. This is multiplied by SUNFP_ONE. This results in a fixed point number with the value of the fractional part - but that is too large by a power of ten! So count the number of characters in the fractional part and divide the intermediate FP by 10^{<number of characters>}. This gives the correct fractional part that is added to the integer part.

However, what about negative numbers? Consider the number "-4.324". In short, the calculation would be roughly as -4 + 324/1000 = -3,676!
To fix it, the integer part is subtracted by 1 and then "1 - fractional part" is added. In the case of "4.32" this would be -4 - 1 + ( 1 - 324/1000) = -5 + 676/100 = -4,324. Yes!
So here a possible implementation:

[Back to top]

5. Basic binary operations

In this chapter, following basic operations are discussed: addition, subtraction, multiplication and division.

[Back to top]

5.1 Addition

Addition is actually ridiculous easy - just adding two fixed-point numbers result in the correct value. This works as we use the same 16.16 format for all numbers. Of course, if the decimal point is not at the same position for both numbers, special care has to be taken. As normal integer addition is used, there is also the problem of overflowing, e.g 32767 + 1 = -32768.

[Back to top]

5.2 Subtraction

The same as for addition also applies for subtraction - just use integer subtraction, nothing interesting here.

[Back to top]

5.3 Multiplication

Let's see if two fixed point numbers can just be multiplied to get the correct value, similar to addition and subtraction:
The numbers 4 and 5 have the values in fixed point format 0x00040000 and 0x00050000. Calculating 0x00040000 * 0x00050000 results in 0x1400000000!
Two problems are obvious here:

• At first, the result is far too big. The expected result of 4 * 5 is 20 (0x14), that in fixed point format 0x00140000.
  So the calculated value 0x1400000000 is too large by factor 0x10000 = 2¹⁶.
• Secondly, this value does also not fit into our 16.16 format.

So first address the result issue itself: Because the result is too large by factor 2¹⁶, divide it by this value which is actually just shifting it by 16 (more specific: by NUM_FRAC_BITS) to the right.

Shifting the result by 16 leads to following possible implementation:

However, there is still the small issue that the intermediate result does not fit into our 16.16 fixed point format, therefore this intermediate result has to be cased to long before shifting.

A possible solution is to treat the integer and fractional parts of both operands separately as following examples demonstrates:

This can be implemented as follows:

[Back to top]

5.4 Division

After having discussed multiplication, it's not surprising that similar issues arise with division. Instead of dividing the result by 2¹⁶, it has to be multiplied by 2¹⁶, thus shifted to the left by 16. However, the shift must now occur before the actual division, otherwise precision is lost.

This means again that we need a long data type for implementation:

The only way that came into my mind to avoid the long data type is follows but unfortunately it introduces a small loss of precision. (If you know a better solution, I would be glad to be informed about it - feel free to contact me!):
Basic idea is to replace the division by multiplication with the reciprocal, i.e. instead of a / b, it's calculated a * (1 / b).
Well, nothing gained, there is still a division left. And the reciprocal 1 / a of value a is in fixed point format (1 << 32) / a.

Again a long is needed to hold the value 1 << 32...

Here comes the step that avoids the long however at which also precision is lost: It is shifted only by 31 bits during reciprocal calculation, then the multiplication is performed, and afterwards the last shift by one to the left is performed.

Let's compare both values with example.

Let's see the implementation without using long data type:

[Back to top]

6. Basic unary operations

In this chapter, following basic unary operations are discussed: floor, ceil, round and absolute.

[Back to top]

6.1 Floor

The floor function applied to a decimal number x gives the greatest integer number that is equal or less than x, e.g. floor(1.2) = 1 and floor (-2.4) = 3.
Due to the fact how this fixed point number is represented, the implementation is very easy - just mask out the fractional part. This also works for negative numbers as following examples make clear:

Decimal numbers that are actually integer values are the trivial case, because in fixed-point format all fractional bits are zero, thus masking them out does not change the actual value, so e.g. floor(-12) = -12.

Here's the implementation:

[Back to top]

6.2 Ceil

The ceil function applied to a decimal number x gives the greatest integer number that is equal or greater than x, e.g. ceil(1.2) = 2 and floor (-2.4) = -2.
Ceil and floor are tightly coupled, so an approach is to implement the ceil function using the floor function: Ceil(x) = Floor(x+1). However this works not for plain integer numbers, as the result is then off by one, so this case has to be treated specially:

• If x has a decimal part (no plain integer), then Ceil(x) = Floor(x+1).
• If x has no decimal part (plain integer), then Ceil(x) = x.

Here's a possible implementation:

[Back to top]

6.3 Round

The round function applied to a decimal number x gives the nearest integer number to x, e.g. round(1.2) = 1, round(1.7) = 2 and round(-2.4) = -2.
Again the floor function can be used to implement the round function as follows:

• If x >= 0, then round(x) = floor(x + 0.5).
• If x < 0, then round(x) = floor(x - 0.5).

Due to the way negative numbers are represented in this fixed-point format, always the first implementation case is required, thus round(x) = floor(x + 0.5).

Note that the rounding of value .5 is ambiguous, e.g. the distance from 2.5 to 2 and from 2.5 to 3 is the same, thus the rounding of values with decimal part .5 has to be defined. This implementation returns the larger value in this case.

Here's the implementation:

[Back to top]

6.4 Absolute

The abs function applied to a decimal number x gives the positive value of x independent of the sign of x, e.g. abs(4.5) = 4.5 and abs(-4.5) = 4.5.
Actually this can be easily implemented by handling two distinct cases:

• If x >= 0, just return x.
• If x < 0, then subtract x from zero, thus abs(x) = 0 - x.

[Back to top]

7. Trigonometric functions

In this chapter, following trigonometric functions are discussed: Sine, cosine and tangent.

[Back to top]

7.1 Sine

To calculate the sine of decimal value x, an approximation using Taylor series is used:

The more terms are used, the better the accurary of the result will be. Take shall be taken of overflow: Already 9! = 362880 does not fit into the 16bit integer part. However x⁹/9! can also be written as x⁹/(7! * 72) and 7! still fits into a 16bit number.
Note that the input value must be in range [-PI, PI] - the transformation of any value to that range is easy, but not part of following implementation:

[Back to top]

7.2 Cosine

To calculate the cosine of decimal value x, again an approximation using Taylor series is used:

Also here the input value must be in range [-PI, PI] . As the implementation is very similar to sine, the listing is left out.

[Back to top]

7.3 Tangent

The easiest calculation of tangent is using the equation tan(x) = sin(x) / cos(x) because the functions for sine and cosine already exist:

[Back to top]

8. Square root, exponentiation and logarithm

[Back to top]

8.1 Square root

To calculate the square root, Newton's method for approximating the root of any function is used. That is for any function f we want to find x so that f(x) = 0 (see also [4]).
This method starts with an initial guess and calculates iteratively a new guess that is closer to the actual value than the previous value. If we have an approximation of the result value x_n, then a better approximation x_n+1 is calculated as:

We want to find the result of sqrt(x). That is the same as to find the root of the function f(x) = x² - a = 0 because the result is sqrt(a).
The derivation of f(x) is f'(x) = 2*x. So let's insert f(x) and f'(x) into Newton's approximation formula:

x_n+1 = x - (x² - a) / 2x
= ((x * 2x) / 2x) - ( (x² - a) / 2x)
= ((2x²) / 2x) - ( (x² - a) / 2x)
= (x² + a) / 2x
= 0.5x + 0.5(a + x)
= 0.5 (x + a/x)

Two small things are missing:

• First we need an initial 'guess value' for x₀. For simplicity we just use 1. However the better the initial start value is, the faster the approximation approaches the target value.
• Secondly, a termination condition is required to abort the iterative approximation. It's unknown how fast the approximation works, thus two conditions are chosen in my sample implementation: First a maximum of ten iterations are allowed, then the algorithm terminates at the latest. Secondly, the algorithm terminates earlier if a good approximation is found. At each step there is a guess x_n that shall be the result of sqrt(x) - so if x_n² is very near to x, the algorithm can be stopped as the required accuracy is reached. Here's a sample implementation:

[Back to top]

8.2 Logarithm

Two different approaches are implemented to calculate the natural logarithm (to base e) - note that there are many more possibilities to approximate the logarithm. The logarithm to any other base than e can be computed by base changing, see e.g. [3].

[Back to top]

8.2.1 Logarithm using Arctanh

The logarithm can be calculated using the inverse hyperbolic tangent according to following equation:

where x must be in range [0, +inf). The inverse hyperbolic tangent itself can be approximated using following Taylor series:

The listing of the implementation is left out, have a look inside the source code.

[Back to top]

8.2.2 Logarithm using Taylor series

The logarithm can be approximated using following Taylor series:

This is only valid for input range (-1, 1). Therefore following steps are performed to handle also the general case:

• The logarithm of x, log(x), shall be calculated instead of log(x+1). Thus the input value is subtracted by 1 - this gives us the input range (0, 2).
• To handle all input values, in range [0, +inf), scaling is used. So the input value is scaled in a way that one part falls in range (0, 2) which can be calculated using above Taylor series, the other part is pre-calculated. The fact log(a * b) = log(a) + log(b) is used.

  Example: Calculate log(25)

  log(25) = log(25/16 * 16) = log(25/16) + log(16).
  log(16) is precalculated, while log(25/16) = log(1.5625) is calculated using the Taylor series.

  The logarithms of all 2^x are pre-calculated and stored in a look-up table. Because the largest integer fixed point value is 2¹⁵ = 32768, the lookup table has entries for x in range [0, 15].

The appropriate divider from the 2^x lookup table is searched, then the input value is subtracted by one, then the logarithm is actually computed. Here's the implementation:

[Back to top]

8.3 Exponentiation

Exponentiation to base e is calculated with following equation:

Here's the implementation:

[Back to top]

9. Conclusion & References

That's it! Those were the notes of my own try to implement a fixed-point class. Hope it was interesting for you and have fun!

Sunshine, April 2017

References

History

• 2017-04-18: First version
• 2021-03-13: Corrected the Taylor series for cosine from cos = x - (...) to cos = 1 - (...).
• 2022-07-03: Added definition of SUNFP_ONE, changed FP_ONE to SUNFP_ONE for consistency, fixed some typos.