关于math数学库

MathWare是一个开源软件

GitHub - DeepMathLLM/DeepMath: 一个开源数学大模型项目,旨在探索大模型是否具有数学创造能力,以及大模型在前沿数学研究中的潜在能力。

GitHub - Kaida-Amethyst/math.mbt  Moonbit 数学库是一个在 Moonbit 编程语言中实现的数学函数集合。该库旨在提供高精度的数学运算,涵盖了三角函数、指数函数、对数函数和特殊函数等,

 

Math.net集合图形库:开源数学库

MathFlow 是什么？开源数学计算库
git clone https://github.com/mathflow/mathflow.git cd mathflow

Math.js:一款 JavaScript 和 Node.js 最强数学库

Math Utility开源数学工具库

MathFu是一款由Google开发的C++数学库，专为游戏开发设计，强调简洁和效率。它的核心目标是提供一套矢量、矩阵和四元数类，以支持基础几何运算，适用于图形库如OpenGL或动画、物理系统的计算需求。MathFu库采用了便携式C++编写

[Android]: http://www.android.com
[Linux]: http://en.m.wikipedia.org/wiki/Linux
[MathFu Google Group]: http://groups.google.com/group/mathfulib
[MathFu Issues Tracker]: http://github.com/google/mathfu/issues
[OS X]: http://www.apple.com/osx/
[OpenGL]: http://www.opengl.org/
[SIMD]: http://en.wikipedia.org/wiki/SIMD
[Windows]: http://windows.microsoft.com/
[geometry]: http://en.wikipedia.org/wiki/Geometry
[landing page]: http://google.github.io/mathfu
[matrix]: http://en.wikipedia.org/wiki/Matrix_(mathematics)
[quaternion]: http://en.wikipedia.org/wiki/Quaternion
[stackoverflow.com]: http://stackoverflow.com/search?q=mathfu
[vector]: http://en.wikipedia.org/wiki/Euclidean_vector
[CONTRIBUTING]: http://github.com/google/mathfu/blob/master/CONTRIBUTING

 

以下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状态，不再更新；

 

https://github.com/FrogGuaGuaGua

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

本文Github地址：https://github.com/FrogGuaGuaGua/CSharpCodePrinciple/blob/master/CSharp-MathOptimization.md

华为公司的C语言编程规范在开头就强调了：

一般情况下，代码的可阅读性高于性能，只有确定性能是瓶颈时，才应该主动优化。

本文讲述的方法没有经过大项目和大公司的检验，所以，请批判性地阅读本文。本文中的大部分结论有测试代码支持，参见SpeedTest.cs. 虽然C#的编译器会在release版本中执行一些优化，C#的运行时也有一些优化，但偶尔会遇到debug版本正常，而release版本异常的问题，比如我在github上fork了已停止维护的屏幕录像软件Captura，debug模式能启动，release版本无法启动，我短时间内解决不了这个问题，如果要发布，只能发布debug版本。所以手工执行一些虽然编译器(在release版本中)会做但也简单易读的优化，还是有意义的。同时建议把未经优化的代码作为注释附在旁边，提高可读性。

1). const, readonly, in 这3个关键词在能用的地方要尽量用。这样可以让编译器执行更激进的优化策略，同时提高代码的安全性、可读性和可维护性。

2). 正整数乘以或除以 2𝑛 (n也是整数)，使用移位来代替。但对乘以非 2𝑛 的整数不要使用此方法，比如不要把 i * 12 改写成 (i << 2 + i << 3). 对于负整数的乘除运算，也不要用移位，否则代码可读性很差，容易出错。 正整数 x 对 2𝑛 求余数，可以直接提取x的低n位，即 x % 2𝑛 = x & ( 2𝑛 -1). 特别地，判断正整数 x 的奇偶性，可以用 if(x & 1 == 0) 来实现。对非 2𝑛 求余数的优化，参见Faster Remainder by Direct Computation. 对于int型整数，除法耗时是乘法的13倍，是移位的17倍。对于long型整数，除法耗时是乘法的1.4倍，是移位的9.8倍。(数据来源)

3). 除以浮点型常数的除法，改写为乘以这个数的倒数。比如把 x / 2.0 改写为 x * 0.5. 对于条件语句if(a/b > c)，可以改写为if( (b > 0 && a > b * c) || (b < 0 && a < b * c) ). 对于某些有很多分数嵌套的数学公式，请先进行恒等变形，只保留最长的一条分数线，其它的分数一律去掉分母，除法变成乘法。例如：

(1𝑎2+𝑘2𝑏2)𝑥2+2𝑘𝑚𝑏2𝑥+𝑚2𝑏2−1=0𝑥1+𝑥2=−2𝑘𝑚𝑏21𝑎2+𝑘2𝑏2=−2𝑎2𝑘𝑚𝑏2+𝑎2𝑘2𝑥1𝑥2=𝑚2𝑏2−11𝑎2+𝑘2𝑏2=𝑚2𝑎2−𝑎2𝑏2𝑏2+𝑎2𝑘2=(𝑚2−𝑏2)𝑎2𝑏2+𝑎2𝑘2
因为浮点除法的耗时是浮点乘法耗时的13倍。(数据来源)

4). 除非测试结果表明使用float型比double型更快，或者因为数据量巨大，float型比double型显著节省空间，否则都应该使用double型浮点数。因为float型的计算速度有时比double更慢，而且float型的精度太差，可能带来一些难以察觉的bug，比如 for(float i = 0.0f; i < 20000000; i++) 就是一个难以察觉的死循环，当 i 达到 224 =16777216 时，会由于float的精度太低，无法区分 16777216 和 16777217，即无法自增1. 另外，使用float型，所有的字面量都要加 f 后缀，所有的Math库函数前面都要加(float)进行强制类型转换(或者使用MathF库中的函数)，写起来麻烦，看起来丑陋，所以要尽量避免。在以下情形(但不限于)可考虑使用float型：
> 训练AI模型，数据量巨大但对计算精度要求不高，float型可显著节省存储空间。
> 程序要在32位处理器上运行，或者要在没有硬件浮点单元的处理器上运行。
> 需要使用SIMD技术加速，使用float型可以在一条指令中处理两倍于double型的数据。

5). 如需计算 𝑥𝑛 ，当 n=2,3 时，不要使用Math.Pow(x,n), 而是直接写成 x * x 和 x * x * x. 当 𝑛=2𝑘(𝑘∈𝑁+) 时，可用y=x, 再执行k次 y*=y 来代替。当 n 取其它值时才可调用Math.Pow.

6). 引入一些额外的变量来存储函数调用的结果，或者复杂运算过程中的子过程的值，避免重复调用和计算。比如计算二维坐标旋转:

    x1 = x * cos(a) - y * sin(a) 
    y1 = x * sin(a) + y * cos(a)

同一个角的正弦和余弦值都要使用两次。一元二次方程求根， Δ2𝑎 会使用两次。二元一次方程求根，系数矩阵的行列式值会使用两次。在循环中如果要以同样的参数调用某个函数，或者有一些不随循环变化的子过程，则应提到循环外部，用变量存储。

7). 对浮点数进行取整操作时，如果确定浮点数的大小不超出int(或long)型的范围，以及不会出现NaN，则可以用强制类型转换结合条件语句替代Floor、Ceiling和Round函数，显著提高速度。使用 (int)(t ± 0.5) 来代替Math.Round(t)则需谨慎，因为当t的小数部分为0.5时，Round(t)的结果取决于中点值舍入模式的设定，默认是MidpointRounding.ToEven，即向最近的偶数舍入。其它模式还有ToZero, AwayFromZero, ToNegativeInfinity, ToPositiveInfinity. 要根据不同的舍入模式选择不同的替代写法。

    // 替代 a = (int)Math.Floor(t)
    a = (t < 0 ? (int)t - 1 : (int)t);

    // 替代 b = (int)Math.Ceiling(t)
    b = (t < 0 ? (int)t : (int)t + 1);

    // 替代 c = (int)Math.Round(t, MidpointRounding.ToZero)
    c = (t < 0 ? (int)(t - 0.5) : (int)(t + 0.5));

8). 对于Array of Struct(AoS)和Struct of Array(SoA)两种数据结构，
内存布局：
AoS：每个结构体实例的所有字段在内存中是连续存储的。
SoA：每个字段的所有值在内存中是连续存储的，但不同字段的值分开存储。
性能：
AoS：在需要频繁访问单个结构体实例的所有字段时性能较好。
SoA：在需要频繁访问所有实例的单个字段时性能较好，特别是在SIMD优化中表现更佳。
对于有限元程序，需要存储大量节点的编号、坐标和位移，需要视情况选择AoS或SoA.

9). 对于较小的结构体，可以考虑用ref struct代替struct，强制结构体存储在栈上(注意防范栈溢出)，避免装箱操作，同时减少垃圾回收的性能损失。

10). 对于局部变量，使用 Span<T>, ReadOnlySpan<T> 和 stackalloc 在栈上分配连续的小段内存(注意防范栈溢出)，比使用数组(存储在堆上)速度更快。

11). 对于分支较多的流程，优先使用模式匹配而不是大量的if else，既能提高程序可读性，又能提高运行速度。比如分段函数就应该使用模式匹配。

static double Foo(double x) => x switch
    {
        < 0 => -x,                // 当 x < 0 时，f(x) = -x
        >= 0 and <= 1 => x * x,   // 当 0 <= x <= 1 时，f(x) = x^2
        > 1 and <= 2 => 2 * x,    // 当 1 < x <= 2 时，f(x) = 2x
        > 2 => x + 1,             // 当 x > 2 时，f(x) = x + 1
        _ => throw new ArgumentOutOfRangeException(nameof(x), "Invalid input")
    };

12). 尽量避免编写含递归调用的函数。比如阶乘函数 n!，递推数列(斐波那契数列、汉诺塔问题等)，二分查找等，均可以用循环替代递归。

13). 对于那些参数的允许范围比较小的函数，优先考虑用查表法实现。比如阶乘函数 n!，因为阶乘函数增长太快，在大多数情况下，阶乘函数允许的参数的范围很小，
13!=6227020800>232−1=4294967295 = uint.MaxValue
21!=5.109×1019>264−1=1.845×1019 = ulong.MaxValue
28!=3.049×1029>296−1=7.923×1028 = decimal.MaxValue
35!=1.033×1040>2128=3.403×1038 = float.MaxValue
171!=1.241×10309>21024=1.798×10308 = double.MaxValue
至多占用171*8 = 1368Byte的存储空间，就能满足double型计算的需求。不仅速度快，而且没有多次浮点乘法带来的累积误差。
特殊情况下，指数函数的自变量如果只能取正整数，那么自变量的范围一般也不会很大，比如 𝑒709<21024<𝑒710 ，那么可以考虑对不超过某一阈值的整数采用查表法，超过该阈值则调用标准库。或者为自变量取等差数列时的函数值建立数表，然后用少量运算就能得到0~709内任意整数的函数值(参见 https://zhuanlan.zhihu.com/p/5221342896 )。
二项式系数(组合数)和阶乘的自然对数ln(n!)也可以采用部分查表法。
CRC(循环冗余校验)也经常使用查表法加速。

14). 小于255的素数(质数)一共有54个，如下：

static readonly byte[] PrimesLessThan255 = [2, 3, 5, 7, 11, 13, 17, 19, 23,
        29,  31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 
        103, 107, 109, 113, 127, 131, 137, 139, 149, 151,157, 163, 167, 173, 
        179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251];

对正整数n进行素性测试时，可以先用上表的素数进行试除(比使用255以内的奇数试除要快)，若都不能整除，就可以继续用从257到 𝑛 之间的奇数进行试除，从257开始是因为253(=11*23)和255都是合数。试除法是最简单但并不高效的素性测试方法，比较高效的方法是 Miller-Rabin测试法。但因为绝大多数正整数都有一个不大的素因子，比如：88％的正整数有一个小于100的素因子，92％的正整数有一个小于1000的素因子(数据来源： 大数因子分解算法综述.刘新星)。大于255但不超过1024的素数有118个。同时，根据素数定理，不超过n的正整数中，素数占的比例大致是1/ln(n). 因此，构建一张比较大的素数表，采用先除素数再除奇数的试除法，对于不太大的整数，是一种勉强能用的素性测试方法，同时也是寻找素因子的方法。

15). 以e为底的指数函数有一种快速近似算法：

public static double FastExp(double x) {  
        long tmp = (long) (1512775 * x + 1072632447);  
        return BitConverter.Int64BitsToDouble(tmp << 32);  
    }

该方法的速度大致是Math.Exp的5倍，原理参见《 A Fast, Compact Approximation of the Exponential Function》. 对于神经网络中的Sigmoid函数中的指数函数，就可以采用这种近似算法。

16). 以e为底的对数函数有一种快速近似算法：

public static double FastLn(double x) // 抛弃对x<=0的检查。
    {
        long longx = BitConverter.DoubleToInt64Bits(x);
        double k = (longx >> 52) - 1022.5; 
        return k * 0.693147180559945309;  
    }

该方法实际上就是Math.Log的算法的前半部分，用位运算提取了IEEE 754浮点数的阶码，而抛弃了尾数的对数，速度大致是Math.Log的4倍，其中，-1022.5 = - 1023 + 0.5，0.693147……就是ln(2)，该算法可以保证绝对误差不超过ln(2)/2=0.346573... 但该算法有一个不可忽视的弊端：设 n 为正整数，则对于区间 [2𝑛−1,2𝑛) 内的任意实数，该算法会返回完全一样的结果。以2为底或以10为底的对数函数也可以使用该方法，把最后一行与k相乘的常数换掉即可，以2为底就是return k，以10为底就是return k*0.301029995663981196.

17). 求正实数的平方根的倒数有快速算法，原理参见平方根倒数快速算法。速度测试参见平方根倒数快速近似算法在现代CPU上还有必要吗？.

float InvSqrt(float x) // 比MathF.ReciprocalSqrtEstimate()更慢，不值得使用
{
 float xHalf = 0.5f * x;
 int i = BitConverter.SingleToInt32Bits(x);
 i = 0x5f3759df - (i >> 1);
 x = BitConverter.Int32BitsToSingle(i);
 x = x * (1.5f - xHalf * x * x);
    //x = x * (1.5f - xHalf * x * x); //多迭代一次可以提高精度
 return x;
}

double InvSqrt(double x) // 大多数个人电脑的CPU不支持AVX512指令，该算法比Math.ReciprocalSqrtEstimate()更快，值得使用
{
 double xHalf = 0.5 * x;
 long i = BitConverter.DoubleToInt64Bits(x);
 i = 0x5fe6ec85e7de30daL - (i >> 1);
 x = BitConverter.Int64BitsToDouble(i);
 x = x * (1.5 - xHalf * x * x);
    //x = x * (1.5 - xHalf * x * x); //多迭代一次可以提高精度
 return x;
}

18). 避免在循环中做以下事情：
> 创建对象。
> 使用try catch.
> 打开和关闭同一个文件、数据库等。
> 创建和断开对同一个URI的链接。

19). 避免不加测试地用Parallel.For代替for循环，因为前者需要创建和管理多个线程，会带来额外的开销。当循环次数太少或者单次循环所做的运算太简单时，使用Parallel.For反而会降低性能，而且很可能出现计算结果不正确的问题。比如函数f(x)在某个区间上做数值积分，有sum+=f(xi)*dx这样的累加运算，需要测试Parallel.For的耗时是否更短以及结果是否正确。又比如一些写不出通项公式的递推数列， 𝑎𝑛+1=𝑎𝑛+𝑎𝑛+1 ，每次循环都依赖上一次的结果，注定无法并行，此时禁止使用Parallel.For.

20). 考虑使用[SkipLocalsInit]属性，省略CLR将方法中声明的所有局部变量初始化为其默认值的操作，提高速度。注意：此属性需要 AllowUnsafeBlocks 编译器选项，同时要重点检查代码中是否存在访问未初始化的变量的行为。

21). 绝大多数时候，矩阵求逆都是非必须的(而且计算代价很大的)，除非就是要得到逆矩阵本身。比如对于线性方程组 𝐴𝑥=𝑏, 𝑥=𝐴−1𝑏 ，使用逆矩阵表达方程组的解只具有形式意义，不可直接用于计算。请使用高斯消去法、LU分解法、Jacobi迭代法、Gauss-Seidel迭代法、Cholesky分解法等。矩阵的逆几乎不会单独出现，几乎总是会和其它矩阵做乘法，总有不求逆的替代方案。

22). 多项式求值优先使用秦九韶算法， 𝑎𝑛𝑥𝑛+𝑎𝑛−1𝑥𝑛−1+⋯+𝑎1𝑥+𝑎0=(⋯((𝑎𝑛𝑥+𝑎𝑛−1)𝑥+𝑎𝑛−2)𝑥+⋯+𝑎1)𝑥+𝑎0
对于阶数不太高的多项式(比如小于10阶)，不要使用循环语句来实现这个算法，而应该手工进行循环展开。还可以使用融合乘加指令Fma.MultiplyAdd进行进一步加速。秦九韶算法是一个串行的算法，无法并行。如果某个n次多项式的全部根均为实数(设为 𝑥1,𝑥2,⋯,𝑥𝑛 ，需要提前计算出来)，那就可以使用SIMD指令进行并行计算： 𝑎𝑛(𝑥−𝑥1)(𝑥−𝑥2)⋯(𝑥−𝑥𝑛)

23). 利用泰勒级数计算double型函数值时，多项式阶数通常不应该超过17阶，太高的阶数没有意义(因为浮点运算的累积误差)。泰勒级数具有局部性，离展开点越远，精度越差。所以如果要提高计算精度，首先应考虑更换展开点，而不是提高多项式的阶数。

24). 除非绝对必要(比如要求很高的精度或比long型更大的范围)，否则不要使用decimal型变量，因为其计算耗时是double的几十倍甚至百倍(decimal除法尤其缓慢)。

25). 免费的数学库推荐ALGLIB免费版，收费的数学库推荐ALGLIB、ILNumerics和Dew.Math. 不推荐 MathNET Numerics，其代码质量低下，原因参见点评10多个C#的数学库.

参考文章：

• Writing Faster Managed Code: Know What Things Cost
• 新版C#高效率编程指南
• C#中那些举手之劳的性能优化