超越.NET极限：我打造的高精度数值计算库
还记得那一天，我大学刚毕业，紧张又兴奋地走进人生第一场.NET工作面试。我还清楚地记得那个房间的气氛，空调呼呼地吹着，面试官的表情严肃而深沉。我们进行了一番交谈：
- 面试官，眼镜后的目光犀利：“你还有其它问题吗？”
- 我，有点颤抖但决心坚定：“可惜
C#只有大整数BigInteger，没有大小数。如果我想计算特别大或者特别精确的数字，C#就无能为力。你看Java那边，有BigFloat，请问面试官，您知道C#这边处理大小数的需求有什么办法吗？” - 面试官，微微一笑：“
C#怎么没有大小数？decimal和双精度符点你没听说过吗？” - 我，坚定地反驳：“不，这些还不够长，比如我要计算1后面有1万个0的数字，
C#就算不了咯。” - 面试官，眼神犀利：“你的问题很有趣，如果真的遇到了
C#解决不了的问题，那可能更多的是我们需要重新审视这个问题……” - 我，心里默默立下决心：“……”
就这样，这场面试，像一颗种子，在我心里种下了对高精度数值计算的追求。我知道大家都会说，谁还不会算个数呢？但其实，还真有很多需要高精度数值计算的情况。
比如你正在计算一个飞向火星的火箭的轨道，一个微小的计算误差，可能就会让火箭偏离数百万公里。或者你正在使用GPS导航，在城市的密集街道中，几十米的误差就可能让你迷路。又或者你是一个气候科学家，正在预测全球变暖的趋势，一点点的计算误差，就可能导致模型的预测结果大相径庭。
从那个时候起，我开始关注一些技术问题的本质，而不是仅仅将需求解决。我也很想知道C#/.NET的极限在哪，但可惜那个时候我还只是一个function caller，能力有限。
但随着时间的推移，我逐渐积累了更多的知识和技能。为实现当年的梦想，今年(2023)年初趁过年放假期间，我把自己关在家里，连续几个晚上熬夜工作，基于GMP和MPFR两个知名的开源项目，最终成功开发了.NET的高精度数值计算库：Sdcb.Arithmetic，现在经过多个版本的迭代，已经相当稳定了。
市面上已有的GMP和MPFR封装库及其问题
在打造我的.NET高精度数值计算库之前，我了解到市面上已经存在一些GMP和MPFR的封装库，如machinecognitis/Math.Gmp.Native和emphasis87/mpfr.NET。然而，我发现它们存在一系列的问题，这也是促使我开发新库的原因。
首先，让我们看看machinecognitis/Math.Gmp.Native。这个项目的主要问题是，尽管它提供了对GMP库的封装，但是这个封装仅限于低级API，换句话说，你需要对GMP库有深入的理解才能有效地使用这个库。这个项目没有提供任何高级API，因此它的易用性相当差。作为开发者，我们自然希望能够尽可能地提升开发效率，而这需要高级API的支持。一个好的封装库应该提供既直观又方便的API，而不是仅仅提供底层函数的封装。
接下来，我们看看emphasis87/mpfr.NET。这个项目虽然提供了对MPFR库的封装，但是这个项目已经有些年头了。更糟糕的是，它是通过C++/CLI来实现的，这使得它对.NET Core的支持并不是很好，同时也限制了它在Linux上的使用。在现今这个跨平台开发日益重要的时代，一个好的库应该能够在多种平台上进行无缝的运行。
此外，上述两个项目都存在一个共同的问题，那就是它们的版本都过于陈旧，且很久没有得到维护和更新。在快速变化的技术世界中，一个库如果连续数年没有更新，那么它可能会失去与最新技术接轨的机会，而这对于用户来说是无法接受的。
当然，除了以上的原因外，还有一个更重要的原因驱使我创造新的库，那就是我想要做出一个更好用的GMP和MPFR封装库。我相信，只有当我们不断地挑战自我，才能创造出更好的产品。在接下来的章节中，我将介绍我是如何实现这一目标的。
NuGet包简介
这个项目分为两部分，GMP和MPFR：
GMP可以支持高精度整数、小数和分数，然而高精度小数的功能有限，例如不支持三角函数Sin/Cos。MPFR主要用于处理高精度小数，功能更为丰富，提供超过300个MPFR库函数。
如果你熟悉我的另一个项目PaddleSharp，你会发现这里同样需要同时安装.NET封装包和动态库包。其中带runtime的包为动态库包（例如runtime.win64表示支持64位Windows）。值得一提的是，MPFR依赖于GMP库，因此如果你安装Sdcb.Arithmetic.Mpfr，系统会自动带上Sdcb.Arithmetic.Gmp。
对于Linux用户，你可能并不需要安装我的Linux动态库包。Linux系统大部分都自带了libgmp.so动态库，你还可以通过系统自带的包管理工具安装libmpfr.so。例如，Ubuntu 22.04用户可以使用以下命令进行安装（这也意味着你不需要安装*.runtime.linux64的NuGet包）：
sudo apt-get install libmpfr-dev 最后，所有的Windows动态库包都是由我自己使用vcpkg编译的，而Linux动态库包则来自Ubuntu 22.04。因此，如果你使用我的动态库包，可能主要只支持Ubuntu 22.04和Debian。
值得一提的是，GMP和MPFR的动态库是GPL或者LGPL协议，因此我的动态库nuget包和.NET封装包的协议不相同，.NET封装包是MIT协议，动态库包是LGPL协议。
libgmp
| Package Id | Version | License | Notes |
|---|---|---|---|
| Sdcb.Arithmetic.Gmp |  |
MIT | .NET binding for libgmp |
| Sdcb.Arithmetic.Gmp.runtime.win64 | |
LGPL | native lib in windows x64 |
| Sdcb.Arithmetic.Gmp.runtime.win32 | |
LGPL | native lib in windows x86 |
| Sdcb.Arithmetic.Gmp.runtime.linux64 | |
LGPL | native lib in Linux x64 |
mpfr
| Package Id | Version | License | Notes |
|---|---|---|---|
| Sdcb.Arithmetic.Mpfr | |
MIT | .NET binding for libmpfr |
| Sdcb.Arithmetic.Mpfr.runtime.win64 | |
LGPL | native lib in windows x64 |
| Sdcb.Arithmetic.Mpfr.runtime.win32 | |
LGPL | native lib in windows x86 |
| Sdcb.Arithmetic.Mpfr.runtime.linux64 | |
LGPL | native lib in linux x64 |
使用示例 - 计算2^65536的最后20位数字
在这个章节中，我们首先展示如何使用.NET自带的大整数类进行数值计算，然后将演示如何使用我们的库Sdcb.Arithmetic进行同样的计算。然后，我们将展示如何通过使用C API来优化我们的库，最后，我们将介绍一种更高级的优化策略——使用InplaceAPI。
原生.NET实现
许多程序员可能已经熟悉.NET Framework 3.5引入的大整数类System.Numeric.BigInteger。我们可以使用这个类来计算2^65536的最后20位数字，代码如下：
Stopwatch sw = Stopwatch.StartNew(); int count = 10; BigInteger b = new BigInteger(); for (int c = 0; c < count; ++c) { b = 1; for (int i = 1; i <= 65536; ++i) { b *= 2; } } Console.WriteLine($"耗时：{sw.Elapsed.TotalMilliseconds / count:F2}ms"); Console.WriteLine($"2^65536最后20位数字：{b.ToString()[^20..]}"); 在我的i9-9880h电脑上，这个程序的运行结果如下：
耗时：94.00ms 2^65536最后20位数字：45587895905719156736 使用Sdcb.Arithmetic库
下面我们将展示如何使用我们的库Sdcb.Arithmetic来进行同样的计算。为了安装这个库，你需要使用以下的NuGet包：Sdcb.Arithmetic.Gmp和Sdcb.Arithmetic.Gmp.runtime.win64（或其它对应环境包）。
Stopwatch sw = Stopwatch.StartNew(); int count = 10; GmpInteger b = new GmpInteger(); for (int c = 0; c < count; ++c) { b = 1; for (int i = 1; i <= 65536; ++i) { b *= 2; } } Console.WriteLine($"耗时：{sw.Elapsed.TotalMilliseconds / count:F2}ms"); Console.WriteLine($"2^65536最后20位数字：{b.ToString()[^20..]}"); 运行结果如下：
耗时：89.52ms 2^65536最后20位数字：45587895905719156736 可以看到，Sdcb.Arithmetic.Gmp的结果与.NET原生实现相匹配，并且计算速度相近。
性能优化
虽然上述Gmp代码已经达到了与原生.NET实现相近的性能，但是我们可以通过使用底层的C API来进一步优化我们的库。以下是优化后的代码：
// 安装NuGet包：Sdcb.Arithmetic.Gmp // 安装NuGet包：Sdcb.Arithmetic.Gmp.runtime.win64 （或其它对应环境包） // 函数需要标注unsafe // 项目需要启用unsafe编译选项 Stopwatch sw = Stopwatch.StartNew(); int count = 10; Mpz_t mpz; GmpLib.__gmpz_init((IntPtr)(&mpz)); for (int c = 0; c < count; ++c) { GmpLib.__gmpz_set_si((IntPtr)(&mpz), 1); for (int i = 1; i <= 65536; ++i) { GmpLib.__gmpz_mul_si((IntPtr)(&mpz), (IntPtr)(&mpz), 2); } } sw.Stop(); IntPtr ret = GmpLib.__gmpz_get_str(IntPtr.Zero, 10, (IntPtr)(&mpz)); string wholeStr = Marshal.PtrToStringUTF8(ret)!; GmpMemory.Free(ret); Console.WriteLine($"耗时：{sw.Elapsed.TotalMilliseconds / count:F2}ms"); Console.WriteLine($"2^65536最后20位数字：{wholeStr[^20..]}"); GmpLib.__gmpz_clear((IntPtr)(&mpz)); 在同一台电脑上，输出结果如下：
耗时：20.87ms 2^65536最后20位数字：45587895905719156736 你可以参考下面的源代码了解我是如何封装PInvoke API的：
- https://github.com/sdcb/Sdcb.Arithmetic/blob/master/Sdcb.Arithmetic.Gmp/GmpLib.cs
- https://github.com/sdcb/Sdcb.Arithmetic/blob/master/Sdcb.Arithmetic.Mpfr/MpfrLib.cs
20.87ms显然比之前基于高级API封装的GmpInteger速度89.52ms快很多，但易用性却非常差，这些代码会很难理解、很难维护、也很容易造成内存泄露问题。
值得注意的是，根据最佳实践，上面的代码涉及内存释放时（__gmpz_clear和GmpMemory.Free），本质需要改为多个try...finally来避免在调用前触发异常而导致内存泄露，但加上try...finally会导致代码简洁性进一步下降。
速度易用两全其美 - 使用InplaceAPI优化
为了做到在性能和可维护性之间找到平衡，我还开发了InplaceAPI，这是一个直接调用底层API，但用起来很方便的函数。以下是使用该API后的代码：
// 安装NuGet包：Sdcb.Arithmetic.Gmp // 安装NuGet包：Sdcb.Arithmetic.Gmp.runtime.win64 （或其它对应环境包） Stopwatch sw = Stopwatch.StartNew(); int count = 10; using GmpInteger b = new GmpInteger(); for (int c = 0; c < count; ++c) { b.Assign(1); for (int i = 1; i <= 65536; ++i) { GmpInteger.MultiplyInplace(b, b, 2); } } Console.WriteLine($"耗时：{sw.Elapsed.TotalMilliseconds / count:F2}ms"); Console.WriteLine($"2^65536最后20位数字：{b.ToString()[^20..]}"); 运行结果如下：
耗时：21.13ms 2^65536最后20位数字：45587895905719156736 可见代码相比最开始的和如下变化：
GmpInteger加上了using，值得一提的是，我同时也写了Finalizer，因此它同时也支持不写using，但内存回收较慢；- 赋值时不使用
=，而是改为调用.Assign()函数，这样可以避免创建临时对象 - 计算乘法时，我使用了
MultiplyInplace()而非运算符重载，这样可以省掉创建大量临时对象；
这个速度21.13ms相比纯粹使用C API的20.87ms稍慢，但这几乎可以认为是误差范围之内，但代码的简洁性、可维护性比C API简单很多。
对比Java
作为参考，当然少不了和Java的性能对比，这是用于对比的Java代码：
import java.math.BigInteger; import java.time.Duration; import java.time.Instant; public class Main { public static void main(String[] args) { Instant start = Instant.now(); int count = 10; BigInteger b = BigInteger.ONE; for (int c = 0; c < count; ++c) { b = BigInteger.ONE; for (int i = 1; i <= 65536; ++i) { b = b.multiply(BigInteger.valueOf(2)); } } Instant finish = Instant.now(); long timeElapsed = Duration.between(start, finish).toMillis(); String str = b.toString(); String last20Digits = str.substring(str.length() - 20); System.out.printf("耗时：%f msn", (double) timeElapsed / count); System.out.println("2^65536最后20位数字：" + last20Digits); } } 我使用的Java版本是OpenJDK version "11.0.16.1" 2022-08-12 LTS，使用相同的电脑，输出结果如下：
耗时：103.100000 ms 2^65536最后20位数字：45587895905719156736 可见速度比.NET原生的BigInteger稍慢。
总结 - 性能比较表格
| 实现方式 | 平均耗时（ms） | 结果 |
|---|---|---|
原生.NET实现 |
94.00 | 45587895905719156736 |
无优化的Sdcb.Arithmetic |
89.52 | 45587895905719156736 |
使用C API优化的Sdcb.Arithmetic |
20.87 | 45587895905719156736 |
使用InplaceAPI优化的Sdcb.Arithmetic |
21.13 | 45587895905719156736 |
Java - BigInteger |
103.10 | 45587895905719156736 |
使用示例 - 计算100万位π
在这个示例中，我将展示Sdcb.Arithmetic.Gmp计算小数点后100万位π的计算方式，这段代码著名的Chudnovsky算法来计算圆周率π的值，该算法由Chudnovsky兄弟在1980年代提出的，它基于Ramanujan的公式，但在计算精度和效率上进行了改进。这个算法的一个显著特点是它的超级线性收敛性，即每增加一个迭代步骤，就可以大幅增加精确到的小数位数。实际上，Chudnovsky算法每次迭代大约能增加14个准确的小数位数。这使得它非常适合计算几百万位甚至更高的π值。
// Install NuGet package: Sdcb.Arithmetic.Gmp // Install NuGet package: Sdcb.Arithmetic.Gmp.runtime.win-x64(for windows) using Sdcb.Arithmetic.Gmp; Stopwatch sw = Stopwatch.StartNew(); using GmpFloat pi = CalcPI(); double elapsed = sw.Elapsed.TotalMilliseconds; Console.WriteLine($"耗时：{elapsed:F2}ms"); Console.WriteLine($"结果：{pi:N1000000}"); GmpFloat CalcPI(int inputDigits = 1_000_000) { const double DIGITS_PER_TERM = 14.1816474627254776555; // = log(53360^3) / log(10) int DIGITS = (int)Math.Max(inputDigits, Math.Ceiling(DIGITS_PER_TERM)); uint PREC = (uint)(DIGITS * Math.Log2(10)); int N = (int)(DIGITS / DIGITS_PER_TERM); const int A = 13591409; const int B = 545140134; const int C = 640320; const int D = 426880; const int E = 10005; const double E3_24 = (double)C * C * C / 24; using PQT pqt = ComputePQT(0, N); GmpFloat pi = new(precision: PREC); // pi = D * sqrt((mpf_class)E) * PQT.Q; pi.Assign(GmpFloat.From(D, PREC) * GmpFloat.Sqrt((GmpFloat)E, PREC) * (GmpFloat)pqt.Q); // pi /= (A * PQT.Q + PQT.T); GmpFloat.DivideInplace(pi, pi, GmpFloat.From(A * pqt.Q + pqt.T, PREC)); return pi; PQT ComputePQT(int n1, int n2) { int m; if (n1 + 1 == n2) { PQT res = new() { P = GmpInteger.From(2 * n2 - 1) }; GmpInteger.MultiplyInplace(res.P, res.P, 6 * n2 - 1); GmpInteger.MultiplyInplace(res.P, res.P, 6 * n2 - 5); GmpInteger q = GmpInteger.From(E3_24); GmpInteger.MultiplyInplace(q, q, n2); GmpInteger.MultiplyInplace(q, q, n2); GmpInteger.MultiplyInplace(q, q, n2); res.Q = q; GmpInteger t = GmpInteger.From(B); GmpInteger.MultiplyInplace(t, t, n2); GmpInteger.AddInplace(t, t, A); GmpInteger.MultiplyInplace(t, t, res.P); // res.T = (A + B * n2) * res.P; if ((n2 & 1) == 1) GmpInteger.NegateInplace(t, t); res.T = t; return res; } else { m = (n1 + n2) / 2; PQT res1 = ComputePQT(n1, m); using PQT res2 = ComputePQT(m, n2); GmpInteger p = res1.P * res2.P; GmpInteger q = res1.Q * res2.Q; // t = res1.T * res2.Q + res1.P * res2.T GmpInteger.MultiplyInplace(res1.T, res1.T, res2.Q); GmpInteger.MultiplyInplace(res1.P, res1.P, res2.T); GmpInteger.AddInplace(res1.T, res1.T, res1.P); res1.P.Dispose(); res1.Q.Dispose(); return new PQT { P = p, Q = q, T = res1.T, }; } } } public ref struct PQT { public GmpInteger P; public GmpInteger Q; public GmpInteger T; public readonly void Dispose() { P?.Dispose(); Q?.Dispose(); T?.Dispose(); } } 在我的i9-9880h电脑中，输出如下（100万位中间有...省略）：
耗时：435.35ms 结果：3.141592653589793238462643383...83996346460422090106105779458151 可见速度是非常快的，100万位π的值可以在这个链接进行参考。
同时也作为比较，我也参考了Github上网友写的这段C++计算100万位π的代码：https://gist.github.com/komasaru/68f209118edbac0700da
我使用VS2022通过Release-x64编译，在同一台电脑中，输出如下：
**** PI Computation ( 1000000 digits ) TIME (COMPUTE): 0.425 seconds. TIME (WRITE) : 0.103 seconds. 我也参考了Github上另一段用C写的同样计算100万位π的代码：https://github.com/natmchugh/pi/blob/master/gmp-chudnovsky.c
我同样使用VS2022通过Release-x64编译，输出如下：
#terms=70513, depth=18 sieve time = 0.003 .................................................. bs time = 0.265 gcd time = 0.000 div time = 0.037 sqrt time = 0.022 mul time = 0.014 total time = 0.343 P size=1455608 digits (1.455608) Q size=1455601 digits (1.455601) 这是3种编程语言计算100万位π耗时的比较表格：
| 耗时（ms） | |
|---|---|
| C# | 435.35 |
| C++ | 425 |
| C | 343 |
可见C#和C++几乎不相上下（在我的本地测试中甚至有时C#更快），使用C确实稍快，但从代码可以看出C的实现方式有许多优化，比如递归求值部分的内存是一次性分配的，速度快的主要原因是算法有优化。
尾声
在这篇文章中，我分享了我从大学毕业的第一场面试中获得的启示，以及我如何把这种启示转化为实践，打造了一个.NET的高精度数值计算库——Sdcb.Arithmetic。这个开源项目弥补了C#在处理大数运算方面的不足，打破了Java的BigFloat优势，使得C#也能轻松处理高精度计算的需求。
这个库包括GMP和MPFR两部分，前者支持高精度的整数、小数和分数的运算，后者则是处理高精度小数的利器，提供了超过300个MPFR库函数。通过我对这个库的介绍和展示，你可以看到这个库在处理大数计算，例如计算2^65536的最后20位数字，或者计算π的百万位小数上的表现非常出色。
我深信开源的力量，因此我把这个项目开源了，并希望能够帮助到所有需要高精度数值计算的.NET开发者。想尝试Sdcb.Arithmetic的朋友，欢迎访问我的Github，我希望你能去我的项目主页上给我一个star🌟，这将对我是莫大的鼓励。我会继续保持对开源的热爱，做出更多有价值的贡献。
喜欢的朋友，也请关注我的微信公众号：【DotNet骚操作】
![洛谷 P11345 [KTSC 2023 R2] 基地简化 题解](http://www.itfaba.com/wp-content/themes/kemi/timthumb.php?src=http://www.itfaba.com/wp-content/themes/kemi/img/random/3.jpg&w=218&h=124&zc=1)
