How to estimate floating-point precision of function?

Question

Let's say I have a function that consists solely of floating-point operations where the last operation rounds the computed value to a predefined number of digits. And I feed this function with a range of floats.

How do I conclude - just by looking at the decimals of the operand(s), and the structure of the code / sequence of floating-point calculations - that the round-up errors may or may not yield an unexpected result?

Do you know of any rule of thumb or mathematical methodology? Anything I can put to general use? I don't want to actually test the code.

Just one example for reference; range of input, required precision, and implementation of calc() are arbitrary:

Input: 0, 0.01, 0.02 ... 999.98, 999.99, 1000 (delta = 0.01)

Required precision: Rounded accurately to 2nd decimal

Pseudo code:

function main(a) {
  res = calc(a);
  return round(res); // Rounds res to desired decimals
}

function calc(a) {
  float b = 1.1;
  return a * b;
}

Update #1:

I found an interesting blogpost on the topic: Introduction to Scientific Computing: Error Propagation

It seems to contain the mathematical tools one needs to evaluate how "error prone" a computation is.

Update #2:

I worked out the following method to estimate the error propagation and will describe how I applied it on my use case.

It would be great if you comment on my thoughts. I'm not sure if I did it right. And thanks for your help so far! :)

Pseudo code:

// Floats conform to IEEE 547

function main() {
  float a,b; // Real value between 199 and 684
  a = some_value;
  b = another_value;
  res_a = calc(a);
  res_b = calc(b);
  res = res_a + res_b;
  return round(res); // Rounds to 2nd decimal
}

function calc(x) {
  float c = yet_another_value; // Real value between 1 and 1.4
  return c * x;
}

Estimation:

1. As I'm using IEEE 547 floats, the machine epsilon ε is 1.1e-16.
2. The relative error δ is 2ε - according to the article linked under update #1.
3. So the absolute error of my maximal input e_max won't exceed 3.8304e-13. *
4. That means the result will be precise up to the 12th decimal.
5. Since I just desire a precision of 2, the calculation does work precisely enough.

*) e_max = 2ε * res_max = 2 * 1.1e-16 * 2 * 684 * 1.4 = 3.8304e-13.

Bonus consideration:

Q: Would one add up the results of calc() infinitely, how many results would be "needed" approximately to contaminate the rounded total?

A: One would have to call calc() at least 4 trillion times. \o/

• e has to be at least 0.001 to force a round-up error
• The maximum result of calc() is x_max = 684 * 1.4 = 957.6
• e = δ * x => x = e / δ = e / 2ε = 1e-3 / (2 * 1.1e-16) = 1e+13 / 2.2
• n = x / x_max = 1e+13 / (2.2 * 957.6) ≈ 4e+9

Answer 1

Assuming that you're rounding from the thousandths place, the code in your example will always be accurate. Floating point numbers remain useful because they keep their imprecisions quite small relative to the most significant digit.

For instance, using the 8-byte IEEE754 standard (doubles in Java), 10.0/3.0 becomes3.3333333333333335.

You also asked for a rule of thumb to quickly determine whether your math will be accurate. Assuming that you are using the 8-byte standard, you can reasonably guess whether the calculation you are doing will be inaccurate by remembering that inaccuracies begin roughly 15-16 digits beyond the largest digit. If you're using the 4-byte standard (like you did in your example), you only get about 7-8 digits beyond the largest digit.

With that understanding, you can easily spot places you might run into trouble. If you're going to use a loop to multiply your number repeatedly, and then check whether the 1s digit is even or odd, you will likely run into inaccuracies.

One slightly surprising way to find trouble is to check simple equality. Consider the following Java code:

double x = 10.0;
double y = 3.0;
System.out.println( (1*x/(y*1)) == (10*x/(y*10)) );

It ultimately prints false, because it is really testing 3.3333333333333335 == 3.333333333333333. So, for equality testing, you will need to make some decisions about levels of precision, and run some rounding algorithm to obtain that.

Answer 2

There's no rule of thumb. There is analysis, and that analysis can be more or less elaborate, depending on what you want.

Assume IEEE 754 double precision arithmetic in default rounding mode (round to nearest even). If the highest bit of the mantissa has a value of $2^k$, that is if $2^k ≤ |x| < 2^{k+1}$ then the lowest bit of the mantissa has a value of $2^{k-52}$, the rounding error in each single operation will be at most $2^{k-53}$, and the relative error is less than $2^{-53}$.

For multiplication and division, the relative error of the result is limited by the relative errors of both inputs, plus an additional $2^{-53}$. For addition and subtraction, the absolute error of the result is limited by the relative errors of both inputs, plus an additional $2^{k-53}$, depending on the size of the result.

For a more precise analysis, consider that x-y is calculated with no rounding error at all if x/2 ≤ y ≤ 2x (the actual range is a bit larger). Relative errors are lower when a result is a little below a power of two. Multiplication / division by powers of two is exact, and so on.