“ULP”代表“单位在最后的位置。”一个就ulp进行阴谋采样一个基本的数学函数，如$\sin{x}$，或一个更深奥的函数，如贝塞尔函数。将样品与从更高精度计算得到的更精确的值进行比较。以ulps测量的精度图揭示了关于底层算法的有价值的信息。

内容

fdlibm

libm是支持C编译器的基本数学函数库。万博1manbetxfdlibm“可自由分发”的来源是什么libm在25年前由K. C. Ng开发并投入公共领域，可能还有Sun Microsystems的其他一些人。我写了fdlibm在我们2002年的通讯中，四角函数与数值工艺．

在数学上fdlibm展示了完美的工艺。我们今天在初等超越函数中仍然使用它。我怀疑所有其他数学软件项目也是如此。如果他们没有，他们应该这么做。

就ulp进行情节

就ulp进行(x)是距离x到下一个更大的浮点数。它和每股收益(x)．

评估计算值的准确性

y = f (x)

将它与从符号数学工具箱中获得的更精确的值进行比较

Y = f(信谊(x,“f”））

的“f”标志说转换x到一个信谊没错，不用去猜它是10的逆次方√6的一些东西。相对误差y，用最后一处的单位来度量，则为

u = (y - y)/eps(abs(y))

因为这是相对误差，它是接近零点的严格测量f (x)．

如果y浮点数是否通过正确四舍五入得到Y那么，要加倍精度

-0.5 <= u <= 0.5

这是我们所能期望的最好结果。的精确的数学值f (x)并做一个舍入误差得到最终结果。半ulp精度在算法上难以获得，而且执行时间太长。MATLAB中的所有函数都是由fdlibm精度优于一ulp。

-1.0 < u < 1.0

下面的每一个图都包含100,000个随机参数x，均匀分布在适当的区间内。

罪

我们从这个样本中看到大约0.8 ulp的准确性。这是典型的。

参数减少是计算$\sin{x}$的第一步。从参数中减去$\pi/2$的整数倍$n$，将其带入区间

$ $ - \压裂{\π}{4}\ le x - n \压裂{\π}{2}\ le \压裂{\π}{4}$ $

然后，根据$n$是奇数还是偶数，对$ sin$或$ cos$的13次多项式近似给出了简化后的论证的近似四舍五入的结果，也因此给出了原始论证的近似四舍五入的结果。ulps图显示了在$\pi/4$的奇数倍处的精度的轻微下降，这是多项式逼近的极端点。

值得注意的是，即使在样本区间的端点$0$和$2\pi$附近，精度也比一个ulp更好。这就是$\sin{x}$趋近$0$的地方，并且必须仔细地逼近，以便相对误差保持有界。

棕褐色

同样，大约0.8 ulp精度。

类似的参数减少会导致类似的行为，接近奇数倍的$\pi/4$。在这两个点之间，在$\pi/2$和$3\pi/2$处，$\tan{x}$有一个极点，近似也必须如此。该算法使用倒数和恒等式

$ tan x = -1/ tan{(x+ frac{\pi}{2})} $

当你接近极点时，这接近于除以零，但是得到的近似仍然比一个ulp要好。

:

好到略高于0.8 ulp。基本的近似是一个分段多项式，断点为1/16的几个倍数，在图中很明显，并标记在坐标轴上。

经验值

精确度和之前的图差不多。

参数缩减涉及键值

$ r = ln{2} \约0.6931

和身份

$ exp{(x-n r)} = 2^n {exp{(x-n r)

结果ulps图显示了误差在奇数倍$r/2$处的极值。

兰伯特W

现在来看看两个不在的函数fdlibm．如果你关注这个博客，你可能已经注意到我是Lambert W函数的忠实粉丝。几年前我在博客上写过，朗伯W函数．维基百科的文章很好，兰伯特W函数．你只需要谷歌“朗伯w函数”就可以找到更多有趣的链接。

Lambert W函数在MATLAB本身中是不可用的。符号数学工具箱有一个@double方法，用于访问类型double参数的符号代码。或者，我上面提到的博客文章中有一些简单(但很优雅，如果我自己这么说的话)的代码。

这两个代码都不是一个ulp精确。函数的一次分支在原点处是零。当我们接近0时，测量的相对误差就ulp进行是无限的。绝对精度是可以的，但相对精度就不行。事实上，您可能会看到十亿ulps的错误。这是一个gigaulp,或狼吞虎咽地吃为短。

作为lambertw (x)经过2的幂，最后一位的单位，每股收益(lambertw (x))是美国的两倍。其中三个点在ulps图的x轴上用蜱标记。

为= (1/8 1/4 1/2) z = fzero (@ (x) lambertw (x) - a, 5)结束

Z = 0.1416 Z = 0.3210 Z = 0.8244

贝塞尔函数

我们的贝塞尔函数代码有相当好的绝对精度，虽然不是一个ulp。但它们在零点附近的相对精度不高。这里是第一类零阶贝塞尔函数$J_0(x)$的前5个零，以及覆盖前2个零的ulps图。

为A = (1:5)* z = f0 (@(x) besselj(0,x)， A)结束

Z = 2.4048 Z = 5.5201 Z = 8.6537 Z = 11.7915 Z = 14.9309

erfinv

这是误差函数的逆函数。看起来很有趣，但我还没有调查过。

代码

我最近更新了克里夫的实验室在MATLAB中心进行文件交换。最新版本包括ulpsapp.m，它在本博客中生成ulps图。

获取MATLAB代码

发布与MATLAB®R2017a