# 浮點數運算和定點數操作
contributor: chenishi (陳奕熹), jesus255221 (謝立郇)

==[解說影片](https://youtu.be/b8QeU_z98PY)==

相較於一般整數的運算，我們似乎把交換律、結合律等等這些運算規則視為理所當然；不過同樣的運算操作搬到浮點數，上述規則還會成立嗎？

## 浮點數表示法

我們都知道在 C 裡面 float 是遵照 IEEE 754 來設計的，以下是 IEEE 754 規定的數字位


與一般的 32 bits integer 不同，這種設計方法的數值範圍更廣大，也具有表示小數點的能力，甚至還有 INF 與 NAN 的概念。但此一方法的設計會造成數值運算變得非常複雜，普通的四則運算就要考慮到三種不同的模式。也就是 float 的三種表示模式:
1. Normalize 模式，普通模式，表 $1.x *2^{exp}$
2. Denormalize 模式，表 $0.x * 2^{-126}$
3. NaN 與 INF 模式

從上途中可以看到，主要透過 Exponent 與 Fraction 來表示不同的取值模式

## 案例: 浮點數乘以二

### 程式邏輯

推導：欲使給定的浮點數乘以 `2`，必須考慮到 3 種情況
1. Normalize 模式的話，必須將指數加一後放回。
2. Denormalize 模式的話，比較有趣，可以直接將原數左移一單位，Mantissa 的 overflow 會自動變成 exp 的 1，之後再將 sign bit 放回即可。
3. 遇到 NaN, INF，就直接 return

換成數學語言即:
$\text {For any input x we have}
\begin{cases} \text{sign = the 31th bit of x}\\
\text{exponent = the 30 to 23 bits of x}\\
\text{mantissa = the 22 to 0 bits of x}
\end{cases}$

$\text{we have output y}=\begin{cases}sign \& (exponent + 1) \& mantissa & \text{for 0 < exponent < 255 (Within range)}\\
x & \text{for exponent = 255 (NaN of INF)}\\
sign \& (mantissa << 1) & \text{for exponent = 0}
\end{cases}$

轉換成 C 語言程式:

/* Return bit-level equivalent of expression 2 * uf
 * for floating point uf
 */
 unsigned float_twice (unsigned uf) {
    unsigned uf_sgn = (uf >> 31) & 0x01; //擷取首位正負號
    unsigned uf_exp = (uf >> 23) & 0xff; //擷取八位指數位
    if (uf_exp == 0) { //檢查為 denormalized 模式
        uf <<= 1; // 將數值左移一位也就是乘以二
        uf += (uf_sgn << 31); //將 sign bit 放回原數
    }
    else if (uf_exp != 0xff) { //檢查不是 INF or NaN （也就是剩下來的 Normalized 模式）
        uf_exp += 1; //讓 exp + 1 來乘以二
        /* 更新 uf */
        uf &= 0x807fffff; //把 exp 的部分遮蔽掉
        uf = uf + (uf_exp << 23); //把處理完的 exp 放回去
    }
  return uf;                                                           
}

### 數值範圍檢查
1. case NORMALIZE:
NORMALIZE 下能表示的數字為 ±(1.F) \* 2 ^ (E - 127)
其中 E 範圍在 1 ~ 254 之間
以最大值 E = 254 為例，經過 `float_twice` 之後 EXP 部份變成 254 + 1 = 255 也就是 infinity 或者是 NaN 的表示法
:::warning
理論上這邊結果不應該是 NaN 而應該是 infinity
不過要怎麼確認 mantissa 部份在這個時候會變成零
是不是要補個條件判斷或 mask
> [name=chenishi]

/* Return bit-level equivalent of expression 2 * uf
 * for floating point uf
 */
 unsigned float_twice(unsigned uf) {
    unsigned uf_sgn = (uf >> 31) & 0x01;
    unsigned uf_exp = (uf >> 23) & 0xff;
    if (uf_exp == 0) {
        uf <<= 1;
        uf += (uf_sgn << 31);
    }
    else if (uf_exp != 0xff) {
        uf_exp += 1;
        /* 如果從 normalize 加到 denormalized */
        if (uf_exp == 0xff) {
            uf &= 0x80000000;// 歸零 mantissa
        /* 正常從 normalize 到 normalize */
        } else {
            uf &= 0x807fffff;// 留下 mantissa
        }
        uf = uf + (uf_exp << 23);
    }
  return uf;                                  
}

:::

2. case DENORMALIZE:
DENORMALIZE 下條件為 EXP 部份全部為零，如果 `float_twice` 前，mantissa 為 0b1......（也就是首位為 1）
> 舉 0-00000000-1111111111111111111111 為例（denormalize 情況下的最大數）
> E = 1 - Bias = -126
> M = 2^23 - 1 / 2^23

`float_twice` 之後
> 0-00000001-1111111111111111111110
> E = e - Bias = -126
> M = (2^23 - 1 / 2^23) * 2 (可以看成原本的 Mantissa 左移一位)

這正好符合 CS:APP 中關於 NORMALIZE/ DENORMALIZE 中 EXP 的解釋 (FLOAT 等級距的性質)

3. case NaN/ infinity:
直接回傳不影響

## 案例: 浮點數的絕對值

考慮以下程式碼

void absf(double *x) {
    *(((int *) &x) + S1) S2= 0x7fffffff;
}

推導：浮點數絕對值
$\text{For any input float number x, we have output y such that}
\begin{cases} y = -x & \text{if y < 0} \\
y = x & \text{if y >= 0}
\end{cases}$
我們知道，浮點數的正負號是由最前面的 MSB 控制的，所以想要取得絕對值，只要

$\text{For any input float number x, we have output y such that }\\
y=f(x) \text{ where f(x) will set the MSB of x to 0}$

那問題就變成要怎麼把 $x$ 的 MSB 設為 0。

#include<stdio.h>
#include<stdlib.h>

void absf(double *x) {
    int a = 0x1;
    char little_end = *((char *) &a);
    *(((int *) x) + little_end) &= 0x7fffffff;
}

int main() {
    double x = -1;
    printf("%f \n", x);
    absf(&x);
    printf("%f \n", x);
}

當此機器為小端機器時，程式碼中的變數`char little_end`會變為`0x1`，而大端機器則會得到`0x0`。

接著再利用此變數來控制 x 必須進行且運算的位元。現在假設最低記憶體位址為 a 最高記憶體位址為 a + 7。對大端機器來說，MSB 即在 a 這個地址。所以必須要與 `0x7FFFFFFF`進行且運算的位址即 a ~ a + 3。反之，小端機器要進行且運算的位址為 a + 4 ~ a + 7。根據此觀察，只有小端機器必須讓地址位移 + 1，此時的 +1 可以用 little_end 變數來取代。來達到無分支絕對值的功能。

*用 QEMU 模擬 MIPS 架構下 big-endian 的行為*

*驗證此機器為 big-endian 架構*

### Big & small endian
> 「一直都以為有兩種 byte orderings，現在才知道，其實只有一種。
> 開玩笑的，但是其中一個的確比較受推崇 :-)」
>
> 出自: Beej's Guide to Network Programming: [03-IP address、結構與資料轉](http://beej-zhtw.netdpi.net/03-ip-address-structure-translation/3-2-byte-order)

相較於現今 little endian 架構在筆電，桌機上繁盛的景象，big endian 架構看似日漸消失在我們的生活中 ([Big Endian 與他們的現狀](/s/BJRaOrGCX))；事實上，big endian 也被稱作 Network Byte Order

### 現實中的浮點數運算
實際上，現代 CPU 遇到浮點數運算時可能不會透過軟體加速，相反的大多會有專門處理浮點數運算的硬體單元 (Float Point Unit)，關於 C 語言的型態轉換可參照 [C 語言強制轉型](/s/rJJfgHEa7)。

### 兩種 NaN
根據 IEEE 754 2008 年版本 6.2 節
有兩種不同的 NAN
`X111 1111 1AXX XXXX XXXX XXXX XXXX XXXX`
如果在 A 的位置為 1 ，則為 quiet NAN
反之如果 A 的位置為 0 ，其他的位置為非 0 ，則為 signal NAN
6.2 節也提到當要傳送這個 NAN 時，應採取 quiet NAN 的形式，當要送給系統信號時，應當採取 signal NAN

為了更好的掌握 float 型態的 bit pattern，我們必須將 float 的位元表示法一位一位的輸出。然而 `printf("%x")` 並不能如 int 型態將 float 型態的完整 bit pattern 印在螢幕上。我們必須尋求其他方法！這時候，就是用 `union` 的好時機了

union ufloat {
    float f;
    unsigned u;
};

如此一來，只要初始話這個結構中的 unsigned 成員，即可將 bit pattern 印出。可惜的是，此方法在大端機器將會失效。

以下是使用此`union`的例子

#include <stdio.h>
#include <math.h>

union ufloat {
    float f;
    unsigned u;
};

int main() {
    float f1, f2;
    *(int *)&f2 = 0;
    *(int *)&f1 = 0;
    union ufloat u1;
    u1.f = f1 / f2;
    printf("%f\n", u1.f);
    printf("%x\n", u1.u);
    return 0;
}

輸出為

可以清楚的觀察到 A 的位置為 1，很明顯的這是個 quiet NAN。

### 如何把 float 設為 nan
在決定什麼運算會產生 NAN 之前，我們必須要先會把 float 型態指定成 NAN。然而

float a = (float) 0x7F800001

不會把 float 設定成 NAN ，只會把 a 指派成 0x7F800001 的整數數值。

必須要使用以下的方法

*(int*) &f2 = 0x7F800001;

我想，看到 nan 的那一刻，才開始真正會寫程式吧。

> 這是不是因為上面那個 `float a = (float)0x7F800001` 實際上會涉及到整數轉浮點數 `cvtsi2ss` 這樣的指令
> 下面就只是把值取出來用 int 讀而已
> [name=EC]

### 會產生 NaN 的運算
根據 IEEE 754 7.2 節的說明，會產生 quiet NAN 的運算如下：
1. any general-computational or signaling-computational operation on a signaling NaN (see 6.2), except for some conversions (see 5.12)
2. multiplication: multiplication(0, ∞) or multiplication(∞, 0) fusedMultiplyAdd: fusedMultiplyAdd(0, ∞, c) or fusedMultiplyAdd(∞, 0, c) unless c is a quiet NaN; if c is a quiet NaN then it is implementation defined whether the invalid operation exception is signaled
3. addition or subtraction or fusedMultiplyAdd: magnitude subtraction of infinities, such as: addition(+∞, −∞)
4. division: division(0, 0) or division(∞, ∞)
5. remainder: remainder(x, y), when y is zero or x is infinite and neither is NaN
6. squareRoot if the operand is less than zero
7. quantize when the result does not fit in the destination format or when one operand is finite and the other is infinite

會產生 signal NAN 的運算如下：
1. conversion of a floating-point number to an integer format, when the source is NaN, infinity, or a value that would convert to an integer outside the range of the result format under the applicable rounding attribute
2. comparison by way of unordered-signaling predicates listed in Table 5.2, when the operands are unordered
3. logB(NaN), logB(∞), or logB(0) when logBFormat is an integer format (see 5.3.3).

## NAN 的 sign bit

不像 INF 在規格書中有詳細規定正負 sign bit 的運算結果，根據 6.3 節的敘述，除了特殊的 bit string 與 totalOrder 函數以外，其他運算並不保證 NAN 的正負。

## 數學函數的測試

針對這些平常數學運算中不會看到涉及 INF，NaN；但是事實上 C 語言有針對這些運算定義其運算結果(像是 NaN 會在數值運算像是~~瘟疫~~一樣傳遞)

我們知道常見的數學函數皆有其定義域，諸如 log, arcsin, 與偶次方根。在不考慮複數的情況下（浮點數無法表示虛數），以下函數沒有定義。

### 除法
除以零並沒有定義
測試結果：
1. INF : 0
2. -INF : -0
3. 0 :
1. if devidend is INF then INF
2. if devidend is -INF then -INF
3. if devidend 0 then NAN
4. if devidend NAN then NAN
5. NAN : NAN

### log
對數函數對負數並沒有定義
測試結果：
1. INF : INF
2. -INF : NAN
3. 0 : -INF
4. NAN : NAN
5. negative number : NAN

### arcsin
反三角函數對大於一以及小於負一的區域沒有定義
測試結果：
1. INF : NAN
2. -INF : NAN
3. positive number greater than 1 : NAN
4. negative number less than -1 : NAN
5. NAN : NAN

### 偶次方根
偶次方根對負數沒有定義
1. INF : INF
2. -INF : NAN
3. -0 : -0
4. negative number : NAN
5. NAN : NAN

## 定點數與浮點數比較

#### 什麼是定點數？
定點數與浮點數的差別最大在於「十進制中小數點位置是否固定」
像我們一般用的 int 整數就是「定點數」的一種，其小數點固定在最右方
這樣的表示方法用在小數上，代表我們只能夠表示固定大小的小數位
舉例來說，假設我們定點數小數點固定放在右三到右四之間，那麼
> 10 / 3 = 3.333

如果是左五到左四之間
> 10 / 3 = 3.3333

也就是說，隨著我們設計小數點位置不同，定點數能表示的範圍以及數值精度都會受到影響

:::info
簡單來說，相較於 IEEE 754 以類似「科學記號」的方式定義下的浮點數小數表示法，定點數就是可以自定義數值範圍的小數表示法 (自行定義你的小數點位置)，不過代價就是犧牲表示極大/小數的能力

關於定點數怎麼決定「小數點的位置」，你可以直接寫死 (舉例直接強制定義小數點在左邊數來第二位之類的)，但是也可以在變數中偷偷用一部分的位元紀錄該數值的小數點位置 (類似 IEEE 754 的 exponential 的概念)，不過需要取值運算的時候就需要自行做位移與解讀
:::

以 32 byte 數值為例，IEEE 754 定義下的浮點數在 Normalize 模式下，可以表示從 $2^{127}$ 到 $2^{-127}$ 的值域，儘管精度只有 Mantissa 部份也就是 23 位（比定點數的 32 位小），但是可表示的值域遠大於定點數的極大值 $2^{31}$ 到 $2^{-31}$

假設我們設計一個定點數，資料型態為`unsigned`，跟 denormalize 模式的浮點數可表示最小位相同。因此，我們假定此定點數的第 22 位為$2^{-127}$，與普通浮點數不同的是我們前面不考慮 exp 以及 sign 的部份，也就是

XXXX XXXX XXXX XXXX XXXX XXXX XXXX XXXX
            |_(2^-127)

測試環境

CPU : Intel E3-1230V2
RAM : 20 GB
MB  : H77 chip set
OS  : Ubuntu 18.04 desktop version

#### 單一浮點數運算測試
- [ ] 定點數測試

int main() {
    unsigned a;
    unsigned b;
    *&a = 0x7FFFFF;
    for (int i = 0; i < 10000; i++) {
        b = a * 2;
    }
}

- [ ] 浮點數測試

int main() {
    float a;
    *(int *) &a = 0x7FFFFF;
    for (int i = 0; i < 10000; i++) {
        b = a * 2;
    }
}


上半部是定點數的測試結果，平均執行時間: **49.35** ns
下半部是浮點數的測試結果，平均執行時間: **50.38** ns

兩套測試程式各執行 10000 次測試，其實 floating point 與 fixed point 的表現是差不多的，反而是定點數需要更多**指令數**與**cycle數**

:::info
如果不是乘以二而是乘以其他數值如 7 或者 5 等奇數會怎樣呢？
:::

以下是 `b = a * 5` 的結果，雖然定點數的組合語言比較多，但是**定點數所花的時間卻是浮點數的一半**，大幅增加了運算速度

上半部是定點數的測試結果
下半部是浮點數的測試結果

定點數的組合語言

浮點數的組合語言

#### [大量隨機浮點數運算測試](https://github.com/chenIshi/floatLab)

針對幾個浮點數區域分析
> N2N: 運算輸入為 Normalized，運算輸出也為 Normalized
> N2D: 運算輸入為 Normalized，運算輸出為 Denormalized
> D2N: 運算輸入為 Denormalized，運算輸出為 Normalized
> D2D: 運算輸入為 Denormalized，運算輸出也為 Denormalized
> N&S: Normalized 數值與特殊值（NaN 與 INF）運算
> D&S: Denormalized 數值與特殊值（NaN 與 INF）運算


D2D，其中使用的浮點數測試運算是 "x \*= 2" ，浮點數平均用時 15.63 奈秒，定點數平均用時為 14.86 奈秒，定點數用時較少的測資佔 57.77%，平手的狀況為 35.08%, 浮點數贏機率為 10.15%


另外一張 D2D，其中運算採用 "x \*= 5"，可以看到兩者用時開始出現巨大差異，浮點數平均用時變成 19.09 奈秒，定點數平均用時變成 14.90 奈秒；定點數贏得機率大大提昇到 99.78%（浮點數機率剩下 0.18%）

在這邊可以看到，定點數的平均用時保持固定，但是浮點數的用時上升了 22%


D2D，同樣的樣本數，運算也是 "x \*= 2"，可以看到散布圖情況也與 N2N 在「乘與二」時差不多，平均浮點用時 15.65 奈秒，定點用時 14.65 奈秒；定點數贏得機率 55.43%，浮點數則為 9.87%

D2N, 十萬筆 "x \*= 2"，浮點數平均用時為 15.96 奈秒，定點數平均為 15.01 奈秒，浮點數用時較高的比例為 92.43%

很令人意外的，原本以為可能是 「乘與二」 可以轉換成 `shift` 加速的原因，這邊把運算換成 "x \*= 3" 之後結果與 「乘與二」 大抵類似，浮點用時 15.68 奈秒，定點 15.00 奈秒
定點數贏得機率依然是 52.68%，浮點數 36.56%

換成 "x \*= 5" 之後結果依然類似

不過在 D2N 裡面又看到用時上懸殊的差異了，定點數贏得機率又上升到 99.73%


由於定點數除非特殊自己定義，不然無法表達 NaN 與 INF 的觀念，
以上這兩章圖是 「隨機 normal/denormal 浮點數與 NaN 相乘」 「隨機 normal/denormal 浮點數與 INF 相乘」

數據上，涉及到 NaN 運算中， Normalized 用時 （15.35 奈秒） 與 Denormalized （15.31 奈秒） 差不多， 不過在 INF 時有極端的差異， Normal時還是很正常的平均用時 15.33 奈秒，但是到了 denormal 時跳到 47.68 奈秒，整整多了 211%

> 需要額外注意的幾點是
> 1. NaN 數值中 Mantissa 部份是可以放非零的數值的，這部份數值跟產生 NaN 的原因有關（要查 IEEE 754 裡面有沒有更進一步的解釋）
> 2. INF 的圖後面有著很詭異的突起，而且這個情況在兩種浮點數中表現是類似的，這可能是 gitter ?
> 3. NaN 的運算中， 標準差好像變小了？
> [name=EC]

## 浮點數最佳化的相關應用

#### 深度學習
身為近最常被提及的 Buzzword，為了加速學習速度，簡化模型運算量是一件非常重要事情，在 2018年 11 月初 Facebook 工程師 Jeff Johnson 發布〈[Making floating point math highly efficient for AI hardware](https://code.fb.com/ai-research/floating-point-math/)〉，提及浮點運算在硬體上的限制，**基本上就是 FPU 需要大量額外硬體支援**:
1. Large word size
運算前的資料搬運，像是 external DRAM (主記憶體) 到 internal SRAM (L1~L3) 到暫存器
因此浮點數需要的位元組越大，對應更多在 IO 的時間開銷
2. General fixed point machinery
浮點數的 Significand (也被稱為 Mantissa) 其實就是定點數，因此需要特殊的硬體(像是定點數加法器)支援運算，如果想要增加有效位元數就必須提高 significand 位元數
3. General floating point machinery
浮點數的運算實作涉及到像是 Leading Zero Count 電路實作，為了達到浮點數精度需求需要大量硬體支援
4. IEEE 754 specific machinery
除了上述浮點數通用機制上的大量硬體需求(像是 LZ counter, 額外位移器)，IEEE 754 的浮點數為了 gradual underflow 的機制有了 denormalized float 的設計，這更需要額外的硬體支援

其實 Facebook 早在 2017 年 12 月就發布論文〈[Applied Machine Learning at Facebook: A Datacenter Infrastructure Perspective](https://research.fb.com/wp-content/uploads/2017/12/hpca-2018-facebook.pdf)〉，其中最後一段 `Future Directions in Co-design: Hardware, Software, And Algorithms` 便提到定點數的使用。

Google 也在 2018 年六月發表一篇 [Quantizing deep convolutional networks for efficient inference: A whitepaper](https://arxiv.org/pdf/1806.08342.pdf) 的論文，則更激進地使用自訂的數值系統來取代浮點運算。

> 目前我們正在嘗試用 C 建立一個機器學習模型 [FPNN](https://github.com/chenIshi/FPNN)，可設定調整權重的實作採用 IEEE 754 float/ 32-bit fixed-point/ fake float using unsigned 來比較浮點數/ 定點數/ 軟體實做的假浮點數差異


Intel x86-64 已非 8087 模式，有新的 FPU，透過 SSE 與 AVX 進行運算。每個核存在特定的 mask 以控制 FPU 行為，包含 rounding 及 denormalize number的 處理。另外可選擇要不要用軟體處理 Float point exception，例如 overflow/underflow，沒有 mask 的話，會在需要軟體介入的時候送 Exception。在 GNU/Linux 上，是由 glibc 跟 Linux 核心一同為新的任務做好設定。詳細可以參考 MXCSR 及 FXSAVE/FXRSTOR