# 2020q3 Homework5 (quiz5)
contributed by < `ccs100203` >
###### tags: `linux2020`
> [第 5 週測驗題](https://hackmd.io/@sysprog/2020-quiz5)
## 測驗 1
浮點數除法程式: (fdiv.c)
假設 divop() 的第二個參數必為大於 0 的整數，而且不超過 int 型態能表達的數值上界。
```cpp=
#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;
}
```
- 利用遞迴的方式做除法。簡單來說程式會不斷的對被除數與除數同除 2，直到符合終止條件就會把 `result` 傳回去。但遇到 `slots` (除數) 是奇數時，因為會在第 8 行視為偶數直接運算，所以藉由第 10 行補上誤差。
- 第 5 行的 if 就是作為遞迴的終止條件，當 `slots` 等於 1 時就會把當前的 `orig` (也就是當前運算的 `result`) `return` 回去上一層。而 `orig == 0` 的判斷是用來處理 `orig` (被除數)是 0，或者是數字已經小到精度不夠，程式判定他為 0 時，也會結束並回傳。
- 第 8 行會直接做 $\div 2$ 並作為參數進入 recursive 的下一層。注意到他有對 `od` 做特別處理，也就是對 odd (奇數) 會特別做處理。當 `slots` 是奇數時會先做 +1 才除 2，也就是無條件進位。
- 但這樣子在 `slots` 是奇數時，數值會比實際答案小，因為除數無條件進位了，所以在第 9~10 行就是用來彌補奇數時的值，也就是會取近似值，不斷的把值透過遞迴呼叫加到 `result`，維持除法精準度。
假設 $10.0 / 3$，程式會把 `slots` 為 1 的結果全部相加起來，得以近似 3.3333...
```cpp
orig: 2.500000, slots: 1
orig: 0.625000, slots: 1
orig: 0.156250, slots: 1
orig: 0.039062, slots: 1
orig: 0.009766, slots: 1
orig: 0.002441, slots: 1
orig: 0.000610, slots: 1
...
```
寫成數學式的話會變成:
$$ \sum\limits_{i = 1}^\infty {\frac{10}{4^i}} = \frac{10}{3} $$
換句話說，就是利用等比級數的方式去算出誤差值，會計算出 $S_n$，將這些結果累加求得最後的 `result`
$$ S_n = \frac{a(1-r^n)}{1-r} $$
數學推導，假設 A 除以 n
$$\frac{A}{n} = \frac{A(n+1)}{n(n+1)} = \frac{An+A}{n(n+1)} \\
= \frac{A}{n+1} + \frac{A}{n(n+1)} \\
= \frac{A}{n+1} + \frac{A}{n} * \frac{1}{n+1} \\
= \frac{A}{n+1} + \frac{1}{n+1}(\frac{A}{n+1} + \frac{A}{n} * \frac{1}{n+1})$$
這段式子可以不斷延伸下去，就會發現符合等比級數的數學式
其中在 n 屬於偶數時會使用 $\div2$ 的方式去處理，式子就會類似下列情況，還是符合推導
$$\frac{A \div 2}{n \div 2} \lor \frac{A \div 2}{(n+1) \div 2}$$
## 測驗 2
假設 float 為 32-bit 寬度，考慮以下符合 IEEE 754 規範的平方根程式
```cpp=
#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;
}
```
### 關於 do...while 的解釋，請參考 [為什麼需要 do {...} while(0)](https://hackmd.io/@cccccs100203/linux2020-quiz6#%E7%82%BA%E4%BB%80%E9%BA%BC%E9%9C%80%E8%A6%81-do--while0)
### 程式運作原理
- `EXTRACT_WORDS` 輸入 float，藉由 union 的方式得到 int
- `INSERT_WORDS` 輸入 int，藉由 union 的方式得到 float
關於 [union](https://en.wikipedia.org/wiki/Union_type) 的說明參考這裡
#### `x` is less than or equal 0
根據 [Division by zero](https://en.wikipedia.org/wiki/Division_by_zero#Computer_arithmetic)
> In IEEE 754 arithmetic, a ÷ +0 is positive infinity when a is positive, negative infinity when a is negative, and NaN when a = ±0. The infinity signs change when dividing by −0 instead.
如果 x 是 -0 的話就直接 `return`，其餘則是 -NaN
```cpp=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 */
}
```
#### `x` is NaN or +INF
根據 IEEE 754 的規範，Exponent 的部分全為 1 時會是 INF 或 NaN。
程式遇到 NaN 或 +INF 時會直接 `return`
```cpp=43
/* take care of +INF and NaN */
if ((ix0 & 0x7f800000) == 0x7f800000) {
/* sqrt(NaN) = NaN, sqrt(+INF) = +INF,*/
return x;
}
```
#### `x` is denormalized
根據 IEEE 754 的規範，Exponent 為 0 時，可能遇到 subnormal 的情況
`m` 代表的就是 Exponent 的部分，如果 `m` 是 0，那程式會強制 normalize，會從 $fraction$ 找最高位的 1，不斷對 $fraction$ 往左 shift，一直到補進了 normalized 的 1，再將 shift 多少位數從 Exponent 的地方扣掉。
```cpp=48
/* normalize x */
m = (ix0 >> 23);
if (m == 0) { /* subnormal x */
for (i = 0; (ix0 & 0x00800000) == 0; i++)
ix0 <<= 1;
m -= i - 1;
}
```
#### `x` is normalized and Greater than 0
這時候的 `x` 一定是 normalized 的數字並且大於 0。
- `m -= 127` 是把 bias 的部分扣除，得到真正的 Exponent
- `ix0` 存放的是 $significand$ 也等同 $fraction + 1$
- 再來會做 `m >>= 1`，意即 $m = \lfloor m/2 \rfloor$，若 `m` 是奇數的話則會將多除掉的部分補到 `ix0` 內，所以才會做 `ix0 += ix0`。
```cpp=55
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] */
```
再來的部分就是做 sqrt，關於直式開 2 次方根，參考 [Binary numeral system (base 2)](https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Binary_numeral_system_(base_2)) 以及 [直式開二次方根(page 2~5)](http://science.hsjh.chc.edu.tw/upload_works/106/44f6ce7960aee6ef868ae3896ced46eb.pdf) 的做法。
```
q (r)
----------
s0 / ix0
(t)
t: 用來協助 s0 的紀錄，以及判斷當前的 ix0 是否大於 t
r: 用來協助 q 的紀錄，為當前正在處理的 bit。
從 0x01000000 開始是因為 ix0 最高位為 1 的 bit 會在第 25 (可能因為 m 為奇數而在第 26)
```
解說 while loop 運作原理:
1. 在直式運算中，會先從右到左，將 `ix0` 每兩個 bit 劃分為一組
2. 每一次會先將 `r` (當前處理的 bit) 與 `s0` (直式左邊的數字) 相加存入 `t`
3. 從直式運算的規則可以得知若 `ix0 >= t` 則從 `ix0` 扣掉 `t`，`q` (商) 則會增加 `r`。而根據直式運算 2 次方根的規則，`s0` 應更新為上次 `s0` + `2 * r`，由於 `t = s0 + r`，故程式內的寫法成立。
4. 若 `ix0 < t` 或執行完上述第 3 點，在直式規則中應一次從 `ix0` 新增兩個位數來處理，但這裡利用做 `r >>= 1` 以及 `ix0 += ix0`，恰巧符合運算規則，能夠接著運算 `ix0` 下兩個 bit。這樣的好處是不需要在 `ix0 < t` 時對 `s0` 做額外補 0 的動作，並且也不用移動 `s0` 來對齊下一步的 `ix0`。在 `q` 上也會因為 `r` 一次只往右 shift 1，符合其移動規則，不須額外調整。
```cpp=62
/* 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;
}
```
#### Rounding
參考 [Rounding](https://en.wikipedia.org/wiki/Rounding)、[Rounding modes](https://en.wikipedia.org/wiki/Floating-point_arithmetic#Rounding_modes)、[Rounding rules](https://en.wikipedia.org/wiki/IEEE_754#Rounding_rules)
我想程式是利用 `one` 與 `tiny` 的運算去判斷現在的系統屬於哪種 Rounding Mode，根據當前的 rounding mode 做出符合規範的 rounding。
- `one - tiny < 1.0`，代表當前的 mode 不會做進位，所以不需要特別處理 `q`，因為之後往右 shift 就會自動捨去
- `one - tiny >= 1.0`，代表當前的 mode 會進位，但要判斷屬於無條件進位或是四捨五入
- 藉由 `one + tiny > 1.0` 可以知道系統屬於無條件進位，所以直接進位，做 `q += 2`
- 藉由 `one + tiny <= 1.0` 可以知道系統屬於四捨五入，那只會在 `q & 1` 為 1 的時候才會進位
```cpp=78
/* 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);
}
}
```
#### 將數值轉換回 float
- 由於 `r` 是從第 25 位開始算，所以 `q` 需要往右 shift 1 才能夠對齊 IEEE 754 規範中 $significand$ 的位置。
- 然後再加上 126 來符合 bias。 (原本是 127，但已經從 $significand$ 得到 1，所以只需 126)
- 把先前的 `m` (Exponent) 往左 shift 23，將他放在 IEEE 754 規範的位置。
```cpp=90
ix0 = (q >> 1) + 0x3f000000;
ix0 += (m << 23);
```
### [LeetCode 69. Sqrt(x)](https://leetcode.com/problems/sqrtx/)
TODO
## 測驗 3
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;
}
```
我看不懂題目上的推導敘述，所以重新從等差數列的方法去理解程式。
程式的作法會在每一次的迴圈做 `x += 1 - i` 然後再去判斷 `x` 能不能整除 `i`
- 用下表說明 (程式內的 `i` 是從 2 開始的)
`i` 即為程式內的 `i`，代表輸入的數字 `N` 可以相等於被 `i` 個連續數字累加
`n` 則是代表連續數值累加中的第一個數字。
|`1-i`| -0 | -1 | -2 | -3 | -4 | -5 |
| :-: | :-: | :-: | :-: | :-: | :-: | :-: |
| `i` | 1 | 2 | 3 | 4 | 5 | 6 |
| `n` | n | n+1 | n+2 | n+3 | n+4 | n+5 |
1. 當 `i` 為 2 的時候，代表 $n+(n+1)=N$，而 $+1$ 的部分藉由 `1 - i` 扣掉了，所以實際上是判斷 ${(N-1)}\mod{2}$ 是否為 0 (可否整除)，這樣便可以判別是否能由 2 個連續數字累加為 $N$。
2. 再來當 `i` 為 3，代表 $n+(n+1)+(n+2) = N$，新的 $+2$ 也在這次的 `1 - i` 扣掉，$+1$ 則是在上一次的迴圈就扣掉了，所以判斷式變成 ${(N-1-2)}\mod{3}$ 是否為 0。從這邊我們可以往後推論每個 `i` 的情況，而迴圈一直執行到 `i >= x` 會結束，因為 `x` 此時無法繼續扣除 `1 - i`，就算兩者相等扣完之後 `x` 也只剩下 1，不可能整除 `i`。
3. 程式內的 `i` 是從 2 開始的，原因在於 `i` 等於 1 是絕對會成立的情況，所以節省了一次計算。
Sign in with Wallet
Connect another wallet