# 2020q3 Homework5 (quiz5)
###### tags: `sysprog2020`
contributed by < `Msiciots` >
## :page_facing_up: 題目
- [2020q3 第 5 週測驗題](https://hackmd.io/@sysprog/2020-quiz5)
- [2020q3 Homework5 (作業區)](https://hackmd.io/@sysprog/2020-homework5)
:::info
目的: 檢驗學員對 [浮點數運算](https://hackmd.io/@sysprog/c-floating-point), [數值系統](https://hackmd.io/@sysprog/c-numerics) 及 [bitwise 操作](https://hackmd.io/@sysprog/c-bitwise) 的認知
:::
## 測驗 `1`
`D1` = ?D
- `(c)` 2
`D2` = ?
- `(d)` 1
:::success
延伸問題:
1. 解釋上述程式碼運作原理;
2. 以編譯器最佳化的角度，推測上述程式碼是否可透過 [tail call optimization](https://en.wikipedia.org/wiki/Tail_call) 進行最佳化，搭配對應的實驗;
搭配閱讀 [C 語言: 編譯器和最佳化原理](https://hackmd.io/s/Hy72937Me) 及 [C 語言: 遞迴呼叫](https://hackmd.io/s/rJ8BOjGGl)
3. 比照 [浮點數運算和定點數操作](https://hackmd.io/@NwaynhhKTK6joWHmoUxc7Q/H1SVhbETQ)，改寫上述程式碼，提出效率更高的實作，同樣避免使用到 FPU 指令 (可用 `gcc -S` 觀察對應的組合語言輸出)，並與使用到 FPU 的實作進行效能比較
:::
### 延伸問題 1. 程式運作原理
```cpp=1
#include <stdio.h>
#include <stdlib.h>
double divop(double orig, int slots) {
if (slots == 1 || orig == 0)
return orig;
int od = slots & 1;
double result = divop(orig / 2, od ? (slots + 1) >> 1 : slots >> 1);
if (od)
result += divop(result, slots);
return result;
}
```
- 假設被除數為浮點數 $A$，除數為整數 $B$：
```c=8
double result = divop(orig / 2, od ? (slots + 1) >> 1 : slots >> 1);
```
1. 若 $B$ 為偶數，$\frac{A}{B}=\frac{\frac{A}{2}}{\frac{B}{2}}$
2. 若 $B$ 為奇數，$\frac{A}{B}=>\frac{\frac{A}{2}}{\frac{B+1}{2}}$，讓 $B+1$ 成為偶數再進行 `1.`，經過加 1 的除數使結果值變小，我們需要補上差值
```c=10
result += divop(result, slots);
```
$差值=\frac{\frac{A}{B+1}}{B}$，我們加回原本的 `result` 可得到 $\frac{A}{B}=\frac{A}{B+1}+\frac{\frac{A}{B+1}}{B}=\frac{AB+A}{B(B+1)}=\frac{A(B+1)}{B(B+1)}=\frac{A}{B}$
- 遞迴運算至被除數 `orig` 接近`0` 到精確度已經不足，或者除數 `slots == 1`（偶數除法完成）
```c=5
if (slots == 1 || orig == 0)
return orig;
```
以下為計算 $\frac{1.0}{3}$ 結果：
```c=
orig: 0.500000 , slot:2 , result:0.250000 , od:0
orig: 0.125000 , slot:2 , result:0.062500 , od:0
orig: 0.031250 , slot:2 , result:0.015625 , od:0
orig: 0.007812 , slot:2 , result:0.003906 , od:0
orig: 0.001953 , slot:2 , result:0.000977 , od:0
orig: 0.000488 , slot:2 , result:0.000244 , od:0
orig: 0.000122 , slot:2 , result:0.000061 , od:0
orig: 0.000031 , slot:2 , result:0.000015 , od:0
orig: 0.000008 , slot:2 , result:0.000004 , od:0
orig: 0.000002 , slot:2 , result:0.000001 , od:0
.
.
.
orig: 0.000001 , slot:3 , result:0.000000 , od:1
orig: 0.000004 , slot:3 , result:0.000001 , od:1
orig: 0.000015 , slot:3 , result:0.000005 , od:1
orig: 0.000061 , slot:3 , result:0.000020 , od:1
orig: 0.000244 , slot:3 , result:0.000081 , od:1
orig: 0.000977 , slot:3 , result:0.000326 , od:1
orig: 0.003906 , slot:3 , result:0.001302 , od:1
orig: 0.015625 , slot:3 , result:0.005208 , od:1
orig: 0.062500 , slot:3 , result:0.020833 , od:1
orig: 0.250000 , slot:3 , result:0.083333 , od:1
orig: 1.000000 , slot:3 , result:0.333333 , od:1
```
## 測驗 `2`
QQ0 = ?
- `(a)` 0x01000000
QQ1 = ?
- `(d)` 0x3f000000
QQ2 = ?
- `(g)` 23
:::success
延伸題目:
1. 解釋上述程式碼何以運作
- 搭配研讀 [以牛頓法求整數開二次方根的近似值](http://science.hsjh.chc.edu.tw/upload_works/106/44f6ce7960aee6ef868ae3896ced46eb.pdf) (第 19 頁)
2. 嘗試改寫為更有效率的實作，並與 FPU 版本進行效能和 precision /accuracy 的比較
3. 練習 LeetCode [69. Sqrt(x)](https://leetcode.com/problems/sqrtx/)，提交不需要 FPU 指令的實作
:::
### 延伸問題 1. 程式運作原理
延續 從 [$\sqrt{2}$ 的運算談浮點數](https://hackmd.io/s/ryOLwrVtz)，假設 float 為 32-bit 寬度，考慮以下符合 IEEE 754 規範的平方根程式
```c=1
#include <stdint.h>
/* A union allowing us to convert between a float and a 32-bit integers.*/
typedef union {
float value;
uint32_t v_int;
} ieee_float_shape_type;
/* Set a float from a 32 bit int. */
#define INSERT_WORDS(d, ix0) \
do { \
ieee_float_shape_type iw_u; \
iw_u.v_int = (ix0); \
(d) = iw_u.value; \
} while (0)
/* Get a 32 bit int from a float. */
#define EXTRACT_WORDS(ix0, d) \
do { \
ieee_float_shape_type ew_u; \
ew_u.value = (d); \
(ix0) = ew_u.v_int; \
} while (0)
static const float one = 1.0, tiny = 1.0e-30;
float ieee754_sqrt(float x)
{
float z;
int32_t sign = 0x80000000;
uint32_t r;
int32_t ix0, s0, q, m, t, i;
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 */
}
/* take care of +INF and NaN */
if ((ix0 & 0x7f800000) == 0x7f800000) {
/* sqrt(NaN) = NaN, sqrt(+INF) = +INF,*/
return x;
}
/* 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] */
/* 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;
if (t <= ix0) {
s0 = t + r;
ix0 -= t;
q += r;
}
ix0 += ix0;
r >>= 1;
}
/* 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);
}
}
ix0 = (q >> 1) + 0x3f000000;
ix0 += (m << 23);
INSERT_WORDS(z, ix0);
return z;
}
```
以下一步步解說程式運行流程：
### 資料形態處理
```c=3
/* A union allowing us to convert between a float and a 32-bit integers.*/
typedef union {
float value;
uint32_t v_int;
} ieee_float_shape_type;
/* Set a float from a 32 bit int. */
#define INSERT_WORDS(d, ix0) \
do { \
ieee_float_shape_type iw_u; \
iw_u.v_int = (ix0); \
(d) = iw_u.value; \
} while (0)
/* Get a 32 bit int from a float. */
#define EXTRACT_WORDS(ix0, d) \
do { \
ieee_float_shape_type ew_u; \
ew_u.value = (d); \
(ix0) = ew_u.v_int; \
} while (0)
```
- 由於 `float` 與 `uint32_t` 皆佔 4 bytes 空間 利用 [union](https://en.wikipedia.org/wiki/Union_type) 可以存取 `float` 資料型態的 32 bit 數值。
- `EXTRACT_WORDS(ix0, d)` 輸入 float `d`，藉由 union 的方式得到 int 資料型態值
- `INSERT_WORDS(d, ix0)` 輸入 int `ix0`，藉由 union 的方式得到 float 資料型態值
- 在 macro 中使用 do {...} while(0) 來避免 [dangling else](https://en.wikipedia.org/wiki/Dangling_else)，參考自 [C multi-line macro: do/while(0) vs scope block ](https://stackoverflow.com/questions/1067226/c-multi-line-macro-do-while0-vs-scope-block)
### 根號內特殊數值處理
```c=36
/* 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 */
}
/* take care of +INF and NaN */
if ((ix0 & 0x7ff00000) == 0x7ff00000) {
/* sqrt(NaN) = NaN, sqrt(+INF) = +INF,*/
return x;
}
```
根據 [IEEE 754](https://zh.wikipedia.org/wiki/IEEE_754)：
> **特殊值**
這里有三個特殊值需要指出：
1. 如果**指數**是 0 並且尾數的**小數部分**是0，這個數±0（和符號位相關）
2. 如果**指數** = $2^{e}-1$並且尾數的**小數部分**是0，這個數是±∞（同樣和符號位相關）
3. 如果**指數** =$2^{e}-1$並且尾數的**小數部分**非0，這個數表示為[非數（NaN）](https://zh.wikipedia.org/wiki/NaN)。
以上規則，總結如下：
>| 形式 | 指數 | 小數部分 |
>| -------- | -------- | -------- |
>| 零 | 0 | 0 |
>| 非正規形式 | 0 |大於0小於1|
>| 正規形式 | 1 到 $2^{e}-2$ | 大於等於1小於2|
>| 無窮 |$2^{e}-1$ | 0 |
>| NaN | $2^{e}-1$ | 非0|
- $\sqrt{\pm{0}} = 0$
- 根號內為負數回傳 NaN, 不處理虛數狀況
- $\sqrt{\pm{∞}} = \pm{∞}$
- $\sqrt{NaN}$ = $NaN$
### 根號內數值正規化
```c=48
/* 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] */
```
- 根據 [IEEE 754](https://zh.wikipedia.org/wiki/IEEE_754)：
> 如果浮點數的指數部分的編碼值是0，分數部分非零，那麼這個浮點數將被稱為非正規形式的浮點數
>
```c=49
m = (ix0 >> 23);
```
`m` 代表浮點數的 $Exponent$ 位元值，如果 `m` 是 `0`，那程式會進行 normalize，從 $fraction$ 位元的最高為找到 `1` 往左位移直到 `m` 被補上 `1` ，再從 $Exponent$ 扣掉位移的數量
```c=55
m -= 127; /* unbias exponent */
```
- 根據 [IEEE 754 Exponent bias](https://en.wikipedia.org/wiki/Exponent_bias)，$Exponent$ 中的 $1\sim254$ 分別表示 $-126\sim127$，所以將 `m` 減去 127 可以得到 $Exponent$ 位元的實際值
```c=56
ix0 = (ix0 & 0x007fffff) | 0x00800000;
```
- `ix0` 存放 $fraction+1$ = $1.fraction$
```c=57
if (m & 1) { /* odd m, double x to make it even */
ix0 += ix0;
}
m >>= 1; /* m = [m/2] */
```
- 若 $m$ 為偶數，$\sqrt{x \times 2^{m}} = \sqrt{x} \times 2^{\frac{m}{2}}$
- 若 $m$ 為奇數，$\sqrt{x \times 2^{m}} = \sqrt{2x\times2^{m-1}}=\sqrt{2x}\times2^{ \lfloor{\frac{m}{2}}\rfloor}$
### 開根號運算
```c=62
/* generate sqrt(x) bit by bit */
ix0 += ix0;
q = s0 = 0; /* [q] = sqrt(x) */
r = QQ0; /* r = moving bit from right to left */
// QQ0 = 0x01000000
while (r != 0) {
t = s0 + r; // t = s_i + r_i
if (t <= ix0) {
s0 = t + r; // s_{i+1} = t + r_i
ix0 -= t; // x_{i+1}/2 = x_i - t
q += r; // q_{i+1} = q_i + r_i
}
ix0 += ix0; // x_{i+1}
r >>= 1;
}
```
- 參考 [Methods of computing square roots: Binary numeral system (base 2)](https://en.wikipedia.org/wiki/Methods_of_computing_square_roots)
- `q` 為均方根值，`ix0` 為根號內值，`r` 為每次處理的位元
- `if (t <= ix0)` 代表目前處理的 `ix0` 位元為 `1` ，`s0 = s0 + 2r` ，`ix0` 減去目前運算值`t`，目前運算位元 `r` 存入 `q`
- `t > ix0`，則目前 `ix0` 處理位元為 `0`
-
```c=75
ix0 += ix0; // x_{i+1}
r >>= 1;
```
位元移動，繼續處理下兩位值
### 處理進位
```c=75
/* 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);
}
}
```
- 參考 [IEEE 754: Rounding rules](https://en.wikipedia.org/wiki/IEEE_754#Rounding_rules)
-
```c=25
static const float one = 1.0, tiny = 1.0e-30;
```
- 透過 `one - tiny` 的進位行為來判斷當前系統使用何種 Rounding rules
- `one - tiny >= one` 代表系統進行進位，需要進一步判斷是何種進位方式
- `one + tiny > one` 得知系統屬於無條件進位，所以直接進位，`q += 2`
- `one + tiny <= one` 得知系統屬於四捨五入，只會在 `q & 1` 為 `1` 的時候進位
- `one - tiny < one`，代表系統不進位，不需要特別處理 `q`
### 運算值轉回 `float` 回傳
```c=90
ix0 = (q >> 1) + 0x3f000000;
ix0 += (m << 23);
INSERT_WORDS(z, ix0);
return z;
}
```
-
```c=63
r = 0x01000000; /* r = moving bit from right to left */
```
第 63 行 `r`從左邊多一位開始算所以 `q>>1`，加上 $bias$ (這裡加 $126$ 因為當初運算使用 $1.fraction$ 已經加上了 $1$)
- `m` 移到原本 $Exponent$ 位數加上 `ix0`
- 最後將運算值轉成 `float` 型態回傳
## 測驗 `3`
ZZZ = ?
- `(d)` 1 - i
:::success
延伸題目:
1. 解釋上述程式碼何以運作
2. 嘗試改寫為更有效率的實作
:::
### 延伸問題 1. 程式運作原理
LeetCode [829. Consecutive Numbers Sum](https://leetcode.com/problems/consecutive-numbers-sum/)
給定一個正整數 $N$，問 $N$ 能寫成多少種連續正整數之和，例如 $9$ 可寫成 $4+5$，或 $2+3+4$。
```cpp=
int consecutiveNumbersSum(int N)
{
if (N < 1)
return 0;
int ret = 1;
int x = N;
for (int i = 2; i < x; i++) {
x += 1 - i;
if (x % i == 0)
ret++;
}
return ret;
}
```
從數學式可看出
$$
kx = N - \frac{k(k-1)}{2}
$$
- 只要找到 $K$ 滿足$N-(1+2+3...+(k-1))$ 整除 $K$ 就有一種連續正整數組合
- `x` 值就是存放 $N-(1+2+3...+(k-1))$
- `i` 對應到 `k`
Sign in with Wallet
Connect another wallet