Sdcb.Arithmetic 2.0：纯 C# 实现四倍精度 Float128，比 MPFR 快 2 到 12 倍

洛伦兹系统精度阶梯：float、double、DoubleDouble、Float128、QuadDouble、BigFloat 依次叛逃

上面是 sdcb.ai/arithmetic/chaos 的精度阶梯：同一个洛伦兹系统、同一个起点、同一个 RK4 积分器，每条轨迹只有舍入的位数不同。24 位的 float 在 t = 21.3 叛逃，53 位的 double 在 t = 42.7，2.0 新增的三个类型分别撑到了 83.0（DoubleDouble）、86.1（Float128）和 164.6（QuadDouble）。

这几条新轨迹，就是今天发布的 Sdcb.Arithmetic（https://github.com/sdcb/Arithmetic） 2.0 的主角：三个零分配的定长浮点类型。其中 Float128 是 IEEE 754 四倍精度，每个基本运算和超越函数都正确舍入，在同样 113 位精度下，基本运算比原生 MPFR 快 2 到 4 倍，exp、log、sin、atan、pow 快 6 到 12 倍。依然是纯 C#，没有原生依赖，能 AOT，能进浏览器。

dotnet add package Sdcb.Arithmetic

距离上次发布

上一篇（https://www.cnblogs.com/sdcb/p/23446780/20261027-arithmetic-pure-dotnet） 写的时候是 1.1（Chudnovsky 版的 Real.Pi）。之后又发了这些版本：

1.2.x 的三个问题都是 inkydragon（https://github.com/inkydragon） 报告的，在这里表示感谢。具体来说：

• TestBit
  、SetBit、ShiftLeft、ModPow2 这些方法接收 long 的位索引，但内部没检查就转成了 int，大索引会悄悄回绕，算出错误结果。现在按 long 比较，放不下时抛 OverflowException；
• Real
   的指数计算改成 checked，超出 long 范围时抛异常，不再回绕。RootN(x, long k) 里超过 uint.MaxValue 的次数以前会被截断，现在按真实的奇偶性处理特殊值；
• FromBigInt(BigInt, int precision)
   这类入口以前不校验精度，传个非法值会算出损坏的结果。现在所有带 precision 参数的公开方法都统一校验，并且有一个反射测试逐个检查。

这几个都是"正常用法碰不到，碰到了就是静默算错"的问题。对数学库来说，静默算错比抛异常糟糕得多。

为什么要定长浮点

1.x 里只有一个浮点类型 Real（现在叫 BigFloat），任意精度，正确舍入。它很强，但有两个代价：

• 每次运算都要分配，基本运算每次 40 到 56 字节；
• 为了支持任意精度，哪怕你只要 113 位，走的也是通用的多 limb 代码。

实际中很多场景只需要"比 double 多一点"：几何计算里的鲁棒判定、长时间积分、补偿求和、检查 double 结果对不对。这时候用任意精度就太重了。于是 2.0 加了三个 readonly struct：

sdcb.ai 上的 2.0 新增类型

怎么选：结果要逐位可复现、或者要和其他 IEEE binary128 实现交换数据，用 Float128；只要快，用 DoubleDouble 或 QuadDouble。

Float128：每一位都是对的

Float128 的核心算术全部是整数代码：拆包，带 guard 和 sticky 位计算，最后走一个共用的舍入打包。结构参照的是 Berkeley SoftFloat 的 f128。加减乘除、开方、FMA、取余、转换和各种取整都正确舍入，支持次正规数、±0、±∞ 和 NaN。

using Sdcb.Arithmetic;

Float128 a = Float128.Parse("0.1"), b = Float128.Parse("0.2");
Console.WriteLine(a + b);                 // 0.30000000000000000000000000000000004
Console.WriteLine(Float128.Sqrt(2));      // 1.414213562373095048801688724209698
Console.WriteLine(Float128.Exp(1));       // 2.7182818284590452353602874713526623

// 定向舍入：一次得到下界和上界
Console.WriteLine(Float128.Divide(1, 3, RoundingMode.TowardNegative)); // 0.3333333333333333333333333333333333
Console.WriteLine(Float128.Divide(1, 3, RoundingMode.TowardPositive)); // 0.33333333333333333333333333333333337

第一行说明它也是二进制浮点，0.1 + 0.2 照样不等于 0.3，只不过错在第 35 位有效数字上。ToString() 输出的是能往返的最短字符串，Parse 从精确的十进制值正确舍入。

超越函数是这次花功夫最多的地方。34 个函数全部正确舍入，每个都有带 RoundingMode 的重载。如果每次都走 BigFloat 算再舍入，结果是对的，但很慢。所以 exp、log、pow、三角、反三角和双曲这 16 个函数，先用 triple-double（约 159 位）的内核算一遍，这个内核有证明过的误差界：只要误差区间不跨过舍入边界，就直接舍入返回，这就是 Ziv 的策略。跨过了（精确结果、中点、次正规或接近溢出、2^30 以上的三角参数），再交给 BigFloat 算。随机输入里，三角函数大约只有 0.02% 到 0.04% 会走到慢路径。

验证方面，用 BigFloat 做 round-to-odd 参照（避免两次舍入），在约 18 万组随机的加减乘除、开方、FMA 用例上逐位比对，覆盖全部 5 种舍入模式，并且专门偏向次正规数、溢出、指数相差很大、大量抵消和恰好在中点这些难例。

DoubleDouble 和 QuadDouble：快

这两个类型用的是 Hida、Li 和 Bailey 的 QD 库里的算法：精确的（"IEEE 风格"）加法，基于 FMA 的乘法，以及除法和开方。exp、log、三角和双曲函数都是原生实现（参数约化 + 查表 + 级数），只有需要极精确参数约化的情况（巨大的三角参数、%、Round(x, digits)）才回退到 BigFloat。

它们不保证正确舍入，所以我用 BigFloat（300 / 500 位）在固定种子的随机输入上测了最坏相对误差：

Pow 是按 exp(y·log x) 算的，所以 |y·log x| 越大误差越大。FusedMultiplyAdd 在两个类型上都保留未舍入的乘积，算 a·b − fl(a·b) 这种残差也是准的；要快的话用 MultiplyAddEstimate。

三个类型都支持泛型数学，同一份代码可以直接换类型跑：

static T Harmonic<T>(int n) where T : IFloatingPoint<T>
{
    T s = T.Zero;
    for (int i = n; i >= 1; i--) s += T.One / T.CreateChecked(i);
    return s;
}

Harmonic<double>(1000);        // 7.485470860550341
Harmonic<DoubleDouble>(1000);  // 7.4854708605503449126565182043338
Harmonic<Float128>(1000);      // 7.485470860550344912656518204333897
Harmonic<QuadDouble>(1000);    // 7.48547086055034491265651820433390017652167916970880366577362675013166

它们之间，以及和 BigFloat、BigInt、BigRational 之间都能互转。DoubleDouble 转 QuadDouble、三者转 BigFloat 都是隐式的、精确的；往低精度转是显式的，转到 Float128 时正确舍入：

QuadDouble qd = QuadDouble.Pi;
BigFloat exact = qd;            // 隐式，精确
Float128 f = (Float128)qd;      // 显式，正确舍入
QuadDouble q2 = DoubleDouble.Pi; // 隐式，精确

性能

测试环境是 AMD Ryzen 7 5800X、.NET 10，BenchmarkDotNet ShortRun，单位纳秒，越小越快。BigFloat 和 MPFR 都用和对应类型相同的精度。三个定长类型都不分配内存（包括 Float128 超越函数的快路径）。

Float128 对比 113 位：

DoubleDouble 对比 106 位，QuadDouble 对比 212 位：

总结一下：

• Float128
   基本运算比 MPFR 快 2 到 4 倍，超越函数快 6 到 12 倍，而且结果和 MPFR 一样是正确舍入的；
• DoubleDouble
   基本运算快 3 到 16 倍，exp、log、sin 快 9 到 12 倍；
• QuadDouble
   的加和乘和 MPFR 持平，除和开方快 2 到 3 倍，exp、log、sin 快 1.4 到 3.5 倍。也有比 MPFR 慢的：大量抵消的减法要做精确合并，29 ns 对 19 ns，慢 1.5 倍。

这里能快，主要因为定长：位数固定，可以把循环全部展开，用寄存器里的 UInt128 或几个 double 直接算，不用分配，也不用处理任意长度的 limb 数组。MPFR 是通用的任意精度库，113 位对它只是一个普通的精度参数，这样比对它并不完全公平。完整数据见 performance.md（https://github.com/sdcb/Arithmetic/blob/main/docs/performance.md#fixed-width-floats）。

Real 改名 BigFloat

2.0 是大版本号，因为有一个不兼容的改名：

改名的原因：现在库里有了四种浮点，BigFloat 和 BigInt、BigRational 放在一起正好是"任意精度"一族，名字一看就知道它和三个定长类型不是一类。舍入模式也不再是 Real 专属的，Float128 一样用它，所以叫 RoundingMode。

除了名字，成员、行为和二进制布局都没变。迁移有两种办法：全词替换（先把 RealRounding 换成 RoundingMode，再把 Real 换成 BigFloat），或者在项目里任意一个文件加两行别名，旧代码原样编译：

global using Real = Sdcb.Arithmetic.BigFloat;
global using RealRounding = Sdcb.Arithmetic.RoundingMode;

详见 迁移指南（https://github.com/sdcb/Arithmetic/blob/main/docs/migration-2.0.md）。

sdcb.ai 也更新了

sdcb.ai/arithmetic（https://sdcb.ai/arithmetic） 上的 WebAssembly 已经换成 2.0.0，页面上的 Real 全都换成了 BigFloat，类型一节多了上面那一行"2.0 新增"。

实验室里，Rump 的例子多了三个定长类型。这个表达式的真值是 -0.827396059946821368…，要 122 位以上精度才算得出来。113 位的 Float128 和约 106 位的 DoubleDouble 照样算错，约 212 位的 QuadDouble 算对了：

Rump 的例子：DoubleDouble 和 Float128 算错，QuadDouble 算对

Muller 递推也一样。它的真值收敛到 5，但每一次舍入误差都会把序列推向错误的不动点 100。double 在第 14 步左右掉下去，DoubleDouble 和 Float128 撑到第 27 步左右，QuadDouble 撑到第 53 步左右，256 位和 1024 位的 BigFloat 在图里的 60 步内都没掉：

Muller 递推：精度越高，坚持得越久

混沌页（https://sdcb.ai/arithmetic/chaos） 的精度阶梯和三体模式，加上了 DoubleDouble（106 位）和 QuadDouble（212 位），113 位那条现在用的是真正的 Float128。所有轨道的 RK4 都是逐个运算对齐的，区别只在类型。默认只显示原来那几条，DoubleDouble 和 QuadDouble 收在"+2 种 2.0 新类型"里，点开才出现。另外两个小改动：点击右侧的行可以隐藏或找回某条轨道；轨道叛逃 2 秒后自动隐藏，免得画面上堆满已经错掉的线。

开头那张动图最后一帧的右侧面板，就是这次跑出来的叛逃时刻：

各精度的叛逃时刻：float 21.3，double 42.7，DoubleDouble 83.0，Float128 86.1，QuadDouble 164.6，BigFloat 256 位 192.4

DoubleDouble 和 Float128 只差 7 位精度，叛逃时间也只差 3；QuadDouble 的 212 位，和 256 位的 BigFloat 也离得不远。洛伦兹系统的李雅普诺夫指数约 0.9，误差每过一个单位时间大约放大 e^0.9 倍，所以每多 1 位精度，大约多 0.76 个时间单位的可预测未来。

最后

• GitHub：sdcb/Arithmetic（https://github.com/sdcb/Arithmetic）
• NuGet：dotnet add package Sdcb.Arithmetic
• 定长浮点文档：floats.md（https://github.com/sdcb/Arithmetic/blob/main/docs/floats.md）
• 在线 Show Case：sdcb.ai/arithmetic（https://sdcb.ai/arithmetic）