Rounding a double to avoid round off in subsequent summation

Question

How do I implement this?

// round to the nearest double so that x + ref doesn't cause round off error
double round(double x, double ref) { }

So that

double x = ....;
double y = ....;

double x_new = round(x, y);
return x_new + y; // NO ROUND OFF! 

In other terms (y + x_new) - x_new is strictly equal to y

Answer 1

Let us assume that x and y are both positive.

Let S be the double-precision sum x + y.

There are two cases:

• if x ≤ y, then S - y is exact by Sterbenz's lemma. It follows that the addition (S - y) + y is exact (it produces exactly S, which is a double-precision number). Therefore, you can pick S - y for x_new. Not only y + x_new is exact, but it produces the same result S as y + x.

• if x > y, then depending on the number of bits set in the significand of y, you may have a problem. If the last bit in the significand of y is set, for instance, then no number z in a binade after the binade of y can have the property that z + y is exact.

---

This answer is vaguely related to that answer.

Answer 2

A possible solution in Python that can be straightforwardly translated to C

import math
from decimal import Decimal


def round2(x, ref):
    x_n, x_exp     = math.frexp(x)
    ref_n, ref_exp = math.frexp(ref)
    assert x_exp <= ref_exp
    diff_exp = ref_exp - x_exp

    factor = 2. ** (53 - diff_exp)

    x_new_as_int = int(round(x_n * factor))
    x_new_as_norm_float = float(x_new_as_int) / factor
    return math.ldexp(x_new_as_norm_float, x_exp)


x = 0.001
y = 1.0

assert (y + x) - x != y

x_new = round2(x, y)

assert (y + x_new) - x_new == y

print "x:", Decimal(x)
print "x_new:", Decimal(x_new)

print "relative difference:", (x_new/x - 1.)