title: Implementing single-precision division as double-precision multiplication
body: <h2>Question</h2>
  
  <p>For a C99 compiler implementing exact IEEE 754 arithmetic, do values of <code>f</code>, <code>divisor</code> of type <code>float</code> exist such that <code>f / divisor != (float)(f * (1.0 / divisor))</code>?</p>
  
  <p>EDIT: By “implementing exact IEEE 754 arithmetic” I mean a compiler that rightfully defines FLT_EVAL_METHOD as 0.</p>
  
  <h2>Context</h2>
  
  <p>A C compiler that provides IEEE 754-compliant floating-point can only replace a single-precision division by a constant by a single-precision multiplication by the inverse if said inverse is itself representable exactly as a <code>float</code>.</p>
  
  <p>In practice, this only happens for powers of two. So a programmer, Alex, may be confident that <code>f / 2.0f</code> will be compiled as if it had been <code>f * 0.5f</code>, but if it is acceptable for Alex to multiply by <code>0.10f</code> instead of dividing by 10, Alex should express it by writing the multiplication in the program, or by using a compiler option such as GCC's <code>-ffast-math</code>.</p>
  
  <p>This question is about transforming a single-precision division into a double-precision multiplication. Does it always produce the correctly rounded result? Is there a chance that it  could be cheaper, and thus be an optimization that compilers might make (even without <code>-ffast-math</code>)?</p>
  
  <p>I have compared <code>(float)(f * 0.10)</code> and <code>f / 10.0f</code> for all single-precision values of <code>f</code> between 1 and 2, without finding any counter-example. This should cover all divisions of normal <code>float</code>s producing a normal result.</p>
  
  <p>Then I generalized the test to all divisors with the program below:</p>
  
  <pre><code>#include &lt;float.h&gt;
  #include &lt;math.h&gt;
  #include &lt;stdio.h&gt;
  
  int main(void){
    for (float divisor = 1.0; divisor != 2.0; divisor = nextafterf(divisor, 2.0))
      {
        double factor = 1.0 / divisor; // double-precision inverse
        for (float f = 1.0; f != 2.0; f = nextafterf(f, 2.0))
          {
            float cr = f / divisor;
            float opt = f * factor; // double-precision multiplication
            if (cr != opt)
              printf("For divisor=%a, f=%a, f/divisor=%a but (float)(f*factor)=%a\n",
                     divisor, f, cr, opt);
          }
      }
  }
  </code></pre>
  
  <p>The search space is just large enough to make this interesting (2<sup>46</sup>). The program is currently running. Can someone tell me whether it will print something, perhaps with an explanation why or why not, before it has finished?</p>