---
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 找出結果。
在此實作下，我們不需要支援浮點數運算，甚至是整數的除法就可以得到結果，可以看到還不錯的執行時間：
#### 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 而導致錯誤，因此需要對初始狀態做額外的調整
得到的結果如下:
上述的兩個作法中，因為 solution `1` 避免了除法的使用，因此我認為時間上比起 solution `2` 應該會有些許的提升。然而在 leetcode 上，我們無法很容易的計算兩者的效能差異，因此在本地進行一些測試。下圖為執行的時間結果，首先考慮 n 較小的情況下，從 0 到 1000 之間都搜尋時間，可以看到 solution `1` 的時間花費較少。
考慮到在數字很大的時候，binary search 的範圍會變大，可能會導致時間變長，因此也觀察了在 `x = 1,000,000` 到 `x = 10,000,000` 間的執行時間差異。可以看到因為搜尋範圍的增加，不同的 n 時間變動比較大，但是整體而言，在n 達 1,000,000 時，solution `1` 執行所需的時間仍低於 solution `2`。
## 測驗 `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，已經不滿足我們的前提。
Sign in with Wallet
Connect another wallet