Умножение с плавающей запятой и повторное добавление

Пусть 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.

Мой вопрос: почему 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, который показывает, что умножение и суммирование Кахана имеют тот же результат, который, я думаю, показывает, что умножение является более точным чем простая сумма.

Ответ 1