# 2020q3 Homework5 (quiz5)
contributed by < `joe-U16` >
###### tags: `sysprog2020`
## 測驗 `1`
考慮到以下浮點數除法程式: (fdiv.c)
```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 / D1, od ? (slots + D2) >> 1 : slots >> 1);
if (od)
result += divop(result, slots);
return result;
}
// D1 = 2
// D2 = 1
```
假設 `divop()` 的第二個參數必為大於 `0` 的整數，而且不超過 `int` 型態能表達的數值上界。請補完程式碼。
**程式碼運作原理:**
* `int od = slots & 1;`
* 如果 slots 為奇數，則 `od = 1`
* `double result = divop(orig / D1, od ? (slots + D2) >> 1 : slots >> 1);`
* 先看 D2 ，如果 `od = 1` 則下一次的 slots 會是 `(slots + D2) >> 1`
* 如果 D2 是 2、3、4 的話，那每一個遞迴函式的 `slots` 都不會變為 1 ，這會讓最前面 `if (slots == 1 || orig == 0)` slot判斷沒有意義，所以 `D2 = 1`
* 再來是這個程式碼做除法的方式
* 把被除數和除數一起除以 2 ，當除數 slots 為 1 時即代表完成，因為任何數除以 1 都是自己。
* 可以合理推得 `D1 = 1`
## 測驗 `2`
先看到程式的兩個巨集：
```c=
/* 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)
```
上面巨集如同註解所述，把 int 轉換成 float ，且不經過換算，是直接轉換。
e.g. 0x40800000，以 int 表示即為 0x40800000 ，轉成以 float 表示 0x40800000 ，需使用 ieee754 標準，會使其所代表的數值為 4。
```c=
/* 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 轉換成 int ，且不經過換算，是直接轉換。
e.g. 數值 4 以 float 表示需使用 32 bit，以浮點數表示成 0x40800000，直接轉成 int 為 0x40800000。
```c=
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 */
}
```
使用巨集 **EXTRACT_WORDS(ix0, x)** 以 整數 `ix0` 儲存 `x` 的浮點數表示法。
* `if ((ix0 & (~sign)) == 0)` 判斷 Exponent 和 Mantissa 是否全為 `0` ，成立則代表 `ix0` 為 `+-0`，即可 return
* `if (ix0 < 0)` return sNaN
```c=
/* take care of +INF and NaN */
if ((ix0 & 0x7f800000) == 0x7f800000) {
/* sqrt(NaN) = NaN, sqrt(+INF) = +INF,*/
return x;
```
* `if ((ix0 & 0x7f800000) == 0x7f800000)` 判斷 Exponent 是否全為 `1`
```c=
/* 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] */
```
* 做 `sqrt(x)` 可將 x 表示成 $x = y \cdot 2^{2k}$ ，那麼 sqrt(x) = sqrt(y) $\cdot 2^k$
* 要確保 `k` 為整數， `m = 2k`， `m` 需為偶數。
* y 為 `1.fraction`，為了確保 m 為偶數，可能做成 $x = (y\cdot2)\cdot2^{2k-1}$
* 故 y 的範圍會在 `(1, 4]`
* `m = ix0 >> 23` 取到前 9 個bit，即為 Sign + Exponent。
* `if (m == 0)` 判斷是否為 Denormalize
* 若是，則 `i` 儲存 `ix0` fraction 部分的 leading zero
* `ix0 = (ix0 & 0x007fffff) | 0x00800000;` 將 `ix0` 取成 `1.fraction` 的形式，用到後面24個bit
* 浮點數公式為： $(-1)^s\cdot2^{e-127}\cdot(1.m)$
* `s: sign bit`, `e = exponent`, `m = fraction`
* `if (m & 1) ix0 += ix0` ，這邊的 `m` 為 exponent ，當 `m` 為奇數，將 `ix0` 乘 2 ，可看成將 exponent 的 2 丟給 `ix0`，再將 `m` 減 1，將 m 變成偶數，當然因為待會會做 right shift ，有沒有減 1 不影響 right shift 的結果，故可將減 1 這個動作省略
```c=
/* generate sqrt(x) bit by bit */
ix0 += ix0;
q = s0 = 0; /* [q] = sqrt(x) */
r = QQ0; /* r = moving bit from right to left */
while (r != 0) {
t = s0 + r; // t = s_i + 2^{-(i+1)}
if (t <= ix0) { // if t <= y_i
s0 = t + r; // s_{i+1} = s_i + 2^(-i)
ix0 -= t; // y_{i+1} = y_i - s_i - 2{-(i+1)}
q += r;
}
ix0 += ix0; //
r >>= 1;
}
// QQ0 = 0x01000000
```
讓 $q_i$ = sqrt(y) ， $q_0$=1 ， e.g. $q_0.q_1q_2q_3q_4q_5...$ ， `1.fraction`
$s_i = 2\cdot{q_i}\ \ \ \ {and} \ \ \ y_i = 2^{i+1}\cdot{(y - q_i^2)}\tag{1}$
To compute $q_{i+1}$ from $q_i$,one checks whether
${q_i+2^{-(i+1)}}^2 \le y\tag{2}$
If (2) is false, then $q_{i+1} = q_i$; otherwise $q_{i+1} = q_i + 2^{-(i+1)}$
With some algebric manipulation, it is not difficult to see that (2) is equivalent to
$s_i+2^{-(i+1)} \le y_i\tag{3}$
The advantage of (3) is that $s_i$ and $y_i$ can be computed by the following recurrence formula:
$s_{i+1} = s_i\ \ ,\ \ \ \ y_{i+1} = y_i \cdot2;\tag{4}$
:::warning
原本作者是寫：$s_{i+1} = s_i\ \ ,\ \ \ \ y_{i+1} = y_i ;\tag{4}$
:::
otherwise,
$s_{i+1} = s_i+2^{-i}\ \ ,\ \ \ \ y_{i+1} = y_i-s_i-2^{-(i+1)}\tag{5}$
One may easily use induction to prove (4) and (5)
* QQ0 = 0x01000000 等於可以做25次判斷，最開始的 `ix0+=ix0` 是為了讓最後一個位元可以拿來做rounding。
* 最下面的 `ix0 += ix0;` 即為上述式子(4)
---
* `ix0 = (q >> 1) + QQ1; //QQ1 = 0x3f000000`
* 為計算好的 q 加上 exponent的偏移值 `127`
* `ix0 += (m << QQ2); QQ2 = 23`
* 把算好的 exponent 放到bit 24-31 ，即浮點數表示法 exponent 的位子
**練習 LeetCode [69. Sqrt(x)](https://leetcode.com/problems/sqrtx/submissions/)，提交不需要 FPU 指令的實作**
使用 binary search 的方式找到 sqrt(x)
```c=
int mySqrt(int x){
int l = 0;
int h = 65535; //2^16 - 1 = 65535
unsigned int m = h;
while(l < h) {
m = (l + h + 1) / 2;
if (m * m == x) {
return m;
} else if(m * m > x){
h = m - 1;
} else {
l = m;
}
}
return h;
}
```
**Submission Result:**
## 測驗 `3`
[LeetCode 829. Consecutive Numbers Sum](https://leetcode.com/problems/consecutive-numbers-sum/) 給定一個正整數 N，問 N 能寫成多少種連續正整數之和，例如 9 = 9 = 4 + 5 = 2 + 3 + 4，所以需 return `3`。
參考實作如下 [1]:
```c=
int consecutiveNumbersSum(int N)
{
if (N < 1)
return 0;
int ret = 1;
int x = N;
for (int i = 2; i < x; i++) {
x += ZZZ;
if (x % i == 0)
ret++;
}
return ret;
}
// ZZZ = 1-i
```
**解釋程式原理：**
這個等差數列不必從 1 開始，假設從 x 開始，且個數共有 k 個，則可寫出該等差數列為：
`x + (x+1) + (x+2)+...+ k`
`kx + k*(k-1)/2 = N` implies
`kx = N - k*(k-1)/2`
那麼當 `k` 可整除 `kx` 即代表可以從 x 開始每次加 `1` ，共 k 個連續整數相加可得 正整數 N。
Leetcode 實作如下：
```c=
int consecutiveNumbersSum(int N)
{
if (N < 1)
return 0;
int ret = 1;
int x = N;
for (int k = 2; k < sqrt(2 * x); k++) {
if ((x - (k * (k - 1) / 2)) % k == 0) ret++;
}
return ret;
}
```
可以觀察到上下兩個程式碼原理相同，差別在於程式碼[1] 在 for 迴圈裡不停的加上一個負值，且此負值呈等差級數下降。
* 觀察到下面程式碼，當 k 值上升時的差值。
${k(k+1) \over 2} - {k(k-1) \over 2} = k$
* 可以得知程式碼 [1] 為每次扣掉差值 k ，才做後續判斷。
**Submission Result:**
## Reference
https://www.netlib.org/fdlibm/e_sqrt.c
https://leetcode.com/problems/consecutive-numbers-sum/discuss/129015/5-lines-C%2B%2B-solution-with-detailed-mathematical-explanation.
Sign in with Wallet
Connect another wallet