浮点乘法和重复加法答案

【问题标题】：floating point multiplication vs repeated addition浮点乘法和重复加法
【发布时间】：2016-01-14 03:24:42
【问题描述】：

让N 是一个编译时无符号整数。

GCC 可以优化

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

简单地a*N。这可以理解，因为模运算说(a%k + b%k)%k = (a+b)%k。

但是 GCC 不会优化

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

到a*(float)N。

但是通过将关联数学与例如-Ofast 我发现 GCC 可以按log2(N) 的步骤减少这个。例如 N=8 它可以在三个加法中求和。

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

虽然在N=16 之后的某个时间点，GCC 又回到了N-1 sums 中。

我的问题是为什么 GCC 不使用 a*(float)N 和 -Ofast？

而不是O(N) 或O(Log(N))，它可能只是O(1)。由于N 在编译时是已知的，因此可以确定N 是否适合浮点数。即使N 对于浮点数来说太大，它也可以做到sum =a*(float)(N & 0x0000ffff) + a*(float)(N & ffff0000)。事实上，我做了一个小测试来检查准确性，a*(float)N 无论如何更准确（见下面的代码和结果）。

//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);
}

结果是61848396.000000 52707136.000000 52707136.000000，这表明乘法和Kahan summation具有相同的结果，我认为这表明乘法比简单的和更准确。

【问题讨论】：

您是否考虑过三个相加操作可能比浮点相乘更快？这与它在 N=16 处切换回 fp 乘法是一致的。
优化无效，不做优化；一般来说，它会产生不同的结果。浮点算术不遵守您期望的通常的正常算术属性，例如结合性或分配律。乘法不是重复加法。
@R..：他（隐含地）使用-ffast-math，确实允许进行此类优化。 stackoverflow.com/questions/7420665/…
如果-Ofast 隐含-ffast-math，那么这个问题混合了-ffast-math 和Kahan summation，这不是一个好的秘诀（Kahan summation 是不能用非编译的代码的原型示例） -兼容的优化）。
@Zboson：这非常不明显（当我告诉他这个事实时，甚至 Kahan 都感到震惊，尽管一旦我告诉他，他立即明白为什么这是真的）。如果你好奇的话，最简单的判断它是否正确的方法可能是彻底检查 x 的所有可能尾随位模式的舍入。如果 2^n-1 是可表示的，则证明 (2^n - 1)x + x = 2^n x 可能会得到更深层次的证明，然后注意 x + x + x = 3x。该属性适用于 2、3、4 和 5； 6 是第一个失败的 n。

标签： c gcc optimization floating-point

【解决方案1】：

两者之间有一些根本区别

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

和

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

当总和接近 FLT_EPSILON * 大于值时，重复的总和趋向于无操作。因此，任何较大的 N 值都不会导致重复加法的总和发生变化。对于乘法选择，结果（值 * N）需要 FLT_EPSILON * 小于 sum 才能使运算无操作。

因此编译器无法进行优化，因为它无法判断您是想要确切的行为（乘法更好），还是执行的行为，其中总和的比例会影响加法的结果。

【讨论】：