
科学计算在今天已经完全可以使用 C 了。我不是从事科学计算这一领域的工作的不过在工作中也多多少少接触到了一些计算相关的库。今天我就给你介绍几个有用的计算库。Armadillo说到计算你可能首先会想到矩阵、矢量这些东西吧这些计算确实就是科学计算中的常见内容了。这些领域的标准即是一些 Fortran 库定下的如BLAS [1]LAPACK [2]ARPACK [3]它们的实现倒不一定用 Fortran尤其是 BLASOpenBLAS [4] 是用汇编和 C 语言写的Intel MKL [5] 有针对 Intel 的特定 CPU 指令集进行优化的汇编代码Mir GLAS [6] 是用 D 语言写的不管实现的方法是哪一种暴露出来的函数名字是这个样子的ddotdgemvdsyrksgemm……这个接口的唯一好处应该就是它是跨语言并且跨实现的。所以使用这些函数时你可以切换不同的实现而不需要更改代码。唯一需要修改的通常就是链接库的名字或位置而已。假设我们需要做一个简单的矩阵运算对一个矢量进行旋转这么一个简单的操作用纯 C 接口的 BLAS 来表达有点痛苦你需要使用的大概是 dgemv_ 函数而这个函数需要 11 个参数我查阅了一下资料之后也就放弃了给你展示一下如何调用 dgemv_ 的企图我们还是老老实实地看一下在现代 C 里的写法吧#include armadillo #include cmath #include iostream using namespace std; int main() { // 代表位置的向量 arma::vec pos{1.0, 0.0}; // 旋转矩阵 auto pi arma::datum::pi; double angle pi / 2; arma::mat rot { {cos(angle), -sin(angle)}, {sin(angle), cos(angle)}}; cout Current position:\n pos; cout Rotating angle * 180 / pi deg\n; arma::vec new_pos rot * pos; cout New position:\n new_pos; }这就是使用 Armadillo [7] 库来实现矢量旋转的代码。这个代码基本就是上面的数学公式的一一对应了。代码相当直白我只需要稍稍说明一下所有的 Armadillo 的类型和函数都定义在 arma 名空间下。Armadillo 在 arma::datum 下定义了包括 pi 和 e 在内的一些数学常量。vec 是矢量类型mat 是矩阵类型这两个类型实际上是 Col 和 Mat 的缩写别名。Armadillo 支持使用 C11 的列表初始化语法来初始化对象。Armadillo 支持使用流来输出对象。上面代码的输出为Current position:1.00000Rotating 90 degNew position:6.1232e-171.0000e00输出里面的 6.1232e-17 是浮点数表示不精确的后果把它理解成 0 就对了。我们上面已经提到了 vec 实际上是 Col双精度浮点数类型的列矢量。自然Armadillo 也有行矢量 rowvec即 Row也可以使用其他的数字类型如 int、 float 和 complex。此外除了大小不确定的线性代数对象之外Armadillo 也提供了固定大小的子类型如 vec::fixed2 和 mat::fixed2, 2为方便使用还提供了不少别名如 imat22 代表 Mat::fixed2, 2 等。固定大小的对象不需要动态内存分配使用上有一定的性能优势。Armadillo 是一个非常复杂的库它的头文件数量超过了 500 个。我们今天不可能、也不必要描述它的所有功能只能稍稍部分列举一下除了目前提到的列矢量、行矢量和矩阵外Armadillo 也支持三维的数据立方体Cube 模板。Armadillo 支持稀疏矩阵SpMat 模板。除了数学上的加、减、乘运算Armadillo 支持按元素的乘法、除法、相等、不等、小于比较等使用 %、/、、!、 等运算结果的大小跟参数相同每个元素是相应运算的结果。某些运算符可能不太直观尤其是 %不是取模和 返回不是单个布尔值而是矩阵。Armadillo 支持对非固定大小的矢量、矩阵和立方体改变其大小.reshape() 和 resize()。Armadillo 可以方便地按行.col()、列.row()、对角线.diag()读写矩阵的内容包括用一个矢量去改写矩阵的对角线。Armadillo 可以方便地对矩阵进行转置.t()、求反.inv()。Armadillo 可以对矩阵进行特征分解eigen_sym()、eigen_gen() 等。Armadillo 支持傅立叶变换fft()、fft2() 等。Armadillo 支持常见的统计计算如平均值、中位值、标准偏差等mean()、median()、stddev() 等。Armadillo 支持多项式方程求根roots。Armadillo 支持 k‐平均聚类k-means clustering算法kmeans。等等。如果你需要用到这些功能你可以自己去查看一下具体的细节我们这儿只提几个与编程有关的细节。对象的输出我们上面已经展示了直接把对象输出到一个流。我们的写法是cout Current position:\n pos;实际上基本等价于调用 print 成员函数pos.print(Current position:);这个写法可能会更简单些。此外在这两种情况输出的格式都是 Armadillo 自动控制的。如果你希望自己控制的话可以使用 raw_print 成员函数。比如对于上面代码里对 new_pos 的输出我们可以写成需要包含 cout fixed setw(9) setprecision(4); new_pos.raw_print( cout, New position:);这种情况下你可以有效地对格式、宽度和精度进行设置能得到New position:0.00001.0000记得我们说过 vec 是 Col 的别名因此输出是多行的。我们要输出成单行的话转置transpose一下就可以了cout fixed setw(9) setprecision(4); new_pos.t().raw_print( cout, New position:);输出为New position:0.0000 1.0000表达式模板如果你奇怪前面 dgemv_ 为什么有 11 个参数这里有个我没有提的细节是它执行的实际上是个复合操作如果你只是简单地做乘法的话就相当于 α 为 1、β 为 0 的特殊情况。那么问题来了如果你真的写了类似于上面这样的公式的话编译器和线性代数库能不能转成合适的调用、而没有额外的开销呢答案是至少在某些情况下是可以的。秘诀就是表达式模板expression template那什么是表达式模板呢我们先回过去看我上面的例子。有没有注意到我写的是arma::vec new_pos rot * pos;而没有使用 auto 来声明其中部分的原因是rot * pos 的类型并不是 vec而是const GlueMatdouble, Coldouble, glue_times换句话说结果是一个表达式而并没有实际进行计算。如果我用 auto 的话行为上似乎一切都正常但我每次输出这个结果时都会重新进行一次矩阵的乘法而我用 arma::vec 接收的话构造时就直接进行了计算存储了表达式的结果。上面的简单例子不能实际触发对 dgemv_ 的调用我用下面的代码实际验证出了表达式模板产生的优化fill::randu 表示对矢量和矩阵的内容进行随机填充#include armadillo #include iostream using namespace std; using namespace arma; int main() { vec x(8, fill::randu); mat r(8, 8, fill::randu); vec result 2.5 * r * x; cout result; }赋值语句右边的类型是const GlueeOpMatdouble, eop_scalar_times, Coldouble, glue_times当使用这个表达式构造 vec 时就会实际发生对 dgemv_ 的调用。我也确实跟踪到了在将要调用 dgemv_ 时标量值 2.5 确实在参数 alpha 指向的位置上这个接口的参数都是指针。从上面的描述可以看到表达式模板是把双刃剑既可以提高代码的性能又能增加代码被误用的可能性。在可能用到表达式模板的地方你需要注意这些问题。平台细节Armadillo 的文档里说明了如何从源代码进行安装但在 Linux 和 macOS 下通过包管理器安装可能是更快的方式。在 CentOS 下可使用 sudo yum install armadillo-devel在 macOS 下可使用 brew install armadillo。使用包管理器一般也会同时安装常见的依赖软件如 ARPACK 和 OpenBLAS。在 Windows 上Armadillo 的安装包里自带了一个基本版本的 64 位 BLAS 和 LAPACK 库。如果需要更高性能或 32 位版本的话就需要自己另外去安装了。除非你只是做一些非常简单的线性代数计算就像我今天的例子那直接告诉 Armadillo 不要使用第三方库也行。cl /EHsc /DARMA_DONT_USE_BLAS /DARMA_DONT_USE_LAPACK …Boost.Multiprecision众所周知C 和 C甚至推而广之到大部分的常用编程语言里的数值类型是有精度限制的。比如上一讲的代码里我们就用到了 INT_MIN最小的整数。很多情况下使用目前这些类型是够用的最高一般是 64 位整数和 80 位浮点数。但也有很多情况这些标准的类型远远不能满足需要。这时你就需要一个高精度的数值类型了。有一次我需要找一个高精度整数类型和计算库最后找到的就是 Boost.Multiprecision [9]。它基本满足我的需求以及一般意义上对库的期望正确实现我需要的功能接口符合直觉、易用有良好的性能正确实现功能这点我就不多讲了。这是一个基本出发点没有太多可讨论的地方。在我上次的需求里对性能其实也没有很高的要求。让我对 Boost.Multiprecision 满意的主要原因就是它的接口了。接口易用性我在提到了 CLN。它对我来讲就是个反面教材。它的整数类型不仅不提供 % 运算符居然还不提供 / 运算符它强迫用户在下面两个方案中做出选择使用 truncate2 函数得到一个商数和余数使用 exquo 函数当且仅当可以整除的时候不管作者的设计原则是什么这简直就是易用性方面的灾难了——不仅这些函数要查文档才能知晓而且有的地方我真的只需要简单的除法呀……哦对了它在 Windows 编译还很不方便而我那时用的正是 Windows。Boost.Multiprecision 的情况则恰恰相反让我当即大为满意使用基本的 cpp_int 对象不需要预先编译库只需要 Boost 的头文件和一个好的编译器。常用运算符 、-、*、/、% 一个不缺全部都有。可以自然地通过整数和字符串来进行构造。提供了用户自定义字面量来高效地进行初始化。在使用 IO 流时输入输出既可以使用十进制也可以通过 hex 来切换到十六进制。下面的代码展示了它的基本功能#include iomanip #include iostream #include boost/multiprecision/cpp_int.hpp using namespace std; int main() { using namespace boost:: multiprecision::literals; using boost::multiprecision:: cpp_int; cpp_int a 0x123456789abcdef0_cppi; cpp_int b 16; cpp_int c{0400}; cpp_int result a * b / c; cout hex result endl; cout dec result endl; }输出是123456789abcdef81985529216486895我们可以看到cpp_int 可以通过自定义字面量后缀 _cppi只能十六进制来初始化可以通过一个普通整数来初始化也可以通过字符串来初始化并可以使用 0x 和 0 前缀来选择十六进制和八进制。拿它可以正常地进行加减乘除操作也可以通过 IO 流来输入输出。性能Boost.Multiprecision 使用了表达式模板和 C11 的移动来避免不必要的拷贝。后者当然是件好事而前者曾经坑了我一下——我第一次使用 Boost.Multiprecision 时非常困惑为什么我使用 half(n - 1) 调用下面的简单函数居然会编译不过template typename N inline N half(N n) { return n / 2; }我的意图当然是 N 应当被推导为 cpp_inthalf 的结果也是 cpp_int。可实际上n - 1 的结果跟上面的 Armadillo 展示的情况类似是另外一个单独的类型。我需要把 half(n - 1) 改写成 half(N(n - 1)) 才能得到期望的结果。我做的计算挺简单并不觉得表达式模板对我的计算有啥帮助所以我最后是禁用了表达式模板typedef boost::multiprecision:: number boost::multiprecision:: cpp_int_backend, boost::multiprecision::et_off int_type;类似于 Armadillo 可以换不同的 BLAS 和 LAPACK 实现Boost.Multiprecision 也可以改换不同的后端。比如如果我们打算使用 GMP [10] 的话我们需要包含利用 GMP 的头文件并把上面的 int_type 的定义修正一下#include boost/multiprecision/gmp.hpp typedef boost::multiprecision:: number boost::multiprecision::gmp_int, boost::multiprecision::et_off int_type;注意我并不是推荐你换用 GMP。如果你真的对性能非常渴求的话应当进行测试来选择合适的后端。否则缺省的后端易用性最好——比如使用 GMP 后端就不能使用自定义字面量了。我当时寻找高精度算术库是为了做 RSA 加解密。计算本身不复杂属于编程几小时、运行几毫秒的情况。如果你有兴趣的话可以看一下我那时的挑选过程和最终代码 [11]。Boost 里好东西很多远远不止这一样。下一讲我们就来专门聊聊 Boost。内容小结性能和易用性往往是有矛盾的。你对性能和易用性有什么样的偏好呢欢迎留言与我分享。参考资料[1] Wikipedia, “Basic Linear Algebra Subprograms”. https://en.wikipedia.org/wiki/Basic_Linear_Algebra_Subprograms[2] Wikipedia, “LAPACK”. https://en.wikipedia.org/wiki/LAPACK[3] Wikipedia, “ARPACK”. https://en.wikipedia.org/wiki/ARPACK[4] Zhang Xianyi et al., OpenBLAS. https://github.com/xianyi/OpenBLAS[5] Intel, Math Kernel Library. https://software.intel.com/mkl[6] Ilya Yaroshenko, mir-glas. https://github.com/libmir/mir-glas[7] Conrad Sanderson and Ryan Curtin, “Armadillo: C library for linear algebra scientific computing”. http://arma.sourceforge.net/[8] Wikipedia, “Expression templates”. https://en.wikipedia.org/wiki/Expression_templates[9] John Maddock, Boost.Multiprecision. https://www.boost.org/doc/libs/release/libs/multiprecision/doc/html/index.html[10] The GNU MP bignum library. https://gmplib.org/[11] 吴咏炜, “Choosing a multi-precision library for C—a critique”. https://yongweiwu.wordpress.com/2016/06/04/choosing-a-mu