# Robust computation of the mean of two numbers in floating-point?

Let x, y be two floating-point numbers. What's the right way to compute their mean?

The naive way (x+y)/2 can result in overflows when x and y are too large. I think 0.5 * x + 0.5 * y maybe better, but it involves two multiplications (which maybe is inefficient), and I am not sure if it is good enough. Is there a better way?

Another idea I've been playing with is (y/2)(1 + x/y) if x<=y. But again, I am not sure how to analyze this and prove that it satisfies my requirements.

Moreover, I need a guarantee that the computed mean will be >= min(x,y) and <= max(x,y). As pointed out in Don Hatch's answer, maybe a better way of posing this question is: What's an implementation of the mean of two numbers that always gives the most possible accurate result? That is, if x and y are floating-point numbers, how to compute the floating-point number closest to (x+y)/2? In this case, the computed mean is automatically >= min(x,y) and <= max(x,y). See Don Hatch's answer for details.

Note: My priority is robust accuracy. Efficiency is expendable. However, if there are many robust and accurate algorithms, I would pick the most efficient.

• (+1) Interesting question, surprisingly non-trivial. – Kirill Aug 4 '15 at 22:23
• In the past, floating point values were computed and held in a higher precision form for intermediate results. If a+b (64-bit doubles) produces a 80 bit intermediate result and this is what is divided by 2, you don't have to worry about overflow. Loss of precision is less obvious. – JDługosz Aug 5 '15 at 0:41
• The solution to this seems relatively simple (I added an answer). The thing is I am a programmer and not a computer sciences expert, so what am I missing that makes this question so much more difficult? – IQAndreas Aug 5 '15 at 9:25
• Don't worry about the cost of multiplications and divisions by two; your compiler will optimize them for you. – Federico Poloni Aug 5 '15 at 12:58

I think Higham's Accuracy and Stability of Numerical Algorithms addresses how one can analyze these types of problems. See Chapter 2, especially exercise 2.8.

In this answer I'd like to point out something that isn't really addressed in Higham's book (it doesn't seem to be very widely known, for that matter). If you are interested in proving properties of simple numerical algorithms such as these, you can use the power of modern SMT solvers (Satisfiability Modulo Theories), such as z3, using a package such as sbv in Haskell. This is somewhat easier than using pencil and paper.

Suppose I'm given that $0\leq x\leq y$, and I'd like to know if $z=(x+y)/2$ satisfies $x\leq z\leq y$. The following Haskell code

import Data.SBV

test1 :: (SFloat -> SFloat -> SFloat) -> Symbolic SBool
test1 fun =
do [x, y] <- sFloats ["x", "y"]
constrain $bnot (isInfiniteFP x) &&& bnot (isInfiniteFP y) constrain$ 0 .<= x &&& x .<= y
let z = fun x y
return $x .<= z &&& z .<= y test2 :: (SFloat -> SFloat -> SFloat) -> Symbolic SBool test2 fun = do [x, y] <- sFloats ["x", "y"] constrain$ bnot (isInfiniteFP x) &&& bnot (isInfiniteFP y)
constrain $x .<= y let z = fun x y return$ x .<= z &&& z .<= y


will let me do this automatically. Here test1 fun is the proposition that $x \leq \mathit{fun}(x,y) \leq y$ for all finite floats $x,y$ with $0\leq x\leq y$.

λ> prove $test1 (\x y -> (x + y) / 2) Falsifiable. Counter-example: x = 2.3089316e36 :: Float y = 3.379786e38 :: Float  It overflows. Suppose I now take your other formula:$z=x/2+y/2$λ> prove$ test1 (\x y -> x/2 + y/2)
Falsifiable. Counter-example:
x = 2.3509886e-38 :: Float
y = 2.3509886e-38 :: Float


Doesn't work (due to gradual underflow: $(x/2)\times2 \neq x$, which might be unintuitive due to all arithmetic being base-2).

Now try $z=x + (y-x)/2$:

λ> prove $test1 (\x y -> x + (y-x)/2) Q.E.D.  Works! The Q.E.D. is a proof that the test1 property holds for all floats as defined above. What about the same, but restricted to$x\leq y$(instead of$0\leq x\leq y$)? λ> prove$ test2 (\x y -> x + (y-x)/2)
Falsifiable. Counter-example:
x = -3.1300826e34 :: Float
y = 3.402721e38 :: Float


Okay, so if $y-x$ overflows, how about $z = x + (y/2-x/2)$?

λ> prove $test2 (\x y -> x + (y/2 - x/2)) Q.E.D.  So it seems that among the formulas I've tried here,$x + (y/2 - x/2)$seems to work (with a proof, too). The SMT solver approach seems to me a much quicker way of answering suspicions about simple floating-point formulas than going through floating-point error analysis with pencil and paper. Finally, the goal of accuracy and stability is often at odds with the goal of performance. For performance, I don't really see how you can do better than$(x+y)/2$, especially since the compiler will still do the heavy lifting of translating this into machine instructions for you. P.S. This is all with single-precision IEEE754 floating-point arithmetic. I checked$x \leq x + (y/2-x/2) \leq y$with double-precision arithmetic (replace SFloat with SDouble), and it works too. P.P.S. One thing to bear in mind when implementing this in code is that compiler flags like -ffast-math (some forms of such flags are sometimes turned on by default in some common compilers) will not result in IEEE754 arithmetic, which will invalidate the above proofs. If you do use flags that enable, e.g., associative addition optimizations, then there's no point in doing anything other than$(x+y)/2$. P.P.P.S. I got carried away a little looking only at simple algebraic expressions without conditionals. Don Hatch's formula is strictly better. • Hold on; did you claim that if x<=y (regardless of whether x>=0 or not) then x+(y/2-x/2) is a good way of doing it? Seems to me that can't be right, since it gives the wrong answer in the following case when the answer is exactly representable: x=-1, y=1+2^-52 (the smallest representable number greater than 1), in which case the answer is 2^-53. Confirmation in python: >>> x = -1.; y = 1.+2.**-52; print 2**-53, (x+y)/2., x+(y/2.-x/2.) – Don Hatch Aug 5 '15 at 4:00 • @DonHatch I aimed to answer the question about "robustness": ensuring that inequalities$x\leq (x+y)/2\leq y$continue to hold and nothing under-/overflows. Importantly, the example you give is unstable in the large-condition-number sense: small relative perturbations in$x,y$cause large relative changes in the result, i.e., catastrophic cancellation. Other than that you're right:$(x+y)/2$will always produce the correctly rounded result with small relative error. I myself would prefer$(x+y)/2$. – Kirill Aug 5 '15 at 4:14 First, observe that if you have a method that gives a most accurate answer in all cases, then it will satisfy your required condition. (Note that I say a most accurate answer rather than the most accurate answer, since there may be two winners.) Proof: If, to the contrary, you have an accurate-as-possible answer that does not satisfy the required condition, that means either answer<min(x,y)<=max(x,y) (in which case min(x,y) is a better answer, a contradiction), or min(x,y)<=max(x,y)<answer (in which case max(x,y) is a better answer, a contradiction). So I think that means your question boils down to finding a most accurate possible answer. Assuming IEEE754 arithmetic throughout, I propose the following: if max(abs(x),abs(y)) >= 1.: return x/2. + y/2. else: return (x+y)/2.  My argument that this gives a most accurate answer is a somewhat tedious case analysis. Here goes: • Case max(abs(x),abs(y)) >= 1.: • Subcase neither x nor y is denormalized: In this case the computed answer x/2.+y/2. manipulates the same mantissas and therefore gives the exact same answer as computation of (x+y)/2 would yield if we assumed extended exponents to prevent overflow. This answer may depend on rounding mode but in any case it's guaranteed by IEEE754 to be a best possible answer (from the fact that the computed x+y is guaranteed to be a best approximation to mathematical x+y, and the division by 2 is exact in this case). • Subcase x is denormalized (and so abs(y)>=1): answer = x/2. + y/2. = y/2. since abs(x/2.) is so tiny compared to abs(y/2.) = the exact mathematical value of y/2 = a best possible answer. • Subcase y is denormalized (and so abs(x)>=1): analogous. • Case max(abs(x),abs(y)) < 1.: • Subcase the computed x+y is either non-denormalized or denormalized-and-"even": Although the computed x+y may not be exact, it is guaranteed by IEEE754 to be a best possible approximation to the mathematical x+y. In this case the subsequent division by 2 in the expression (x+y)/2. is exact, so the computed answer (x+y)/2. is a best possible approximation to the mathematical (x+y)/2. • Subcase the computed x+y is denormalized and "odd": In this case exactly one of x,y must also be denormalized-and-"odd", which means the other of x,y is denormalized with the opposite sign, and so the computed x+y is exactly the mathematical x+y, and so the computed (x+y)/2. is guaranteed by IEEE754 to be a best possible approximation to the mathematical (x+y)/2. • I'm realizing when I said "denormalized" I really meant something else-- that is, numbers that are as close to each other as numbers get, i.e. the range of numbers that is roughly twice as big as the range of denormalized numbers, i.e. the first 8 ticks or so in the diagram at en.wikipedia.org/wiki/Denormal_number. The point is, the "odd" ones of these are the only numbers for which dividing them by two isn't exact. I need to re-phrase this part of the answer to make this clear. – Don Hatch Aug 5 '15 at 14:48 • You can simplify your analysis by noting that, absent over-/underflow, it always holds that$\mathit{fl}(\mathit{op}(x,y)) = \mathit{op}(x,y)(1+\delta)$where$|\delta|\leq u$, and that divisions by 2 are exact for non-denormal numbers. Since that means that both$x/2+y/2$and$(x+y)/2\$ are always correctly rounded, absent over-/underflow, all that's left is to show nothing over-/underflows, which is easy. – Kirill Aug 5 '15 at 16:10
• @Kirill I'm a bit lost... where did u come from? Also I don't think it's quite true that "divisions by 2 are exact for non-denormal numbers"... this is the same thing I tripped over, and it seems to be a bit awkward to try to get it right. The precise statement is something more like "x/2 is exact as long as abs(x) is at least twice the greatest subnormal number"... argh, awkward! – Don Hatch Aug 6 '15 at 22:03

For IEEE-754 binary floating-point formats, exemplified by binary64 (double precision) computation, S. Boldo formally proved that the simple algorithm shown below delivers the correctly rounded average.

Sylvie Boldo, "Formal verification of programs computing the floating-point average." In International Conference on Formal Engineering Methods, pp. 17-32. Springer, Cham, 2015. (draft online)

Since division by two is accurate in binary floating-point arithmetic, provided underflow does not occur, it seems intuitively clear that by choosing one of the two formulas $$(x+y)/2$$ and $$x/2 + y/2$$ as appropriate (based on the magnitude of the inputs) once should achieve an accurately-rounded average. Boldo's paper shows that for IEEE-754 binary64 any switchover point $$C \in [2^{-967}, 2^{970}]$$ will suffice. One might pick $$C$$ so as to provide the best performance for a particular use case.

This yields the following exemplary ISO-C99 code:

double average (double x, double y)
{
const double C = 1; /* 0x1p-967 <= C <= 0x1p970 */
return (C <= fabs (x)) ? (x / 2 + y / 2) : ((x + y) / 2);
}


In recent follow-up work, S. Boldo and co-authors showed how to achieve best possible results for the IEEE-754 decimal floating-point formats by making use of fused multiply-add (FMA) operations and a well-known precision-doubling building block (TwoSum):

Sylvie Boldo, Florian Faissole, and Vincent Tourneur, "A Formally-Proved Algorithm to Compute the Correct Average of Decimal Floating-Point Numbers." In 25th IEEE Symposium on Computer Arithmetic (ARITH 25), June 2018, pp. 69-75. (draft online)

Although it may not be super-efficient performance-wise, there is a very simple way to (1) make sure none of the numbers are greater than either x or y (no overflows) and (2) keep the floating point as "accurate" as possible (and (3), as an added bonus, even though subtraction is being used, no values will ever be stored as negative numbers.

float difference = max(x, y) - min(x, y);
return min(x, y) + (difference / 2.0);


In fact, if you really want to go for accuracy, you don't even need to perform the division on the spot; just return the values of min(x, y) and difference which you can use to simplify logically or manipulate later.

• What I'm trying to figure out now is how to make this same answer work with more than two items, while keeping all variables staying lower than the greatest of the numbers, and using only one division operation to preserve accuracy. – IQAndreas Aug 5 '15 at 9:28
• @becko Yup, you would be doing division at least twice. Also, the example you gave would make the answer come out wrong. Imagine the mean of 2,4,9, it is not the same as the mean of 3,9. – IQAndreas Aug 5 '15 at 12:23
• You're right, my recursion was wrong. I'm not sure how to fix it right now, without losing precision. – becko Aug 5 '15 at 12:33
• Can you prove that this gives the most accurate possible result? That is, if x and y are floating-point, your computation produces a floating-point closest to (x+y)/2? – becko Aug 5 '15 at 13:05
• Won't this overflow when x,y are the least and greatest expressible numbers? – Don Hatch Aug 5 '15 at 14:57

Convert to higher precission, add the values there and convert back.

There should be no overflow in the higher precission and if both are in the valid floating point range, the calculated number should be inside too.

And it should be in between them, worst case just half of the larger number if the precission is not sufficient.

• This is the brute force approach. It probably works, but I was looking for an analysis that didn't require intermediate higher precision. Also, can you estimate how much intermediate higher precision is required? In any case, don't delete this answer (+1), I just won't accept it as the answer. – becko Aug 5 '15 at 12:28

Theoretically, x/2 can be computed by subtracting 1 from the mantissa.

However, actually implementing bitwise operations like this is not necessarily straightforward, particularly if you don't know the format of your floating point numbers.

If you can do this, the entire operation is reduced to 3 addition/subtractions, which should be a significant improvement.

I was thinking along the same line as @Roland Heath but can't comment yet, here's my take:

x/2 can be computed by subtracting 1 from the exponent (not the mantissa, subtracting 1 from the mantissa is subtracting 2^(value_of_exponent-length_of_mantissa) from the overall value).

Without restriction of the general case, let's assume x < y. (If x > y, relabel the variables. If x = y, (x+y) / 2 is trivial.)

• Transform (x+y) / 2 into x/2 + y/2, which can be performed by two integer subtractions (by one from the exponent)
• However there's a lower bound on the exponent depending on your representation. If your exponent already is minimal before subtracting 1, this method will require special case handling. A minimal exponent on x will make x/2 smaller than representable (assuming mantissa is represented with an implicit leading 1).
• Instead of subtracting 1 from the exponent of x, shift x's mantissa to the right by one (and add the implicit leading 1, if any).
• Subtract 1 from the exponent of y, if it's not minimal. If it's minimal (y is bigger than x, because of the mantissa), shift the mantissa to the right by one (add implicit leading 1, if any).
• Shift the new mantissa of x to the right according to the exponent of y.
• Perform integer addition on the mantissae, unless the mantissa of x has been shifted out completely. If both exponents were minimal, the leading ones will overflow, which is ok, because that overflow is supposed to become an implicit leading one again.
• and one floating point addition.
• Can't think of any special case here; except for rounding, which also applies to shifting described above.