# 2020q3 Homework5 (quiz5)
contributed by < `zzzxxx00019` >
>[2020q3 第 5 週測驗題](https://hackmd.io/@sysprog/2020-quiz5)
## 測驗 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;
}
```
假設 `divop()` 的第二個參數必為大於 0 的整數，而且不超過 int 型態能表達的數值上界。請補完程式碼。
### 原理：
```cpp=5
if (slots == 1 || orig == 0)
return orig;
```
* 判斷除數是否為 1 或是被除數是否為 0 ，如果成立，答案即為被除數 orig
```cpp=7
int od = slots & 1;
double result = divop(orig / D1, od ? (slots + D2) >> 1 : slots >> 1);
```
* 判斷 slots 是否為偶數
* slots 為偶數，可視為 $\dfrac{orig}{slots} = \dfrac{orig/2}{slots/2}$ 並不影響結果，推得 `D1 = 2`
* slots 為奇數，可將 slots 加 1 使其可被 2 整除，故 `D2 = 1`
```cpp=9
if (od)
result += divop(result, slots);
```
* 考慮到 slot 為奇數，若 slot 加 1 ，要維持原結果，則可寫成：
$\dfrac{orig}{slots}=\dfrac{orig}{slots+1}\times\dfrac{slots+1}{slots}＝\dfrac{orig}{slots+1}\times(1＋\dfrac{1}{slots})=\dfrac{orig/2^n}{(slots+1)/2^n}+\dfrac{orig/2^n}{(slots+1)/2^n}\times\dfrac{1}{slots}$
* $\dfrac{orig/2^n}{(slot+1)/2^n}\times\dfrac{1}{slot}$ 的部分即為 `divop(result, slots)`
---
## 測驗 2
### 題目：
延續 [從 $\sqrt{2}$ 的運算談浮點數](https://hackmd.io/s/ryOLwrVtz)，假設 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 = QQ0; /* 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) + QQ1;
ix0 += (m << QQ2);
INSERT_WORDS(z, ix0);
return z;
}
```
請補完程式碼，使上述程式碼得以運作。
### 原理：
* [單精度浮點數](https://zh.wikipedia.org/wiki/%E5%96%AE%E7%B2%BE%E5%BA%A6%E6%B5%AE%E9%BB%9E%E6%95%B8) ，第 1 位表示正負，中間 8 位表示指數，後 23 位儲存有效數位
```
sign exponent(8 bits) fraction(23 bits)
0 01111100 01000000000000000000000
```
以上述為例：
$sign=+1$
$exponent=(-127)+124=-3$
$fraction=1+2^{-2}=1.25$
$value=sign\times 2^{exp}\times fraction=(+1)\times 2^{-3}\times 1.25=0.15625$
因此 `float` 在透過單精度浮點數的變換後，可以用 `unit32_t` 的形式來表達
接著逐步拆解程式的每個步驟：
```cpp=10
/* 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)
```
* `INSERT_WORDS(d, ix0)` 將 `ix0` 轉換成浮點數 `d`
```cpp=17
/* 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)
```
* `EXTRACT_WORDS(ix0, x)` 將浮點數 `d` 的編碼轉換成 `ix0`
```cpp=34
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` 表示，若 `ix0 = 0` ，則 `ix0 & (~sign) = 0` ，而 $\sqrt{0}=0$ ，因此直接 `return 0` ，反之則代表 `x` 為負數浮點數，根據定理， $\sqrt{-N}$ 並不存在，因此 `return (x - x) / (x - x)`
```cpp=43
/* take care of +INF and NaN */
if ((ix0 & 0x7f800000) == 0x7f800000) {
/* sqrt(NaN) = NaN, sqrt(+INF) = +INF,*/
return x;
}
```
* `0x7f800000` 以二進制呈現為以下結果，代表若 `ix0` 的 `exp` 全為 `1` ，表示 `x` 為 `+INF` ，而 `ix0` 的 `exp` 全為 `1` ， `fraction` 為任意值 ，則表示 `x` 為 `NaN` ，直接 `return x`
|sign|exponent|fraction |
| - | - | - |
| 0 | 11111111| 00000000000000000000000|
```cpp=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] */
```
以下參考自 < `RinHizakura` > 文章說明：
* 若 $x=2^{n}\times y$ ，且 $1≤y<4$，則 $\sqrt{x}=\sqrt{2^n\times y}=2^{n/2}\times \sqrt{y}$
* 在第 `48` 行， `ix0 >> 23` 後，留下的為 `sign` 與 `exponent` ，而 `sign` 為 `1` 的情況，在前面已過濾掉，因此 `m` 即為 `exponent`
* 當 `m = 0` ，表示 `float x` 的 `exponent` 為 `0` ，在 `denormalized` 情況下，將浮點數表達為 $x=2^{-126}\times0.fraction$
* 為了使 $x$ 得以表達為 $2^{n}\times y$ 的形式，將 `ix0` 左移直到 `exponent bit` 出現 `1` 為止
* 而原本在 `denormalized` 情況下 `exp = -126` ， `m = 0` ，即 `exp = -127 - m + 1` ，而考慮到 `ix0` 左移後，需將左移量補上，因此 `m = m + 1 - i` ，即 `m -= i - 1`
$float=1\times 2^{-126} \times (2^{-1}+2^{-2})=2^{-127}+2^{-128}$
|sign|exponent|fraction|
| - | -- | --- |
| 0 |00..00 | 11..00 |
**normalized float ，left shift 1 bit ：**
$m=1-i-127=-127$
$float = 1\times 2^{-127} \times (1+2^{-1})=2^{-127}+2^{-128}$
|sign|exponent|fraction|
| - | -- | --- |
| 0 |00..01 | 10..00 |
* 再透過 `(ix0 & 0x007fffff) | 0x00800000` 這樣一連串的邏輯運算，將 `bit 24` 設為 `1` ，如此可保留原本 `1.fraction` 的部分 ( 二進制分數表示法 $2^0 + 2^{-1} + 2^{-2}...=1.fraction$ )
* 檢查 `m` ，若為偶數， $\sqrt{2^m}=2^{m/2}$ ；反之，若為奇數，將多出來的 `m` 補到 `fraction` ，即 `ix0 += ix0` 的部分，而 `m >>= 1` 會自動略去奇數多出來的 `1`
$m=even，\sqrt{2^m\times y}=2^{m/2}\times \sqrt{y}$
$m=odd，\sqrt{2^m\times y}=2^{m/2}\times \sqrt{2y}$
```cpp=62
/* 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;
if (t <= ix0) {
s0 = t + r;
ix0 -= t;
q += r;
}
ix0 += ix0;
r >>= 1;
}
```
* $\sqrt{y}$ 的部分，若 $1≤y<4$ ，則 $1≤\sqrt{y}<2$ ，可從 $1$ 逐漸向 $\sqrt{y}$ 逼近
* 令 $f(x)=x^2-N$ ，透過尋找 $f(x)=0$ 的解，找到 $x=\sqrt{N}$ 的最近似解
:::info
數學推導：
* 令 $q_i$ 為 $\sqrt{y}$ 到小數點第 i 位的精確值，且 $i≥0$ ，則 $q_0=1$
* 透過 $(q_{i}+2^{-(i+1)})^2$ 逐漸逼近 $y$ :
* 若 $(q_i+2^{-(i+1)})^2≤y$，則 $q_{i+1}=q_i+2^{-(i+1)}$ ，即第 `i+1` 位補上 `1`
* 若 $(q_i+2^{-(i+1)})^2>y$，則 $q_{i+1}=q_i$ ，即第 `i+1` 位補上 `0`
* 定義以下變數來簡化算式 :
* $s_i = 2\times q_i$
* $y_i = 2^{i+1}\times (y-q_i^2)$
* 接著將上述 $(q_i+2^{-(i+1)})^2≤y$ 整理：
$\implies q_i^2+q_i\times 2^{-i}+2^{-2(i+1)}≤y$
$\implies q_i\times 2^{-i}+2^{-2(i+1)}≤y-q_i^2$
$\implies 2^{-(i+1)}(2\times q_i^2 + 2^{-(i+1)})≤y-q_i^2$
$\implies q_i\times 2^{-i}+2^{-(i+1)}≤(y-q_i^2)\times 2^{i+1}$
$\implies s_i+2^{-(i+1)}≤y_i$
* 若 $s_i+2^{-(i+1)}≤y_i$ 成立：
* $s_{i+1}=s_i+2^{-i}$
* $y_{i+1}=2\times(y_i-s_i-2^{-(i+1)})$
* 若 $s_i+2^{-(i+1)}≤y_i$ 不成立：
* $s_{i+1}=s_i$
* $y_{i+1}=2\times y_i$
:::
* 回顧程式碼， `ix0 += ix0` 的部分是為了增加精確度，使小數部分有 `24 bits` ，故將整數部分左移，而 `QQ0 = 0x01000000` ，則為了開始從 $2^0$ 開始向小數點迫近
* 回顧上面所定義的變數：
* $s_i = 2\times q_i$
* $y_i = 2^{i+1}\times (y-q_i^2)$
* 而後面 `while` 迴圈就是在執行以下計算：
:::info
$t=s_i+2^{-(i+1)}$
$if( t ≤ y_i):$ ，即 $(s_i+2^{-(i+1)})≤y_i$
$s_{i+1}=t+2^{-(i+1)}=s_i+2^{-(i+1)}+2^{-(i+1)}$ ，即 $s_{i+1}=s_i+2^{-i}$
$y_{i+1}=y_i-t$ ， 即 $y_{i+1}=y_i-s_i-2^{-(i+1)}$
$q_{i+1}=q_i+r$ ， 即 $q_{i+1}=q_i+2^{-(i+1)}$
$y_{i+1}=2\times y_{i+1}$ ，若 $if$ 成立，$y_{i+1}=2\times (y_i-s_i-2^{-(i+1)})$ ；反之， $y_{i+1}=2\times y_i$
:::
```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);
}
}
```
* 四捨五入所算出來的近似值 `q`
:::warning
TO DO
詳細解釋 rounding 的原理。
:::
```cpp=90
ix0 = (q >> 1) + QQ1;
ix0 += (m << QQ2);
INSERT_WORDS(z, ix0);
```
* 第 `90` 行的部分，將近似值 `q` 以精準浮點數表達
* 最後 `q` 是將小數點點在 `24 bit` 與 `25 bit` 間的二進位小數，因此 `q` 必須往右退一位，才能符合 IEEE 754 對浮點數的規範
* 加上 `QQ1` ，使得 `ix0` 得以正確轉換為浮點數，即 $ix0 = 2^0\times 1.fraction$， 因此在 `exponent` 這邊必須為 `0` ，故 `QQ1=0x3f000000`
* 再乘上一開始的 $2^{n/2}$ ，即 `ix0` 的 `exponent` 部分加上本來的 `m` ，而 `exponent` 部分從第 `24 bit` 開始，因此 `QQ2=23` ，將 `m` 左移至 `exponent` 位置
* 最後透過 `INSERT_WORDS` ，將精準浮點數編碼 `ix0` 轉換為 `float` 型態
### 延伸問題：
* 練習 LeetCode [69. Sqrt(x)](https://leetcode.com/problems/sqrtx/)，提交不需要 FPU 指令的實作
* 透過 `left` 與 `right` 兩變數，兩邊同時向 `Sqrt(x)` 迫近
* 根據 [算術-幾何平均值不等式](https://en.wikipedia.org/wiki/Inequality_of_arithmetic_and_geometric_means) 證明 $\dfrac{x+y}{2}≥ \sqrt{xy}$ ，因此將 right 設為 $\dfrac{x}{2}$ 逐漸遞減
* 若 $x=2^n\times y$ ，則 $\sqrt{x}=2^{n/2}\times \sqrt{y}$ ，因此 left 從 $2^{n/2}$ 開始遞增
* 由 `mid` 直接取 `left` 與 `right` 中間值做計算，達到時間複雜度 $O(log$ $n)$的效果
```graphviz
digraph
{
node[shape = record]
Example[label = "0 | 1 | ... | left | → | Sqrt(x) | ← | right | ... | x"]
}
```
```cpp=
int mySqrt(int x) {
if ( x <= 1 ) return x ;
uint32_t m = 32 - __builtin_clz(x) ;
m -= 1 ;
m >>= 1 ;
int l = 1 << m ;
int r = x/2 ;
while(l <= r) {
uint64_t mid = (l + r)>>1 ;
uint64_t check = mid*mid ;
if( check == x) return mid ;
else if ( check > x ) r = mid - 1 ;
else l = mid + 1 ;
}
return r ;
}
```
---
## 測驗 3
### 題目：
LeetCode [829. Consecutive Numbers Sum](https://leetcode.com/problems/consecutive-numbers-sum/)
```cpp
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;
}
```
請補完程式碼，使上述程式碼得以運作。
### 原理：
根據解題想法說明，若要寫成連續正整數之和，可寫成以下表示式：
$x+(x+1)+(x+2)+...+x+(k-1)=kx+\dfrac{(k-1)k}{2}=N$
整理後得到：
$kx=N-\dfrac{(k-1)k}{2}$ ，其中 $\dfrac{(k-1)k}{2}$ 即為 $1+2+3+...+k$ 之和
而 `x += ZZZ` 就是在做 $N-\dfrac{(k-1)k}{2}$ 的部分，而 $\dfrac{(k-1)k}{2}$ 為從 `1` 慢慢往上加到 `k` 個結果，因此 `ZZZ = 1-i` ，從 $N-1$ 慢慢遞往下減
如果 $kx=N-\dfrac{(k-1)k}{2}$ 等式成立，且 $k$ 、 $x$ 為整數，則 $\dfrac{N-\dfrac{(k-1)k}{2}}{k}$ 必也為整數，那麼就可用 `x % k = 0` 來檢查條件是否成立
Sign in with Wallet
Connect another wallet