C#的数学库: ALGLIB, ILNumerics, Dew.Math, Math.NET Numerics等

猿柒

6个月前更新

100

https://zhuanlan.zhihu.com/p/12783824787

以下4个为商业库，性能都不错，发布的dll文件均未加密，容易逆向。(等我以后有空了，再来详细对比一下4个商业库的性能)

ALGLIB(C++/C#/Java numerical analysis library): 有免费版和商业版，支持C++/C#/Java/Delphi/CPython五种语言，涵盖数据分析(分类与回归、统计学)、优化与非线性问题求解、插值与最小二乘拟合、线性代数、FFT等领域。免费版是开源的，官网下载的压缩包中就包含源代码，相对商业版缺少多线程、native HPC Kernels和Intel MKL，不想花钱的首选ALGLIB免费版；
ILNumerics(ILNumerics – Technical Computing): 只有收费版，支持多核并行计算，绘图功能比较强大，涵盖线性代数、FFT等领域，支持HDF5文件；
Dew.Math(Math Library Component Delphi Component C++ Builder Component Net): 只有收费版，涵盖矩阵和向量运算、DSP、统计学、数据挖掘、FFT等领域，利用AVX, AVX2, AVX512以及多核心加速。
Extreme.Numerics(Build financial, engineering and scientific applications faster): 只有收费版，涵盖统计学、概率论、回归、机器学习、线性、代数、矩阵、向量、求解、优化、曲线拟合、fft、积分、微分、插值等领域。

以下均为开源库，性能方面大多不如商业库，

SdcbArithmetic(https://github.com/sdcb/Sdcb.Arithmetic): 开源，作者的知乎账号是 @周杰，P/Invoke调用GMP(GNU Multiple Precision Arithmetic Library)和MPFR(multiple-precision floating-point computations)库，性能不错，可以进行高精度计算。
MathNET Numerics(GitHub – mathnet/mathnet-numerics): 开源，知名度比较高，但项目更新频率低，性能很差，bug很多(下面有几个具体的例子)，只能做做玩具级的应用，或者学生写作业用。
AngouriMath(GitHub asc-community/AngouriMath): 符号运算库，可以解方程(组)，求符号导数，解决集合问题等。
CSharpMath(GitHub – verybadcat/CSharpMath): 根据Latex表达式渲染公式和符号。
StringMath(GitHub – miroiu/string-math): 对字符串形式的数学表达式求值，比如double result = “1 * (2 – 3) ^ 2”.Eval(); // 1.
ncalc(GitHub – ncalc/ncalc): 对字符串形式的数学表达式求值。
DecimalMath(GitHub – nathanpjones/DecimalMath): 包含Decimal版的Sqrt, Pow, Exp, Log, in, Cos, Tan, ASin, ACos, ATan, ATan2. 性能很差(下面有例子)。

以下为停止维护的开源库。

MetaNumerics: github仓库最后一次commit的时间是2020-08-22，基本等于停止维护；
AForgeNET: github仓库最后一次commit的时间是2020-06-19，基本等于停止维护；
AccordNET: github仓库已于2020-11-19转为archived状态，不再更新；

这里推荐一下本人的另一篇博客：C#数学相关开发性能优化方法

免费版ALGLIB 4.04.0的源代码中，diffequations.cs文件的第571~605行，用除法定义浮点常数时加了大量的(double)进行强制类型转换，显得十分繁琐和丑陋，不知为何要这么写，而不是写成1.0 / 5.0, 3.0 / 10.0等。

// Cask-Karp solver // Prepare coefficients table. // Check it for errors // state.rka = new double[6]; state.rka[0] = 0; state.rka[1] = (double)1/(double)5; state.rka[2] = (double)3/(double)10; state.rka[3] = (double)3/(double)5; state.rka[4] = 1; state.rka[5] = (double)7/(double)8; state.rkb = new double[6, 5]; state.rkb[1,0] = (double)1/(double)5; state.rkb[2,0] = (double)3/(double)40; state.rkb[2,1] = (double)9/(double)40; state.rkb[3,0] = (double)3/(double)10; state.rkb[3,1] = -((double)9/(double)10); state.rkb[3,2] = (double)6/(double)5; state.rkb[4,0] = -((double)11/(double)54); state.rkb[4,1] = (double)5/(double)2; state.rkb[4,2] = -((double)70/(double)27); state.rkb[4,3] = (double)35/(double)27; state.rkb[5,0] = (double)1631/(double)55296; state.rkb[5,1] = (double)175/(double)512; state.rkb[5,2] = (double)575/(double)13824; state.rkb[5,3] = (double)44275/(double)110592; state.rkb[5,4] = (double)253/(double)4096; state.rkc = new double[6]; state.rkc[0] = (double)37/(double)378; state.rkc[1] = 0; state.rkc[2] = (double)250/(double)621; state.rkc[3] = (double)125/(double)594; state.rkc[4] = 0; state.rkc[5] = (double)512/(double)1771; state.rkcs = new double[6]; state.rkcs[0] = (double)2825/(double)27648; state.rkcs[1] = 0; state.rkcs[2] = (double)18575/(double)48384; state.rkcs[3] = (double)13525/(double)55296; state.rkcs[4] = (double)277/(double)14336; state.rkcs[5] = (double)1/(double)4; state.rkk = new double[6, n];

在ALGLIB的statistics.cs文件中，7704~7781行，写了19条if语句，没有用switch，也没有用数组，这实在是太过丑陋和低效了。如果在每个if中都写上return xxx，那也尚可接受，可是它没有return，做了大量的无效比较操作，令人唏嘘。

r = 0; w = (int)Math.Round(-(7.141428e+00*s)+1.800000e+01); if( w>=18 ) { r = -6.399e-01; } if( w==17 ) { r = -7.494e-01; } if( w==16 ) { r = -8.630e-01; } if( w==15 ) { r = -9.913e-01; } if( w==14 ) { r = -1.138e+00; } if( w==13 ) { r = -1.297e+00; } if( w==12 ) { r = -1.468e+00; } if( w==11 ) { r = -1.653e+00; } if( w==10 ) { r = -1.856e+00; } if( w==9 ) { r = -2.079e+00; } if( w==8 ) { r = -2.326e+00; } if( w==7 ) { r = -2.601e+00; } if( w==6 ) { r = -2.906e+00; } if( w==5 ) { r = -3.243e+00; } if( w==4 ) { r = -3.599e+00; } if( w==3 ) { r = -3.936e+00; } if( w==2 ) { r = -4.447e+00; } if( w==1 ) { r = -4.852e+00; } if( w<=0 ) { r = -5.545e+00; } result = r; return result;

MathNET Numerics

优点: 包含了线性代数(OpenBLAS)、特殊函数、概率论、回归分析、插值与拟合、积分、优化、解常微分方程等领域的函数。同时还提供了对F#语言的支持。支持导入matlab的.mat文件。支持CUDA和intel MKL.

缺点: 项目更新频率低，维护人员并不是很积极，几年前提的Pull requests和Issues都没有响应，其中某些代码的水平差得令人发指，比如位于src\Numerics\FindRoots.cs中的Quadratic函数，用于求解实系数的一元二次方程ax^2+bx+c=0，

public static (Complex, Complex) Quadratic(double c, double b, double a) { if (b == 0d) { var t = new Complex(-c/a, 0d).SquareRoot(); return (t, -t); } var q = b > 0d ? -0.5*(b + new Complex(b*b - 4*a*c, 0d).SquareRoot()) : -0.5*(b - new Complex(b*b - 4*a*c, 0d).SquareRoot()); return (q/a, c/q); } public static Complex SquareRoot(this Complex complex) { if (complex.IsRealNonNegative()) { return new Complex(Math.Sqrt(complex.Real), 0.0); } Complex result; var absReal = Math.Abs(complex.Real); var absImag = Math.Abs(complex.Imaginary); double w; if (absReal >= absImag) { var ratio = complex.Imaginary / complex.Real; w = Math.Sqrt(absReal) * Math.Sqrt(0.5 * (1.0 + Math.Sqrt(1.0 + (ratio * ratio)))); } else { var ratio = complex.Real / complex.Imaginary; w = Math.Sqrt(absImag) * Math.Sqrt(0.5 * (Math.Abs(ratio) + Math.Sqrt(1.0 + (ratio * ratio)))); } if (complex.Real >= 0.0) { result = new Complex(w, complex.Imaginary / (2.0 * w)); } else if (complex.Imaginary >= 0.0) { result = new Complex(absImag / (2.0 * w), w); } else { result = new Complex(absImag / (2.0 * w), -w); } return result; }

首先，没有处理 a=0 的情况，反而是考虑 b=0 这种并不需要单独考虑的情况；其次，根本没有必要调用对复数求平方根( Complex( , ).SquareRoot())这种计算代价很大的函数，Math.Sqrt(±Delta)足以解决问题；最后，没有必要使用代价较大的复数除法c/q. 我觉得哪怕是编程新手都很难写出这么蠢的代码。我提供了一个写法，在release版本中我的写法的速度是它的3.7倍，如果把我的写法中处理a=0的部分去掉，则速度是它的4.8倍。

public static (Complex, Complex) Quadratic(double c, double b, double a) { Complex x1, x2; if (a == 0.0) { if (b == 0.0) { x1 = Complex.Zero / Complex.Zero; // Complex.NaN;  x2 = x1; } else { x1 = new Complex(-c / b, 0.0); x2 = x1; } } else { a = 1.0 / a; b = -0.5 * b * a; c = c * a; double delta = b * b - c; double sqrtDelta; if (delta < 0.0) { sqrtDelta = Math.Sqrt(-delta); x1 = new Complex(b, sqrtDelta); x2 = new Complex(b, -sqrtDelta); } else { sqrtDelta = Math.Sqrt(delta); x1 = new Complex(b + sqrtDelta, 0.0); x2 = new Complex(b - sqrtDelta, 0.0); } } return (x1, x2); }

再举一个证明其水平低下的例子，src\Numerics\SpecialFunctions\Factorial.cs中的Binomial函数，作用是计算组合数C^n_k=n!/(k!*(n-k)!)=n*(n-1)···(n-k+1)/(k*(k-1)···2*1),

public static double Binomial(int n, int k) { if (k < 0 || n < 0 || k > n) { return 0.0; } return Math.Floor(0.5 + Math.Exp(FactorialLn(n) - FactorialLn(k) - FactorialLn(n - k))); }

这里面调用了3次FactorialLn()，就是对阶乘取自然对数，即ln(n!)，它的阶乘函数是用查表法实现的(有一个含171个double型元素的数组)，当n<171时，FactorialLn(n)的耗时约等于对数函数的耗时。n更大时会调用GammaLn函数，耗时更长。然后再调用指数函数和向下取整函数，这每一个函数的的耗时都不少。它这种写法在k比较大时才合适(我初步估计要大于100)，下面是我用循环实现的写法:

static double MyBinomial(int n, int k) { if (k < 0 || n < 0 || k > n) { return double.NaN; } else if (k == 0 || k == n) { return 1.0; } else { if (k > (n >> 1)) { k = n - k; } double dblN = n; double dblK = k; double Cnk = dblN / dblK; for (double i = 1.0; i < dblK; i++) { Cnk = Cnk * (dblN - i) / (dblK - i); } return Math.Round(Cnk); } }

当k比较小时，我的写法的速度是它的几十倍，浮点误差也远小于它。它如果把FactorialLn在n比较小时也用查表法实现，速度能提高，误差也能减少。所以，更合理的做法是，k比较小时，用循环或者查阶乘表的方法实现，k比较大时，用查阶乘的对数表(n比较小)与GammaLn函数(n比较大)，再结合指数函数Exp的方法。

DecimalMath的水平也不行，下面是它的Exp函数的实现，看一眼就知道没什么性能可言了。对于输入参数d，作者将其分成整数部分t和小数部分，竟然还是递归调用的，整数部分采用了快速幂算法(竟然不是查表)。不知道要把泰勒级数的计算进行循环展开，循环的终止条件竟然是if (nextAdd == 0)，也不知道要把除以常数的除法改成乘法。decimal型的计算本来就很慢了(我做过两个简单的测试，decimal型的耗时是double型耗时的60~80倍)，如果还不努力优化，基本上可以说慢得没法用。decimal型的Exp函数，允许的参数的最大值是 96*ln(2)=66.54，应该直接用一个含66个decimal型元素的数组保存Exp(1), Exp(2), Exp(3),……Exp(66)，占用存储空间66*16=1056 Byte，但是速度是快速幂的几百倍。

public static decimal Exp(decimal d)
{
	decimal result;
	decimal nextAdd;
	int iteration;
	bool reciprocal;
	decimal t;

	reciprocal = d < 0;
	d = Math.Abs(d);

	t = decimal.Truncate(d);

	if (d == 0)
	{
		result = 1;
	}
	else if (d == 1)
	{
		result = E;
	}
	else if (Math.Abs(d) > 1 && t != d)
	{
		// Split up into integer and fractional
		result = Exp(t) * Exp(d - t);
	}
	else if (d == t)   // Integer power
	{
		result = ExpBySquaring(E, d);
	}
	else                // Fractional power < 1
	{
		iteration = 0;
		nextAdd = 0;
		result = 0;

		while (true)
		{
			if (iteration == 0)
			{
				nextAdd = 1;               // == Pow(d, 0) / Factorial(0) == 1 / 1 == 1
			}
			else
			{
				nextAdd *= d / iteration;  // == Pow(d, iteration) / Factorial(iteration)
			}

			if (nextAdd == 0) break;

			result += nextAdd;

			iteration += 1;
		}
	}

	if (reciprocal) result = 1 / result;
	return result;
}

下面是DecimalMath的阶乘函数的实现，根本就不应该用循环，应该直接查表实现，因为 28!=3.049×1029>296−1=7.923×1028 =decimal.MaxValue，只需要存28个decimal值就可以了。用循环就慢到可悲了。

public static decimal Factorial(decimal n) { if (n < 0) throw new ArgumentException("Values less than zero are not supoprted!", "n"); if (Decimal.Truncate(n) != n) throw new ArgumentException("Fractional values are not supoprted!", "n"); var ret = 1m; for (var i = n; i >= 2; i += -1) { ret *= i; } return ret; }