标量代码与SSE代码的数组求和结果差异问题排查

我正在实现数组求和归约操作，但SSE代码的计算结果与标量代码不一致。请问我是否错误使用了SSE intrinsics？或是对SIMD编程的某些基础概念存在误解？以下是我的代码及结果差异图，恳请各位提供帮助！

编译命令 #

gcc -Wall -Wextra -pedantic -std=c11 -march=native -o sum sumSSE.c

代码实现 #

//sumSSE.c
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <assert.h>

#include <immintrin.h>

//Scalar sum
float sc_sum(float *vec, uint32_t n)
{
    float sum = 0.0;
    for (uint32_t i = 0; i < n; i++)
        sum += vec[i];
    return sum;
}
//SSE sum
float li_sum(float *vec, uint32_t n)
{
    float sum = 0.0;
    float t[4];
    uint32_t i, n4 = n>>2<<2;
    __m128 v1;
    __m128 v2 = _mm_setzero_ps();
    for (i = 0; i < n4; i+=4) {
        v1 = _mm_loadu_ps(&vec[i]);
        v2 = _mm_add_ps(v2, v1);
    }
    for (; i < n; i++) sum += vec[i];
    _mm_storeu_ps(t, v2);
    sum += t[0] + t[1] + t[2] + t[3];
    return sum;
}

int main(int argc, char **argv)
{
    srand(time(NULL));
    if (argc < 2) return -1;

    uint32_t n = strtol(argv[1], NULL, 10);
    float *vec = malloc(n * sizeof(float));
    if (!vec) return -1;
    //Fill array with random numbers in range 0-1
    for (uint32_t i = 0; i < n; i++)
        vec[i] = (float)rand()/RAND_MAX;

    float r1 = sc_sum(vec, n);
    fprintf(stderr, "R1: %.5f\n", r1);

    float r2 = li_sum(vec, n);
    fprintf(stderr, "R2: %.5f\n", r2);
    
    fprintf(stdout, "r1-r2 == %.8f\n", r1-r2);
    free(vec);
    assert(r1 == r2);
}

结果差异图 #

问题解答 #

核心原因：浮点数加法的非交换性 #

你没有错误使用SSE intrinsics，结果差异是浮点数精度特性导致的正常现象：

• 标量求和按顺序逐个累加所有元素，每次加法的舍入误差会沿固定路径累积。
• SSE实现采用分组并行累加（4个元素一组分别求和），最后合并各组结果。这种累加顺序和标量完全不同，舍入误差的累积方式自然存在差异，最终出现微小数值差。

单精度浮点数（float）仅保留24位有效数字，当累加大量小数值时，这种误差差异会被放大，这就是你看到的差值来源。

修正建议 #

1. 替换严格相等判断：直接用assert(r1 == r2);必然失败，浮点数比较应判断差值是否在可接受阈值内：

   #include <math.h>
   // ...
   assert(fabs(r1 - r2) < 1e-6f);
   

2. 接受误差或调整累加顺序：大多数数值计算场景中，这种微小误差是可接受的。如果必须追求完全一致的结果，可以修改SSE实现的累加逻辑，让其贴近标量的顺序，但这会损失SIMD的并行性能。

内容的提问来源于stack exchange，提问作者jregalad