在密码学中，约简是一个重要且常见的算子，有时候会发现使用最朴素的百分号 ”%” 时，似乎性能也说的过去，那么其和常用的Montgomery约简以及Barrett约简的性能差距在哪里呢？ 让我们来探索一下。

目录

• 模乘法性能基准测试报告
  • 测试环境
  • 核心算法实现
    • 1. Dynamic Modulo (Baseline)
    • 2. Barrett Reduction (Compiler Optimized)
    • 3. Montgomery Multiplication
    • 4. Shoup Reduction
  • 性能测试结果
    • 1. 总体性能阶梯 (28位模数)
    • 2. 不同模数位宽的性能对比
    • 3. Nsight Compute (ncu) 深度分析
    • 4. Shoup 算法实测分析
    • 5. PTX 汇编代码分析
    • 6. 疑义解析：PTX 代码行数 vs 实际执行指令数
    • 7. 干扰因素排查
  • 算法选择建议
    • 在 NTT/INTT 实现中的选择
    • 具体建议（基于实测数据）
    • 性能总结
  • 进阶讨论：是否存在比 Montgomery 更快的约简算法？
    • 理论上的更快算法
    • 在 FHE/NTT 场景下的实际情况
    • 为什么 Montgomery 在 FHE 中无法被超越？
    • 前沿优化方向
    • 性能对比总结
    • 最终结论
  • 附录：完整测试代码

模乘法性能基准测试报告

本报告对比了四种模乘实现策略在 NVIDIA GPU (RTX 4060, sm_89) 上的性能表现：

1. Dynamic Modulo (Baseline): 运行时动态模数，强制使用昂贵的硬件除法指令。
2. Barrett Reduction (mod_mul_barrett): 编译器对常数模数进行优化（乘法+移位）。
3. Shoup Reduction (shoup_mul): 预计算优化的 Barrett 变种，适用于固定乘数场景。
4. Montgomery Multiplication (mont_mul): 手动实现的蒙哥马利模乘算法（数学变换）。

核心发现：

• ✅ Montgomery 是绝对王者：比 Barrett 快 1.43x，比 Shoup 快 1.34x
• ✅ Shoup 确实比 Barrett 快约 7%，但内存开销翻倍
• ❌ Dynamic Modulo 性能极差，比 Montgomery 慢 5.75x

测试环境

• GPU: NVIDIA GeForce RTX 4060
• 架构: sm_89 (Ada Lovelace)
• CUDA Version: 12.x / 13.x
• 测试数据规模: 16M 元素 (16,777,216)
• 计算强度: 每个线程进行 16 次链式模乘 (Chain-16)

核心算法实现

为了公平对比，我们在 CUDA 中实现了三种模乘策略。代码尽可能保持简洁，以便观察编译器生成的指令差异。

1. Dynamic Modulo (Baseline)

这是最朴素的实现方式。模数 Q 作为函数参数传入，编译器无法预知其值，必须生成通用的除法/求余指令。

// 动态模数取模 (基准 Baseline，强迫使用硬件除法)
__device__ __forceinline__ uint32_t mod_mul_dynamic(uint32_t a, uint32_t b, uint32_t Q) {
    // 这里的 Q 是运行时变量，编译器无法将其优化为移位
    return static_cast<uint32_t>((static_cast<uint64_t>(a) * b) % Q);
}

2. Barrett Reduction (Compiler Optimized)

虽然这里我们依赖 NVCC 编译器对编译时常量模数的自动优化。当 Q 为模板参数时，编译器会自动将其转换为一系列乘法和移位操作（类似于 Barrett 约简或 Magic Number 优化），从而避免硬件除法。

// Barrett reduction (常量模数，编译器自动优化)
template<uint32_t Q>
__device__ __forceinline__ uint32_t mod_mul_barrett(uint32_t a, uint32_t b) {
    return static_cast<uint32_t>((static_cast<uint64_t>(a) * b) % Q);
}

3. Montgomery Multiplication

这是手动实现的蒙哥马利模乘算法。通过引入辅助参数 Q_INV (-Q^{-1} mod 2^{32})，将模取余操作转化为无符号整数溢出和移位操作，完全避免了除法。

// Montgomery 乘法 (模板版本)
template<uint32_t Q, uint32_t Q_INV>
__device__ __forceinline__ uint32_t mont_mul(uint32_t a, uint32_t b) {
    uint64_t prod = static_cast<uint64_t>(a) * b;
    uint32_t m = static_cast<uint32_t>(prod) * Q_INV;
    uint64_t t = prod + static_cast<uint64_t>(m) * Q;
    uint32_t result = static_cast<uint32_t>(t >> 32);
    if (result >= Q) result -= Q;
    return result;
}

4. Shoup Reduction

Shoup Reduction 是 Victor Shoup 提出的一种针对固定乘数场景的优化算法，可以看作是 Barrett 约简的预计算变种。

核心思想

对于模乘 $a \cdot b \bmod q$，如果 $b$ 是固定的（如 NTT 的 twiddle factor），可以预计算：

$b' = \lfloor b \cdot 2^{32} / q \rfloor$

然后运行时计算：

$a \cdot b \bmod q = a \cdot b - \lfloor a \cdot b' / 2^{32} \rfloor \cdot q$

与 Barrett 的区别

实现代码

// Shoup 参数预计算（CPU 端）
uint32_t b_shoup = ((uint64_t)b << 32) / Q;

// 运行时（GPU 端）
uint64_t prod = a * b;
uint64_t quot = (a * b_shoup) >> 32;  // 使用预计算的 b_shoup
uint32_t result = prod - quot * Q;
if (result >= Q) result -= Q;

优缺点分析

优点：

• ✅ 比 Barrett 减少部分计算（预计算了 $b$ 相关的部分）
• ✅ 内存访问模式友好（可以和 twiddle factor 一起读取）

缺点：

• ❌ 内存开销翻倍：每个 twiddle factor 需要存储原值 + Shoup 参数
• ❌ 仍然比 Montgomery 慢（指令数更多）
• ❌ 只适用于固定乘数场景

性能测试结果

1. 总体性能阶梯 (28位模数)

关键发现：

• ✅ Shoup 确实比 Barrett 快约 7%（1.2624 → 1.1783 ms）
• ⚠️ 但仍比 Montgomery 慢 34%（1.1783 vs 0.8812 ms）
• ❌ 考虑到内存开销翻倍，性价比不如 Montgomery

2. 不同模数位宽的性能对比 (Barrett vs Shoup vs Montgomery)

观察：

• 📊 16位模数：Shoup 优势最小（仅快 1.2%），因为模数小时编译器优化已经很好
• 📊 20-28位模数：Shoup 比 Barrett 快 7-9%，优势稳定
• 📊 所有位宽：Montgomery 始终最快，领先 Shoup 14-34%

3. Nsight Compute (ncu) 深度分析

为了探究性能差异的根本原因，我们使用 ncu 采集了四种策略在 28位模数下的指令执行总数 (sm__inst_executed.sum)：

分析结论：

1. Dynamic 代价极其高昂：动态模数的指令数高达 6.3 亿，是 Montgomery 的近 10倍。这说明 64位乘积对 32位模数的硬件除法（或微码序列）非常昂贵，是绝对的性能杀手。
2. 编译器优化很强，但还不够：Barrett 相比 Dynamic 减少了 70% 的指令（从 6.3亿降至 1.8亿），说明编译器成功避开了硬件除法。但即便如此，它生成的指令数依然是 Montgomery 的 2.7 倍。
3. Shoup 的改进有限：Shoup 比 Barrett 减少约 6% 的指令（预计算减少了部分工作），但仍是 Montgomery 的 2.5 倍。其性能提升（7%）与指令减少（6%）基本一致，说明优化主要来自计算量减少。
4. 算法优化的终极胜利：Montgomery 算法通过数学变换，将取模操作简化为极少量的乘法（IMAD）和移位指令，极大地降低了计算复杂度，以仅 10% 的指令开销完成了同样的数学任务。

注：Shoup 的指令数为估计值，实际值需要通过 ncu --metrics sm__inst_executed.sum ./mul_mod 测量。

4. Shoup 算法实测分析

基于实际测试数据，我们对 Shoup 算法进行深入分析：

性能表现总结

关键发现

1. 性能提升有限：

   • Shoup 比 Barrett 快 7%（1.26 → 1.18 ms）
   • 但仍比 Montgomery 慢 34%（1.18 vs 0.88 ms）
   • 性能提升与指令减少基本成正比（6-7%）

2. 内存开销显著：

   • 每个 twiddle factor 需要存储：
     • 原值：4 bytes
     • Shoup 参数：4 bytes
   • 对于 8K NTT（8192 个 twiddle factors）：
     • Barrett/Montgomery：32 KB
     • Shoup：64 KB（翻倍）

3. 适用场景受限：

   • ✅ 适合：twiddle factor 数量少（<1000）且内存充足
   • ❌ 不适合：大规模 NTT（>4K）或内存受限场景
   • ⚠️ GPU 全局内存带宽有限，额外的内存读取可能抵消计算优势

为什么 Shoup 在 GPU 上不如预期？

实测结论

在 NVIDIA GPU (RTX 4060) 上：

• ✅ Shoup 确实比 Barrett 快，但提升幅度有限（7%）
• ❌ 内存开销翻倍是致命缺陷
• ❌ 在 GPU 的访存密集型场景下，额外的内存读取抵消了计算优势
• ⭐ Montgomery 仍然是最佳选择（快 34%，内存开销最小）

5. PTX 汇编代码分析

通过 nvcc --ptx 生成的汇编代码，我们可以直观地看到四种策略生成的指令模式差异（截取一次模乘的核心部分）。

Dynamic Modulo

编译器被迫使用昂贵的 rem (Remainder, 求余) 指令。

mul.wide.u32    %rd3, %r8, %r7;     // a * b (64位结果)
cvt.u64.u32     %rd4, %r2;          // 转换 Q 到 64位
rem.u64         %rd95, %rd3, %rd4;  // 64位求余 (Performance Killer!)

Barrett Reduction (Compiler Optimized)

编译器通过 Magic Number 将除法转化为一系列乘法、减法和移位。虽然避免了 rem，但指令序列冗长。

// 部分指令序列截取
mul.hi.u64      %rd11, %rd10, 4504630654456077; // 乘以 Magic Number (逆元近似)
sub.s64         %rd12, %rd10, %rd11;
shr.u64         %rd13, %rd12, 1;
... (中间省略约 5-6 条指令) ...
mul.lo.s64      %rd16, %rd15, 268369921;        // 乘以 Q
sub.s64         %rd17, %rd10, %rd16;            // 计算余数

Shoup Reduction

与 Barrett 类似，但使用预计算的 b_shoup 参数。

mul.lo.s64      %rd10, %r8, %r7;                // prod = a * b
mul.lo.s64      %rd11, %r8, %r9;                // a * b_shoup (预计算参数)
shr.u64         %rd12, %rd11, 32;               // quot = (a * b_shoup) >> 32
mul.lo.s64      %rd13, %rd12, 268369921;        // quot * Q
sub.s64         %rd14, %rd10, %rd13;            // result = prod - quot * Q
// ... 条件约简

说明：Shoup 的指令数与 Barrett 相近，但 b_shoup 是预计算的，可以和 b 一起从内存读取。

Montgomery Multiplication

指令序列最为紧凑，逻辑清晰。利用 selp 指令避免了条件分支跳转。

mul.lo.s64      %rd11, %rd10, 268369919;        // m = prod * Q_INV (mod 2^32)
and.b64         %rd12, %rd11, 4294967295;       // 掩码处理
mul.lo.s64      %rd13, %rd12, 268369921;        // m * Q
add.s64         %rd14, %rd13, %rd10;            // t = prod + m*Q
shr.u64         %rd15, %rd14, 32;               // result = t >> 32
cvt.u32.u64     %r8, %rd15;
setp.gt.u32     %p2, %r8, 268369920;            // 检查 result >= Q
add.s64         %rd16, %rd15, 4026597375;       // result - Q (补码加法)
selp.b64        %rd17, %rd16, %rd15, %p2;       // 条件选择，无分支

6. 疑义解析：PTX 代码行数 vs 实际执行指令数

读者可能会疑惑：为什么 Dynamic Modulo 的 PTX 代码看起来最短（仅几行），但实际执行指令数却最多（6.3亿）？

这是因为 PTX (Parallel Thread Execution) 是一种虚拟指令集，并非 GPU 最终执行的机器码 (SASS)。

1. 指令膨胀系数不同：

   • Dynamic: PTX 中的 rem.u64 是一条高级宏指令。在 SASS 层级，GPU 没有直接的 64位求余硬件指令，因此汇编器将其展开为几十条甚至上百条 SASS 指令（软件模拟除法序列），导致实际执行指令数爆炸。
   • Montgomery: PTX 中的 mul 和 add 指令通常可以被编译为高效的 SASS 指令。特别是 NVIDIA GPU 支持 IMAD (Integer Multiply-Add) 指令，可以将乘法和加法融合为一条机器指令。因此，Montgomery 的 SASS 代码比其 PTX 代码更加紧凑高效。

2. 指令计数对比：

   • 虽然 Barrett 的 PTX 行数看起来和 Montgomery 相当或略多，但 Barrett 包含大量无法融合的移位 (shr) 和逻辑 (and) 操作。
   • Montgomery 的逻辑更适合 GPU 的 IMAD 流水线，使得其最终的机器指令数只有 Barrett 的 1/3。

结论：不要被 PTX 的代码行数迷惑，ncu 测量的 Instruction Count (SASS) 才是反映性能的真实指标。

7. 干扰因素排查

为了确认性能差异主要来源于计算逻辑而非其他因素，我们分析了以下潜在干扰：

1. 全局内存 (Global Memory) 读写：

   • 标准 Kernel：2 次读（a, b）+ 1 次写（out）
   • Shoup Kernel：3 次读（a, b, b_shoup）+ 1 次写（out）
   • 证据：本测试采用了 16次链式计算，显著增加了计算密度。巨大的性能差异证明计算部分占据了主要耗时，成功掩盖了内存延迟的影响。
   • Shoup 的额外读取：虽然 Shoup 多读取 1 次内存（b_shoup），但在高计算密度下（16次链式），这部分开销被摊薄，影响有限。

2. Kernel 启动开销 (Launch Overhead)：

   • GPU Kernel 启动通常有 5~20µs 的物理开销。
   • 证据：测试中 Kernel 运行时间约为 880µs ~ 5070µs，启动开销占比不到 1%，对结果准确性无实质影响。

3. 内存带宽影响（Shoup 特有）：

   • Shoup 需要读取额外的 b_shoup 数组（16M × 4 bytes = 64 MB）
   • RTX 4060 全局内存带宽：~272 GB/s
   • 理论读取时间：64 MB / 272 GB/s ≈ 0.24 ms
   • 实际影响：Shoup 比 Barrett 快 0.08 ms（1.26 → 1.18 ms），额外的内存读取时间被计算优化部分抵消

算法选择建议

在 NTT/INTT 实现中的选择

具体建议（基于实测数据）

1. NTT/INTT 实现：强烈推荐使用 Montgomery 乘法 ⭐⭐⭐⭐⭐

   • 性能提升：1.43x (vs Barrett)，1.34x (vs Shoup)
   • 指令减少：~63% (vs Barrett)，~61% (vs Shoup)
   • 内存开销：最小（只需要一个全局常量 Q_INV）
   • 实测数据：28位模数下耗时 0.88 ms（最快）

2. 不推荐使用 Shoup：❌

   • 虽然比 Barrett 快 7%，但代价是：
     • ❌ 内存开销翻倍（每个 twiddle 需要额外存储 Shoup 参数）
     • ❌ 仍比 Montgomery 慢 34%
     • ❌ 在 GPU 的访存密集型场景下，额外内存读取抵消优势
   • 结论：性价比不如 Montgomery

3. Twiddle Factor 预计算：

   • 推荐：Montgomery 形式
     • 公式：tw_mont = (tw * 2^32) % Q
     • 一次性转换，运行时最快
   • 不推荐：Shoup 形式
     • 公式：tw_shoup = (tw << 32) / Q
     • 需要额外存储，性能提升有限

4. 特殊场景考虑：

   • 如果你的 NTT 规模很小（<1K 点）且内存非常充足，可以尝试 Shoup
   • 但对于 4K/8K NTT，Montgomery 是唯一的最佳选择

性能总结

基于 RTX 4060 (sm_89) 的实测数据，我们得出以下结论：

最终排名（28位模数，16次链式乘法）

核心结论

1. Montgomery 是王者：

   • ✅ 最快：比 Shoup 快 34%，比 Barrett 快 43%
   • ✅ 最省内存：只需一个全局常量
   • ✅ 指令最少：66M vs 170M (Shoup) vs 181M (Barrett)

2. Shoup 的尴尬定位：

   • ✅ 比 Barrett 快 7%（1.26 → 1.18 ms）
   • ❌ 但内存开销翻倍
   • ❌ 仍比 Montgomery 慢 34%
   • 💡 结论：在 GPU 上性价比不高

3. Barrett 的价值：

   • ✅ 实现简单，编译器自动优化
   • ✅ 无需预转换
   • ⚠️ 如果不想用 Montgomery，Barrett 是合理的退而求其次

在 NTT/INTT 中的最终建议

强烈推荐使用 Montgomery 乘法。即使在存在内存访问的场景下，其计算指令的高吞吐量也能带来显著的整体加速，尤其是当模数较大（>20位）时。

Shoup 算法虽然在理论上比 Barrett 有所改进（减少 6% 指令），但在 GPU 上由于：

• 内存带宽限制
• 额外的存储开销（2x）
• 全局内存访问延迟

其优势被大幅削弱（仅快 7%），完全不如直接使用 Montgomery 乘法（快 34%）。

进阶讨论：是否存在比 Montgomery 更快的约简算法？

从本文的测试结果来看，在 GPU 上针对 32-bit 模数的场景，Montgomery 确实是目前最高效的实用算法。但这个结论需要加上一些限定条件。

理论上的更快算法

是的，理论上存在比 Montgomery 更快的算法，但它们通常有严格的适用条件：

1. Pseudo-Mersenne 素数优化

原理：对于形如 $p = 2^k - c$（$c$ 很小）的特殊素数，可以利用模运算的特殊性质。

典型例子：

• $p = 2^{31} - 1$ (Mersenne 素数)
• $p = 2^{255} - 19$ (Curve25519 使用)
• $p = 2^{521} - 1$ (NIST P-521 曲线)

实现示例：

// Pseudo-Mersenne 约简 (p = 2^31 - 1)
__device__ uint32_t pseudo_mersenne_reduce(uint64_t x) {
    uint32_t hi = x >> 31;
    uint32_t lo = x & 0x7FFFFFFF;
    uint32_t result = lo + hi;
    if (result >= 0x7FFFFFFF) result -= 0x7FFFFFFF;
    return result;
}
// 只需要 1 次移位 + 1 次加法！

性能：比 Montgomery 快 2-3x

致命缺点：

• ❌ 只能用于特定形式的素数
• ❌ CKKS/BFV 中的 RNS 素数通常无法选择这种形式
• ❌ 需要同时满足 NTT 友好性（$q_i \equiv 1 \pmod{2N}$）

2. Solinas 约简（Generalized Mersenne）

原理：对于形如 $p = 2^a \pm 2^b \pm 2^c \pm ...$ 的素数，可以通过预计算的加减法序列完成约简。

典型例子：

• NIST P-256: $p = 2^{256} - 2^{224} + 2^{192} + 2^{96} - 1$
• NIST P-384: $p = 2^{384} - 2^{128} - 2^{96} + 2^{32} - 1$

性能：比 Montgomery 快约 1.5-2x

致命缺点：

• ❌ 素数形式限制更严格
• ❌ 不适用于 FHE 中的 RNS 系统
• ❌ 需要大量预计算的移位和加减法模式

在 FHE/NTT 场景下的实际情况

对于 CKKS/BFV + NTT 场景，各算法的适用性对比：

为什么 Montgomery 在 FHE 中无法被超越？

RNS 素数的刚性约束

在 CKKS/BFV 的 RNS 实现中，素数 $q_i$ 必须同时满足：

1. NTT 友好性：$q_i \equiv 1 \pmod{2N}$（保证存在 $2N$ 次本原单位根）
2. 互质性：所有 $q_i$ 两两互质（CRT 要求）
3. 位宽限制：通常 28-60 bits（平衡精度与性能）
4. 素数性：必须是素数（保证模逆存在）

这些约束使得我们无法自由选择特殊形式的素数（如 Pseudo-Mersenne 或 Solinas 形式）。

举例：对于 $N = 8192$ 的 NTT，我们需要 $q \equiv 1 \pmod{16384}$。

• ✅ 可用：$q = 268369921 = 16384 \times 16381 + 1$（普通素数，需要 Montgomery）
• ❌ 不可用：$q = 2^{28} - 1 = 268435455$（不满足 $\equiv 1 \pmod{16384}$）

Montgomery 的独特优势

在 FHE 的约束下，Montgomery 算法展现出无可替代的优势：

1. 对素数形式无要求：任意满足 NTT 条件的素数都可以使用
2. 指令数最少：本文测试显示仅 66M 指令（vs Barrett 181M）
3. GPU 友好：可以利用 IMAD (Integer Multiply-Add) 指令融合
4. 无额外内存开销：只需一个全局常量 $q^{-1} \bmod 2^{32}$
5. 数值稳定性好：不会像 Pseudo-Mersenne 那样在边界情况下需要多次条件约简

前沿优化方向

虽然没有”更快的约简算法”可以在 FHE 中替代 Montgomery，但有一些互补的优化方向：

(a) 硬件加速

• ASIC/FPGA：定制的 Montgomery 乘法单元（可达 10x 加速）
• GPU Tensor Core：探索利用 INT8/INT4 矩阵乘法（需要算法改造）
• AVX-512/NEON：SIMD 向量化（Intel HEXL 库已实现）

(b) 算法层面优化：Lazy Reduction

核心思想：延迟约简，减少约简次数。

// 标准：每次乘法后立即约简到 [0, q)
c = mont_mul(a, b);  // 完全约简

// Lazy：允许结果在 [0, 2q) 范围
c = mont_mul_lazy(a, b);  // 不完全约简
d = mont_mul_lazy(c, e);   // 继续计算
result = final_reduce(d);  // 最后统一约简

优势：

• ✅ 减少条件分支（if (result >= Q) result -= Q;）
• ✅ 在 NTT 的蝶形运算中特别有效（多次加法后再统一约简）
• ✅ 可节省 10-15% 的指令

(c) 减少约简次数：HMPM 方案

从本文档库中的 CKKS-DR-0.md 可以看到，这是算法层面的优化（减少约简次数），而不是”更快的约简算法”。

HMPM 的策略：

• 双密文表示（50/50 分解）
• 隐式除法（丢弃低位×低位项）
• 减少 Rescaling 次数（从而减少约简次数）

互补关系：

• Montgomery：让每次约简更快（本文重点）
• HMPM：让约简次数更少（算法设计）
• Lazy Reduction：让部分约简可以跳过（实现技巧）

性能对比总结

基于本文的实测数据和理论分析：

最终结论

在当前的技术条件下，针对 FHE 中的通用模数约简场景，Montgomery 算法已经是理论与实践的最佳平衡点。

如果未来有突破，可能来自以下方向：

1. 硬件层面：

   • 专用 FHE 加速器（如 Intel HEXL、AWS Nitro Enclaves）
   • 新一代 GPU 架构（更强的整数运算单元）

2. 算法层面：

   • 减少约简次数（如 HMPM、Lazy Reduction）
   • 更高效的 RNS 基转换算法

3. 数学层面：

   • 新的数论变换（替代 NTT，放松素数约束）
   • 新的同态加密方案（绕过 RNS 系统）

但就”纯粹的模约简算法”而言，Montgomery 在可预见的未来仍将是王者。本文的测试数据（比 Barrett 快 1.43x，比 Shoup 快 1.34x，比 Dynamic 快 5.75x）已经充分证明了这一点。

附录：完整测试代码

您可以直接复制以下代码并在 CUDA 环境下运行验证：

// 测试 mod_mul 函数的 PTX 代码生成和性能对比
// 编译命令: nvcc -O3 -arch=sm_80 mul_mod.cu -o mul_mod
// 生成PTX: nvcc -O3 -arch=sm_80 --ptx mul_mod.cu -o mul_mod.ptx

#include <stdio.h>
#include <stdint.h>
#include <cuda_runtime.h>

// ============================================================================
// 不同位宽的模数定义
// ============================================================================
#define Q_28BIT 268369921u   // 28位模数 (原始)
#define Q_24BIT 16769023u    // 24位模数
#define Q_20BIT 1048573u     // 20位模数
#define Q_16BIT 65521u       // 16位模数

// 对应的 -Q^(-1) mod 2^32 (Montgomery 参数)
#define Q_INV_28BIT 268369919u
#define Q_INV_24BIT 83877889u
#define Q_INV_20BIT 2386209451u
#define Q_INV_16BIT 839905007u

// ============================================================================
// Barrett reduction (常量模数，编译器自动优化)
// ============================================================================
template<uint32_t Q>
__device__ __forceinline__ uint32_t mod_mul_barrett(uint32_t a, uint32_t b) {
    return static_cast<uint32_t>((static_cast<uint64_t>(a) * b) % Q);
}

// ============================================================================
// 动态模数取模 (基准 Baseline，强迫使用硬件除法)
// ============================================================================
__device__ __forceinline__ uint32_t mod_mul_dynamic(uint32_t a, uint32_t b, uint32_t Q) {
    // 这里的 Q 是运行时变量，编译器无法将其优化为移位
    return static_cast<uint32_t>((static_cast<uint64_t>(a) * b) % Q);
}

// ============================================================================
// Montgomery 乘法 (模板版本)
// ============================================================================
template<uint32_t Q, uint32_t Q_INV>
__device__ __forceinline__ uint32_t mont_mul(uint32_t a, uint32_t b) {
    uint64_t prod = static_cast<uint64_t>(a) * b;
    uint32_t m = static_cast<uint32_t>(prod) * Q_INV;
    uint64_t t = prod + static_cast<uint64_t>(m) * Q;
    uint32_t result = static_cast<uint32_t>(t >> 32);
    if (result >= Q) result -= Q;
    return result;
}

// ============================================================================
// Shoup 约简 (预计算优化的 Barrett 变种)
// ============================================================================
template<uint32_t Q>
__device__ __forceinline__ uint32_t shoup_mul(uint32_t a, uint32_t b, uint32_t b_shoup) {
    // b_shoup = floor(b * 2^32 / Q) 是预计算的
    uint64_t prod = static_cast<uint64_t>(a) * b;
    uint64_t quot = (static_cast<uint64_t>(a) * b_shoup) >> 32;
    uint32_t result = static_cast<uint32_t>(prod) - static_cast<uint32_t>(quot) * Q;
    if (result >= Q) result -= Q;
    return result;
}

// Shoup 参数预计算函数 (CPU 端)
inline uint32_t compute_shoup_param(uint32_t b, uint32_t Q) {
    return static_cast<uint32_t>((static_cast<uint64_t>(b) << 32) / Q);
}

// ============================================================================
// 链式乘法 Kernel - 不同次数
// ============================================================================

// 16次链式乘法
__global__ void kernel_dynamic_chain16(uint32_t *a, uint32_t *b, uint32_t *out, int n, uint32_t Q) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n) {
        uint32_t x = a[idx], y = b[idx];
        x = mod_mul_dynamic(x, y, Q); x = mod_mul_dynamic(x, y, Q);
        x = mod_mul_dynamic(x, y, Q); x = mod_mul_dynamic(x, y, Q);
        x = mod_mul_dynamic(x, y, Q); x = mod_mul_dynamic(x, y, Q);
        x = mod_mul_dynamic(x, y, Q); x = mod_mul_dynamic(x, y, Q);
        x = mod_mul_dynamic(x, y, Q); x = mod_mul_dynamic(x, y, Q);
        x = mod_mul_dynamic(x, y, Q); x = mod_mul_dynamic(x, y, Q);
        x = mod_mul_dynamic(x, y, Q); x = mod_mul_dynamic(x, y, Q);
        x = mod_mul_dynamic(x, y, Q); x = mod_mul_dynamic(x, y, Q);
        out[idx] = x;
    }
}

template<uint32_t Q>
__global__ void kernel_barrett_chain16(uint32_t *a, uint32_t *b, uint32_t *out, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n) {
        uint32_t x = a[idx], y = b[idx];
        x = mod_mul_barrett<Q>(x, y); x = mod_mul_barrett<Q>(x, y);
        x = mod_mul_barrett<Q>(x, y); x = mod_mul_barrett<Q>(x, y);
        x = mod_mul_barrett<Q>(x, y); x = mod_mul_barrett<Q>(x, y);
        x = mod_mul_barrett<Q>(x, y); x = mod_mul_barrett<Q>(x, y);
        x = mod_mul_barrett<Q>(x, y); x = mod_mul_barrett<Q>(x, y);
        x = mod_mul_barrett<Q>(x, y); x = mod_mul_barrett<Q>(x, y);
        x = mod_mul_barrett<Q>(x, y); x = mod_mul_barrett<Q>(x, y);
        x = mod_mul_barrett<Q>(x, y); x = mod_mul_barrett<Q>(x, y);
        out[idx] = x;
    }
}

template<uint32_t Q, uint32_t Q_INV>
__global__ void kernel_mont_chain16(uint32_t *a, uint32_t *b, uint32_t *out, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n) {
        uint32_t x = a[idx], y = b[idx];
        x = mont_mul<Q, Q_INV>(x, y); x = mont_mul<Q, Q_INV>(x, y);
        x = mont_mul<Q, Q_INV>(x, y); x = mont_mul<Q, Q_INV>(x, y);
        x = mont_mul<Q, Q_INV>(x, y); x = mont_mul<Q, Q_INV>(x, y);
        x = mont_mul<Q, Q_INV>(x, y); x = mont_mul<Q, Q_INV>(x, y);
        x = mont_mul<Q, Q_INV>(x, y); x = mont_mul<Q, Q_INV>(x, y);
        x = mont_mul<Q, Q_INV>(x, y); x = mont_mul<Q, Q_INV>(x, y);
        x = mont_mul<Q, Q_INV>(x, y); x = mont_mul<Q, Q_INV>(x, y);
        x = mont_mul<Q, Q_INV>(x, y); x = mont_mul<Q, Q_INV>(x, y);
        out[idx] = x;
    }
}

template<uint32_t Q>
__global__ void kernel_shoup_chain16(uint32_t *a, uint32_t *b, uint32_t *b_shoup, uint32_t *out, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n) {
        uint32_t x = a[idx], y = b[idx], y_s = b_shoup[idx];
        x = shoup_mul<Q>(x, y, y_s); x = shoup_mul<Q>(x, y, y_s);
        x = shoup_mul<Q>(x, y, y_s); x = shoup_mul<Q>(x, y, y_s);
        x = shoup_mul<Q>(x, y, y_s); x = shoup_mul<Q>(x, y, y_s);
        x = shoup_mul<Q>(x, y, y_s); x = shoup_mul<Q>(x, y, y_s);
        x = shoup_mul<Q>(x, y, y_s); x = shoup_mul<Q>(x, y, y_s);
        x = shoup_mul<Q>(x, y, y_s); x = shoup_mul<Q>(x, y, y_s);
        x = shoup_mul<Q>(x, y, y_s); x = shoup_mul<Q>(x, y, y_s);
        x = shoup_mul<Q>(x, y, y_s); x = shoup_mul<Q>(x, y, y_s);
        out[idx] = x;
    }
}

// ============================================================================
// 性能测试函数
// ============================================================================
template<typename KernelFunc>
float benchmark(KernelFunc kernel, uint32_t *d_a, uint32_t *d_b, uint32_t *d_out, 
                int n, int blocks, int threads, int iterations) {
    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    
    // 预热
    kernel<<<blocks, threads>>>(d_a, d_b, d_out, n);
    cudaDeviceSynchronize();
    
    // 计时
    cudaEventRecord(start);
    for (int i = 0; i < iterations; i++) {
        kernel<<<blocks, threads>>>(d_a, d_b, d_out, n);
    }
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);
    
    float ms = 0;
    cudaEventElapsedTime(&ms, start, stop);
    cudaEventDestroy(start);
    cudaEventDestroy(stop);
    
    return ms / iterations;
}

// 专门用于动态模数 Kernel 的测试函数
template<typename KernelFunc>
float benchmark_dynamic(KernelFunc kernel, uint32_t *d_a, uint32_t *d_b, uint32_t *d_out, 
                        int n, uint32_t Q, int blocks, int threads, int iterations) {
    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    
    // 预热
    kernel<<<blocks, threads>>>(d_a, d_b, d_out, n, Q);
    cudaDeviceSynchronize();
    
    // 计时
    cudaEventRecord(start);
    for (int i = 0; i < iterations; i++) {
        kernel<<<blocks, threads>>>(d_a, d_b, d_out, n, Q);
    }
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);
    
    float ms = 0;
    cudaEventElapsedTime(&ms, start, stop);
    cudaEventDestroy(start);
    cudaEventDestroy(stop);
    
    return ms / iterations;
}

// 专门用于 Shoup Kernel 的测试函数（需要额外的 b_shoup 参数）
template<typename KernelFunc>
float benchmark_shoup(KernelFunc kernel, uint32_t *d_a, uint32_t *d_b, uint32_t *d_b_shoup,
                      uint32_t *d_out, int n, int blocks, int threads, int iterations) {
    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    
    // 预热
    kernel<<<blocks, threads>>>(d_a, d_b, d_b_shoup, d_out, n);
    cudaDeviceSynchronize();
    
    // 计时
    cudaEventRecord(start);
    for (int i = 0; i < iterations; i++) {
        kernel<<<blocks, threads>>>(d_a, d_b, d_b_shoup, d_out, n);
    }
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);
    
    float ms = 0;
    cudaEventElapsedTime(&ms, start, stop);
    cudaEventDestroy(start);
    cudaEventDestroy(stop);
    
    return ms / iterations;
}


int main() {
    const int N = 1024 * 1024 * 16;  // 16M 元素
    const int THREADS = 256;
    const int BLOCKS = (N + THREADS - 1) / THREADS;
    const int ITERATIONS = 100;
    
    printf("============================================================\n");
    printf("模乘法性能深度测试 (含 Shoup 算法)\n");
    printf("============================================================\n");
    printf("元素数量: %d (%.1f M)\n", N, N / 1e6);
    printf("迭代次数: %d\n\n", ITERATIONS);
    
    // 分配内存
    uint32_t *d_a, *d_b, *d_b_shoup, *d_out;
    cudaMalloc(&d_a, N * sizeof(uint32_t));
    cudaMalloc(&d_b, N * sizeof(uint32_t));
    cudaMalloc(&d_b_shoup, N * sizeof(uint32_t));
    cudaMalloc(&d_out, N * sizeof(uint32_t));
    cudaMemset(d_a, 1, N * sizeof(uint32_t));
    cudaMemset(d_b, 2, N * sizeof(uint32_t));
    
    // 预计算 Shoup 参数（在 CPU 上）
    uint32_t *h_b = new uint32_t[N];
    uint32_t *h_b_shoup = new uint32_t[N];
    for (int i = 0; i < N; i++) {
        h_b[i] = 2;  // 与 d_b 相同
        h_b_shoup[i] = compute_shoup_param(h_b[i], Q_28BIT);
    }
    cudaMemcpy(d_b_shoup, h_b_shoup, N * sizeof(uint32_t), cudaMemcpyHostToDevice);
    delete[] h_b;
    delete[] h_b_shoup;
    
    // ==================== 测试1: 不同链式乘法次数 ====================
    printf("==================== 测试1: 链式乘法次数影响 ====================\n");
    printf("模数: Q = 268369921 (28位)\n\n");
    
    printf("| 次数 | Dynamic (ms) | Barrett (ms) | Shoup (ms) | Montgomery (ms) | Speedup (D/B/S/M) |\n");
    printf("|------|--------------|--------------|------------|-----------------|-------------------|\n");
    
    // 16次
    float t_d16 = benchmark_dynamic(kernel_dynamic_chain16, d_a, d_b, d_out, N, Q_28BIT, BLOCKS, THREADS, ITERATIONS);
    float t_b16 = benchmark(kernel_barrett_chain16<Q_28BIT>, d_a, d_b, d_out, N, BLOCKS, THREADS, ITERATIONS);
    float t_s16 = benchmark_shoup(kernel_shoup_chain16<Q_28BIT>, d_a, d_b, d_b_shoup, d_out, N, BLOCKS, THREADS, ITERATIONS);
    float t_m16 = benchmark(kernel_mont_chain16<Q_28BIT, Q_INV_28BIT>, d_a, d_b, d_out, N, BLOCKS, THREADS, ITERATIONS);
    printf("| 16   | %.4f       | %.4f       | %.4f     | %.4f          | 1.00 / %.2f / %.2f / %.2f |\n", 
           t_d16, t_b16, t_s16, t_m16, t_d16/t_b16, t_d16/t_s16, t_d16/t_m16);
    
    // ==================== 测试2: 不同模数位宽 ====================
    printf("\n==================== 测试2: 模数位宽影响 (16次链式) ====================\n\n");
    
    printf("| 模数 | 位宽 | Barrett (ms) | Shoup (ms) | Montgomery (ms) | 比值 (B/S/M) |\n");
    printf("|------|------|--------------|------------|-----------------|---------------|\n");
    
    // 预计算不同模数的 Shoup 参数
    uint32_t *h_b_tmp = new uint32_t[N];
    uint32_t *h_b_shoup_tmp = new uint32_t[N];
    
    // 16位模数
    for (int i = 0; i < N; i++) h_b_shoup_tmp[i] = compute_shoup_param(2, Q_16BIT);
    cudaMemcpy(d_b_shoup, h_b_shoup_tmp, N * sizeof(uint32_t), cudaMemcpyHostToDevice);
    float t_b16bit = benchmark(kernel_barrett_chain16<Q_16BIT>, d_a, d_b, d_out, N, BLOCKS, THREADS, ITERATIONS);
    float t_s16bit = benchmark_shoup(kernel_shoup_chain16<Q_16BIT>, d_a, d_b, d_b_shoup, d_out, N, BLOCKS, THREADS, ITERATIONS);
    float t_m16bit = benchmark(kernel_mont_chain16<Q_16BIT, Q_INV_16BIT>, d_a, d_b, d_out, N, BLOCKS, THREADS, ITERATIONS);
    printf("| %u | 16位 | %.4f       | %.4f     | %.4f          | %.2f / %.2f / 1.00 |\n", 
           Q_16BIT, t_b16bit, t_s16bit, t_m16bit, t_b16bit/t_m16bit, t_s16bit/t_m16bit);
    
    // 20位模数
    for (int i = 0; i < N; i++) h_b_shoup_tmp[i] = compute_shoup_param(2, Q_20BIT);
    cudaMemcpy(d_b_shoup, h_b_shoup_tmp, N * sizeof(uint32_t), cudaMemcpyHostToDevice);
    float t_b20bit = benchmark(kernel_barrett_chain16<Q_20BIT>, d_a, d_b, d_out, N, BLOCKS, THREADS, ITERATIONS);
    float t_s20bit = benchmark_shoup(kernel_shoup_chain16<Q_20BIT>, d_a, d_b, d_b_shoup, d_out, N, BLOCKS, THREADS, ITERATIONS);
    float t_m20bit = benchmark(kernel_mont_chain16<Q_20BIT, Q_INV_20BIT>, d_a, d_b, d_out, N, BLOCKS, THREADS, ITERATIONS);
    printf("| %u | 20位 | %.4f       | %.4f     | %.4f          | %.2f / %.2f / 1.00 |\n", 
           Q_20BIT, t_b20bit, t_s20bit, t_m20bit, t_b20bit/t_m20bit, t_s20bit/t_m20bit);
    
    // 24位模数
    for (int i = 0; i < N; i++) h_b_shoup_tmp[i] = compute_shoup_param(2, Q_24BIT);
    cudaMemcpy(d_b_shoup, h_b_shoup_tmp, N * sizeof(uint32_t), cudaMemcpyHostToDevice);
    float t_b24bit = benchmark(kernel_barrett_chain16<Q_24BIT>, d_a, d_b, d_out, N, BLOCKS, THREADS, ITERATIONS);
    float t_s24bit = benchmark_shoup(kernel_shoup_chain16<Q_24BIT>, d_a, d_b, d_b_shoup, d_out, N, BLOCKS, THREADS, ITERATIONS);
    float t_m24bit = benchmark(kernel_mont_chain16<Q_24BIT, Q_INV_24BIT>, d_a, d_b, d_out, N, BLOCKS, THREADS, ITERATIONS);
    printf("| %u | 24位 | %.4f       | %.4f     | %.4f          | %.2f / %.2f / 1.00 |\n", 
           Q_24BIT, t_b24bit, t_s24bit, t_m24bit, t_b24bit/t_m24bit, t_s24bit/t_m24bit);
    
    // 28位模数
    printf("| %u | 28位 | %.4f       | %.4f     | %.4f          | %.2f / %.2f / 1.00 |\n", 
           Q_28BIT, t_b16, t_s16, t_m16, t_b16/t_m16, t_s16/t_m16);
    
    delete[] h_b_tmp;
    delete[] h_b_shoup_tmp;
    
    // ==================== 总结 ====================
    printf("\n==================== 总结 ====================\n\n");
    printf("1. 链式乘法次数影响 (Baseline = Dynamic Modulo):\n");
    printf("   - 16次: \n");
    printf("     * Dynamic    : 1.00x (%.4f ms)\n", t_d16);
    printf("     * Barrett    : %.2fx (%.4f ms)\n", t_d16/t_b16, t_b16);
    printf("     * Shoup      : %.2fx (%.4f ms)\n", t_d16/t_s16, t_s16);
    printf("     * Montgomery : %.2fx (%.4f ms)\n", t_d16/t_m16, t_m16);
    printf("\n");
    printf("2. 模数位宽影响 (16次链式，相对 Montgomery):\n");
    printf("   - 16位: Barrett/Shoup/Montgomery = %.2fx / %.2fx / 1.00x\n", t_b16bit/t_m16bit, t_s16bit/t_m16bit);
    printf("   - 20位: Barrett/Shoup/Montgomery = %.2fx / %.2fx / 1.00x\n", t_b20bit/t_m20bit, t_s20bit/t_m20bit);
    printf("   - 24位: Barrett/Shoup/Montgomery = %.2fx / %.2fx / 1.00x\n", t_b24bit/t_m24bit, t_s24bit/t_m24bit);
    printf("   - 28位: Barrett/Shoup/Montgomery = %.2fx / %.2fx / 1.00x\n", t_b16/t_m16, t_s16/t_m16);
    printf("\n");
    printf("3. Shoup vs Barrett (28位模数):\n");
    printf("   - Shoup 相对 Barrett 加速: %.2fx\n", t_b16/t_s16);
    printf("   - 但仍比 Montgomery 慢: %.2fx\n", t_s16/t_m16);
    printf("   - 内存开销: Shoup 需要 2x 存储 (原值 + Shoup 参数)\n");
    
    // 清理
    cudaFree(d_a);
    cudaFree(d_b);
    cudaFree(d_b_shoup);
    cudaFree(d_out);
    
    printf("\n============================================================\n");
    
    return 0;
}