学习LabVIEW（十）——关于Matlab的eps函数（十三）

关于Matlab的eps函数的实现原理，我以前写过十几篇短文了。其中最接近官方的版本可以参考《关于Matlab的eps函数（十）——MATLAB Coder生成的C代码》。那时使用了MATLAB的Coder直接生成了一份eps的C语言版本，应该就是官方的实现方式了。
为什么这么关注eps函数？因为对于数值分析而言，eps是非常重要的，它展示了基于IEEE754浮点数的一个基本属性：相对浮点精度。然而，尽管eps如此重要，能够说清楚eps运算方法的人寥寥无几。现代社会，人人都要会编程，然而如果只是接受了普通的编程入门学习，没有像计算机专业科班出身的那群人那样花了好多好多的学时学习数值分析的话，就会倾向于将计算机中的浮点数与数学中的实数这个概念等价起来。可事实却是：浮点数和数学概念中的实数有着巨大的差异。比如，概念中的实数，是可以连续取值的。而基于754标准的浮点数，则和整数类似，取值是离散的！而eps正是衡量一个浮点数取值到下一个浮点数取值之间的距离的工具。与整数不同的是，浮点数之间的间隔不是一成不变，而是随着浮点数数值的变化而变化的；所以eps不是一个常数而是一个函数。

• 在《关于Matlab的eps函数》中，我们讨论了eps的本质，并使用Matlab的typecast函数（用来实现reinterpret cast）实现了eps的计算；
• 在《关于Matlab的eps函数（续）》，《关于Matlab的eps函数（又续）》，《关于Matlab的eps函数（六）》和《关于Matlab的eps函数（七）——“又续”的续》中，我们利用位段结构体实现了eps计算，这种版本的可读性比最初用typecast实现的版本要好多了；
• 在《关于Matlab的eps函数（再续）》中，我们根据mathworks上的一个帖子，讨论了利用纯粹的数学运算而不是位运算实现eps的方法；
• 在《关于Matlab的eps函数（五）》和《Article 4 in 1: 关于Matlab的eps函数（八） &读Matlab7.7的rank函数 &读Matlab7.7的orth》中，讨论了eps函数的几个实际应用；
• 在《关于Matlab的eps函数（九）——Java也有"eps"函数》中，我们展示了Java中，与eps功能类似的方法ulp的用法，并讨论了ulp和eps的不同；
• 在《关于Matlab的eps函数（十）——MATLAB Coder生成的C代码》中，我们通过MATLAB Coder窥见了官方的eps实现方法；
• 在《GPU Powered Matlab(三)——关于Matlab的eps函数（十一）》和《GPU Powered Matlab（三点一）——关于Matlab的eps函数（十二）》中，我们试着将eps搬到CUDA GPU上执行。

今天要做的事情，是在LabVIEW中用纯粹的G语言实现eps运算。其实我现在并不需要在LabVIEW中使用eps，只是想通过实现eps来了解LabVIEW到底有多强的位运算能力。
首先，LabVIEW在“编程->数值”中提供了名为“计算机?”的节点，这个节点是一个常数节点（而非函数），其数值等于MATLAB中的eps(1)：


上图中数值显示控件“数值3”输出的数值为0，说明LabVIEW节点“计算机?”和Matlab的eps(1)相等。实际上，由于MATLAB的函数支持可变个数的参数，且函数调用的时候可以省去括号，因此eps(1)也可以写成eps这种形式，这就是导致很多人认为Matlab中的eps是一个常数而非函数的原因。如果勤快的话，在Matlab的Command Window中敲上一行“doc eps”就能看到关于eps的更多信息。
既然“计算机?”节点不是函数，我们就需要自己添加额外的程序来实现eps的功能了。实际上，C++中也提供了eps(1)的数值，而我们在《关于Matlab的eps函数（六）》中的实现方式，就是用位运算的方式获取给定的浮点数R的指数部分E，然后使用表达式
eps(1) * 2 ^ E
得到了eps(R)的值。这里我们也可以使用这种手法。
为了提取浮点数的指数部分，我们首先测试一下LabVIEW的位运算能力。想要对浮点数进行位操作，就需要语言提供一种reinterpret cast的机制。在Matlab中是typecast，而LabVIEW也提供了typecast，根据文档，LabVIEW的typecast节点可以实现：
*(type *) &R
的变换，正是我们所需要的。
浮点数本质上和整数没有区别，就是内存中的几个字节的集合，里面的位的取值包含了信息。和整数唯一的不同就是解读这些位的方式。然而没有学过低阶语言的编程者，受到了高阶语言编译器/解释器的照顾，浮点数的存储和计算这些细节完全被语言的机制所掩盖。这是一件好事，因为大家编程的时候就不用关心太多的细节。这也是一件坏事，让大家不用去了解浮点数的底层机制。时间长了，大家就形成了错误的印象：既然语言不让我对浮点数做位操作，那浮点数大概就是不能做位操作的吧。我以前在大学的时候，隔壁有一个项目组，因为不知道如何在C语言中将浮点数拆成字节数组以放到串口缓冲区中，差点连硬件电路都改了。
LabVIEW的typecast的用法几乎和Matlab的typecast函数一样。我们先复习一下Matlab的typecast函数。
首先将输出格式调成16进制：
>> format hex
取一个双精度浮点数15，转换成字节数组：
>> bytes = typecast(15, uint8)

bytes =

0000000000002e40

我们知道，浮点数是由符号，指数，尾数三个部分组成的。这里试着将字节数组最后一个字节的最高位改成1，对应的浮点数就是符号位变成1，会变成负数：
>> bytes(end) = hex2dec(c0);
>> format; dbl = typecast(bytes, double)

dbl =

-15

再试试对浮点数的指数部分进行操作。双精度浮点数有8个字节，64位。其中符号位1位，指数位11位，尾数位52位，如果想把指数加上1（等效于浮点数乘以2），只需要执行下面的操作：
>> typecast(typecast(15, uint64) bitshift(uint64(1), 52), double)

ans =

30

首先将双精度浮点数15 reinterpret成unsigned int64，然后加上1左移52位，再reinterpret成double。就是这么简单。
上面两个实例我们在LabVIEW用G语言是这样实现的：


非常重要的一点，上面这些蓝色的数值常量，需要手工指定其具体的整数类型。在常量的蓝色方框内点击右键，在弹出的菜单中选择“表示法->UINT64”显式指定数据的具体整数类型，否则运算会出错，无法得到期望的结果。


通过上面的实验，我们终于认识到，LabVIEW拥有类似于C的位操作能力。下一步是不是用typecast以及按位与，左移右移来提取浮点数的指数呢？其实不必这么麻烦，LabVIEW提供了一个专门的提取浮点数指数和尾数的节点。试着用“尾数与指数”节点获取浮点数15的指数：


这正是我们需要的数值。15位于8和16的中间，8是2的3次方，16是2的4次方，所以15的指数部分是3。现在将“计算机?”乘以2的三次方（利用“按2的幂缩放”节点），就得到了eps(15)的值：


我们将计算得到的数值与Matlab的eps(15)进行了比较。“数值5”控件显示计算的结果为1.77636E-15；“数值6”控件显示LabVIEW计算的结果与Matlab计算的结果的差，输出为0，说明和Matlab是一致的。利用这么简单的G语言框图，我们就实现了eps的计算。