面向AVX/AVX2/AVX512的双精度正态分布SIMD随机生成器C++实现咨询

核心思路 #

我们优先选择逻辑最简单的实现路线，降低学习门槛：

• 基础均匀分布：并行多通道LCG（线性同余生成器），每个SIMD lane对应独立的生成器实例，参数选择glibc标准64位LCG参数，统计质量满足学习需求
• 正态分布转换：使用Box-Muller变换，公式简单易向量化，不需要复杂的查表操作
• 所有操作直接对应SIMD intrinsic调用，没有隐藏逻辑，方便逐行学习

1. 向量化均匀分布实现（AVX2示例，一次输出4个double） #

LCG的核心更新公式为 state[i] = a * state[i] + c，其中i对应SIMD的每个通道，我们直接用SIMD整数指令并行计算所有通道的状态更新，再将64位整数状态转换为[0,1)区间的双精度浮点数。
核心转换逻辑：取64位状态的高52位作为双精度浮点数的尾数，拼接上1.0对应的指数位，再减去1.0即可得到[0,1)区间的均匀数，比浮点除法效率更高。

2. 向量化正态分布实现 #

基于Box-Muller变换，输入两个独立的[0,1)均匀随机数u1、u2，输出两个独立的标准正态分布数：

z1 = sqrt(-2 * ln(u1)) * cos(2π * u2)
z2 = sqrt(-2 * ln(u1)) * sin(2π * u2)
因为AVX2的__m256d一次处理4个double，我们可以并行计算4组变换，一次输出4个或8个正态分布随机数。

完整可运行AVX2示例代码 #

#include <immintrin.h>
#include <cstdint>
#include <cmath>
#include <array>

// 向量化均匀随机数生成器（AVX2版本，一次生成4个[0,1)双精度数）
class SIMDUniformGen {
private:
    __m256i state;
    static constexpr uint64_t LCG_A = 0x5851f42d4c957f2d;
    static constexpr uint64_t LCG_C = 0x14057b7ef767814f;
public:
    // 初始化需要传入4个不同的64位种子，避免通道序列重复
    explicit SIMDUniformGen(const std::array<uint64_t, 4>& seeds) {
        state = _mm256_set_epi64x(seeds[3], seeds[2], seeds[1], seeds[0]);
    }

    __m256d next() {
        // 更新所有通道的LCG状态
        const __m256i a = _mm256_set1_epi64x(LCG_A);
        const __m256i c = _mm256_set1_epi64x(LCG_C);
        state = _mm256_add_epi64(_mm256_mul_epi64(state, a), c);

        // 64位整数转[0,1)双精度
        const __m256i mantissa = _mm256_srli_epi64(state, 12); // 取高52位作为尾数
        const __m256i exp = _mm256_set1_epi64x(0x3FF0000000000000); // 对应1.0的指数位
        const __m256d one = _mm256_set1_pd(1.0);
        return _mm256_sub_pd(_mm256_castsi256_pd(_mm256_or_si256(mantissa, exp)), one);
    }
};

// 向量化正态分布生成器（标准正态分布，均值0方差1）
class SIMDNormalGen {
private:
    SIMDUniformGen uniform_gen;
public:
    explicit SIMDNormalGen(const std::array<uint64_t, 4>& seeds) : uniform_gen(seeds) {}

    // 一次生成4个标准正态分布双精度数
    __m256d next() {
        // 生成两组均匀数，u1限制最小为1e-10避免ln(0)出错
        __m256d u1 = _mm256_max_pd(uniform_gen.next(), _mm256_set1_pd(1e-10));
        __m256d u2 = uniform_gen.next();

        const __m256d neg_two = _mm256_set1_pd(-2.0);
        const __m256d two_pi = _mm256_set1_pd(2.0 * M_PI);

        __m256d r = _mm256_sqrt_pd(_mm256_mul_pd(neg_two, _mm256_log_pd(u1)));
        __m256d theta = _mm256_mul_pd(two_pi, u2);

        // 取cos计算的结果返回，也可以同时返回sin的结果得到8个随机数
        return _mm256_mul_pd(r, _mm256_cos_pd(theta));
    }
};

// 使用示例
int main() {
    // 初始化4个不同的种子
    std::array<uint64_t, 4> seeds = {12345, 67890, 13579, 24680};
    SIMDNormalGen gen(seeds);

    // 生成1组4个正态分布随机数
    __m256d res = gen.next();

    // 结果存到数组里打印验证
    double arr[4];
    _mm256_storeu_pd(arr, res);
    // 这里可以加打印逻辑输出arr的4个值
    return 0;
}

注：代码中的_mm256_log_pd、_mm256_cos_pd属于SVML intrinsic，MSVC、ICC默认支持，GCC下可以通过添加编译选项-mavx2 -ffast-math -lm启用，或者手动替换为标量函数循环实现，不影响核心逻辑学习。

扩展到其他指令集的方法 #

• AVX（__m128d，一次2个double）：将代码中所有__m256前缀的类型换成__m128，_mm256前缀的intrinsic换成_mm，种子数改为2个即可，核心逻辑完全不变
• AVX512（__m512d，一次8个double）：将类型换成__m512d/__m512i，intrinsic前缀换成_mm512，种子数改为8个即可，不需要修改算法逻辑

学习拓展建议 #

如果掌握了基础实现，可以尝试以下优化方向：

• 替换LCG为向量化xoroshiro128+生成器，提升随机数统计质量
• 使用Box-Muller极坐标版本替代三角函数版本，消除cos/sin调用提升性能
• 添加参数支持自定义均值和标准差

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