C++ 与 ‘SIMD’ (单指令多数据流)：利用 AVX-512 指令集加速高性能物理模拟与加解密

各位来宾，各位对高性能计算和底层优化技术抱有浓厚兴趣的朋友们，大家好。

今天，我们将深入探讨一个在现代计算领域至关重要的技术——SIMD（单指令多数据流），特别是围绕 Intel 的 AVX-512 指令集，如何在 C++ 环境下，为高性能物理模拟和加解密算法带来革命性的加速。随着处理器核心数的不断增加，以及对数据处理吞吐量要求的日益严苛，传统的标量计算已经难以满足需求。SIMD 技术应运而生，它允许单个指令同时操作多个数据元素，从而极大地提升了并行处理能力。而 AVX-512 作为目前 x86 架构上最先进、最强大的 SIMD 指令集之一，其 512 位宽的向量寄存器和丰富的功能集，为计算密集型任务的优化打开了全新的大门。

本讲座将从 SIMD 的基本概念出发，逐步深入到 AVX-512 的特性，并通过具体的 C++ 代码示例，演示如何利用编译器内在函数（Intrinsics）和适当的数据结构设计，将这些强大的指令集应用于实际问题中。我们将聚焦于物理模拟中的力计算和流体网格操作，以及加解密领域的 AES 和 SHA-256 算法，剖析其优化策略和潜在挑战。

开篇：SIMD 与高性能计算的时代脉动

1.1 SIMD 的核心思想与演进

SIMD，即 Single Instruction, Multiple Data，其核心思想是让一个指令能够同时处理多个数据元素。想象一下，你有一队工人需要搬运砖块，传统的标量（Scalar）处理就像是每个工人一次只能搬一块砖。而 SIMD 则像是一个工人能够一次性搬运一堆（比如 8 块、16 块）砖块。在处理器层面，这意味着使用一个指令周期完成多个独立但相同的数据操作，例如同时对四个浮点数进行加法，或者对八个整数进行乘法。

SIMD 技术在 x86 架构上经历了漫长的演进：

• MMX (MultiMedia eXtensions)：Intel 于 1997 年推出，提供了 64 位寄存器，主要用于整数运算，面向多媒体应用。
• SSE (Streaming SIMD Extensions)：1999 年随 Pentium III 引入，提供了 128 位寄存器，支持浮点运算，是现代 SIMD 的基石。后续的 SSE2, SSE3, SSSE3, SSE4.1, SSE4.2 不断扩展功能。
• AVX (Advanced Vector Extensions)：2011 年随 Sandy Bridge 架构发布，将 SIMD 寄存器宽度扩展到 256 位（YMM 寄存器），引入了三操作数指令和非破坏性源操作数，极大地提高了编程灵活性和表达能力。
• AVX2 (Advanced Vector Extensions 2)：2013 年随 Haswell 架构发布，将整数运算也扩展到 256 位，并增加了 Gather 指令，允许从非连续内存地址加载数据。
• AVX-512 (Advanced Vector Extensions 512)：2015 年随 Knights Landing（KNL）架构首次亮相，并在 2017 年的 Skylake-X 处理器上进一步普及。这是我们今天关注的重点，它将 SIMD 寄存器宽度再次翻倍至 512 位（ZMM 寄存器），引入了多达 32 个 ZMM 寄存器，并带来了丰富的全新指令集，如掩码操作、融合乘加（FMA）指令、完善的 Gather/Scatter 指令、位操作指令等。

1.2 AVX-512：性能突破的利器

AVX-512 带来了前所未有的并行度。对于单精度浮点数（float），512 位寄存器可以同时处理 16 个元素；对于双精度浮点数（double），可以处理 8 个元素；对于 32 位整数（int），可以处理 16 个元素；对于 64 位整数（long long），可以处理 8 个元素。这意味着相比 AVX2，理论吞吐量翻倍。

其主要特性包括：

• ZMM 寄存器：32 个 512 位寄存器（zmm0 – zmm31），相比 AVX/AVX2 的 16 个 YMM 寄存器，数量和宽度都翻倍。
• 掩码寄存器 (K 寄存器)：8 个 16 位掩码寄存器（k0 – k7），用于对向量操作进行元素级的条件控制。这极大地简化了条件分支的向量化，提高了代码效率，避免了昂贵的 if/else 语句在循环内部的开销。
• 新的指令前缀 (EVEX)：支持更灵活的指令编码，允许三操作数、掩码、广播、内存操作数尺寸控制等功能。
• 丰富的指令集：
  • Foundation (AVX512F)：基础指令集，涵盖所有数据类型和基本算术逻辑操作。
  • Conflict Detection Instructions (AVX512CD)：用于检测向量中的冲突，在某些算法（如稀疏矩阵）中很有用。
  • Exponential and Reciprocal (AVX512ER)：用于快速计算超越函数，主要用于 KNL。
  • Prefetch (AVX512PF)：预取指令，主要用于 KNL。
  • Vector Bit Manipulation (AVX512VBMI)：更强大的位操作，如字节置换。
  • Vector Neural Network Instructions (AVX512_VNNI)：专为深度学习推理优化，支持 8 位/16 位整数乘累加。
  • Galois Field New Instructions (GFNI)：用于伽罗瓦域运算，对密码学有益。
  • 还有其他如 VPCLMULQDQ (Carry-less Multiplication), VPOPCNTDQ (Popcount) 等，极大地扩展了 AVX-512 的应用范围。

AVX-512 的强大功能带来了显著的性能提升潜力，尤其是在数据并行度高、计算密集型的应用中，如科学计算、金融建模、图像/视频处理、机器学习和密码学。

C++ 编程范式下的 SIMD 接口

在 C++ 中利用 SIMD 指令集，主要有三种方式：编译器内在函数（Intrinsics）、向量化库和自动向量化。

2.1 编译器内在函数 (Compiler Intrinsics)

内在函数是 C++ 编译器提供的一种特殊函数，它们直接映射到 CPU 的 SIMD 指令。使用内在函数的好处是能够精确控制生成的汇编代码，最大限度地发挥硬件性能。缺点是代码可读性较差，且不具备跨平台可移植性（不同的 CPU 架构有不同的内在函数集）。

对于 AVX-512，相关的头文件是 <immintrin.h>。它定义了诸如 __m512 (单精度浮点数向量)、__m512d (双精度浮点数向量)、__m512i (整数向量) 等数据类型，以及大量的 _mm512_* 前缀的函数。

示例：AVX-512 单精度浮点数向量加法

#include <iostream>
#include <immintrin.h> // 包含 AVX-512 intrinsics
#include <vector>
#include <numeric>
#include <chrono>

// 确保内存对齐函数
void* aligned_malloc(size_t size, size_t alignment) {
    void* ptr;
    #ifdef _MSC_VER
        ptr = _aligned_malloc(size, alignment);
    #else
        if (posix_memalign(&ptr, alignment, size) != 0) {
            ptr = nullptr;
        }
    #endif
    return ptr;
}

void aligned_free(void* ptr) {
    #ifdef _MSC_VER
        _aligned_free(ptr);
    #else
        free(ptr);
    #endif
}

int main() {
    const int N = 1024 * 1024; // 大量数据
    const int ALIGNMENT = 64; // AVX-512 要求 64 字节对齐

    // 分配对齐内存
    float* a = static_cast<float*>(aligned_malloc(N * sizeof(float), ALIGNMENT));
    float* b = static_cast<float*>(aligned_malloc(N * sizeof(float), ALIGNMENT));
    float* c_scalar = static_cast<float*>(aligned_malloc(N * sizeof(float), ALIGNMENT));
    float* c_simd = static_cast<float*>(aligned_malloc(N * sizeof(float), ALIGNMENT));

    if (!a || !b || !c_scalar || !c_simd) {
        std::cerr << "Failed to allocate aligned memory." << std::endl;
        return 1;
    }

    // 初始化数据
    for (int i = 0; i < N; ++i) {
        a[i] = static_cast<float>(i);
        b[i] = static_cast<float>(i * 2);
    }

    // 标量加法
    auto start_scalar = std::chrono::high_resolution_clock::now();
    for (int i = 0; i < N; ++i) {
        c_scalar[i] = a[i] + b[i];
    }
    auto end_scalar = std::chrono::high_resolution_clock::now();
    std::chrono::duration<double> diff_scalar = end_scalar - start_scalar;
    std::cout << "Scalar addition time: " << diff_scalar.count() << " s" << std::endl;

    // AVX-512 SIMD 加法
    auto start_simd = std::chrono::high_resolution_clock::now();
    // AVX-512 可以同时处理 16 个 float
    const int VEC_SIZE = 16;
    for (int i = 0; i < N; i += VEC_SIZE) {
        // 从内存加载 16 个 float 到 __m512 向量
        __m512 va = _mm512_load_ps(&a[i]);
        __m512 vb = _mm512_load_ps(&b[i]);

        // 执行向量加法
        __m512 vc = _mm512_add_ps(va, vb);

        // 将结果存储回内存
        _mm512_store_ps(&c_simd[i], vc);
    }
    auto end_simd = std::chrono::high_resolution_clock::now();
    std::chrono::duration<double> diff_simd = end_simd - start_simd;
    std::cout << "AVX-512 SIMD addition time: " << diff_simd.count() << " s" << std::endl;

    // 验证结果（部分）
    bool correct = true;
    for (int i = 0; i < N; ++i) {
        if (c_scalar[i] != c_simd[i]) {
            correct = false;
            break;
        }
    }
    std::cout << "Results are " << (correct ? "correct" : "incorrect") << std::endl;

    // 释放内存
    aligned_free(a);
    aligned_free(b);
    aligned_free(c_scalar);
    aligned_free(c_simd);

    return 0;
}

编译命令示例:
g++ -O3 -march=skylake-avx512 -mavx512f -Wall -std=c++17 simd_add.cpp -o simd_add
或者
cl /O2 /arch:AVX512 /EHsc simd_add.cpp /Fe:simd_add.exe

上述代码展示了如何使用 _mm512_load_ps 加载 512 位单精度浮点数，使用 _mm512_add_ps 执行向量加法，然后使用 _mm512_store_ps 存储结果。注意，这里使用了 _mm512_load_ps 和 _mm512_store_ps，它们是要求内存地址 64 字节对齐的。如果不确定对齐，可以使用 _mm512_loadu_ps 和 _mm512_storeu_ps，但它们可能性能稍差。

2.2 向量化库 (Vectorization Libraries)

为了提高可移植性和编程便利性，许多向量化库应运而生。它们通常提供 C++ 类或模板，封装了底层 SIMD 指令，并提供类似数学向量操作的接口。开发者无需直接接触内在函数，就可以编写更易读、更易维护的向量化代码。

• VCL (Vector Class Library by Agner Fog)：一个非常流行且高效的 C++ 模板库，支持多种 SIMD 指令集（SSE, AVX, AVX2, AVX-512），并能在不同平台上编译。
• Eigen：一个 C++ 模板库，用于线性代数，它内部也大量使用 SIMD 指令进行优化。
• xsimd：一个现代 C++ SIMD 抽象库，旨在提供统一的接口来利用各种 SIMD 指令集。

这些库通常会根据编译器的目标架构自动选择最合适的 SIMD 指令集。使用它们可以显著降低 SIMD 编程的复杂性。

2.3 自动向量化 (Auto-Vectorization)

现代 C++ 编译器（如 GCC, Clang, MSVC）都具备强大的自动向量化能力。它们会分析循环结构和数据访问模式，尝试将标量代码自动转换为 SIMD 指令。这对于那些符合特定模式（如数据独立、循环计数已知、内存访问连续）的循环非常有效。

优点：无需手动编写 SIMD 代码，易于维护。
缺点：编译器并非总是能成功向量化，或者生成的代码不如手动优化的高效。对于复杂的算法或带有分支的代码，自动向量化可能失败。

为了帮助编译器更好地进行自动向量化，开发者可以使用一些技巧：

• #pragma omp simd (OpenMP 5.0 及更高版本)：明确告诉编译器尝试向量化某个循环。
• __restrict__ 关键字：告诉编译器指针之间没有别名（即不指向同一块内存），这有助于编译器进行更积极的优化。
• 循环展开 (Loop Unrolling)：有时手动或通过编译器选项进行循环展开，可以为 SIMD 提供更多并行机会。
• 数据对齐：使用 alignas 关键字或 _aligned_malloc 确保数据对齐，这对 SIMD 性能至关重要。

在实际项目中，通常是自动向量化作为第一步，对于性能瓶颈，再考虑使用内在函数或向量化库进行精细优化。

案例分析一：高性能物理模拟的 AVX-512 加速

物理模拟是高性能计算的典型应用场景，它通常涉及大量重复的数学运算，具有高度的数据并行性，非常适合 SIMD 优化。

3.1 粒子系统模拟：N-body 问题的力计算

N-body 问题（如万有引力模拟、分子动力学）需要计算系统中每对粒子之间的相互作用力。对于 $N$ 个粒子，这通常涉及 $O(N^2)$ 次力计算。

假设我们有一个简化的 N-body 模拟，每个粒子有位置 $(x, y, z)$ 和质量 $m$。我们需要计算每个粒子受到的合力。

数据结构的选择：SoA vs. AoSoA

在 SIMD 编程中，数据结构的选择至关重要。

• AoA (Array of Arrays) / SoA (Structure of Arrays)：将不同分量的数据分别存储在独立的数组中。
  • float x[N], y[N], z[N], mass[N];
  • 优点：内存访问连续，缓存命中率高，非常适合 SIMD 加载和存储。一次加载 x 数组的 16 个元素，就得到了 16 个粒子的 x 坐标。
  • 缺点：访问单个粒子的所有属性需要从不同内存区域读取。
• AoS (Array of Structures)：将每个粒子的所有属性打包成一个结构体，然后存储结构体数组。
  • struct Particle { float x, y, z, mass; }; Particle particles[N];
  • 优点：访问单个粒子的所有属性方便。
  • 缺点：内存访问不连续，对于 SIMD 而言，一次加载一个 Particle 结构体，可能只用到了其中的一部分数据，其余部分是浪费的。
• AoSoA (Array of Structures of Arrays)：结合 AoS 和 SoA 的优点。将数据分为多个小块，每块内部使用 SoA 布局。
  • struct ParticleBlock { float x[VEC_SIZE], y[VEC_SIZE], z[VEC_SIZE], mass[VEC_SIZE]; }; ParticleBlock blocks[N/VEC_SIZE];
  • 这种结构对于 SIMD 更加灵活，可以在处理每个块时进行高效的向量操作。

对于 N-body 模拟，SoA 或 AoSoA 是更好的选择。这里我们使用 SoA 简化演示。

力计算公式：
对于粒子 $i$ 和粒子 $j$，它们之间的力是：
$F_{ij} = G frac{m_i mj}{r{ij}^2} cdot frac{mathbf{r}{ij}}{r{ij}}$
其中 $mathbf{r}_{ij} = mathbf{p}_j – mathbf{p}_i = (x_j-x_i, y_j-y_i, z_j-zi)$，$r{ij} = |mathbf{r}_{ij}|$ 是距离。

C++ 示例：N-body 力计算（简化版）

#include <iostream>
#include <vector>
#include <cmath>
#include <numeric>
#include <chrono>
#include <immintrin.h> // AVX-512 intrinsics

// G 常量
const float G = 6.674e-11f;

// 确保内存对齐函数 (同上)
void* aligned_malloc(size_t size, size_t alignment) { /* ... */ return nullptr; }
void aligned_free(void* ptr) { /* ... */ }

// SoA 数据结构
struct SoAParticles {
    float* x;
    float* y;
    float* z;
    float* mass;
    float* fx;
    float* fy;
    float* fz;
    int N;

    SoAParticles(int n, int alignment = 64) : N(n) {
        x = static_cast<float*>(aligned_malloc(N * sizeof(float), alignment));
        y = static_cast<float*>(aligned_malloc(N * sizeof(float), alignment));
        z = static_cast<float*>(aligned_malloc(N * sizeof(float), alignment));
        mass = static_cast<float*>(aligned_malloc(N * sizeof(float), alignment));
        fx = static_cast<float*>(aligned_malloc(N * sizeof(float), alignment));
        fy = static_cast<float*>(aligned_malloc(N * sizeof(float), alignment));
        fz = static_cast<float*>(aligned_malloc(N * sizeof(float), alignment));

        // 检查分配是否成功，这里省略了详细检查
        if (!x || !y || !z || !mass || !fx || !fy || !fz) {
            std::cerr << "Failed to allocate particle memory." << std::endl;
            // 实际应用中需要更严谨的错误处理
        }
    }

    ~SoAParticles() {
        aligned_free(x); aligned_free(y); aligned_free(z); aligned_free(mass);
        aligned_free(fx); aligned_free(fy); aligned_free(fz);
    }

    void init_random() {
        for (int i = 0; i < N; ++i) {
            x[i] = static_cast<float>(rand()) / RAND_MAX * 100.0f;
            y[i] = static_cast<float>(rand()) / RAND_MAX * 100.0f;
            z[i] = static_cast<float>(rand()) / RAND_MAX * 100.0f;
            mass[i] = static_cast<float>(rand()) / RAND_MAX * 10.0f + 1.0f;
            fx[i] = fy[i] = fz[i] = 0.0f;
        }
    }
};

// 标量计算函数
void compute_forces_scalar(SoAParticles& p) {
    for (int i = 0; i < p.N; ++i) {
        p.fx[i] = p.fy[i] = p.fz[i] = 0.0f; // 重置力
        for (int j = 0; j < p.N; ++j) {
            if (i == j) continue;

            float dx = p.x[j] - p.x[i];
            float dy = p.y[j] - p.y[i];
            float dz = p.z[j] - p.z[i];

            float dist_sq = dx * dx + dy * dy + dz * dz;
            // 避免除零和超近距离引力过大，引入一个小的平滑因子 epsilon
            const float epsilon = 1e-6f;
            dist_sq += epsilon;
            float dist = std::sqrt(dist_sq);

            float inv_dist_cubed = 1.0f / (dist_sq * dist);
            float magnitude = G * p.mass[i] * p.mass[j] * inv_dist_cubed;

            p.fx[i] += dx * magnitude;
            p.fy[i] += dy * magnitude;
            p.fz[i] += dz * magnitude;
        }
    }
}

// AVX-512 SIMD 计算函数
void compute_forces_simd(SoAParticles& p) {
    const int VEC_SIZE = 16; // 512 bits / 32 bits per float = 16 floats
    const __m512 v_G = _mm512_set1_ps(G);
    const __m512 v_epsilon = _mm512_set1_ps(1e-6f);

    for (int i = 0; i < p.N; ++i) {
        // 重置力向量
        __m512 v_fx_i = _mm512_setzero_ps();
        __m512 v_fy_i = _mm512_setzero_ps();
        __m512 v_fz_i = _mm512_setzero_ps();

        __m512 v_xi = _mm512_set1_ps(p.x[i]);
        __m512 v_yi = _mm512_set1_ps(p.y[i]);
        __m512 v_zi = _mm512_set1_ps(p.z[i]);
        __m512 v_mass_i = _mm512_set1_ps(p.mass[i]);

        for (int j = 0; j < p.N; j += VEC_SIZE) {
            // 处理剩余元素，这里假设 N 是 VEC_SIZE 的倍数，实际应用需要处理余数
            // 或使用掩码指令 _mm512_mask_load_ps 等
            if (j + VEC_SIZE > p.N) { // 处理尾部，简化起见，这里跳过
                // 实际应用中，这里需要一个标量循环处理剩余的元素，或者使用掩码加载/存储
                // _mm512_mask_load_ps, _mm512_mask_store_ps
                continue;
            }

            // 加载 16 个粒子的位置和质量
            __m512 v_xj = _mm512_load_ps(&p.x[j]);
            __m512 v_yj = _mm512_load_ps(&p.y[j]);
            __m512 v_zj = _mm512_load_ps(&p.z[j]);
            __m512 v_mass_j = _mm512_load_ps(&p.mass[j]);

            // 计算相对位移
            __m512 v_dx = _mm512_sub_ps(v_xj, v_xi);
            __m512 v_dy = _mm512_sub_ps(v_yj, v_yi);
            __m512 v_dz = _mm512_sub_ps(v_zj, v_zi);

            // 计算距离平方: dist_sq = dx*dx + dy*dy + dz*dz
            __m512 v_dist_sq = _mm512_fmadd_ps(v_dx, v_dx,
                               _mm512_fmadd_ps(v_dy, v_dy,
                               _mm512_mul_ps(v_dz, v_dz))); // FMA 指令在此处非常高效

            // 加入平滑因子
            v_dist_sq = _mm512_add_ps(v_dist_sq, v_epsilon);

            // 计算距离倒数立方: 1.0f / (dist_sq * dist)
            // rsqrtps 计算近似倒数平方根，精度可能不够，这里使用精确的
            __m512 v_dist = _mm512_sqrt_ps(v_dist_sq);
            __m512 v_inv_dist = _mm512_rcp_ps(v_dist); // 近似倒数
            // 迭代改进倒数精度，或者直接使用 _mm512_div_ps(v_one, v_dist);
            // 这里为了简化，直接用近似倒数，或假设编译器优化_mm512_div_ps

            __m512 v_inv_dist_cubed = _mm512_mul_ps(_mm512_mul_ps(v_inv_dist, v_inv_dist), v_inv_dist);

            // 计算力的大小: G * mass_i * mass_j * inv_dist_cubed
            __m512 v_magnitude = _mm512_mul_ps(v_G, v_mass_i);
            v_magnitude = _mm512_mul_ps(v_magnitude, v_mass_j);
            v_magnitude = _mm512_mul_ps(v_magnitude, v_inv_dist_cubed);

            // 计算分量力并累加
            v_fx_i = _mm512_fmadd_ps(v_dx, v_magnitude, v_fx_i);
            v_fy_i = _mm512_fmadd_ps(v_dy, v_magnitude, v_fy_i);
            v_fz_i = _mm512_fmadd_ps(v_dz, v_magnitude, v_fz_i);
        }

        // 将累加的力向量横向求和并存储
        // _mm512_reduce_add_ps 是一个方便的内在函数，可以进行横向求和
        p.fx[i] = _mm512_reduce_add_ps(v_fx_i);
        p.fy[i] = _mm512_reduce_add_ps(v_fy_i);
        p.fz[i] = _mm512_reduce_add_ps(v_fz_i);
    }
}

int main() {
    const int N_PARTICLES = 1024; // 粒子数量
    SoAParticles particles(N_PARTICLES);
    particles.init_random();

    // 标量计算
    auto start_scalar = std::chrono::high_resolution_clock::now();
    compute_forces_scalar(particles);
    auto end_scalar = std::chrono::high_resolution_clock::now();
    std::chrono::duration<double> diff_scalar = end_scalar - start_scalar;
    std::cout << "Scalar N-body force calculation time: " << diff_scalar.count() << " s" << std::endl;

    // 清空力，重新计算
    for(int i=0; i<N_PARTICLES; ++i) particles.fx[i] = particles.fy[i] = particles.fz[i] = 0.0f;

    // AVX-512 SIMD 计算
    auto start_simd = std::chrono::high_resolution_clock::now();
    compute_forces_simd(particles);
    auto end_simd = std::chrono::high_resolution_clock::now();
    std::chrono::duration<double> diff_simd = end_simd - start_simd;
    std::cout << "AVX-512 SIMD N-body force calculation time: " << diff_simd.count() << " s" << std::endl;

    // 验证结果 (简化，只检查第一个粒子)
    // 实际验证需要比较浮点数，考虑精度误差
    // std::cout << "Particle 0 force X (scalar): " << particles.fx[0] << std::endl;
    // std::cout << "Particle 0 force X (SIMD): " << particles.fx[0] << std::endl;

    return 0;
}

关键优化点：

• SoA 数据布局：确保内存访问连续，利于 _mm512_load_ps。
• 广播操作：_mm512_set1_ps(value) 将单个标量值广播到向量的所有 16 个元素中，以便与向量进行操作。
• 融合乘加 (FMA)：_mm512_fmadd_ps(a, b, c) 计算 a*b + c，在一个指令周期内完成乘法和加法，减少延迟，提高吞吐量。在距离平方和力分量累加中非常有用。
• 近似函数：_mm512_rcp_ps (倒数) 和 _mm512_rsqrt_ps (倒数平方根) 可以提供非常快的近似值。如果需要更高精度，可能需要通过牛顿迭代法对结果进行 refinement，或者直接使用 _mm512_div_ps 和 _mm512_sqrt_ps。
• 横向求和：_mm512_reduce_add_ps 用于将 16 个浮点数相加得到一个标量结果，这在累加每个粒子受到的总力时是必要的。
• 处理剩余元素：在 compute_forces_simd 函数中，为了简化，我们假设 N 是 VEC_SIZE 的倍数。在实际应用中，循环末尾剩余的不足 VEC_SIZE 的元素需要单独用标量代码处理，或者利用 AVX-512 的掩码指令 _mm512_mask_load_ps 和 _mm512_mask_store_ps 来只处理有效元素。

3.2 流体动力学：基于格点的数值计算

流体动力学模拟，特别是基于格点的数值方法（如有限差分法、格子玻尔兹曼方法），通常涉及对网格中的每个点应用一个“模板”操作。例如，计算一个点的新值可能取决于它周围邻居点的值。这类操作天然具有局部性和并行性。

示例：2D 雅可比迭代（Jacobi Iteration）

雅可比迭代是求解线性方程组的一种方法，在物理模拟中常用于扩散过程或泊松方程的求解。在 2D 网格上，每个点的新值是其四个邻居（上、下、左、右）的平均值。

$U{i,j}^{k+1} = frac{1}{4} (U{i-1,j}^k + U{i+1,j}^k + U{i,j-1}^k + U_{i,j+1}^k)$

我们需要两个网格：一个存储旧值 (U_old)，一个存储新值 (U_new)。

#include <iostream>
#include <vector>
#include <chrono>
#include <immintrin.h> // AVX-512 intrinsics

// 确保内存对齐函数 (同上)
void* aligned_malloc(size_t size, size_t alignment) { /* ... */ return nullptr; }
void aligned_free(void* ptr) { /* ... */ }

// 辅助函数：初始化网格
void init_grid(float* grid, int M, int N) {
    for (int i = 0; i < M; ++i) {
        for (int j = 0; j < N; ++j) {
            grid[i * N + j] = static_cast<float>(i + j); // 简单初始化
        }
    }
}

// 标量雅可比迭代
void jacobi_scalar(float* U_new, const float* U_old, int M, int N, int iterations) {
    for (int iter = 0; iter < iterations; ++iter) {
        for (int i = 1; i < M - 1; ++i) { // 避免边界
            for (int j = 1; j < N - 1; ++j) {
                U_new[i * N + j] = 0.25f * (U_old[(i - 1) * N + j] + // 上
                                            U_old[(i + 1) * N + j] + // 下
                                            U_old[i * N + (j - 1)] + // 左
                                            U_old[i * N + (j + 1)]); // 右
            }
        }
        // 交换 U_new 和 U_old，以便下一轮迭代
        std::swap(U_new, const_cast<float*&>(U_old));
    }
}

// AVX-512 SIMD 雅可比迭代
void jacobi_simd(float* U_new, const float* U_old, int M, int N, int iterations) {
    const int VEC_SIZE = 16;
    const __m512 v_quarter = _mm512_set1_ps(0.25f);

    for (int iter = 0; iter < iterations; ++iter) {
        for (int i = 1; i < M - 1; ++i) {
            // 注意这里 j 从 1 开始，每次处理 16 个元素，所以 j_start 应是 VEC_SIZE 的倍数
            // 如果 N-1 不是 VEC_SIZE 的倍数，需要处理剩余元素。
            // 这里为了简化，假设 (N-2) 是 VEC_SIZE 的倍数，且 N 足够大。
            for (int j = 1; j < N - 1; j += VEC_SIZE) {
                // 加载周围邻居的值
                __m512 v_up = _mm512_load_ps(&U_old[(i - 1) * N + j]);
                __m512 v_down = _mm512_load_ps(&U_old[(i + 1) * N + j]);
                __m512 v_left = _mm512_load_ps(&U_old[i * N + (j - 1)]);
                __m512 v_right = _mm512_load_ps(&U_old[i * N + (j + 1)]);

                // 累加
                __m512 v_sum = _mm512_add_ps(v_up, v_down);
                v_sum = _mm512_add_ps(v_sum, v_left);
                v_sum = _mm512_add_ps(v_sum, v_right);

                // 乘以 0.25
                __m512 v_result = _mm512_mul_ps(v_sum, v_quarter);

                // 存储结果
                _mm512_store_ps(&U_new[i * N + j], v_result);
            }
        }
        // 交换指针
        std::swap(U_new, const_cast<float*&>(U_old));
    }
}

int main() {
    const int M = 1024;
    const int N = 1024;
    const int ALIGNMENT = 64;
    const int ITERATIONS = 100;

    float* U_scalar_1 = static_cast<float*>(aligned_malloc(M * N * sizeof(float), ALIGNMENT));
    float* U_scalar_2 = static_cast<float*>(aligned_malloc(M * N * sizeof(float), ALIGNMENT));
    float* U_simd_1 = static_cast<float*>(aligned_malloc(M * N * sizeof(float), ALIGNMENT));
    float* U_simd_2 = static_cast<float*>(aligned_malloc(M * N * sizeof(float), ALIGNMENT));

    if (!U_scalar_1 || !U_scalar_2 || !U_simd_1 || !U_simd_2) {
        std::cerr << "Failed to allocate memory." << std::endl;
        return 1;
    }

    init_grid(U_scalar_1, M, N);
    init_grid(U_simd_1, M, N);
    // 复制初始状态到另一个网格，用于交换
    std::copy(U_scalar_1, U_scalar_1 + M*N, U_scalar_2);
    std::copy(U_simd_1, U_simd_1 + M*N, U_simd_2);

    auto start_scalar = std::chrono::high_resolution_clock::now();
    jacobi_scalar(U_scalar_2, U_scalar_1, M, N, ITERATIONS);
    auto end_scalar = std::chrono::high_resolution_clock::now();
    std::chrono::duration<double> diff_scalar = end_scalar - start_scalar;
    std::cout << "Scalar Jacobi iteration time: " << diff_scalar.count() << " s" << std::endl;

    auto start_simd = std::chrono::high_resolution_clock::now();
    jacobi_simd(U_simd_2, U_simd_1, M, N, ITERATIONS);
    auto end_simd = std::chrono::high_resolution_clock::now();
    std::chrono::duration<double> diff_simd = end_simd - start_simd;
    std::cout << "AVX-512 SIMD Jacobi iteration time: " << diff_simd.count() << " s" << std::endl;

    // 验证结果 (简化，只检查部分)
    // 最终的 U_scalar_1/2 和 U_simd_1/2 会因为迭代次数是偶数或奇数而不同
    // 这里假设 U_scalar_1 和 U_simd_1 是最终结果（如果 ITERATIONS 是偶数）
    // 否则应该是 U_scalar_2 和 U_simd_2
    bool correct = true;
    float* final_scalar_grid = (ITERATIONS % 2 == 0) ? U_scalar_1 : U_scalar_2;
    float* final_simd_grid = (ITERATIONS % 2 == 0) ? U_simd_1 : U_simd_2;

    for (int i = 1; i < M - 1; ++i) {
        for (int j = 1; j < N - 1; ++j) {
            // 浮点数比较需要容忍误差
            if (std::abs(final_scalar_grid[i * N + j] - final_simd_grid[i * N + j]) > 1e-4f) {
                correct = false;
                break;
            }
        }
        if (!correct) break;
    }
    std::cout << "Results are " << (correct ? "correct (within tolerance)" : "incorrect") << std::endl;

    aligned_free(U_scalar_1); aligned_free(U_scalar_2);
    aligned_free(U_simd_1); aligned_free(U_simd_2);

    return 0;
}

关键优化点：

• 内存访问模式：为了最大化 SIMD 效率，内层循环应该尽可能地访问连续的内存。在 2D 数组 grid[i * N + j] 中，改变 j 值是连续的，改变 i 值是跨行的。因此，将 j 作为内层循环的变量，然后进行向量化，可以实现高效的 _mm512_load_ps 和 _mm512_store_ps。
• 边界处理：模板操作通常不适用于网格边界。在循环中，我们从 i=1 到 M-2 和 j=1 到 N-2，跳过边界，因为边界条件通常是预设的，或者需要特殊处理，不适合通用模板。
• 双缓冲：使用 U_old 和 U_new 两个网格，并在每次迭代后交换它们，是处理依赖关系的常见模式，可以避免数据竞争和读写冲突。

案例分析二：加解密算法的 SIMD 优化

加解密算法通常涉及大量的位操作、查表、算术运算等，这些操作在 SIMD 指令集下也能得到显著加速，尤其是当需要并行处理多个数据块时。

4.1 AES 加密算法的核心循环

AES（Advanced Encryption Standard）是一种对称密钥分组密码，以 128 位数据块为单位进行加密。其核心操作包括：SubBytes（字节替换）、ShiftRows（行移位）、MixColumns（列混淆）和 AddRoundKey（轮密钥加）。其中，MixColumns 和 SubBytes 是主要的计算瓶颈。

MixColumns 优化

MixColumns 是一个矩阵乘法操作，它在伽罗瓦域 $GF(2^8)$ 上对 4×4 状态矩阵的每一列进行操作。
一个字节 $b$ 乘以常数 $c$ 在 $GF(2^8)$ 上的乘法可以表示为一系列移位和异或操作，或者通过查表实现。SIMD 在这里可以并行处理多个列或多个数据块的 MixColumns 操作。

AVX-512 的 VPCLMULQDQ（Carry-less Multiplication）和 GFNI (Galois Field New Instructions) 指令集可以直接加速伽罗瓦域乘法，但它们不是所有 AVX-512 处理器都支持的。更通用的方法是利用查表或位操作。对于查表，AVX-512 的 _mm512_i32gather_epi32 (32位整数的 gather 操作) 可以高效地从内存中收集非连续的查表结果。

这里我们演示一个简化的 AES 轮函数中的 MixColumns，假设我们已经准备好了查表数据。

C++ 示例：简化版 MixColumns (查表方式)

#include <iostream>
#include <vector>
#include <chrono>
#include <immintrin.h> // AVX-512 intrinsics

// 假设我们有预计算好的 MixColumns 查找表
// aes_mul2_table, aes_mul3_table, aes_mul1_table, aes_mul_inv_table 等
// 它们是 256 字节的数组，存储 GF(2^8) 乘法结果
// 为了简化，这里仅模拟一个简单的查表操作
// 实际的 AES MixColumns 涉及多个表的组合和异或

// 简化版查表函数，模拟 GF(2^8) 乘法
unsigned char gf_mul_table[256][4]; // 假设有四个乘法表，针对 0x01, 0x02, 0x03, 0x09, 0x0B, 0x0D, 0x0E

void init_gf_mul_table() {
    // 实际的 GF(2^8) 乘法表生成非常复杂，这里只是一个占位符
    // 假设 gf_mul_table[byte_val][factor_idx] = byte_val * factor_val in GF(2^8)
    for (int i = 0; i < 256; ++i) {
        // 伪代码，实际要计算 GF(2^8) 乘法
        gf_mul_table[i][0] = i; // 乘以 1
        gf_mul_table[i][1] = (i << 1) ^ ((i & 0x80) ? 0x1B : 0); // 乘以 2
        gf_mul_table[i][2] = gf_mul_table[i][1] ^ i; // 乘以 3 = 2x ^ x
        gf_mul_table[i][3] = gf_mul_table[gf_mul_table[gf_mul_table[i][1]][1]][1] ^ i; // 乘以 9 = 8x ^ x
    }
}

// 模拟一个 AES 状态块 (16 字节)
using AESBlock = unsigned char[16];

// 标量 MixColumns 对单个 4x4 状态矩阵的列操作
void mix_column_scalar(AESBlock& state) {
    // 简化：只处理一列，实际处理四列
    unsigned char s0 = state[0], s1 = state[1], s2 = state[2], s3 = state[3];

    // 假设 MixColumns 矩阵的第一行是 [02 03 01 01]
    // s0' = (02 * s0) ^ (03 * s1) ^ (01 * s2) ^ (01 * s3)
    state[0] = gf_mul_table[s0][1] ^ gf_mul_table[s1][2] ^ gf_mul_table[s2][0] ^ gf_mul_table[s3][0];
    // ... 其他行类似
    state[1] = gf_mul_table[s0][0] ^ gf_mul_table[s1][1] ^ gf_mul_table[s2][2] ^ gf_mul_table[s3][0];
    state[2] = gf_mul_table[s0][0] ^ gf_mul_table[s1][0] ^ gf_mul_table[s2][1] ^ gf_mul_table[s3][2];
    state[3] = gf_mul_table[s0][2] ^ gf_mul_table[s1][0] ^ gf_mul_table[s2][0] ^ gf_mul_table[s3][1];
}

// AVX-512 SIMD MixColumns (并行处理多个 AES 块)
// 这里需要处理多个 AES 块，每个块是 16 字节
// AVX-512 可以同时处理 512 / 8 = 64 字节，即 4 个 AES 块
void mix_columns_simd(AESBlock* states, int num_blocks) {
    const int BLOCKS_PER_VEC = 4; // 512 bits / (16 bytes/block * 8 bits/byte) = 4 blocks
    // 假设查表数据已经转换为 __m512i 友好的格式，或者使用 gather 指令

    // 我们需要将每个块的列数据提取出来，以便 SIMD 处理
    // 例如，处理所有块的第一列，然后所有块的第二列，以此类推
    // 或者，将每个块的四个字节打包成一个 32 位整数，然后操作

    // 由于 AES 的 MixColumns 跨字节操作，这里 SIMD 化需要更精细的数据重排
    // 对于 AVX-512，一个典型的优化是并行处理多个 AES 块。
    // 将 4 个 AES 块 (4 * 16 = 64 字节) 放入一个 __m512i 寄存器
    // 然后对每个字节执行查表和异或。
    // 但是，查表操作如果直接使用 gather，需要 32 位索引。
    // 字节级的 shuffle 和 permute 指令在 AVX-512 上非常强大。

    // 这里我们演示一种利用 gather 的思路，假设我们将查表结果打包成 32 位整数
    // 且每个字节的输入，通过查表得到 4 个字节的输出
    // 实际的 MixColumns 是一个 4x4 矩阵乘法，每个输出字节是 4 个输入字节的函数。

    // 更实际的 SIMD AES 实现通常使用 `AESENC` 指令（AES-NI），或通过位操作和 shuffle 模拟。
    // 如果没有 AES-NI，或者需要更底层的控制，可以利用 AVX-512 的通用指令。

    // 简化演示：假设我们并行处理 4 个块的某个字节的 MixColumns 变换
    // 例如，同时处理 states[0].col0.row0, states[1].col0.row0, states[2].col0.row0, states[3].col0.row0

    for (int i = 0; i < num_blocks; i += BLOCKS_PER_VEC) {
        if (i + BLOCKS_PER_VEC > num_blocks) { /* 处理剩余 */ continue; }

        // 加载 4 个 AES 块，共 64 字节
        __m512i v_state = _mm512_loadu_si512(states[i]); // _mm512_loadu_si512 允许非对齐加载

        // 复杂的 MixColumns 变换需要很多 shuffle 和 GF(2^8) 乘法。
        // 例如，对于状态矩阵的第一列 (s0, s1, s2, s3)
        // new_s0 = mul2(s0) ^ mul3(s1) ^ mul1(s2) ^ mul1(s3)
        // 这里需要将 s0, s1, s2, s3 从 v_state 中提取出来，并对 4 个块并行操作。

        // 提取并打包 4 个块的 s0 (第一个字节)
        __m512i v_s0 = _mm512_and_si512(v_state, _mm512_set1_epi8(0xFF)); // 假设 s0 在每个块的第一个字节
        // ... 这需要复杂的字节提取和重排

        // 使用 GFNI 指令会更直接，例如 _mm512_gf2p8affineinv_epi64_epi8

        // 由于 MixColumns 涉及 4x4 矩阵乘法，每个输出字节依赖于 4 个输入字节。
        // 最直接的 SIMD 优化是：
        // 1. 将 4 个 AES 块（共 64 字节）加载到 ZMM 寄存器
        // 2. 利用 `_mm512_shuffle_epi8` 进行字节重排，将同一“行”或“列”的字节汇聚
        // 3. 利用 `_mm512_permute_pd` (或者更通用的 _mm512_permutexvar_epi32)
        // 4. 使用 `_mm512_xor_si512` 进行异或操作
        // 5. 如果没有 GFNI，则使用查表。查表通过 `_mm512_i32gather_epi32` (或 8/16 位版本)
        //    需要将字节索引扩展为 32 位整数，然后用它作为 `gather` 的索引。
        //    _mm512_i32gather_epi32 需要 16 个 32 位索引，因此可以同时查 16 个表项。
        //    对于 4 个 AES 块，共有 64 个字节。需要多次 gather。

        // 这是一个高度复杂的 SIMD 优化示例，超出了本讲座直接展示完整代码的范围。
        // 主要思想是：
        // - 将多个 AES 块打包进 512 位寄存器。
        // - 利用 AVX-512 的 shuffle/permute 指令实现字节和字级别的数据重排，以适应 MixColumns 的矩阵操作。
        // - 利用 gather 指令高效执行 GF(2^8) 乘法查表，或利用 GFNI 指令集。
        // - 大量使用 _mm512_xor_si512 进行异或操作。

        // 简化存储：这里只是一个占位符，表示最终结果存储回内存
        _mm512_storeu_si512(states[i], v_state);
    }
}

int main() {
    init_gf_mul_table();
    const int NUM_BLOCKS = 1024;
    std::vector<AESBlock> scalar_states(NUM_BLOCKS);
    std::vector<AESBlock> simd_states(NUM_BLOCKS);

    // 初始化随机数据
    for (int i = 0; i < NUM_BLOCKS; ++i) {
        for (int j = 0; j < 16; ++j) {
            scalar_states[i][j] = static_cast<unsigned char>(rand() % 256);
            simd_states[i][j] = scalar_states[i][j];
        }
    }

    // 标量计算
    auto start_scalar = std::chrono::high_resolution_clock::now();
    for (int i = 0; i < NUM_BLOCKS; ++i) {
        mix_column_scalar(scalar_states[i]);
    }
    auto end_scalar = std::chrono::high_resolution_clock::now();
    std::chrono::duration<double> diff_scalar = end_scalar - start_scalar;
    std::cout << "Scalar MixColumns time (per column for each block): " << diff_scalar.count() << " s" << std::endl;

    // SIMD 计算
    auto start_simd = std::chrono::high_resolution_clock::now();
    mix_columns_simd(simd_states.data(), NUM_BLOCKS);
    auto end_simd = std::chrono::high_resolution_clock::now();
    std::chrono::duration<double> diff_simd = end_simd - start_simd;
    std::cout << "AVX-512 SIMD MixColumns time (parallel blocks): " << diff_simd.count() << " s" << std::endl;

    // 验证结果 (简化，因为 SIMD 部分是占位符)
    // 实际应比较 scalar_states 和 simd_states 的内容

    return 0;
}

4.2 哈希函数：以 SHA-256 为例

SHA-256 是一种密码学哈希函数，它将任意长度的数据映射为固定长度（256 位）的哈希值。其核心算法是一个消息调度器和 64 轮迭代。每轮迭代涉及大量的位操作（AND, OR, XOR, NOT）、右移、循环右移和加法。

SIMD 加速 SHA-256 的主要策略是并行处理多个消息块。即，同时计算 8 个或 16 个独立消息的 SHA-256 哈希值。

SHA-256 核心操作：

• 消息调度：将 16 个 32 位字扩展到 64 个 32 位字。涉及循环右移 (ROTR) 和异或。
• 压缩函数：64 轮迭代，每轮更新 8 个工作变量 (a-h)，涉及：
  • Ch (选择函数): (e & f) ^ (~e & g)
  • Maj (多数函数): (a & b) ^ (a & c) ^ (b & c)
  • $Sigma_0, Sigma_1, sigma_0, sigma_1$：循环右移和异或的组合。
  • 常数加法和模 $2^{32}$ 加法。

AVX-512 优化思路：

• 并行处理 16 个独立的 SHA-256 实例：将 16 个 32 位 SHA-256 状态变量（a, b, c, d, e, f, g, h）加载到 16 个 __m512i 寄存器中。
• 位操作：_mm512_and_si512, _mm512_or_si512, _mm512_xor_si512, _mm512_andnot_si512 可以并行执行 16 个 32 位整数的位操作。
• 循环右移：AVX-512 提供了 _mm512_rorv_epi32 (Vector Rotate Right) 指令，可以对向量中的每个 32 位元素执行独立的可变位移。这对于 SHA-256 的 $Sigma$ 和 $sigma$ 函数至关重要。
• 加法：_mm512_add_epi32 执行 16 个 32 位整数的并行加法。由于 SHA-256 的加法是模 $2^{32}$ 的，这直接对应于 CPU 的默认整数加法行为。
• 常量加载：_mm512_set1_epi32 用于广播 SHA-256 常量。
• 消息调度：需要将 16 个消息块的 32 位字排列成适合向量处理的格式。可以通过 _mm512_shuffle_epi32, _mm512_permutexvar_epi32 等指令进行高效重排。

C++ 示例：SHA-256 压缩函数中的部分 SIMD 步骤

#include <iostream>
#include <vector>
#include <array>
#include <chrono>
#include <immintrin.h> // AVX-512 intrinsics

// SHA-256 常量 (前 8 个素数的平方根的小数部分的前 32 位)
const std::array<uint32_t, 64> K = {
    0x428a2f98, 0x71374491, 0xb5c0fbcf, 0xe9b5dba5, 0x3956c25b, 0x59f111f1, 0x923f82a4, 0xab1c5ed5,
    0xd807aa98, 0x12835b01, 0x243185be, 0x550c7dc3, 0x72be5d74, 0x80deb1fe, 0x9bdc06a7, 0xc19bf174,
    // ... 更多常量，共 64 个
};

// 32 位整数的循环右移
inline uint32_t rotr(uint32_t val, int shift) {
    return (val >> shift) | (val << (32 - shift));
}

// SHA-256 的 Ch 和 Maj 函数 (标量)
inline uint32_t Ch(uint32_t e, uint32_t f, uint32_t g) { return (e & f) ^ (~e & g); }
inline uint32_t Maj(uint32_t a, uint32_t b, uint32_t c) { return (a & b) ^ (a & c) ^ (b & c); }

// SHA-256 的 Sigma0, Sigma1, sigma0, sigma1 函数 (标量)
inline uint32_t Sigma0(uint32_t a) { return rotr(a, 2) ^ rotr(a, 13) ^ rotr(a, 22); }
inline uint32_t Sigma1(uint32_t e) { return rotr(e, 6) ^ rotr(e, 11) ^ rotr(e, 25); }
inline uint32_t sigma0(uint32_t w) { return rotr(w, 7) ^ rotr(w, 18) ^ (w >> 3); }
inline uint32_t sigma1(uint32_t w) { return rotr(w, 17) ^ rotr(w, 19) ^ (w >> 10); }

// 标量 SHA-256 压缩函数的一个轮次（简化）
void sha256_round_scalar(std::array<uint32_t, 8>& H, uint32_t Wt, uint32_t Kt) {
    uint32_t a = H[0], b = H[1], c = H[2], d = H[3], e = H[4], f = H[5], g = H[6], h = H[7];

    uint32_t T1 = h + Sigma1(e) + Ch(e, f, g) + Kt + Wt;
    uint32_t T2 = Sigma0(a) + Maj(a, b, c);

    h = g;
    g = f;
    f = e;
    e = d + T1;
    d = c;
    c = b;
    b = a;
    a = T1 + T2;

    H[0] = a; H[1] = b; H[2] = c; H[3] = d; H[4] = e; H[5] = f; H[6] = g; H[7] = h;
}

// AVX-512 SHA-256 压缩函数的一个轮次（并行处理 16 个哈希状态）
void sha256_round_simd(__m512i H[8], __m512i v_Wt, __m512i v_Kt) {
    __m512i v_a = H[0], v_b = H[1], v_c = H[2], v_d = H[3], v_e = H[4], v_f = H[5], v_g = H[6], v_h = H[7];

    // Sigma1(e)
    __m512i v_sigma1_e = _mm512_rorv_epi32(v_e, _mm512_set1_epi32(6));
    v_sigma1_e = _mm512_xor_si512(v_sigma1_e, _mm512_rorv_epi32(v_e, _mm512_set1_epi32(11)));
    v_sigma1_e = _mm512_xor_si512(v_sigma1_e, _mm512_rorv_epi32(v_e, _mm512_set1_epi32(25)));

    // Ch(e, f, g) = (e & f) ^ (~e & g)
    __m512i v_Ch_efg = _mm512_and_si512(v_e, v_f);
    __m512i v_not_e = _mm512_xor_si512(v_e, _mm512_set1_epi32(0xFFFFFFFF)); // ~e
    v_Ch_efg = _mm512_xor_si512(v_Ch_efg, _mm512_and_si512(v_not_e, v_g));

    // T1 = h + Sigma1(e) + Ch(e, f, g) + Kt + Wt
    __m512i v_T1 = _mm512_add_epi32(v_h, v_sigma1_e);
    v_T1 = _mm512_add_epi32(v_T1, v_Ch_efg);
    v_T1 = _mm512_add_epi32(v_T1, v_Kt);
    v_T1 = _mm512_add_epi32(v_T1, v_Wt);

    // Sigma0(a)
    __m512i v_sigma0_a = _mm512_rorv_epi32(v_a, _mm512_set1_epi32(2));
    v_sigma0_a = _mm512_xor_si512(v_sigma0_a, _mm512_rorv_epi32(v_a, _mm512_set1_epi32(13)));
    v_sigma0_a = _mm512_xor_si512(v_sigma0_a, _mm512_rorv_epi32(v_a, _mm512_set1_epi32(22)));

    // Maj(a, b, c) = (a & b) ^ (a & c) ^ (b & c)
    __m512i v_Maj_abc = _mm512_and_si512(v_a, v_b);
    v_Maj_abc = _mm512_xor_si512(v_Maj_abc, _mm512_and_si512(v_a, v_c));
    v_Maj_abc = _mm512_xor_si512(v_Maj_abc, _mm512_and_si512(v_b, v_c));

    // T2 = Sigma0(a) + Maj(a, b, c)
    __m512i v_T2 = _mm512_add_epi32(v_sigma0_a, v_Maj_abc);

    // 更新 H 寄存器
    H[7] = v_g;
    H[6] = v_f;
    H[5] = v_e;
    H[4] = _mm512_add_epi32(v_d, v_T1);
    H[3] = v_c;
    H[2] = v_b;
    H[1] = v_a;
    H[0] = _mm512_add_epi32(v_T1, v_T2);
}

int main() {
    const int NUM_PARALLEL_HASHES = 16; // 同时处理 16 个哈希
    const int MESSAGE_BLOCK_SIZE_WORDS = 16; // 16 * 32-bit words = 64 bytes
    const int NUM_ROUNDS = 64;

    // 标量哈希状态
    std::vector<std::array<uint32_t, 8>> H_scalar_states(NUM_PARALLEL_HASHES);
    std::vector<std::array<uint32_t, MESSAGE_BLOCK_SIZE_WORDS>> W_scalar_messages(NUM_PARALLEL_HASHES);

    // SIMD 哈希状态 (8 个 ZMM 寄存器，每个存储 16 个 32 位整数)
    __m512i H_simd_states[8];
    // 假设消息调度后的 W 也在 SIMD 寄存器中
    __m512i W_simd_messages[MESSAGE_BLOCK_SIZE_WORDS];

    // 初始化哈希状态和消息 (简化)
    for (int i = 0; i < NUM_PARALLEL_HASHES; ++i) {
        // SHA-256 初始哈希值
        H_scalar_states[i] = {0x6a09e667, 0xbb67ae85, 0x3c6ef372, 0xa54ff53a,
                              0x510e527f, 0x9b05688c, 0x1f83d9ab, 0x5be0cd19};
        for (int j = 0; j < MESSAGE_BLOCK_SIZE_WORDS; ++j) {
            W_scalar_messages[i][j] = static_cast<uint32_t>(rand());
        }
    }

    // 将标量初始状态打包到 SIMD 寄存器
    for (int j = 0; j < 8; ++j) { // H[0]...H[7]
        alignas(64) uint32_t temp_H_values[NUM_PARALLEL_HASHES];
        for (int i = 0; i < NUM_PARALLEL_HASHES; ++i) {
            temp_H_values[i] = H_scalar_states[i][j];
        }
        H_simd_states[j] = _mm512_load_epi32(temp_H_values);
    }
    // 假设 W_simd_messages 也已类似方式打包和调度

    // 标量计算（一个消息块的完整压缩过程）
    auto start_scalar = std::chrono::high_resolution_clock::now();
    for (int i = 0; i < NUM_PARALLEL_HASHES; ++i) { // 模拟并行计算 NUM_PARALLEL_HASHES 个哈希
        // 实际的 SHA-256 消息调度在这里
        // std::array<uint32_t, 64> W_expanded;
        // for (int t = 0; t < 16; ++t) W_expanded[t] = W_scalar_messages[i][t];
        // for (int t = 16; t < 64; ++t) {
        //     W_expanded[t] = sigma1(W_expanded[t-2]) + W_expanded[t-7] + sigma0(W_expanded[t-15]) + W_expanded[t-16];
        // }

        std::array<uint32_t, 8> current_H = H_scalar_states[i];
        for (int t = 0; t < NUM_ROUNDS; ++t) {
            sha256_round_scalar(current_H, W_scalar_messages[i][t % MESSAGE_BLOCK_SIZE_WORDS], K[t]); // 简化 Wt
        }
        H_scalar_states[i] = current_H; // 更新最终哈希值
    }
    auto end_scalar = std::chrono::high_resolution_clock::now();
    std::chrono::duration<double> diff_scalar = end_scalar - start_scalar;
    std::cout << "Scalar SHA-256 rounds time (for " << NUM_PARALLEL_HASHES << " hashes): " << diff_scalar.count() << " s" << std::endl;

    // SIMD 计算
    auto start_simd = std::chrono::high_resolution_clock::now();
    for (int t = 0; t < NUM_ROUNDS; ++t) {
        // 假设 v_Wt 和 v_Kt 已经准备好，包含 16 个消息块的 Wt 和 K[t]
        // 这里只是一个简化，实际需要复杂的 SIMD 消息调度来生成 v_Wt
        __m512i v_Wt = _mm512_set1_epi32(W_scalar_messages[0][t % MESSAGE_BLOCK_SIZE_WORDS]); // 简化：用第一个消息的 Wt
        __m512i v_Kt = _mm512_set1_epi32(K[t]);

        sha256_round_simd(H_simd_states, v_Wt, v_Kt);
    }
    auto end_simd = std::chrono::high_resolution_clock::now();
    std::chrono::duration<double> diff_simd = end_simd - start_simd;
    std::cout << "AVX-512 SIMD SHA-256 rounds time (for " << NUM_PARALLEL_HASHES << " hashes): " << diff_simd.count() << " s" << std::endl;

    // 验证结果 (将 SIMD 结果解包并与标量结果比较)
    // ...

    return 0;
}

关键优化点：

• 并行化粒度：SHA-256 的 SIMD 优化通常在处理多个独立消息块上实现。
• 位操作指令：AVX-512 提供了完整的位操作指令，包括 AND, OR, XOR, NOT。
• 循环移位：_mm512_rorv_epi32 是一个非常关键的指令，它允许对向量中每个 32 位元素执行独立的循环右移操作，这在 SHA-256 的各种 $Sigma$ 和 $sigma$ 函数中大量使用。
• 加法：_mm512_add_epi32 执行并行 32 位整数加法，自然符合 SHA-256 的模 $2^{32}$ 加法。
• 数据打包与解包：将 16 个 SHA-256 状态变量（a-h）的数据打包到 8 个 __m512i 寄存器中，每个寄存器存储 16 个不同哈希实例的同一状态变量。

AVX-512 性能优化的关键考量与最佳实践

利用 AVX-512 实现高性能需要周密的规划和细致的优化。

5.1 数据对齐与内存访问模式

AVX-512 寄存器是 512 位宽，即 64 字节。对于 _mm512_load_ps、_mm512_store_ps 等指令，它们要求内存地址 64 字节对齐。未对齐访问会严重降低性能，甚至可能引发崩溃（如果硬件不允许）。

• 使用 alignas(64) 关键字定义变量或数组。
• 使用 _aligned_malloc / aligned_alloc 等函数分配对齐内存。
• 如果无法保证对齐，使用 _mm512_loadu_ps (u for unaligned) 等非对齐加载/存储指令，但性能可能受损。
• 内存访问模式：尽量保持内存访问连续性和局部性，以最大化缓存命中率。SoA 和 AoSoA 结构在这方面表现优异。

5.2 缓存层级与局部性

SIMD 优化往往会增加内存带宽需求。如果数据无法很好地适应 CPU 缓存（L1, L2, L3），处理器可能会频繁地从主内存读取数据，导致性能瓶颈从计算转移到内存带宽。

• 数据平铺 (Tiling)：将大问题分解为小块，使每个小块的数据都能放入缓存，处理完一个块再处理下一个。
• 算法选择：选择那些具有良好数据局部性的算法。

5.3 分支预测与指令吞吐

• 条件分支：SIMD 单元通常不喜欢条件分支。在循环内部的 if/else 语句会中断 SIMD 流水线，导致性能下降。AVX-512 的掩码寄存器 (K 寄存器) 是解决这个问题的一个强大工具。它允许你对向量操作的每个元素进行条件执行，而无需实际的分支跳转。例如 _mm512_mask_add_ps 只对掩码为真的元素执行加法。
• 指令吞吐：关注 CPU 的微架构，了解不同指令的延迟和吞吐量。FMA 指令（融合乘加）是一个典型的例子，它能在一个周期内完成乘法和加法，显著提高浮点运算吞吐。

5.4 处理剩余元素 (Remainder Handling)

当数据总量不是向量宽度（例如 16 个 float）的整数倍时，循环末尾会有一些剩余元素。

• 标量处理：最简单的方法是在 SIMD 循环结束后，用一个普通的标量循环处理剩余的元素。
• 掩码处理：使用 AVX-512 的掩码加载/存储指令 (_mm512_mask_load_ps, _mm512_mask_store_ps) 或掩码算术指令 (_mm512_mask_add_ps)，可以只对有效元素进行操作，而对剩余的“无效”元素不产生影响，从而保持 SIMD 循环的统一性。

5.5 编译器选项与构建环境

• 启用 AVX-512 指令集：需要向编译器传递特定的标志。
  • GCC/Clang: -march=skylake-avx512 (或更具体的如 -mavx512f -mavx512dq -mavx512cd 等，skylake-avx512 包含了基础的 AVX-512 指令集)
  • MSVC: /arch:AVX512
• 优化级别：使用 -O3 (GCC/Clang) 或 /O2 (MSVC) 等高优化级别，以启用更积极的编译优化和自动