---
# System prepended metadata

title: 2020q3 Homework5 (quiz5)
tags: [進階電腦系統理論與實作]

---

---
tags: 進階電腦系統理論與實作
---
# 2020q3 Homework5 (quiz5)
contributed by < `RinHizakura` >
> [第 5 週測驗題](https://hackmd.io/@sysprog/2020-quiz5)
> [GitHub](https://github.com/RinHizakura/NCKU-System-Software-Quiz/tree/master/quiz5)

## 測驗 `1`
### 程式原理
假設被除數為浮點數 $A$，除數為整數 $B$，我們可以把 $A / B$ 拆解：
$$
\begin{split}
\frac{A}{B} = \frac{A(B+1)}{B(B+1)} &= \frac{AB+A}{B(B+1)}  \\
&= \frac{AB}{B(B+1)}  + \frac{A}{B(B+1)}\\
&= \frac{A}{B+1} + \frac{A}{B(B+1)} \\
&= A/(B+1) + \frac{A/(B+1)}{B} \\
&= \frac{A}{2}/\frac{(B+1)}{2} + \frac{\frac{A}{2}/\frac{(B+1)}{2}}{B}
\end{split}
$$

這有甚麼意義呢?根據上面的式子，我們可以巧秒的把除法做拆解，從而使得真正的操作只有**對 2 的除法**(準確的說是對 1 和 2 的除法，不過我們小學的時候就知道任何數除以 1 就等於他自己...)。對於 $A \space divop \space B$，$A$ 為浮點數，而 B 為 整數：

Step 1. 判斷
1. 如果 B 為 2 的倍數(偶數)，則 $A/B$ 等同於 $\frac{A}{2}  \space divop \space \frac{B}{2}$，然後回到 step 1 重新判斷
2. 如果 B 為奇數，則把 $A  \space divop \space B$ 拆解成 
$$
\frac{A}{2} \space divop \space\frac{(B+1)}{2} + (\frac{A}{2} \space divop \space \frac{(B+1)}{2}) \space divop \space {B}
$$ 
    * 因為 B 是奇數，所以 $\frac{(B+1)}{2}$ 仍是整數 
    * 則 $\frac{A}{2} \space divop \space\frac{(B+1)}{2}$ 可以回到 step 1 向下拆解
    * 你可能會想問，$(\frac{A}{2} \space divop \space \frac{(B+1)}{2}) \space divop \space {B}$ 這一項不是有可能除到 2 以外的數嗎? 這是不會的，因為  $\frac{A}{2} \space divop \space\frac{(B+1)}{2}$ (假設叫 $C$ 好了) 一定比 A 還小，則 $C \space divop \space B$，最後一定可以拆成 "某個數" $+ \space (0 \space divop \space B)$，而 "某個數" 是拆解出來不需除以 B 的那些項

Step 2. 則我們可以把除法不斷向下拆解，直到 $A=0$ 或者 $B=1$

:::warning
表達的可能不是很好懂，你可以照著這個規則嘗試幾個數字來追蹤看看
:::

因此 `D1 = (c) 2`、`D2 = (d) 1`。



## 測驗 `2`

### 程式原理
首先來理解程式中的資料結構與巨集的作用:
* `ieee_float_shape_type` 這個結構使得我們可以把一段 IEEE 754 編碼的 float，轉換成 uint32_t 來表達
* `INSERT_WORDS(d, ix0)` 可以將整數 `ix0` 的編碼轉換成浮點數 `d`
* `EXTRACT_WORDS(ix0, x)` 反之將浮點數 `d` 的編碼轉換成表達整數 `ix0` 

舉例來說：
```
SIGN  EXPONENT  FRACTION
0     00000010  01000000000000000000000
```
* 如果用 `uint32_t` 的角度來解釋，就是 $2^{24}+2^{21}$
* 如果從 `float`(IEEE 754)的角度來解釋，則是 $(1 + 0.25) \times 2^{(2 - 127)}$

接著，我們逐步來理解程式是如何執行的:
```cpp
    int32_t sign = 0x80000000;
    EXTRACT_WORDS(ix0, x);

    /* take care of zero */
    if (ix0 <= 0) {
        if ((ix0 & (~sign)) == 0)
            return x; /* sqrt(+-0) = +-0 */
        if (ix0 < 0)
            return (x - x) / (x - x); /* sqrt(-ve) = sNaN */
    }
```
* 將浮點數 `x` 的編碼存入 `ix0`
* `if (ix0 <= 0)` 分成兩種情形處理
    * 如果 `(ix0 & (~sign))`，表示除了 sign bit 以外都是 0，則此編碼對應 0 或 -0，開根號得到自己
    * 若上述條件不滿足，且滿足 `ix0 < 0`，表示 `x` 的 sign bit 為 1，是負浮點數(或者 -INF 等不能開根號的數字)，則回傳 0.0 / 0.0 = NaN 
:::warning
粗略看了一下 IEEE 754 規範或者 C 語言的規範沒有看到 0.0 / 0.0 必須為 NaN 的說法，需要進一步確認哪裡有對應的描述。
:::
```cpp
    /* take care of +INF and NaN */
    if ((ix0 & 0x7ff00000) == 0x7ff00000) {
        /* sqrt(NaN) = NaN, sqrt(+INF) = +INF,*/
        return x;
    }
```
* 對於 +INF:
    * sign bit 為 0
    * exponent bit 的 8 bit 會是 `11111111`
    * fraction bit 則全零
* 對於 NaN:
    * sign bit 為 0 或 1 都可以(事實上為 1 的 case 在前一個 if 已經排除掉)
    * exponent bit 的 8 bit 會是 `11111111`
    * fraction bit 為非全零

對應的位置為 1 的 bitmask 是 0x7F800000，對於這兩個數字，開根號也是等於自己本身。
```cpp
    /* normalize x */
    m = (ix0 >> 23);
    if (m == 0) { /* subnormal x */
        for (i = 0; (ix0 & 0x00800000) == 0; i++)
            ix0 <<= 1;
        m -= i - 1;
    }
    m -= 127; /* unbias exponent */
    ix0 = (ix0 & 0x007fffff) | 0x00800000;
    if (m & 1) { /* odd m, double x to make it even */
        ix0 += ix0;
    }
    m >>= 1; /* m = [m/2] */
```
先理解一些前提可以更容易思考：
* 對於我們要開根號的數字 x，我們想把它表示成 $2^{2k} \times y$，則 `sqrt(x)` (k 為整數、而 y 要滿足 $1 \leq y < 4$，這個範圍是刻意為之，後面會解釋)，那麼求 $\sqrt x = 2^k \times \sqrt y$
* 先想像一下我們把小數點點在 `ix0` 的第 23 位和第 24 位之間，使得 `ix0` 變成小數點的二進位表達式來表示 y

接著再回顧程式碼：
* 取出 `ix0` 的前 9 個 bit，如果等於 0，表示 x 是 denormalized number(注意這是因為我們此前已經排除 0)
    * 此時的編碼表示的數字是 $0.FRACTION \times 2^{-126}$ (注意此時的 m 是 0，因此轉成 EXPONENT 是 -126 而不是後面的 -127)
    * `0x00800000` 是第 24 位，也就是說，為使得可以表達成 $y \times 2^{2k}$ 的形式，因此我們要位移使得 `ix0` 表示的 `y` 符合範圍
    * 則用來表示 exponent 的 m 也要減掉對應的位移量，應該要 `m -= i`，但注意到第 1 點的說明，因此調整成 `m = m - i + 1` 也就是 `m -= (i - 1)`
* m 代表 exponent 對應的數字，因此根據 IEEE 的規範，減掉 127 後可以得到對應數字 $(1.FRACTION) \times 2^{EXPONENT}$ 中的 ${EXPONENT}$
* `ix0 = (ix0 & 0x007fffff) | 0x00800000` 留下編碼後的 23 位，並在第 24 位加上 1，則符合我們前提中對 `ix0` 的想像
* 對於 $\sqrt{2^{EXPONENT}}$，如果 ${EXPONENT}$ 是偶數的話，那就是 $2^{EXPONENT/2}$，為此，如果是奇數 E(`if(m & 1)`)，則可以先把多出來的 2 乘去 $1.FRACTION$ 那邊，這是 `m >>= 1` 的目的 
    * 則我們只需要求對 $1.FRACTION$ 或 $2 \times 1.FRACTION$ 開根號的結果即可，這也是為甚麼前提中假設 $y < 4$，對應 $2 \times 1.FRACTION$ 的範圍
```cpp
    /* generate sqrt(x) bit by bit */
    ix0 += ix0;
    q = s0 = 0;      /* [q] = sqrt(x) */
    r = 0x01000000;  /* r = moving bit from right to left */

    while (r != 0) {
        t = s0 + r;      // t = s_i + 2^(-(i+1))
        if (t <= ix0) {  // t <= y_i ?
            s0 = t + r;  // s_{i+1} = s_i + 2^(-i)
            ix0 -= t;    // y_{i+1} = yi - t
            q += r;      // q_{i+1} = q_{i}+ 2^(-i-1)
        }
        ix0 += ix0;
        r >>= 1;
    }
```

因為 `y < 4`，綜合前面我們說想像成把小數點點在第 23 位與第 24 位之間，因此要從第 25 位開始判斷，所以 `QQ0 = (a) 0x01000000`。下面更仔細解釋此段程式：

如果 y 滿足 $1 \leq y < 4$，則 $\sqrt y$ 滿足 $1 \leq y < 2$，這使得我們可以從 1 開始逼近 $\sqrt y$ 的正確答案。假設 $q_i$ 代表 $\sqrt y$ 取小數點後 $i$ 位的精確度表示數值，則 $q_0 = 1$。

接著，我們就可以遞迴判斷小數點後要寫 0 或 1，補足至可以接受的精度。轉換成數學式的說法就是判斷
$$
式(1): {(q_{i} + 2^{-(i + 1)})}^2 \leq y  	
$$
如果成立則補 1($q_i = q_{i-1} + 2^{-(i + 1)}$)，否則補 0($q_i = q_{i-1}$)。

現在，我們定義：
$$
s_i = 2 \times q_i \\
y_i = 2^{i+1} \times (y - q_i^2)
$$


接著，把式 1 的展開為：
$$
\begin{split}
&q_i^2 + 2*qi*2^{-(i+1)} + 2^{-1} \leq y \\
&\implies 2^{-i}*q_i+2^{-2(i+1)} \leq y - q_i^2 \\
&\implies 2*q_i +2^{-(i+1)} \leq (y-q_i^2) *2^{i+1} \\
&\implies 式(2): s_i + 2^{-(i+1)} \leq y_i
\end{split}
$$

為甚麼要這麼做呢？這是因為 $y_i$ 和 $s_i$ 可以很容易的透過下列遞迴式運算：
* 如果式(2)不成立：
$$
s_{i+1} = s_i   \\
y_{i+1} = y_i
$$
這個很容易理解，因為此情況下 $q_{i+1}$ = $q_{i}$

* 如果式(2)成立：
$$
s_{i+1} = s_i + 2^{-i} \\
y_{i+1} = y_i - s_i - 2^{-(i+1)}
$$

此外要注意 `ix0 += ix0` 是把 `ix0` 向左位移，這是為了使 `while` 回圈中的 `if` 是在判斷 $s_i + 2^{-(i+1)} \leq y_i$。綜合這些概念後，再回頭看程式和我特別補充的標註應該就可以更好理解。

:::info
我是參考 [e_sqrt.c](https://www.netlib.org/fdlibm/e_sqrt.c) 的註解說明的，如有整理的不詳盡之處請直接回顧原說明
:::

```cpp
    /* use floating add to find out rounding direction */
    if (ix0 != 0) {
        z = one - tiny; /* trigger inexact flag */
        if (z >= one) {
            z = one + tiny;
            if (z > one)
                q += 2;
            else
                q += (q & 1);
        }
    }
```
這裡的目的是對精確度不足表示的部份做四捨五入，但如此進行的理由尚不理解...

* 如果 1 - 很小的數字

```cpp
    ix0 = (q >> 1) + 0x3f000000;
    ix0 += (m << 23);

    INSERT_WORDS(z, ix0);
    return z;
```
最後，要把計算出的結果還原成浮點數的表示
*  `ix0 = (q >> 1) + 0x3f000000`
    *  從程式中可以想像成 `q` 表示的是小數點擺在第 25 和 24 bit 之間的 $sqrt(y)$ 的二進位表示
    * 因此，要右移一位讓 $FRACTION$ 落在正確的位置
    * `QQ1 = 0x3f000000`: 照裡來說，我們要加上 bias 對應的 mask 是 `0x3f800000`，不過因為 `q` 特行如第一點所述且表示的 `1.FRACTION`，所以最右邊的 1 其實已經包含在 `(q >> 1)` 中被我們加進去了！
* `ix0 += (m << 23)`: 是算出的 $EXPONENT$ 部份，所以左移 23 位到正確的地方，所以 `QQ2 = (g) 23`

### [69. Sqrt(x)](https://leetcode.com/problems/sqrtx/submissions/)

#### solution `1`

根據上面的理解，接下來嘗試練習 LeetCode [69. Sqrt(x)](https://leetcode.com/problems/sqrtx/submissions/)，提交不需要 FPU 指令的實作。從題目可以看到我們只需要求到整數位的精準度就可以了，由此我的初步想法為：
```cpp
int mySqrt(int n){ 
    int32_t sign = 0x80000000;
    int32_t ix0, m, i;

    float x = n;
    EXTRACT_WORDS(ix0, x);
    
    /* take care of zero */
    if (ix0 <= 0) {
        if ((ix0 & (~sign)) == 0)
            return x; /* sqrt(+-0) = +-0 */
        if (ix0 < 0)
            return (x - x) / (x - x); /* sqrt(-ve) = sNaN */
    }   
    /* normalize x */
    m = (ix0 >> 23);
    if (m == 0) { /* subnormal x */
        for (i = 0; (ix0 & 0x00800000) == 0; i++)
            ix0 <<= 1;
        m -= i - 1;
    }    
    m -= 127; /* unbias exponent */
    ix0 = (ix0 & 0x007fffff) | 0x00800000;
    if (m & 1) { /* odd m, double x to make it even */
        ix0 += ix0;
    }   
    m >>= 1; /* m = [m/2] */
     
    // binary search 'm'     
    unsigned int head = 1 << m;       // 2^m
    unsigned int tail = 1 << (m + 1);  // 2^(m+1)
    unsigned int mid = (head + tail) >> 1; // same as (2^m + 2^(m+1)) / 2
    while(1){
        if(n > (mid + 1) * (mid + 1)){
            head = mid; 
            mid = (head + tail) >> 1;
        }
        else if(n < mid * mid){
            tail = mid; 
            mid = (head + tail) >> 1;
        }
        else
            break;
    }

    return mid;
}
```
延伸自原本程式的想法，減少了對一些測資中不會出現 input 做額外的判斷，並只是先求出整數 `n` 開根號後結果的 $1.FRACTION \times 2^{EXPONENT}$ 的 $EXPONENT$ 部份，則我們知道 n 開根號後的結果 $k$ 應該滿足 $2^{EXPONENT} \leq k < 2^{EXPONENT+1}$，因此後續可以用 binary search 找出結果。

在此實作下，我們不需要支援浮點數運算，甚至是整數的除法就可以得到結果，可以看到還不錯的執行時間：
![](https://i.imgur.com/KZGkbFO.png)

#### solution `2`
對於求開根號的近似值，我們也可以使用[牛頓法](https://en.wikipedia.org/wiki/Newton%27s_method)來求解。

```cpp
int mySqrt(int n){
    if(n == 0)
        return 0;
    if(n < 4)
        return 1;

    unsigned int ans = n / 2;
    
    if(ans > 65535)  // 65535 = 2^16 - 1
        ans = 65535;
    
    while( ans * ans > n  || (ans + 1) * (ans + 1) <= n){
        ans = (ans + n / ans) / 2;
    }   
    return ans;
}
```
不談論其數學的原理的話，程式基本上相當精簡(數學原理請參照維基百科或者原作業的 reference)，只要注意一些細節即可:
* `ans` 可以初始為 `n / 2`，是因為假設 $\sqrt n = x$，當 $x >= 2$ 時，$n = x^2 >= 2x$ 衡成立，也就是 $n / 2 >= x$，是符合的初始狀態
* `if(ans > 65535) ans = 65535;` 是為了避免 `ans * ans` 可能會 overflow 而導致錯誤，因此需要對初始狀態做額外的調整

得到的結果如下:
![](https://i.imgur.com/xGM7S2G.png)


上述的兩個作法中，因為 solution `1` 避免了除法的使用，因此我認為時間上比起 solution `2` 應該會有些許的提升。然而在 leetcode 上，我們無法很容易的計算兩者的效能差異，因此在本地進行一些測試。下圖為執行的時間結果，首先考慮 n 較小的情況下，從 0 到 1000 之間都搜尋時間，可以看到 solution `1` 的時間花費較少。

![](https://i.imgur.com/nH5nDMf.png)

考慮到在數字很大的時候，binary search 的範圍會變大，可能會導致時間變長，因此也觀察了在 `x = 1,000,000` 到 `x = 10,000,000` 間的執行時間差異。可以看到因為搜尋範圍的增加，不同的 n 時間變動比較大，但是整體而言，在n 達 1,000,000 時，solution `1` 執行所需的時間仍低於 solution `2`。

![](https://i.imgur.com/xx9lv8Z.png)

## 測驗 `3`
### 程式原理
延續原題目的說明，我們可以理解此題想求的是，在 k 和 x 皆為正整數的前提下
$$
N - \frac{(k-1)k}{2} = kx
$$
中，使得等號成立的整數 x 和 k 組合有多少個，如果用更好理解的形式表達，那就是:
$$
N - \{從 \space 0 \space 加到\space  (k - 1)\space 之和\} = kx
$$

因此，`ZZZ = (d) 1 - i`。

回顧原程式(下面 `i` 改成 `k`，與前式子做對應)，假設 N = 10，舉實際例子來解說：
```c=
int consecutiveNumbersSum(int N)
{
    if (N < 1)
        return 0;
    int ret = 1;
    int x = N;
    for (int k = 2; k < x; k++) {
        x += (1 - k); // = x -= (k - 1)
        if (x % k == 0)
            ret++;
    }
    return ret;
} 
```
`x` 初始為 10，則第一個 iteration 中:
* `k` 初始為 2
* `x -= (k - 1)` 變成 9，這是 N - {0 加到 1 之和}
* 因此檢查 `x % k` 是否為 0 就知道是否存在可行的整數解 x

第二個 iteration 中:
* `k++` 變成 3 
* `x -= (k - 1)` 變成 7，這是 `N - {0 加到 2 之和}`
    * 因為 `N - {0 加到 2 之和}` 就是 `N - {0 加到 1 之和} - 2`
* 同理，檢查 `x % k` 是否為 0 就知道是否存在可行的整數解

以此類推，直到 `k >= x`，這是因為如果 `k >= x`，下一個 `x -= (k - 1)` 會使得 x <= 0，已經不滿足我們的前提。 