float 在 CUDA

这篇文章主要为了说明float的相关背景和float在CUDA中的特性。

先说说 float 相关背景：

关于 float 在 MSDN 中的说明可以参考：

https://msdn.microsoft.com/zh-cn/library/hd7199ke.aspx

浮点数使用 IEEE 75（电气和电子工程师协会）格式。 浮点类型的单精度值具有 4 个字节，包括一个符号位、一个 8 位 excess-127 二进制指数和一个 23 位尾数。 尾数表示一个介于 1.0 和 2.0 之间的数。 由于尾数的高顺序位始终为 1，因此它不是以数字形式存储的。 此表示形式为 float 类型提供了一个大约在 3.4E–38 和 3.4E+38 之间的范围。

这里的有效位就是数据表示的精度，6~7位有效数字。

float总结：
共计 32 位，4 字节
由最高到最低位分别是第 31 、 30 、 29 、……、 0 位
31 位是符号位， 1 表示该数为负， 0 反之。
30-23 位，一共 8 位是指数位。
22-0 位，一共 23 位是尾数位。

更详细表示参见博文：浮点数在计算机中存储方式(转) 

注意点：

1. MSDN： “由于指数是以无符号形式存储的，因此指数的偏差为其可能值的一半。 对于 float 类型，偏差(offset 偏移)为 127；对于 double 类型，偏差(offset 偏移)为 1023。 您可以通过将指数值减去偏差值来计算实际指数值。 ”

这句话中提到指数是以无符号形式存储的，但是我们知道指数是有正有负的，这里该如何理解？

在这篇博文中提到：

再回来看指数，一共 8 位，可以表示范围是 0 - 255 的无符号整数，也可以表示 -128 - 127 的有符号整数。但因为指数是可以为负的，所以为了统一把十进制的整数化为二进制时，都先加上 127 ，在这里，我们的 16 加上 127 后就变成了 143 ，二进制表示为： 10001111

2. MSDN：“存储为二进制分数的尾数大于或等于 1 且小于 2。 对于 float 和 double 类型，最高有效位位置的尾数中有一个隐含的前导 1，这样，尾数实际上分别为 24 和 53 位长，即使最高有效位从未存储在内存中也是如此。 ”

这句话中提到最高有效位位置的尾数中有一个隐含的前导1，“120.5用二进制表示为：1110110.1，1110110.1可以科学记数法表示为1.1101101*2^6，任何一个十进制数的二进制科学计数法表示都为1.xxx*2^n(xxx是0，1表示的数字), 尾数部分就可以表示为xxxx,第一位都是1嘛，干嘛还要表示呀？可以将小数点前面的1省略，所以23bit的尾数部分，可以表示的精度却变成了 24bit，道理就是在这里。

那24bit能精确到小数点后几位呢，我们知道9的二进制表示为1001，所以4bit能精确十进制中的1位小数点， 24bit就能使float能精确到小数点后6位，而对于指数部分，因为指数可正可负，8位的指数位能表示的指数范围就应该为:-127-128了，所以指数部分的存储采用移位存储，存储的数据为原数据+127”。

指数有正有负，存储上却是用无符号来存储，那就再存入之前入一些偏移：在这篇博文中也提到：“再回来看指数，一共 8 位，可以表示范围是 0 - 255 的无符号整数，也可以表示 -128 - 127 的有符号整数。但因为指数是可以为负的，所以为了统一把十进制的整数化为二进制时，都先加上 127 ，在这里，我们的 16 加上 127 后就变成了 143 ，二进制表示为： 10001111”。

部分摘自：浮点数在计算机中存储方式(转)  

3. 浮点数对无穷大和非数也进行了表示：

Also, encodings to represent infinity and not-a-number (NaN) data are reserved. The IEEE 754 Standard [2] describes floating point encodings in full.

4. 浮点数精度问题：并非所有的实数都能被浮点数精确的表示，就是十进制小数在向二进制小数转换的过程中，会发现尾数的最大长度 23 位也除不尽，这就产生了浮点数精度问题。

例如分数2/3的二进制表示为 0.10101010... 是一个无穷的二进制小数，那么2/3就必须先进行近似以便用有限精度的浮点数来表示。而近似的规则和模式在 IEEE 754 中也做了详细的说明。这里不做详细展开，最常用的就是 round-to-nearest 方式。

2/3在这种方式下表示为：

指数-1被存为正整数 (-1) + offset，其中对于 float 类型，offset 为 127；对于 double 类型，offset 为 1023。

再说说 float 在 CUDA 中的特性。

在 CUDA 中浮点数的详细说明可以参见官方文档：Floating Point and IEEE 754

该文档将 CUDA 中的浮点数涉及到的问题描述得及其细致，就其中几点总结如下：

1. IEEE 754 标准要求对一系列操作（加减乘除开方等）的支持，这些操作对给定的格式和近似模式都保证在不同的操作实现方式下都能产生一样的结果，即与硬件平台无关。

2. 数学上数字运算规则和性质对浮点数来说并非是一定成立的，这是由于浮点数精度限制造成的。

例如数学上(A + B) + C = A + (B + C)是成立的，如果用 rn(x) 表示依据 IEEE 754 标准对 x 的近似操作，那么 rn(rn(A + B) + C) 和 rn(A + rn(B + C)) 与数学上的结果都不相同。

这表明相同的计算会产生不同的结果，即使他们所有的基本操作都是依据 IEEE 754 标准执行的。而结果的不同并不取决于执行主体，任何一个支持单精浮点数的微处理器 CPU，GPU 都会生成相同的结果。

3. FMA

在2008年 IEEE 754 标准引入FMA，而FMA 计算 rn(X * Y +Z) 只需要一步近似处理，因此较之前的 rn(rn(X * Y) + Z)更加精确，官方文档中示例很好的说明了这一点。而选择是否使用FMA是取决于你的平台和编译器的。

通过查看自己平台的处理器支持的指令集： AVX，AVX2，FMA 来知道是否支持 FMA

4. 精度示例：向量点乘

恰好最近在进行 CUDA 点乘运算，这个示例很好的解释了我遇到的现象。

在数学上，这个运算很简单也没有歧义，但是如果我们采用不同的策略，采用加，乘，FMA间不同的组织方式的话，就会产生不同的结果。以上的讨论都是在 IEEE 754 下。

方法1：serial，通过一个循环，串行计算，将乘和加分开

方法2：FMA 方法，将方法1中的加乘放在一起

显然方法2要比方法1精确。

方法3：parallel，分治的思想

这是一种递归的策略，由于各个部分是相互独立的，因此也可以成为并行方法，虽然实现上也可以采用串行。

三种方法在都使用 IEEE 754 标准下的比较：

• a = [1.907607, -.7862027, 1.148311, .9604002]
• b = [-.9355000, -.6915108, 1.724470, -.7097529]

Read more at: http://docs.nvidia.com/cuda/floating-point/index.html#ixzz3YQmGKgkh
Follow us: @GPUComputing on Twitter | NVIDIA on Facebook

可以看到每种方法产生的结果都与正确结果不太一样。一般来说精确度上有：FMA > parallel > serial。

5.  CUDA 对浮点数的支持

Current generations of the NVIDIA architecture such as Tesla Kxx, GTX 8xx, and GTX 9xx, support both single and double precision with IEEE 754 precision and include hardware support for fused multiply-add in both single and double precision.
CUDA 中 GPU 对特性的支持反映在 capability 号上。2.0及其以上均支持以上特性。

IEEE 754 支持四种近似模式：round-to-nearest(default), round towards positive, round towards negative, and round towards zero
编译器保留字可以指定近似模式：__fadd_[rn | rz | ru | rd] (x, y) etc. 具体参见手册。
编译也可以选择也可以相关参数，具体参见手册。

6. 一些建议

1.）使用 FMA

2.）仔细检查结构

3.）了解自己 GPU 的capability

4.）尽量使用 CUDA 数学库函数

至此我们对 CUDA 中的浮点数略知一二了。

再次说明，以上均来自官方文档：Floating Point and IEEE 754

在CUDA中，对双精度数支持还有一定要求：

设置GPU架构 sm_13 (如果你的GPU支持)，命令行: -arch=sm_13

细节可以参考： Double precision floating point in CUDA