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

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;
}

 

编辑于 2025-01-03 21:29・IP 属地广东