Rounding error in computing average

Question

I have some problems with rounding errors in c++. If I have to compute the average of two floats a and b, then why is it better to do a+0.5*(b-a) than (a+b)/2? I can't understand why there should be any difference in the two ways of computing it.

Answer 1

Your formula is correct in case if you are computing average of many numbers. In this case you can do the following:

μ_n = 1/nΣx_i

but here while adding 101st number you will need to add x₁₀₁ to μ₁₀₀, where μ₁₀₀ could be pretty big compared to x₁₀₁, and so you will lose some precision. To avoid this problem you can to like this:

μ₁₀₁ = μ₁₀₀ + 1/n(x₁₀₁ - μ₁₀₀)

This formula is much better if you x_i is of the same order of magnitude, as you avoid dealing with arithmetic operations between two big number and x_i.

You might want to read article Numerically stable computation of arithmetic means

Let's look into how numbers are represented in IEEE floating point. Consider C++ float:

Interval [1,2] goes with step 2^-23, so you can represent numbers 1+n*2^-23, where n belongs to {0, ..., 2²³}.

Interval [2^j, 2^j+1] goes with like [1,2] but multiplied by 2^j.

To, see how precision is lost you can run this program:

#include <iostream>
#include <iomanip>
int main() {
    float d = pow(2,-23);
    std::cout << d << std::endl;
    std::cout << std::setprecision(8) << d + 1 << std::endl;
    std::cout << std::setprecision(8) << d + 2 << std::endl; // the precision has been lost
    system("pause");
}

The output is

1.19209e-07
1.0000001
2