# 2020q3 Homework5 (quiz5)
contributed by < `ptzling310` >
>[2020q3 第 5 週測驗題](https://hackmd.io/@sysprog/2020-quiz5)
## 測驗 `1`：浮點數除法
```c=
#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;
}
```
- [x] D1 = `2`
- [X] D2 = `1`
### 理解
#### 參數
- `double orig`： 被除數
- `int slots`： 除數
#### 程式
```c=4
if (slots == 1 || orig == 0)
return orig;
```
當以下 case 時，結果即為 `orig`：
- $\dfrac{orig}{slots} = \dfrac{orig}{1} = orig$
- $\dfrac{orig}{slots} = \dfrac{0}{1} = 0 =orig$
```c=7
int od = slots & 1;
```
`od` 用來記錄 `slots` 的 LSB 是否為 bit 1 。
若 `od` == 1 ，表示 `slots` 為奇數。
```c=8
double result = divop(orig / 2, od ? (slots + 1) >> 1 : slots >> 1);
```
- `od ? (slots + 1) >> 1 : slots >> 1`
- `od == 0` 表示 `slots` 為偶數，則運算 $\dfrac{orig}{slots} = \dfrac{\dfrac{orig}{2}}{\dfrac{slots}{2}}$，不會影響 `result` 。 故 D1 = `2`。
- `od == 1` 表示 `slots` 為奇數，則先利用 `+1` 把 `slots` 變成偶數 ($slots'$)再右移。但此舉會造成 `result` 比正確結果來的小 (因為分母變大了)。
```c=9
if (od)
result += divop(result, slots);
```
在 line 8 時，將 slots 做了 `+1` 的動作，所以回傳的 `result` 並不是真正的答案，這裡就是在進行修正。
舉例：
$\dfrac{1}{3}$ (已知答案為 $0.\overline{33}$)， 由於 3 不為偶數，故先將 `3+1`，
再求 $\dfrac{1}{4}=\dfrac{0.5}{2}=\dfrac{0.25}{1}=0.25$。(依照 line 8 )。
目前的 result = 0.25。
$result\times slot = 0.25 \times 3 = 0.75 \neq 1$。
故用 `result` 距離 $1$ 還差了 $1-0.75 = 0.25$ ，正好為目前 `result` 的值。
所以只要再補上 $\dfrac{0.25}{3}=0.125$。則 result = 0.25+0.125 = 0.375。
-----------------------------------------
## 測驗 `2`：浮點數平方根程式 (利用牛頓法找近似值)
### 程式
```c=
#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;
}
```
### 理解
#### IEEE 754 表示
float 表示式：$S\times1.Frac\times2^{127+E}$
```graphviz
digraph structs {
node[shape=record]
NaN[label="float",shape=plaintext];
struct1 [label="S|E|E|E|E|E|E|E|E|F|...|F"];
{rank= NaN; struct1}
}
```
:::spoiler 浮點數轉 IEEE格式 過程
6.25 = 110.1 x $2^0$
= 1.101 x $2^2$
所以 6.25 儲存的格式：
Sign(S) = (0)~2~ (因為正數)
Exponent(E) = 2 + 127 = 129 = (10000001)~2~
Frac(F) = (1010..0)~2~
:::
#### 概念
將 `float x` 表示為 $x=y\times2^{m}$，當要算 $\sqrt{x}$ 時，就計算$\sqrt{y}\times2^{m/2}$。其中 $2^{m/2}$ 用位移即可達成；而要處理的部分就剩下 $\sqrt{y}$ 的估算。
#### struct
```c=4
typedef union {
float value;
uint32_t v_int;
} ieee_float_shape_type;
```
`value` 是 float 型態，而 `v_int` 是以 `int` 型態保存。
故同一二進制的數值，如果是 value 則代表一個浮點數，如果是 v_int 用來表示整數。
#### function
##### INSERT_WORDS(d, ix0)
```c=9
/* 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` 型式
##### EXTRACT_WORDS(ix0, d)
```c=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)
```
將 `float` 轉成 `int` 型式。
##### ieee754_sqrt(float x)
```c=36
EXTRACT_WORDS(ix0, x);
```
`x` 為 `float` ，藉由 `EXTRACT_WORD` 轉為 `int` 型式的 `ix0`。
```c=43
/* 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 */
}
```
- `ix0 <= 0` ，也就是 `ix0` 的 MSB (b~31~) = 1：或者 `ix0` = (00..0)~2~。
- `(ix0 & (~sign))==0)`
`sign = 0x80000000` = (1000 ... 0000)~2~ , 則 `~sign = 0x80000000` = (0111 ... 1111)~2~。 = `0x7FFFFFFF`。取 `AND` 代表只取 Exponent(E)、Fraction(F) 的部分。
所以這裡條件若成立，表示`ix0` 為 (10...0)~2~ = +0 或 (0...0)~2~ = -0 。
又因 0 開根號還是 0，故 `retrun x`。
- `ix0 < 0`，若成立代表 `ix0` 的 MSB(b~31~) 為 1，表示 `x` 是負數的，但負數不能開根號。 而 `(x - x) / (x - x)` 會為 `0/0` ，歸類在 `Not a Numbe`。
```c=43
/* take care of +INF and NaN */
if ((ix0 & 0x7f800000) == 0x7f800000) {
/* sqrt(NaN) = NaN, sqrt(+INF) = +INF,*/
return x;
```
```graphviz
digraph structs {
node[shape=record]
NaN[label="NaN",shape=plaintext];
struct1 [label="0/1|11111111|nonZero"];
{rank= NaN; struct1}
INF[label="INF",shape=plaintext];
struct2 [label="0/1 | 11111111| 0..0"];
{rank= INF; struct2}
}
```
若 Exponent(E) 的部分皆為 1 ，則為 NaN 或 INF 的情況，故 `retrun x`及可。
```c=48
/* normalize x */
m = (ix0 >> 23);
if (m == 0) { /* subnormal = denormal x */
for (i = 0; (ix0 & 0x00800000) == 0; i++)
ix0 <<= 1;
m -= i - 1;
}
```
`m` 是保留 `ix0` 的 Sign(S) 、 Exponent(E)。
加上在前面已經過濾掉 MSB(b~31~) = 1 的 `ix0`，這裡如果要符合 denormal 的條件的話，就會為：
```graphviz
digraph structs {
node[shape=record]
NaN[label="denormal",shape=plaintext];
struct1 [label="0|00000000|nonZero"];
{rank= NaN; struct1}
}
```
- 若 `m = 0` 則符合 denormal 的條件 (S、E 皆 0 )。
若為 denormal，則 `float x` = $0.Fraciton\times2^{-126}$。
- `& 0x00800000` 就是在取第9個 bit (b~23~)，正為 Exponent 的最後一個 bit。
故 `ix0 & 0x00800000` 的目的是為了讓 Exponent 的位置產生 1。但這樣的移動會使得這個值變得不正確。
:::info
舉例：
| 0 | 00000000 | 0110...0 |
| -------- | -------- | -------- |
此為一 denormal float point，值為：$(2^{-2}+2^{-3})\times2^{-126}=2^{-128}+2^{-129}$
經由 line `50`~`52`，會將 Fraction 的最左 1 ，轉為以下：
| 0 | 00000001 | 1000...0 |
| -------- | -------- | -------- |
此為一 normal float point，值為：$(1+2^{-1})\times2^{-126}=2^{-126}\times2^{-127}$。
所以必定要進行修正，將值修為正確的值。
修正後值的計算：$(1+2^{-1})\times2^{1-127-3+1}=(1+2^{-1})\times2^{-128}=2^{-128}+2^{-129}$
※ 因為 for 的 `i++` ，所以會記錄到 `i=3`。
:::
- denormal float point 仍可開根號，為了讓這類型的數字可以一同套用同一個估算開根號的程式，先將該值轉為 normal 的形式 (就是 exponent 非全0)。
```c=55
m -= 127; /* unbias exponent */
```
由於儲存為 binary 的 expoenet 的部分是經由 `+127` 後所得，為得到表示式中的二的次方之值，將 `m - 127` 還原。此舉可以取的浮點數真正的指數值。
```c=56
ix0 = (ix0 & 0x007fffff) | 0x00800000;
```
取 `ix0` 的 Frac 以及將 **Exponent(E) 的** LSB 設 1。
如此可得$1.fraction$
```c=57
if (m & 1) { /* odd m, double x to make it even */
ix0 += ix0;
}
m >>= 1; /* m = [m/2] */
```
`m` 有以下兩種情況：
- $m$ 為偶數： $x = y\times2^{m},\sqrt{x}=\sqrt{y\times2^{m}}=\sqrt{y}\times2^{m/2}$
- $m$ 為奇數：$x=y\times 2^{m} , \sqrt{x} = \sqrt{y\times2\times2^{m-1}}=\sqrt{2y\times2^{m-1}}=\sqrt{2y}\times2^{m-1/2}$
```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 */
while (r != 0) {
t = s0 + r;
if (t <= ix0) {
s0 = t + r;
ix0 -= t;
q += r;
}
ix0 += ix0;
r >>= 1;
}
```
> 參考 <[RinHizakura](https://hackmd.io/@RinHizakura/SJnXiSewD#%E6%B8%AC%E9%A9%97-2)>、<[zzzxxx00019](https://hackmd.io/@jT29vQ4-SxmlBU5UMj1TGA/quiz5#%E6%B8%AC%E9%A9%97-2)>
此段程式就是在處理 $\sqrt{y}$ 的估算。
假設：$1\leq y<4 \Rightarrow 1 \leq \sqrt{y}<2$
:::info
為何 $y$ 的範圍介於 $[ 1 , 4 )$?
這裡目標是將 $\sqrt{x} = \sqrt{y} \times 2^{m/2}$ 求出，
目前已在 line `55` 找出 $m$ ， line `56` 將 `ix0` = $y$ 找出。
又 $y$ 為 $1.fracion$ ，必定在區間 $[ 1 , 2)$，
在 line `58`，若 $m$ 為奇數，則會提出一個 2 乘入，
所以 $y$ 的範圍會改為： $[ 1 , 4 )$。
:::
:::spoiler 如何逼近
假設現在要從 (1)~2~ 逼近 (1.011)~2~，
- (1.1)~2~ > (1.011)~2~，故小數點後一位 bit 應設 0 。
- (1.01)~2~ < (1.011)~2~，故小數點後兩位 bit 應設 1 。
- (1.011)~2~= (1.011)~2~，已相等，完成。
:::info
- $q_i$：逼近 $\sqrt{y}$ 至小數點後 $i$ 位。$q_0 = 1$ ，我們從 1 開始逼近。
- $(q_i+2^{-(i+1)})^{2}$ 來逼近 $y$。目前已經判斷到小數點後第 $i$ 位 ($q_i$)，則就在往後判斷一位 ($2^{-(i+1)}$) 要為 0 或 1。
- 規則：從上點 "如何逼近" 便可理解，要做的事情就是判斷該位子該填 bit 0 或 bit 1。
- $(q_i+2^{-(i+1)})^{2} \leq y$，代表小數點後第 `i+1`個位子要補上 bit 1，則 $q_{i+1} = q_i +2^{-(i+1)}$。
- $(q_i+2^{-(i+1)})^{2} > y$，代表小數點後第 `i+1` 個位子要補上 bit 0，則 $q_{i+1} = q_i$。
- 算式整理
$(q_i+2^{-(i+1)})^{2} \leq y$
展開上式 $\Rightarrow (q_i)^{2}+2\times q_i \times 2^{-(i+1)}+2^{-2(i+1)} \leq y$
移項$(q_i)^2$ $\Rightarrow 2\times q_i \times 2^{-(i+1)}+2^{-2(i+1)} \leq y-(q_i)^{2}$
提出 $2^{-(i+1)}$ $\Rightarrow 2^{-(i+1)}\times(2\times q_i+2^{-(i+1)}) \leq y-(q_i)^2$
同乘$2^{i+1}$ $\Rightarrow (2\times q_i+2^{-(i+1)}) \leq (y-(q_i)^2)\times 2^{i+1}$
令 $s_i=2\times q_i, y_i=(y-(q_i)^2)\times 2^{i+1}$
$\Rightarrow s_i+2^{-(i+1)} \leq y_i$
而汰換成此型式可以利用遞迴的方式求得 $s_i$、$y_i$
- 求 $s_i$、$y_i$
- $s_i+2^{-(i+1)} \leq 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$
:::
:::spoiler s_i、y_i推導
- $s_i+2^{-(i+1)} \leq y_i$
$\begin{aligned}s_{i+1} &= 2 \times q_{i+1} \\
&=2\times(q_i+2^{-(i+1)}) \\
&=2\times q_i + 2\times 2^{-(i+1)}\\
&=s_i+2^{-i}
\end{aligned}$
$\begin{aligned} y_{i+1} &= (y-({q_{i+1})}^{2})\times 2^{(i+1)+1}\\
&=(y-q_i^{2}-2^{-(i+1)}\times s_i-2^{-2(i+1)})\times 2^{i+2}\\
&=(y-q_i^{2})\times2^{i+2}-2^{i+2}\times2^{-(i+1)}\times s_i-2^{i+2}\times2^{-2(i+1)}\\
&=(y-q_i^{2})2^{i+1}\times2-2s_i-2^{-i}\\
&=2y_i-2s_i-2^{-i}\\
&=2\times(y_i-s_i-2^{-(i+1)})
\end{aligned}$
- $s_i+2^{-(i+1)} > y_i$
$\begin{aligned}s_{i+1} &= 2\times q_{i+1}\\
&=2\times q_i\\
&=s_i
\end{aligned}$
$\begin{aligned} y_{i+1} &= (y-({q_{i+1})}^{2})\times 2^{(i+1)+1}\\
&=(y-q_i^{2})\times 2^{(i+2)}\\
&=(y-q_i^{2})2^{i+1}\times 2\\
&= 2y_{i}
\end{aligned}$
:::
整理目前所知資訊：
- 定義 $s_i=2q_i$，且可使用遞迴求解
$s_{i+1}= \left\{\begin{array}{lr}
s_i+2^{-i}, s_i+2^{-(i+1)} \leq y_i\\
s_i , s_i+2^{-(i+1)} > y_i\\
\end{array} \right.$
- 定義 $y_i=(y-(q_i)^2)\times 2^{i+1}$，且可使用遞迴求解
$y_{i+1}= \left\{\begin{array}{lr}
2\times(y_i-s_i-2^{-(i+1)}), s_i+2^{-(i+1)} \leq y_i\\
2\times y_i , s_i+2^{-(i+1)} > y_i\\
\end{array} \right.$
Line `63` ：在 line `78` 開始有進行進位判斷，所以這裡將 `ix0` 左移一位，以使小數部分可以多一位來進行進位(現在小數部分有 24 bits了)。
Line `65` ： `r` 的功能是用來一個 bit 一個 bit 的檢查，首先小數部分要先保留 24 bits (b~0~-b~23~)，第 25 bit (b~24~)才為整數之部分，所以從第 25 bit (b~24~) 開始檢查，故 QQ0 = `0x01000000` ，從 $2^0 = 1$ 開始逼近。
Line `68` ： `t = s0 + r`。即：`t` = $s_i+2^{-(i+1)}$。
從 line `69`~`72` 處理 $= s_i+2^{-(i+1)} \leq y_i$ 之情況。
Line `69` ： `if (t <= ix0)`。即：若`t` $= s_i+2^{-(i+1)} \leq y_i$
Line `70` ： `s0 = t + r;`。即： $s_{i+1}=s_i+2^{-(i+1)}+2^{-(i+1)}=s_i+2^{-i}$
Line `71` ： `ix0 -= t;`。即： $y_{i+1}=y_i-(s_i+2^{-(i+1)})=y_i-s_i-2^{-(i+1)}$
Line `72` ： `q += r;`。即：$q_{i+1}=q_i+2^{-(i+1)}$
Line `74` ： `ix0 += ix0;`。即 $y_i=2\times2y_i$。在$y_{i+1}$定義中，不論條件為何，皆有一個$\times 2$這個動作，就在此處完成。如此對$y_{i+1}$的處理才算完整。
Line `75` ： `r >>= 1;`。即查看下一個 bit (從$2^{-(i+1)}$ 變 $2^{-(i+2)}$)。
```c=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);
}
}
```
Line `80` ：照理來說 `z` 是一個比 1 還略小的數字(0.9....99)。在 line `81` 會判斷說這個 `z` 是否有小於 1，若無則在表這個系統有進行進位。
Line `82` ： 照理來說 `z` 是一個比 1 還略大的數字(1.0...01)。在 line `83` 會判斷說這個 `z` 是否有大於 1，若有則表示此系統採無條件進位。
目前第 1 bit (b~0~) 是用來判斷進位用的，由於這裡是無條件進位，藉由 `+2` = +(10)~2~ 將 1 加到 b~1~。
Line `86` ：此處表示 `z` = 1，也就是此系統是採用四捨五入，唯有第 1 個 bit(b~0~) = 1 才會進位。
```c=90
ix0 = (q >> 1) + QQ1;
ix0 += (m << QQ2);
INSERT_WORDS(z, ix0);
```
- `(q >> 1)` 的目的是要將第 1 個 bit(b~0~) 消去，因為這個 bit 是多出來做進位判斷用的。
- 由於浮點數的 Exponent(E) 的必須 `+127`，但在前面我們已將 Exponent(E) 取 1，所以這裡在補上 126，故 QQ1 = `0x3F000000`。
- 將 `m` 補回正確位置，故 QQ2 = `23`。
- 最後在 line `93` 轉回 `float` 型式。
------------------
## 測驗`3`： [LeetCode 829. Consecutive Numbers Sum ](https://leetcode.com/problems/consecutive-numbers-sum/)
### 程式
```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;
}
```
- [x] ZZZ=`1-i`
### 理解
```c=5
int ret = 1;
```
`ret` 用來記錄組合數。
```c=7
for (int i = 2; i < x; i++) {
x += ZZZ;
```
$\begin{aligned}
令N
&= x + (x+1) + ... + [x+(k-1)] \\
&= kx\times(1+2+...+(k-1)) \\
&= kx + \dfrac{k\times(k-1)}{2} \\
整理後得： \\
kx &= N - \dfrac{k\times(k-1)}{2} …(1)
\end{aligned}$
由於$k$是給定任意值，所以是會對應到變數`i`。
而在 line `6` 中將 `N` 指派給 `x`，所以 `x` 應用來紀錄 $N-\dfrac{k\times(k-1)}{2}$。
`x += ZZZ;` 就是在處理 $-\dfrac{k\times(k-1)}{2}$。
$\begin{aligned}
-\dfrac{k\times(k-1)}{2}
&= - [1+2+...+(k-1)]\\
&= [-1-2-...-(k-1)]\\
\end{aligned}$
每次都在處理 $-(k-1)=-k+1=1-k$， 又這裡的 $k$ 就是 `i`，所以 ZZZ = `1-i`。
- 為何不從 i = 1 開始？
由於在 line `5`：`ret=1`，若這裡 `i=1` ，則 $x=N$，會在將 `ret` 又多算一次。
```c=9
if (x % i == 0)
ret++;
```
回顧 $算式 (1)$，如果可得一$x$為整數解，則代表可以對應到一個整數組合。
$kx = N - \dfrac{k\times(k-1)}{2}\Rightarrow x= \dfrac{N - \dfrac{k\times(k-1)}{2}}{k} 為整數$。
所以只要確認 $x$ 是否為 $k$ 之倍數就可判斷這組整數組合是否存在。
###### tags: `sysprog2020`
Sign in with Wallet
Connect another wallet