Floating point multiplication and re-addition

Question

Let N be an unsigned compilation integer.

GCC can optimize

 unsigned sum = 0; for(unsigned i=0; i<N; i++) sum += a; // a is an unsigned integer 

just a*N This can be understood since modular arithmetic says (a%k + b%k)%k = (a+b)%k .

However gcc will not optimize

 float sum = 0; for(unsigned i=0; i<N; i++) sum += a; // a is a float 

up to a*(float)N

But using associative math, for example, -Ofast I found that GCC can reduce this in order of log2(N) . For example, for N=8 he can make the sum in three additions.

 sum = a + a sum = sum + sum // (a + a) + (a + a) sum = sum + sum // ((a + a) + (a + a)) + ((a + a) + (a + a)) 

Although some point after N=16 GCC returns to the execution of the sums of N-1 .

My question is: why does GCC not make a*(float)N with -Ofast ?

Instead of O(N) or O(Log(N)) it could just be O(1) . Since N known at compile time, you can determine if N is suitable in the float. And even if N too large for the float, it can do sum =a*(float)(N & 0x0000ffff) + a*(float)(N & ffff0000) . In fact, I did a little test to check the accuracy, and a*(float)N more accurate anyway (see Code and results below).

 //gcc -O3 foo.c //don't use -Ofast or -ffast-math or -fassociative-math #include <stdio.h> float sumf(float a, int n) { float sum = 0; for(int i=0; i<n; i++) sum += a; return sum; } float sumf_kahan(float a, int n) { float sum = 0; float c = 0; for(int i=0; i<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<<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); } 

The result is 61848396.000000 52707136.000000 52707136.000000 , which shows that the Kahan multiplication and summation have the same result, which, I think, shows that the multiplication is more accurate than a simple sum.

Answer 1

There is some fundamental difference between

  float funct( int N, float sum ) { float value = 10.0; for( i = 0; i < N ;i ++ ) { sum += value; } return sum; } 

and

 float funct( int N, float sum ) { float value = 10.0; sum += value * N; return sum; } 

When the sum approaches FLT_EPSILON * more than the value, the repeated sum tends to no-op. Thus, any large value of N will lead to a change in the amount for re-addition. To select a multiplication, the result (value * N) must be FLT_EPSILON * less than the sum for an operation without an operation.

Thus, the compiler cannot do the optimization because it cannot determine whether you want the exact behavior (where multiplication is better) or the implemented behavior, where the scale of the sum affects the result of the addition.