Precisional Predicament: Floating point inaccuracies and Latent Bugs

Alright, alright! I know precisional is not a real word, but I couldn’t find a better alliteration with the same rhythm, so I decided to go with it!

Recently at work, I tracked down a quite frustrating and hard-to-catch bug caused by precision errors. It was, despite the frustation however, actually quite enjoyable. To dissect such a bug gives a tremendous amount of satisfaction when the cause is identified. Precision errors are well known, but this specific instance illustrates how elusive they can be when everything comes together with just the right values and rounding in the right places.

For quite a while, customers would very sporadically report a bug where an order was paid in full but didn’t close, and we could not figure out why. It took me almost ten hours of stepping through the code line by line, over and over, until I had identified the problem.

Floating point precision errors

Before going into details about the bug, it is important to understand why precision errors arise in the first place. Not unlike the literacy skills of certain male models, computers can’t math particularly good. The cause of this is how floating point numbers are represented in binary. While IEEE 754 ensures deterministic results for basic operations on the same architecture, cumulative rounding and differing order of operations due to library implementations or compiler optimizations can produce small divergences between systems.

Most are familiar with the trouble of representing a fraction such as $\frac{1}{3}$ in decimal numbers, where it will result in a repeating decimal. The same essentially happens with numbers represented in binary, but with different fractions than those that present a problem in decimal. For instance, simple decimal numbers like 0.1 or 0.2 cannot be represented exactly in binary floating point, which is why 0.1 + 0.2 != 0.3.

Binary representation of floating point values

To understand why this is, we have to know how floating point numbers are represented in binary. I will just use 32-bit single precision to illustrate the problem here, but the principle applies similarly to 64-bit double precision. The IEEE 754 standard requires the representation of a 32-bit signed float to have three components: A sign bit to indicate positive or negative values, eight bits representing the exponent and the mantissa¹ consisting of 23 bits. The exponent is stored with a bias of 127 meaning we will add this to the exponent².

To calculate the binary representation of a decimal number, we have to do the following steps.

1. First we convert the fractional part to binary by repeatedly multiplying by 2. If the result is ≥ 1, record a 1 and subtract 1 from the result. If < 1, record a 0 and repeat with the result. We continue until the fractional values start to repeat.
2. Next we normalize this by moving the decimal point until right after the first 1 in the sequence. The number of places we moved it is our exponent, which we then add 127 to (the bias) and convert to 8-bit binary.
3. The mantissa stores the 23 bits after the leading 1, which is implicit, and thus not stored. We now have our 32-bit representation by putting together the sign bit, 8-bit exponent, and the 23-bit mantissa.

Conversion example

As an example we can calculate 0.1 to be represented in binary as a float. We convert the fractional part to binary by starting with 0.1 and repeatedly multiplying by 2:

I continued further than necessary here, since it starts to repeat after six values have been computed. We can now take the integer values which are our binary representation $0.00011001100…_2$ and normalize it by moving the floating point and dropping the first 1, which gives us $10011001…_2$. This is our mantissa. The exponent is then -4, since we had to move four places. Adding the bias we get $-4+127 = 123$ which we then also convert to binary $01111011_2$. We can now put the three parts together

If we now calculate this back to base 10, we will see the problem. The formula for a floating point number is:

$$(-1)^{sign} \times 2^{exponent - bias} \times (1 + mantissa)$$

First we convert the exponent and mantissa back to base 10. The exponent is the same, but here we already notice a problem with the mantissa as we sum the fractions after the implicit 1. Each bit in the mantissa ($10011001100110011001100_2$) corresponds to a power of 2 in the fraction, where the first bit is the coefficient of $\frac{1}{2}$, the second bit is the coefficient of $\frac{1}{4}$, and so on: $$ \frac{1}{2} + \frac{0}{4} + \frac{0}{8} +\frac{1}{16} + \frac{1}{32} + \frac{0}{64} + \frac{0}{128} + \frac{1}{256} + … = 0.60000002384185791015625$$

We put it together with the rest of the formula

Putting it all together: $$(-1)^0 \times 2^{-4} \times (1+ 0.60000002384185791015625) = 0.100000001490116119384765625$$

So instead of 0.1, we get 0.100000001490116119384765625 and thus the problem arises.

Bug Hunt

Now that we understand how floating point precision errors can appear, let’s look at how they manifested in our system. As mentioned earlier, we had customers sporadically reporting that orders wouldn’t close despite being paid in full. I quickly suspected precision errors, since it sounded exactly like something it could cause. The code is somewhat convoluted and there are quite a lot of calculations being done to account for taxes, discounts of different types as well as currency conversions and whatnot with sometimes overzealous rounding added in. This made it hard to reason about, but after a while I had boiled it down to two values, which should be the same but on rare occasions aren’t.

The production system is written in Dart, but for illustrating the underlying problem, I have written a simplified example in C which mimics the behavior of the actual system, with everything stripped away except the core functionality.

Example code

#include <stdio.h>
#include <math.h>

float round_lowest_denomination(float amount)  {
    return rintf(amount * 100.0f) / 100.0f;
}

float get_remainder(float amount_paid, float order_total) {
    return order_total - amount_paid;
}

float amount_to_pay(float order_total) {
    return round_lowest_denomination(order_total);
}

int close_order(float remainder) {
    return (round_lowest_denomination(remainder) <= 0);
}

void run_example(char* info, float order_subtotal, float discount) {
    float order_total = order_subtotal - (order_subtotal * discount);
    float order_total_rounded = round_lowest_denomination(order_total);
    float paid_amount = amount_to_pay(order_total);
    float remainder = get_remainder(paid_amount, order_total);
    float rounded_remainder = round_lowest_denomination(remainder);
    int should_close_order =  close_order(remainder);

    printf("=== %s ===\n", info);
    printf("--- On screen values ---\n");
    printf("Order subtotal: $%.2f\n", order_subtotal);
    printf("Discount: %.0f%%\n", discount * 100);
    printf("Total: $%.2f\n", order_total_rounded);
    printf("Customer pays: $%.2f\n", paid_amount);

    printf("\n--- System internal values---\n");
    printf("Order total (stored value): %.8f\n", order_total);
    printf("Paid amount (stored value): %.8f\n", paid_amount);

    printf("\n--- Closing order ---\n");
    printf("Remainder: %.8f\n", remainder);
    printf("Remainder rounds to: %.2f\n", rounded_remainder);
    printf("Order closed: %s\n", should_close_order ? "true" : "false");
    printf("\n\n\n");

}

int main() {
    run_example("Example 1 - No precision error", 20.00, 0.0);
    run_example("Example 2 - Precision error but no bug ", 19.4, 0.0);
    run_example("Example 3 - Precision error causing bug", 19.89, 0.5);
    return 0;
}

There are essentially just four functions imitating the system, which are fairly straight forward to understand. The first is round_lowest_denomination(), which simply rounds the 8-decimal number to the lowest denomination possible in the given currency. This is two decimals for most currencies and also what I used in this example. Note that the rounding here uses rintf() for convenience³. The function get_remainder() just returns the difference between what has been paid and what the order total is. close_order() takes the remainder and returns true (1), if the remainder is equal to or lower than 0. Lastly, amount_to_pay() rounds the order total to two decimals, so the customer can pay. In reality, the value given to this function could be a partial payment, so not necessarily the full order total, but for simplicity, we assume the order is paid in full at once.

The run_example() function uses these functions to essentially calculate how much is to be paid and applies discounts if there are any. It then rounds the eight-decimal value to a payable amount of two decimals and stores this paid amount as eight decimals again. Then it checks if there is any remainder and if not, it closes the order.

Example output

Following are the outputs of the three examples, which illustrates the problem. In the first example there is no problem:

=== Example 1 - No precision error ===
--- On screen values ---
Order subtotal: $20.00
Discount: 0%
Total: $20.00
Customer pays: $20.00

--- System internal values---
Order total (stored value): 20.00000000
Paid amount (stored value): 20.00000000

--- Closing order ---
Remainder: 0.00000000
Remainder rounds to: 0.00
Order closed: true

Here there is no problem since 20.0 can be represented exactly in binary, so the order would close exactly as expected. In the second example there is a precision error present:

=== Example 2 - Precision error but no bug  ===
--- On screen values ---
Order subtotal: $19.40
Discount: 0%
Total: $19.40
Customer pays: $19.40

--- System internal values---
Order total (stored value): 19.39999962
Paid amount (stored value): 19.39999962

--- Closing order ---
Remainder: 0.00000000
Remainder rounds to: 0.00
Order closed: true

The total is a value, which cannot be represented exactly in binary, yet this does not cause any problems here, since the stored values of both total and paid amount are the same. The next example illustrates where the problem arises.

=== Example 3 - Precision error causing bug ===
--- On screen values ---
Order subtotal: $19.89
Discount: 50%
Total: $9.94
Customer pays: $9.94

--- System internal values---
Order total (stored value): 9.94499969
Paid amount (stored value): 9.93999958

--- Closing order ---
Remainder: 0.00500011
Remainder rounds to: 0.01
Order closed: false

The applied discount brings the total from \$19.89 to \$9.945 $(19.89 \times 0.5 = 9.945)$. This value cannot be represented precisely in binary, so the stored value ends up being slightly off at 9.94499969, but rounds to \$9.94, which the customer pays.

The critical issue is that when amount_to_pay() rounds 9.94499969 to display \$9.94, it creates a new 8-decimal float value, and 9.94 cannot be represented precisely in binary either, so this new value is stored as 9.93999958. Now we have two different binary approximations of what should be the same number: the original order total (9.94499969) and the rounded payment amount (9.93999958).

When we subtract these two imprecise numbers, we get a remainder of 0.00500011. This rounds up to \$0.01, so the system believes the customer still owes a penny despite having paid the exact amount displayed. The order then fails to close leaving the customer understandably confused.

So essentially, for the bug to trigger, the order total has to be a value that produces a precision error AND the rounded value has to also produce a precision error AND the difference between these two imprecise values has to be greater than or equal to 0.005. In example 3 above, the second value for the paid amount is slightly smaller than the precise value. Had it been slightly larger, the bug would not trigger.

All these conditions have to be in place for the bug, which of course happens very rarely, which is why tracking it down was not exactly easy.

Fixing it

Quick fix

I made a quite straight forward fix, which only required a single line changed. I simply made the remainder take the difference of the rounded numbers as such.

float get_remainder(float amount_paid, float order_total) {
    float rounded_order_total = round_lowest_denomination(order_total);
    float rounded_amount_paid = round_lowest_denomination(amount_paid);
    return rounded_order_total - rounded_amount_paid ;
}

This ensures that the order will close because the precision errors will never be so large, that the rounded values would be different. Theoretically I think it could happen if precision errors accumulated, but that would require an exceptionally large number of items being summed on a single order for the accumulated error to be large enough.

A Better Solution

A far better solution, but also far more extensive in terms of work, would be to not use floats at all to represent the small denominations. Instead everything should be counted as the smallest possible denomination in integers. So instead of having, say, \$1.25, we would have ¢125. This is a fairly common way to avoid these types of errors in financial systems, but it would basically require a rewrite of a large part of the code base to accommodate such a change. We are planning to implement it in a future iteration, where most of the code will be rewritten anyways, but for now, the quick fix will have to suffice.

If there is a lesson here, it is to be careful with floating point numbers, even if you believe you don’t need to worry because your use case doesn’t require high accuracy. These types of errors can still cause unwanted and hard to predict behavior as shown here.