title: floating point multiplication vs repeated addition
body: <p>Let <code>N</code> be an a compile time unsigned integer.</p>
  
  <p>GCC can optimize</p>
  
  <pre><code>unsigned sum = 0;
  for(unsigned i=0; i&lt;N; i++) sum += a; // a is an unsigned integer   
  </code></pre>
  
  <p>to simply  <code>a*N</code>. This can be understood since modular arithmetic says <code>(a%k + b%k)%k = (a+b)%k</code>.</p>
  
  <p>However GCC will not optimize</p>
  
  <pre><code>float sum = 0;
  for(unsigned i=0; i&lt;N; i++) sum += a;  // a is a float
  </code></pre>
  
  <p>to <code>a*(float)N</code>.</p>
  
  <p>But by using associative math with e.g. <code>-Ofast</code> I discovered that GCC can reduce this in order <code>log2(N)</code> steps.  E.g for <code>N=8</code> it can do the sum in three additions.</p>
  
  <pre><code>sum = a + a
  sum = sum + sum // (a + a) + (a + a)
  sum = sum + sum // ((a + a) + (a + a)) + ((a + a) + (a + a))
  </code></pre>
  
  <p>Though some point after <code>N=16</code> GCC goes back to doing <code>N-1</code> sums. </p>
  
  <p><strong>My question is why does GCC not do <code>a*(float)N</code> with <code>-Ofast</code>?</strong>  </p>
  
  <p>Instead of being <code>O(N)</code> or <code>O(Log(N))</code> it could be simply <code>O(1)</code>. Since <code>N</code> is known at compile time it's possible to determine if <code>N</code> fits in a float. And even if <code>N</code> is too large for a float it could do <code>sum =a*(float)(N &amp; 0x0000ffff) + a*(float)(N &amp; ffff0000)</code>.  In fact, I did a little test to check the accuracy and <code>a*(float)N</code> is more accurate anyway (see the code and results below).</p>
  
  <pre><code>//gcc -O3 foo.c
  //don't use -Ofast or -ffast-math or -fassociative-math
  #include &lt;stdio.h&gt;   
  float sumf(float a, int n)
  {
    float sum = 0;
    for(int i=0; i&lt;n; i++) sum += a;
    return sum;
  }
  
  float sumf_kahan(float a, int n)
  {
    float sum = 0;
    float c = 0;
    for(int i=0; i&lt;n; i++) {
      float y = a - c;
      float t = sum + y;
      c = (t -sum) - y;
      sum = t;
    }
    return sum;
  }  
  
  float mulf(float a, int n)
  {
    return a*n;
  }  
  
  int main(void)
  {
    int n = 1&lt;&lt;24;
    float a = 3.14159;
    float t1 = sumf(a,n);
    float t2 = sumf_kahan(a,n);
    float t3 = mulf(a,n);
    printf("%f %f %f\n",t1,t2,t3);
  }
  </code></pre>
  
  <p>The result is <code>61848396.000000 52707136.000000 52707136.000000</code> which shows that multiplication and the <a href="https://en.wikipedia.org/wiki/Kahan_summation_algorithm" rel="noreferrer">Kahan summation</a> have the same result which I think shows that the multiplication is more accurate than the simple sum.</p>