pmarreck/blip_mp
GitHub: pmarreck/blip_mp
纯 Zig 实现的任意精度整数库，采用 BLIP 自描述可变长度编码替代传统基数组存储，在 Apple Silicon 上常见密码学运算大小中击败 GMP，同时规避了 LGPL 许可限制。
Stars: 0 | Forks: 0
# blip_mp
**一个纯 Zig 的多精度整数库，将值存储为 BLIP 编码的字节（可变长度、自描述），而不是 GMP 的固定宽度基数组表示。** 在 Apple Silicon 上的最常见大数运算中击败了 GMP。纯 Zig 实现，无内联汇编，无 LGPL 链接限制。
## 为什么会存在这个项目？
GMP 的 `mpz_t` 是严肃数值计算软件中事实上的大数表示。它的弱点在于：**每个值的开销是固定和结构性的，无论实际整数有多大。**
```
mpz_t = { int _mp_alloc, int _mp_size, mp_limb_t *_mp_d }
≈ 16 bytes of struct
+ heap allocation for _mp_d
+ alignment padding
+ allocator bookkeeping
Storing the value 5 takes ~24 bytes spread across two cache lines,
mediated by a malloc round-trip.
```
对于由中小型数字主导的工作负载（累加器循环、EC 标量、多项式系数、哈希派生整数），GMP 的结构性开销会导致内存占用比实际信息内容高出 5–25 倍，并且在每次创建或销毁值时都需要经历一次分配器交互。
**BLIP**（[规范](https://github.com/pmarreck/BLIP)）是一种自描述的可变长度整数编码：
| 值 | BLIP 字节 | GMP 存储 |
|---|---|---|
| `5` | 1 字节 (立即数) | 24+ 字节 (结构体 + 堆) |
| `2^32` | 5 字节 (头部 + 4 字节小端序) | 24+ 字节 |
| `2^256` | 33 字节 | 16 字节结构体 + 32 字节堆 |
| `2^4096` | 513 字节连续 | 16 + 512 + 额外开销 |
`blip_mp` 直接在 BLIP 存储之上构建了一个任意精度整数库。没有单独的长度字段，小数值没有堆指针间接寻址，累加器循环没有分配器交互。**字节即数值。**
## 我们测量了什么？
Apple Silicon (M 系列)，aarch64-darwin，Zig 0.16.0 ReleaseFast，libc malloc。
### 核心胜利
- **所有适合 `i64` 的值：比 GMP 快 1.95–2.66 倍**
- **所有三种标准 RSA 乘法大小（1024、2048、3072 位）以 1.19–1.45 倍的优势击败 GMP** —— 在 iter-33 Karatsuba 叶子清理后增加了 4096、6144 和 8192 位的新胜利，目前总共有 12 种乘法大小击败了 GMP。
- **RSA-2048 大满贯**：blip_mp 在全球部署最广泛的加密操作数大小上，于 `mul` (1.16 倍) + `divMod` (1.31 倍) + `powm` (1.11 倍) 中击败了 GMP。
- **768 位及以上的加法：比 GMP 快 1.04–1.28 倍**（在 M5/M9 内务处理清理后，现有 8 种大小击败了 GMP；之前仅在 4096+ 位时达到 1.03 倍）
- **256 位及以上的减法：比 GMP 快 1.03–2.11 倍** —— 在分块 `subPayloads` 修复镜像了长期存在的 addPayloads 优化后，10 种 Mp.sub 大小击败了 GMP。核心亮点：6144 位减法从 606 ns 降至 48 ns (12.7 倍)。
- **大数乘法（16384+ 位，通过 Toom-3）：快 1.03 倍**（幅度不大）
- **模幂运算 (RSA-2048)：比 GMP 快 13%**（带有 Montgomery 的 Mp.powm，M7-4.3 + Möller-Granlund 改进的内部除法）。在 1024 位时我们以 3% 的优势胜出；在 3072 位时为 8%。这是任何严肃加密工作负载（RSA 加密/解密/签名、DH 密钥交换、ECC 标量乘法）的核心关注点。
- **长除法（2048 位 / 1024 位）：比 GMP 快 24%**（带有基于 u64 的 Knuth 算法 D 的 Mp.divMod）—— 比最初落后 28.8 倍的基于字节的实现快了 36 倍。
- **正确性：12029/12029 随机 GMP 交叉验证测试通过**，涵盖 add、sub、mul、div、mod、divMod、powm、invMod —— 完整的模块化算术 API。
### 坦诚的落败
我们在以下方面比 GMP 慢：
- **128 位加法和减法**（约为 GMP 的 0.56–0.65 倍）。在最小尺寸上剩余的差距来自于 blip_mp 较重结构体的 load+store ABI 开销（一个缓存行 + 8B，而 GMP 是 24 字节的 `mpz_t`）—— 这是基础性的，而非算法层面的。在优化落地后，其他较小的大小（192-512 位）现在已达到 GMP 的 70-90%。
- **128–256 位乘法**（0.40–0.67 倍）。相同的每次操作开销。
- **16384+ 位乘法**（iter-33 后为 0.55–0.87 倍）。8192 位在 iter-33 前落后 GMP 0.72 倍；Karatsuba 叶子分块将其反转为 1.04 倍的胜利。16384 位从 0.68 倍缩小至 0.87 倍。32768 位仍处于 FFT 领域（0.55 倍）。GMP 在约 16K 位以上使用 Schönhage-Strassen FFT 乘法。我们有一个完整的纯 Zig FFT 栈（单素数 NTT + 双素数 CRT + NEON-SIMD 蝶形运算）已经发布并通过了正确性验证，但在生产环境中被禁用 —— 即使矢量化带来了 1.40 倍的加速（32K 位 FFT 路径：191K -> 135K ns），Toom-3 在 iter-33 后以 100K ns 的成绩依然获胜。PLAN.md 中的 M6-4-E 阶梯目标是实现翻转所需的分配消除 + Stockham + 内联汇编杠杆。
### 惊喜：GMP 手写优化的 aarch64 汇编在 Apple Silicon 上的优势约为 0%
我们构建了一个带有 `--disable-assembly` 的第二个 GMP 变体，并运行了相同的基准测试。在我们的所有测试大小中，汇编的优势基本为零。在 4096 位和 32768 位的加法中，**C 参考代码实际上比汇编更快。** 现代 clang `-O3` 从 C `__builtin_add_overflow` 生成的 ADCS 链近乎最优，手写汇编无法超越 —— 并且汇编成为了一个不透明的调用边界，破坏了内联。
**含义：** 与 GMP 的每一项差距纯粹是算法层面的。我们不需要内联汇编。我们需要 FFT 乘法和更紧凑的内务处理。这两者都可以通过纯 Zig 实现。
## 快速开始
需要 [Nix](https://nixos.org/)（处理 Zig 0.16.0 工具链和用于基准测试的 GMP 构建）。
```
git clone https://github.com/pmarreck/blip_mp
cd blip_mp
./build # native ReleaseFast build via nix
./test # 76+ unit tests + 8240-check GMP cross-validation
./result/bin/blip_mp_bench # run the bench (after nix build .#packages..bench)
```
在你自己的 Zig 代码中：
```
const std = @import("std");
const blip_mp = @import("blip_mp");
const Mp = blip_mp.Mp;
pub fn main() !void {
const allocator = std.heap.c_allocator;
var a = Mp.init(allocator);
defer a.deinit();
var b = Mp.init(allocator);
defer b.deinit();
var r = Mp.init(allocator);
defer r.deinit();
try a.setI64(123_456_789);
try b.setI64(987_654_321);
try r.mul(&a, &b);
// r.bytes() is the canonical signed-BLIP encoding — also the wire form.
// No mpz_export round-trip needed.
std.debug.print("Product encoded as {d} bytes\n", .{r.bytes().len});
const value = try r.getI64(); // works because product fits in i64
std.debug.print("Value: {d}\n", .{value});
}
```
## 一段话架构概述
`Mp` 是一个 72 字节的结构体（一个缓存行 + 8B）。内联模式将最多 24 字节的编码值存储在 `inline_buf` 中；堆模式使用带有 `heap_offset` 字段的 `heap_buf`，该字段允许写入结果而无需移位（头部直接放置在有效载荷之前）。对于内联的长度前缀值，一个内部不变量将 `inline_buf[1..9]` 维护为完整的符号扩展 i64 —— 算术热路径通过单个 `LDR` u64 加载来读取它。缓存的 `(payload_offset, sign, payload_len)` 字段跳过了每次操作的头部解析。层级分派是自动的：适合 i64 的值使用本机算术；较大的值通过 BLIP 有效载荷字节使用直接字节的两补码加/减法或 Karatsuba/Toom-3 乘法。**不存在任何辅助的“基数组”数据结构** —— 我们通过 `readInt`/`writeInt` 直接从字节有效载荷中读取 u64/u128/u256/u512 块。字节即数值，贯穿始终。
完整细节见 [`CODE_MINIMAP.md`](CODE_MINIMAP.md)，基准测试历史见 [`BENCHMARK_RESULTS.md`](BENCHMARK_RESULTS.md)，综合结果撰写见 [`RESULTS.md`](RESULTS.md)。
## 权衡与局限性
**本库是什么：** 一个研究级的任意精度整数库，验证了 BLIP 存储范式，并在 Apple Silicon 上的常见大小中击败了 GMP。纯 Zig，无汇编，无 LGPL 限制。
**它（还）不是什么：**
- **FFT 乘法已发布并确保正确性，但被禁用** —— 完整的单素数 NTT + 双素数 CRT + NEON-SIMD 矢量化蝶形运算位于 `src/fft.zig` 中，在高达 256K 位的大小下，所有 8240/8240 交叉检查均与 GMP 位对位一致。但常量因子使得 Toom-3 在我们支持的范围内的每个操作数大小上都保持领先（M 系列特有发现：纯 NEON Montgomery 的集成速度慢于现有的标量内置矢量形式，因为它挤占了 NEON 管道并使 M4 的双标量乘法管线处于饥饿状态）。要缩小剩下的 13–15% 差距需要分配消除 + 内联汇编，已在 M6-4-E 中计划。
- **模逆在 1024-2048 位上落后 GMP 约 4-6 倍**（在 M9 Lehmer 之前为 29-38 倍；在 M10 宽窗口 Lehmer 之后进一步降至 7-8 倍）。核心亮点：2048 位 invMod 现在落后 GMP 3.96 倍（在 M10 之前为 7.85 倍）。缩小剩余差距需要真正的递归半 GCD，已计划作为 M11。
- **仅验证了单一平台** —— 上述数字均为 aarch64-darwin (Apple M 系列)。x86_64 可能会改变这一局面，尤其是关于 asm 与 clang 的结果。
- **尚无 C FFI** —— 公共接口仅限 Zig。添加 `include/blip_mp.h` 是一个明确的扩展方向。
- **未针对非对齐操作数大小进行优化** —— `tier3Op` 适用于任何大小，但在有效载荷长度为 8 字节倍数时速度最快（大多数加密大小正是如此）。
**它不打算成为什么：**
- 不是 GMP 的完整替代品。没有 `mpf_t`（浮点数），没有 `mpq_t`（有理数），没有 `mpfr`（扩展精度浮点数）。
- 不追求超大数记录。GMP 针对每个架构的汇编调优是数十年的工作；在 FFT 落地之前，我们在 64K+ 位操作数上不具备竞争力。
- 未进行汇编调优。M5-5 对照实验表明汇编在 M 系列上优势约为 0%。如果情况不同，我们可能需要在 x86_64 上重新审视。
## 路线图
**按优先级顺序：**
1. **完成 FFT 与 Toom-3 的翻转**（PLAN.md 中的 M6-4-E）。FFT 原语、CRT 扩展和 NEON-SIMD 蝶形运算均已发布并验证了正确性；将 Toom-3 的差距从 1.93 倍缩小至 1.15 倍。剩下的 13–15% 需要调用方提供的暂存空间（消除了 4 次每次调用的分配，约 6–9K ns），将 Stockham 接入生产环境，并可能需要针对蝶形运算内循环进行手动调度的 aarch64 内联汇编。
2. ~~**更紧凑的 `tier3Op` 内务处理**~~ —— 已于 2026-05-02 完成。新增五项加法胜利（768/1536/2048/3072/4096 位）；128 位差距缩小了约 40%；1024 位达到持平。实现：针对固定有效载荷大小的快速路径特化，编译为单个 uN +% ADC 链，以及一条适合内联的栈路径。
3. **在 x86_64 Linux + Windows 上进行跨平台验证。** 两项 M 系列特有的发现需要在 x86_64 上进行验证： GMP 汇编在 M 系列上优势约为 0%，AVX-512 可能会改变这一点"（M5-5）； 纯 NEON Montgomery 输给了标量内置矢量形式，因为 M4 的双标量乘法管线"（M6-4-A.6） —— 不同的调度器可能会此结果。
4. **用于 `Mp.invMod` 的真正递归半 GCD (M11)** —— 缩小与 GMP 之间剩余的 4-6 倍差距。M10（宽窗口 Lehmer）刚刚落地，在 2048 位时达到 GMP 的 3.96 倍；带有多个精度矩阵条目的真正递归 HGCD 是亚二次方的，将缩小剩下的大部分差距。
5. **C FFI 头文件**（`include/blip_mp.h`），供下游消费者使用。
6. **Toom-Cook 4-way**，用于 4K-16K 位乘法。已被降级 —— 其区区 15-30% 的收益不值得付出实现成本，因为 FFT 依然是重头戏。M6-2.1 + M6-2.2 辅助函数（`divExactBy5`、`mulSmallSignedConst`）作为未来的构建模块已位于 tier3.zig 中。
7. **在 `./bm` 中集成 `hyperfine`**，以进行适当的统计基准聚合。当前数字为 3 次运行的手动中位数。
## 许可意图
MIT 或 Apache-2.0（在标记发布前待定）。两者均允许与 libgmp (LGPLv3) 链接及下游商业使用。
## 致谢
- Peter Marreck 的 [BLIP 编码规范](https://github.com/pmarreck/BLIP)。
- GMP 团队提供了我们用于交叉验证的参考实现。
- Marco Bodrato 和 Alberto Zanoni 提供了 Toom-3 插值公式。
*有关最初的设计假设，请参阅 [SPEC.md](SPEC.md)；有关全面的技术说明，请参阅 [RESULTS.md](RESULTS.md)；有关每次运行的历史记录，请参阅 [BENCHMARK_RESULTS.md](BENCHMARK_RESULTS.md)；有关里程碑检查清单，请参阅 [PLAN.md](PLAN.md)；有关按文件索引，请参阅 [CODE_MINIMAP.md](CODE_MINIMAP.md)。*
标签：Apple Silicon, bignum, BLIP编码, GMP替代, M系列芯片, RSA, Zig语言, 任意精度, 内存优化, 加密算法, 变长编码, 基础组件, 多精度算法, 大数运算, 密码学, 底层开发, 性能优化, 手动系统调用, 数值计算, 无内联汇编, 检测绕过, 算法库, 纯Zig, 自描述编码, 高性能计算, 高精度整数