c - Rounding a double to avoid round off in subsequent summation -

how implement this?

// round nearest double 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 strictly equal y

let assume x , y both positive.

let s double-precision sum x + y.

there 2 cases:

if x ≤ y, s - y exact sterbenz's lemma. follows addition (s - y) + y exact (it produces s, double-precision number). therefore, can pick s - y x_new. not y + x_new exact, produces same result s y + x.
if x > y, depending on number of bits set in significand of y, may have problem. if last bit in significand of y set, instance, no number z in binade after binade of y can have property z + y exact.

this answer vaguely related that answer.

Test