Is it always true that x * y = ((x * y) / y) * y under IEEE 754 semantics?

Question

Given two nonzero, finite, double-precision floating point numbers x and y, is it always true that the equality

x * y == ((x * y) / y) * y

holds under default IEEE 754 semantics?

I've searched programmatically over billions of possibilities (including in the subnormal range) and am unable to find a counterexample, yet am also unsure how to go about proving the assertion to be true.

(What I do know is that the simpler assertions x == (x / y) * y and x == (x * y) / y are both false, as I can easily find counterexamples for them.)

Answer 1

In double precision (binary64), i.e. a precision of p = 53 bits, one has the following counter-example:

#include <stdio.h>
int main (int argc, char *argv[])
{
  double x = 0x1.0000004000001p+0;
  double y = 0x1.ffffff8000001p+0;
  double r1 = x * y;
  double r2 = ((x * y) / y) * y;
  printf ("r1 = %a\n", r1);
  printf ("r2 = %a\n", r2);
  return 0;
}

which outputs

r1 = 0x1p+1
r2 = 0x1.fffffffffffffp+0

I found this counter-example by doing exhaustive tests (the property does not depend on the exponent of x and y, which can be set to a constant) in small precisions with GNU MPFR, and I found a pattern: in precision p = 2k+3, the following pair written in base 2 is a counter-example:

x = 1.0{k}10{k}1
y = 1.1{k}00{k}1

where {k} means that the digit is repeated k times.

Note: This problem is a bit similar to what yielded a bug in GNU Emacs + Cairo, where the expression ((1/s)*b)*s) is involved. See my article.