数值近似初等函数

原创
飞起转圈
2023-08-20
703
C++

Image from X (formerly Twitter) And I really don't know who drew this.

数值近似初等函数

引入

上回说到，我在写一个小小的玩具3D数学库。
既然是数学库，那就一定得有初等函数的计算啊\ (￣︶￣*\ ))

......好吧其实我一开始没想做这方面的。
但是当我复习到了泰勒展开的时候，我在想.....如果我把函数泰勒展开了，再配上simd指令，是不是就能比标准库快了捏.....
然后我试了试，发现不对--误差太大啦 {{{(>_<)}}}
然后就生气了，我要和初等函数决一死战╰（‵□′）╯啊啊啊啊啊啊

什么是初等函数

反函数，对数函数，幂函数，指数函数，三角函数

为什么要自己造轮子

在任何情况下都不该造这种轮子（除非吃饱了撑的
我们有可靠的工具来计算这些函数，比如Glibc
甚至芯片厂商都在设计计算初等函数的指令，fyl2x,fcos,fsin,fsqrt.....

under x86

当然，我也确实有理由说明为什么要自己造轮子，理由如下:

  1. 吃饱了撑的
  2. 标准库是性能与精度的平衡，但我可以放弃精度要速度
  3. 硬件实现的浮点运算在不同厂商下结果不同，而且精度不行

不管怎么说，这篇博客的标题都太大了，我必须声明一下，我只是在记录自己的最佳实践。

数值近似这个话题也是个高深数学的分支，至于对初等函数的计算，有一本书叫 “Elementary Functions: Algorithms and Implementation - Second Edition”，讲了很多相关的算法，有兴趣可以去找来看看。

误差的计算

对于误差有两个硬指标，绝对误差和ULP
绝对误差很好理解，就是$|近似的数值-真正的数值|$。
那么ULP是什么呢？Unit In The Last Place
这个概念十分抽象，我们都知道浮点数在电脑中无法精确的表达。事实上在现实生活中很多浮点数也无法被精确表达，只能靠一些约定俗成的符号来记录，比如说$\pi,\frac{1}{3}$。
电脑是用二进制来储存浮点数的，其浮点数是用2的不同幂次按位权加和得到。假如无法用一定的的位数组合出这个浮点数，那么计算机就会进行截断，ULP就是进行截断的最低位有效数字。
更精准的定义是：
$ulp(x)$ is the distance between the two closest straddling floating-point numbers，assuming that the exponent range is not upper-bounded.

好吧确实太抽象了，举个小小的例子
你上幼儿园了老师教了你除法，你会做$6\div3=2$了，但是有个上了五年级的坏批给你出了道题叫$2\div3=?$，鉴于你浅薄的知识你只能写成这样：$1\div3=0.....2$，于是你把一个在0-1之间的小数进行了截断，此时ULP就是1.
当你也上了五年级你终于会做小数除法了，你找到了尘封的那道题并且进行了计算：$2\div3=0.667$，很显然你不想写下去了，进行了四舍五入（当然，也是截断的方法），于是此时的ULP就是0.001

那么在电脑中如何计算ULP呢，这需要对于浮点数格式的了解了...

想快速了解我建议去找王道408计算机组成原理网课第三章

#include <cmath>   
    #include <math.h> 

    double ulp(double x)
    {
      double ulp;                /* ulp of x */
      int    exp;                /* exponent of x, where exp = E+1 */

      frexp(x, &exp);            /* extract exponent of x */
      ulp = ldexp(0.5, exp-52);  /* calculate ulp = 0.5*2^(exp-52) */

      return ulp                 /* return ulp */
    }

frexp，用来提取浮点数的指数和尾数。但神奇的是，在IEEE754的约定里，尾数应该在$[1,2)$中，但是frexp返回值是在$[0.5,1)$中，因此，其返回值是真正的指数+1。
于是ULP也就是$0.5*2^{exp-52}$，其中52是浮点数尾数的位数。小于这个数的任何数都会被截断。

这里我写的比较详细，因为资料太少了

对于之前提到的标准库与硬件实现的方法，其都有ULP和绝对误差的规定来保证结果的正确性。

在我的误差测试中，我懒得算ulp了，只对最大绝对误差上限进行了计算。

几种简单的近似思路

泰勒展开

这种方法的缺点很明显，只有在展开处附近的精度可以，一旦离开了一点，精度就会掉掉掉....

投影到基函数

这种方法比泰勒好在：其精度在整个函数域上都挺高的(其实也高不到哪去)
但坏也坏在这，泰勒可以通过在某点不断展开来取得更高的精度；基函数方法当然也可以投影到更高维度，但在固定点的收敛速度明显不及泰勒展开。
常用勒让德多项式基和切比雪夫多项式基

查表

就是放张表嗯查。
缺点嘛..没有周期性的函数你得准备张多大的表呢？

其他神奇的多项式逼近方法

依问题而论，总之非常神奇

三角函数的近似

Sin具有周期性，完全可以用诱导公式把它化到$[0,\frac{\pi}{4}]$内，然后泰勒展开来计算。
在我一开始实现的某版本中，我是把sin转化到了$[0,\pi]$中，然后直接开始查一个1k元素的表。
这个版本的实现我还留着，最大绝对误差上界大概在$10^{-7}$次方左右，速度在MinGW-w64 O3下比std::sin快一丢丢，到了linux G++下就完全没优势了。

那么投影到基函数上呢？我尝试过把Sin投影到勒让德多项式上，结果是投影到7阶的勒让德上其误差上界在$10^-4$次方左右。完全不够看

对于Sin，Cos这种有周期性的函数投影是很好计算的，会存在奇数项系数为0，或者偶数项系数为0的情况，完全可以多计算几项。但是鉴于这可怜的收敛速度.....比泰勒好点吧

我个人的最佳实践是利用Sin(a+b) = SinACosB+SinBCosA来计算。
比如说,如果计算Sin(X)，X=1.234，那么可以把X分开位1和0.234，其中Sin1和Cos1查表获得，Sin0.234和Cos0.234通过泰勒展开获得。
这个函数的实现在这里，其在Linux g++12.1下能比标准库的Sin快三倍，误差上界在$10^{-13}$左右。
对于Cos也可以以同样的思想得到
至于Tan，Sin/Cos就行啦

当然，如果是想了解通用的，高精度的三角函数求法，可以去了解下Cordric算法。

反三角函数

我只实现了Arctan，但是Arcsin和Arccos都能用其计算出来.......好吧我有点懒得查资料了
我这有两篇文章，首先是 - New Fast Arctangent Approximation Algorithm for Generic Real-Time Embedded Applications
这篇文章中提到的精度最高的算法，其最大误差上界大概是$10^{-2}$左右，太拉跨了。

另一篇文章A NEW METHOD FOR OBTAINING HIGHLY ACCURATE APPROXIMATIONS TO THE ARCTAN FUNCTION 则好用的多。 其思想简单来说就是对于计算$arctan(\frac{1}{a})$，就拿偶数项勒让德多项式除以$(a^2+x^2)$并对其在0-1上做积分再乘以a，得到一个$a\int^1_0\frac{P^{2n}}{a^2+x^2}dx$的形式。
然后利用平方差公式把这个式子化为两部分，第一部分是关于a的多项式，第二部分是a的多项式乘以$\frac{1}{a}*arctan(\frac{1}{a})$
注意到当勒让德多项式的次数足够高，a足够大，那么这个式子会趋向于0。移项，得到arctan(1/a)
以4阶勒让德为例

$$ \int^1_0\frac{P_{4}}{a^2+x^2}dx $$ $$ \Longrightarrow\frac{1}{8}\int^1_0\frac{35x^4-30x^2+3}{a^2+x^2}dx $$
$$ \Longrightarrow\frac{1}{8}\int^1_0\frac{35x^4-35a^4+35a^4-30x^2+30a^2-30a^2+3}{a^2+x^2}dx $$

$$ \Longrightarrow-\frac{5}{24}(11+21a^2)+\frac{1}{8}(3+30a^2+35a^4)*\frac{1}{a}*arctan(\frac{1}{a}) $$
令等式=0则有
$$ \Longrightarrow\frac{5}{3}\frac{a(11+21a^2)}{(3+30a^2+35a^4)} = \frac{1}{a}arctan\frac{1}{a} $$
又因为$arctan(\frac{1}{a})+arctan(a)=\frac{\pi}{2}$
那么就有arctan(a)力
至于他们咋想到的，我也不知道。毕竟原文真写了个 “The idea is to recognize that....”
可能是我注意力不集中吧 QWQ

实测这种方法再MinGW-w64下是std::atan速度一倍多，linux下快了0.5倍，误差上界是$10^{-6}$左右。

对数函数

对数函数的实现乍一看没有头绪，事实上可以对Ln函数用换底公式得到。
那么问题来了...ln(x)怎么实现...
通过点简单的高中数学我们可知
$ln(x)=ln(2^k*u)=kln(2)+ln(u)$
在前文我说过，计算机中浮点数的储存事实上是尾数和指数，那么只要想办法把指数，尾数提取出来。其中尾数的范围是$[1,2)$
对于double，可以通过位运算达到这点

FM_FORCE_INLINE int64_t FM_CALL mfrexp(double t, double *rem) {
  union {
    double s;
    uint64_t i;
  };
  s = t;
  uint64_t _rem = i;
  _rem = _rem & 0x000fffffffffffff; // set exponent to 0
  _rem = _rem | 0x3FF0000000000000; // set exponent to 1023( - 1023 = 0 as its real value)
  uint64_t exp = i;
  exp = ((exp & 0x7fffffffffffffff) >> 52) - 1023;
  *rem = *(double *)&_rem;
  return exp;
}

接着只要预计算$ln2$就可以了
对于项$ln(u)$该如何处理呢？一个Naive的想法是把其化为$ln(1+k)$的形式然后再0.5处泰勒展开
然而我试过了，精度很差....
然后我找到了一个很天才的想法，在一个国人的开源库Morn中实现，它把u化为了n*m，其中n是形如1.xx的形式，那么对于m，其必然为1于一个极小小数的和，这样我们可以对ln(n)进行查表，对ln(m)进行泰勒展开。
我的实现在这里
实测，其速度在MinGW-w64下是std::log的6倍快，而在Linux下比std::log慢了一倍(?)

我估计是GCC直接调用指令计算ln了

但总而言之，赢！

我的log误差大概在$10^{-12}$次方，但是加一个前提，特定情况下....

浮点数还有非规格数的特例，非规格数大概是$2^{-1023-53}$大，但我又懒得写特例.....所以对于小于最大非规格数的一切数字我是直接返回-inf的
(。・ω・。)

指数函数和幂函数

先来实现指数函数exp吧
对于exp(x)
令$k = \lfloor\frac{x}{ln2}\rfloor$
再令$x-kln2=r$
其中$0< r<ln2$
那么$e^x = e^r2^k$
又可以把r拆为a+b则有
那么$e^x = e^ae^b2^k$ 其中a可以查表，b可以泰勒
这种实现在MinGW下比标准库快一倍，但是在Linux下和标准库速度相当
至于精度..可以保证前10位数字的精准


对于幂函数可以用$x^a=e^{a*lnx}$得到
这里由于标准库和我的代码都是如此实现，所以根据ln和exp的速度就能猜到这个的速度了吧.....

后记

有没有发现我一直没有提到MSVC上的速度？
因为MSVC对我的代码负优化严重
MinGW的标准库在win上用的是和MSVC同款的，所以两边的数学函数速度一摸一样
但是MSVC编译出来的代码，比MinGW编译出来的慢了4-5倍不止

而在WSL的G++下，我的代码运行速度比win+MinGW下还快，标准库也比win上的快

谁是答辩，我不好评价




还是放松一下吧

END