# 2018q1 Homework2 (assessment)
contributed by < `TingL7` >
## 第 1 週 測驗 `一`
考慮以下程式碼:
```clike
#include <stdlib.h>
int isMultN(unsigned int n) {
int odd_c = 0, even_c = 0; /* variables to count odd and even SET bits */
if (n == 0) // return true if difference is 0.
return 1;
if (n == 1) // return false if the difference is not 0.
return 0;
while (n) {
if (n & 1) // odd bit is SET, increment odd_C
odd_c++;
n >>= 1;
if (n & 1) // even bit is SET, increment even_c
even_c++;
n = n >> 1;
}
/* Recursive call till you get 0/1 */
return(isMultN(abs(odd_c - even_c)));
}
```
其作用為檢查輸入整數是否為 N 的倍數，那麼 N 為多少？
---
### 解題想法&思考
* 程式碼解讀
* 由這兩段程式碼，可以得知這是一個遞迴程式，結束條件是 odd_c 和 even_c 相等或相差 1 ，且相等時，為N的倍數。
```clike
if (n == 0) // return true if difference is 0.
return 1;
if (n == 1) // return false if the difference is not 0.
return 0;
```
`return(isMultN(abs(odd_c - even_c)));`
* 而由這段程式碼可以知道 odd_c 及 even_c 分別在計算奇數位及偶數位為 1 的個數
```clike
while (n) {
if (n & 1) // odd bit is SET, increment odd_C
odd_c++;
n >>= 1;
if (n & 1) // even bit is SET, increment even_c
even_c++;
n = n >> 1;
}
```
* 綜合上述，並簡單化這個題目，在程式只有進行一次遞迴時，奇數位及偶數位總合相等時，即為 N 的倍數。將選項帶入有此特性的只有 3(11~2~) 。
### 延伸問題: 為什麼上述程式碼可運作呢？
>數字與文字間也要空白喔！
>[name=課程助教][color=red]
>>好的，已修改
>>[name=TingL7]
* 程式碼可以運作，是依據 3 的倍數，在二進位的情況下的特性，即奇數位數字和及偶數位數和相等或相差 3 的倍數。因此，只要證明 3 的倍數有此特性，即可確認上述程式碼是運作良好的。
* 這個特性與十進位中的 11 很像。參考[ 100 以內的質數倍數判別法及其延伸](http://www.shs.edu.tw/works/essay/2015/11/2015111502175190.pdf)中， 11 的倍數判別法，可寫出以下證明:
設有一 (n+1) 位數 N = (a~n~a~n-1~..a~1~a~0~)~2~ ，其中 a~n~$\neq$ 0
令 k 為 $\le n$ 的正整數
\begin{align}
a_n &=\sum_{i=0}^na_i\times2^i\\
& = [a_0+a_2(3+1)+a_4(15+1)+...+a_{2k}(2^{2k}-1+1)+...]\\
& +[a_1(3-1)+a_3(9-1)+a_5(33-1)+...+a_{2k+1}(2^{2k+1}+1-1)+...]\\
& = (a_0+a_2+a_4+...+a_{2k}+...)\\
& +[a_2(3)+a_4(15)+...+a_{2k}(2^{2k}+1)+...]\\
& - (a_1+a_3+a_5+...+a_{2k+1})\\
& + [a_1(3)+a_3(9)+a_5(33)+...+a_{2k+1}(2^{2k+1}+1)+...]...第一式
\end{align}
利用歸納證明法，可證明 2^2n^-1 及 2^2n+1^+1 為 3 的倍數:
n=1:則 2^2^-1=3 為 3 的倍數
令 n=k 成立，即 2^2k^-1=3x ， x 可為任一正整數
n=k+1:
2^2(k+1)^-1=2^2k+2^-1=4\*2^2k^-1=4(2^2k^-1)+3=4(3x)+3=3(4x+1)
得證: 2^2n^-1 為 3 的倍數
因此， [a~2~(3)+a~4~(15)+...+a~2k~(2^2k^+1)+...] 為 3 的倍數...**第二式**
n=0:則 2^1^-1=3 為 3 的倍數
令 n=k 成立，即 2^2k+1^+1=3x ， x 可為任一正整數
n=k+1:
2^2(k+1)+1^+1=2^2k+3^+1=4\*2^2k+1^+1=4(2^2k+1^+1)-3=4(3x)-3=3(4x-1)
得證: 2^2n+1^-1 為 3 的倍數
因此， [a~1~(3)+a~3~(9)+a~5~(33)+...+a~2k+1~(2^2k+1^+1)+...] 為 3 的倍數...**第三式**
由第一、二、三式可得知:
如果 (a~0~+a~2~+a~4~+...+a~2k~+...)-(a~1~+a~3~+a~5~+...+a~2k+1~+...) 為 3 的倍數，則 N 必為 3 的倍數。
因此，要判定二進位中，是否為 3 的倍數，可由 奇數位數字和及偶數位數和相等或相差 3 的倍數 的特性去判斷。
### 延伸問題: 將 N 改為其他質數再改寫為對應的程式碼，一樣使用遞迴
* 要將 N 改為其他質數，要先找出其他質數是否類似 3 的特殊判斷特性。
* N=5
* 參考[從二進位判斷數字是否被 5 整除](https://www.ptt.cc/bbs/Prob_Solve/M.1223098109.A.99F.html)， N 是否為 5 的倍數判定，可由 (N>>2)-(N&3) 是否為 5 的倍數決定。
* 證明:
(N>>2)-(N&3) 為 N 除以 4 的商數 (a) 減去 N 除以 4 的餘數 (b)。寫成算式為
$N\div4=a...b\\
=>N =a\times4+b=a\times5-a+b\\
=>如果 b-a 為 5 的倍數，則 N 為 5 的倍數$
* 程式碼
```clike
#include <stdlib.h>
int isMultN(unsigned int n) {
if (n == 0|| n == 5) // return true if n is 0 or 5.
return 1;
if (n < 5) // return false n is smaller than 5
return 0;
n = (n >> 2) + (n & 3);
/* Recursive call till n <= 5 */
return(isMultN(n));
}
```
* N=7
* 由 N=5 的證明，是推廣到 N=7 ，甚至可以說是任何數字都可以用，只要找出與 2 的次方的關係。
* $N =a\times8+b=a\times7+(a+b)\\
=>如果 a+b 為 7 的倍數，則 N 為 7 的倍數$
* 程式碼
```clike
#include <stdlib.h>
int isMultN(unsigned int n) {
if (n == 0|| n == 7) // return true if n is 0 or 7.
return 1;
if (n < 7) // return false n is smaller than 7.
return 0;
n = (n >> 3) + (n & 7);
/* Recursive call till n <= 5 */
return(isMultN(n));
}
```
* N=3
* 由上述推得的方法，可以將 N=3 改寫。
* 改寫一:
$N =a\times4+b=a\times3+(a+b)\\
=>如果 a+b 為 3 的倍數，則 N 為 3 的倍數$
* 程式碼
```c
#include <stdlib.h>
int isMultN(unsigned int n) {
if (n == 0|| n == 3) // return true if n is 0 or 3.
return 1;
if (n < 3) // return false n is smaller than 3.
return 0;
n = (n >> 2) + (n & 1);
/* Recursive call till n <= 5 */
return(isMultN(n));
}
```
* 改寫二:
$N =a\times2+b=a\times3+(-a+b)\\
=>如果 b-a 為 3 的倍數，則 N 為 3 的倍數$
* 程式碼
```c
#include <stdlib.h>
int isMultN(unsigned int n) {
if (n == 0) // return true if n is 0.
return 1;
if (n == 1) // return false n is 1.
return 0;
n = (n >> 1) - (n & 1);
/* Recursive call till n <= 5 */
return(isMultN(n));
}
```
* 通用
* 其實，上面的方法是可以通用的，不論是否為質數，只要算清楚與 2 的次方差多少，再利用 shift 及 + 或 - 就可以判斷倍數。
## 第 2 週 測驗 `一`
請完成下方程式碼，依循 IEEE 754 單精度規範，輸出 $2^x$ 的浮點數表達法，而且考慮兩種狀況：
- 當 $x$ 過小的時候，回傳 $0.0$
- 當 $x$ 過大的時候，回傳 $+∞$
注意：這裡假設 `u2f` 函式返回的浮點數值與其無號數輸入有著相同的位數，也就是至少 32-bit
```c
#include <math.h>
static inline float u2f(unsigned int x) { return *(float *) &x; }
float exp2_fp(int x) {
unsigned int exp /* exponent */, frac /* fraction */;
/* too small */
if (x < 2 - pow(2, Y0) - 23) {
exp = 0;
frac = 0;
/* denormalize */
} else if (x < Y1 - pow(2, Y2)) {
exp = 0;
frac = 1 << (unsigned) (x - (2 - pow(2, Y3) - Y4));
/* normalized */
} else if (x < pow(2, Y5)) {
exp = pow(2, Y6) - 1 + x;
frac = 0;
/* too large */
} else {
exp = Y7;
frac = 0;
}
/* pack exp and frac into 32 bits */
return u2f((unsigned) exp << 23 | frac);
}
```
---
### 解題想法&思考
* 要解這題，需要充分了解浮點數的邊界判定。
* 浮點數 = (-1)^s^ \* 2^E^ \* M
* bias = 2^8-1^-1 = 127
* denormalized:(0~1)
* s: +/-，本題中沒有負的
* E: 1-bias= -126 = 2-2^7^
* M: 0.f~22~f~21~f~20~...f~0~
* normalized:(1~max)
* s: +/-，本題中沒有負的
* E: e-bias = e+1-2^7^
* M: 1.f~22~f~21~f~20~...f~0~
* too small
```
if (x < 2 - pow(2, Y0) - 23) {
exp = 0;
frac = 0;
```
* 判斷d enormalized 可以表達的最小數值，在單精度的情況下為
`0 00000000 00000000000000000000001` 即為 $(-1)^0 \times 2^{2-2^{7}} \times 2^{-23}=2^{2-2^{7}-23}$
因此， ++Y0=7++
* x 太小，回傳 0.0 。
* denormalized
```
} else if (x < Y1 - pow(2, Y2)) {
exp = 0;
frac = 1 << (unsigned) (x - (2 - pow(2, Y3) - Y4));
```
* 判斷 denormalized 可以表達的最大數值為
`0 00000000 11111111111111111111111` 即為 $(-1)^0 \times 2^{2-2^{7}} \times \sum_{i=1}^{23}2^{-i}=2^{2-2^{7}} \times (1-2^{-23})$
* (1-2^-23^) 取近似值為 1 ，因此 ++Y1=2++ ， ++Y2=7++ 。
* $2^x = 2^{2-2^7}\times(frac\times2^{-23})$ ，則 $frac=2^{x-(2-2^7-23)}$
因此， ++Y3=7++ ， ++Y4=23++ 。
* normalized
```
} else if (x < pow(2, Y5)) {
exp = pow(2, Y6) - 1 + x;
frac = 0;
```
* 判斷 normalized 可以表達的最大數值為
`0 11111110 11111111111111111111111` 即為 $(-1)^0 \times 2^{2^8-1-1-(2^7-1)} \times (1+1-2^{-23})=2^{2^7-1}\times (2-2^{-23})=2^{2^7}-2^{104}$
* 最大可表示的值不會超過 $2^x$ 為 $2^{2^7}$ ，因此 ++Y5=7++ 。
* 由於要表達 2 的次方，不需要小數點後的表示，因此 frac=0 。
* $2^x = 2^{exp-bias}=2^{exp-(2^7-1)}$ ，則 $exp=x+(2^7-1)$
因此， ++Y6=7++ 。
* too large
```
} else {
exp = Y7;
frac = 0;
}
```
* x 太大，回傳 $+∞$， INF 表示 exp 為 `11111111` = `0xFF` ， frac=0 。
因此， ++Y7=0xFF++ 。
### 延伸問題: 給定一個單精度的浮點數值，輸出較大且最接近的 $2^x$ 值，需要充分考慮到邊界
* 邊界
* denormalized 最小值:
`0 00000000 00000000000000000000001` = $(-1)^0 \times 2^{2-2^{7}} \times 2^{-23}=2^{2-2^{7}-23}$
* denormalized 最大值:
`0 00000000 11111111111111111111111` = $(-1)^0 \times 2^{2-2^{7}} \times \sum_{i=1}^{23}2^{-i}=2^{2-2^{7}} \times (1-2^{-23}) = 2^{log_2(2^{-126}-2^{-149})}$
* normalized 最大值:
`0 11111110 11111111111111111111111` = $(-1)^0 \times 2^{2^8-1-1-(2^7-1)} \times (1+1-2^{-23})=2^{2^7-1}\times (2-2^{-23})=2^{log_2(2^{128}-2^{104})}$
* exp&frac
* too small
* 為 0
* exp = 0
* frac = 0
* denormalized
* exp = 0
* $2^x = 2^{2-2^7}\times(frac\times2^{-23})$，則$frac=2^{x-(2-2^7-23)}$
* normalized
* $2^x = 2^{exp-bias}\times(1+frac\times 2^{-23})=2^{exp-127}\times (1+frac\times 2^{-23})$
* $exp=127+x(整數部分)$
* $frac=(2^{x(小數部分)}-1)\times 2^{23}$
* too large
* INF
* exp = 0xFF
* frac = 0
* 程式碼
```c
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
static inline float u2f(unsigned int x) { return *(float *) &x; }
float exp2_fp(float x) {
unsigned int exp /* exponent */, frac /* fraction */;
float frac_f;
/* too small */
if (x < 2 - pow(2, 7) - 23) {
exp = 0;
frac = 0;
/* denormalize */
} else if (x < log(pow(2, -126) - pow(2, -149))/log(2)) {
exp = 0;
frac_f = pow(2.0, (x + 149));
frac = ((int)(frac_f*2)&1)?(int)frac_f+1:(int)frac_f;//向上捨入
/* normalized */
} else if (x < log(pow(2, 128) - pow(2, -104))/log(2)) {
exp = 127 + (int)((x<0)?x-0.99:x);//x為負數時，無條件進位
frac_f = ((pow(2.0,x-(float)((int)((x<0)?x-0.99:x)))-1.0)*(1<<23));
frac = ((int)(frac_f*2)&1)?(int)frac_f+1:(int)frac_f;//向上捨入
/* too large */
} else {
exp = 0xFF;
frac = 0;
}
/* pack exp and frac into 32 bits */
return u2f((unsigned) exp << 23 | frac);
}
```
* 執行結果
```
2^-150
exp2_fp():0.000000000000000000000000000000000000000000000000000000000000
powf():0.000000000000000000000000000000000000000000000000000000000000
2^-149.5
exp2_fp():0.000000000000000000000000000000000000000000000000000000000000
powf():0.000000000000000000000000000000000000000000001401298464324817
2^-127
exp2_fp():0.000000000000000000000000000000000000005877471754111437500000
powf():0.000000000000000000000000000000000000005877471754111437500000
2^-126.3
exp2_fp():0.000000000000000000000000000000000000009547961485342182800000
powf():0.000000000000000000000000000000000000009547961485342182800000
2^-126
exp2_fp():0.000000000000000000000000000000000000011754943508222875000000
powf():0.000000000000000000000000000000000000011754943508222875000000
2^127
exp2_fp():170141183460469230000000000000000000000.000000000000000000000000000000000000000000000000000000000000
powf():170141183460469230000000000000000000000.000000000000000000000000000000000000000000000000000000000000
2^127.5
exp2_fp():240615965050037600000000000000000000000.000000000000000000000000000000000000000000000000000000000000
powf():240615965050037600000000000000000000000.000000000000000000000000000000000000000000000000000000000000
2^128
exp2_fp():1.#INF00000000000000000000000000000000000000000000000000000000
powf():1.#INF00000000000000000000000000000000000000000000000000000000
2^0.5
exp2_fp():1.414213538169860800000000000000000000000000000000000000000000
powf():1.414213538169860800000000000000000000000000000000000000000000
```
* 結果探討
:::info
- 可以從結果看出，基本上與 powf() 的結果相同，除了 -149.5 ，為何 -149~-150 之間仍可以表示?
- (2^0.5^)^2^=2，但 1.41421353816986080^2^=1.99999993154 ，比實際值略小，如何輸出較大且最接近的 $2^x$ 值?
:::
:::warning
這發現很重要！你應該要能從數學的角度計算出上述程式碼和 powf(x, 0.5) 以及理論值的差異，在 IEEE 754 中，已有明確規範。
另外，請改善程式碼縮排，一率用 4 個 space，而非 tab，務必讓版面簡潔清爽。
:notes: jserv
:::
## 第 3 週 測驗 `一`
延續 [從 2 – √ 2 的運算談浮點數](https://hackmd.io/s/ryOLwrVtz)，假設 double 為 64-bit 寬度，考慮以下符合 IEEE 754 規範的平方根程式，請補上對應數值，使其得以運作：
```clike=
#include <stdint.h>
/* A union allowing us to convert between a double and two 32-bit integers.
* Little-endian representation
*/
typedef union {
double value;
struct {
uint32_t lsw;
uint32_t msw;
} parts;
} ieee_double_shape_type;
/* Set a double from two 32 bit ints. */
#define INSERT_WORDS(d, ix0, ix1) \
do { \
ieee_double_shape_type iw_u = { \
.parts.msw = ix0, \
.parts.lsw = ix1, \
}; \
(d) = iw_u.value; \
} while (0)
/* Get two 32 bit ints from a double. */
#define EXTRACT_WORDS(ix0, ix1, d) \
do { \
ieee_double_shape_type ew_u; \
ew_u.value = (d); \
(ix0) = ew_u.parts.msw; \
(ix1) = ew_u.parts.lsw; \
} while (0)
static const double one = 1.0, tiny = 1.0e-300;
double ieee754_sqrt(double x)
{
double z;
int32_t sign = 0x80000000;
uint32_t r, t1, s1, ix1, q1;
int32_t ix0, s0, q, m, t, i;
EXTRACT_WORDS(ix0, ix1, x);
/* take care of INF and NaN */
if ((ix0 & KK1) == KK2) {
/* sqrt(NaN) = NaN, sqrt(+INF) = +INF, sqrt(-INF) = sNaN */
return x * x + x;
}
/* take care of zero */
if (ix0 <= 0) {
if (((ix0 & (~sign)) | ix1) == 0)
return x; /* sqrt(+-0) = +-0 */
if (ix0 < 0)
return (x - x) / (x - x); /* sqrt(-ve) = sNaN */
}
/* normalize x */
m = (ix0 >> 20);
if (m == 0) { /* subnormal x */
while (ix0 == 0) {
m -= 21;
ix0 |= (ix1 >> 11);
ix1 <<= 21;
}
for (i = 0; (ix0 & 0x00100000) == 0; i++)
ix0 <<= 1;
m -= i - 1;
ix0 |= (ix1 >> (32 - i));
ix1 <<= i;
}
m -= KK3; /* unbias exponent */
ix0 = (ix0 & 0x000fffff) | 0x00100000;
if (m & 1) { /* odd m, double x to make it even */
ix0 += ix0 + ((ix1 & sign) >> 31);
ix1 += ix1;
}
m >>= 1; /* m = [m/2] */
/* generate sqrt(x) bit by bit */
ix0 += ix0 + ((ix1 & sign) >> 31);
ix1 += ix1;
q = q1 = s0 = s1 = 0; /* [q,q1] = sqrt(x) */
r = 0x00200000; /* 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 + ((ix1 & sign) >> 31);
ix1 += ix1;
r >>= 1;
}
r = sign;
while (r != 0) {
t1 = s1 + r;
t = s0;
if ((t < ix0) || ((t == ix0) && (t1 <= ix1))) {
s1 = t1 + r;
if (((t1 & sign) == sign) && (s1 & sign) == 0)
s0 += 1;
ix0 -= t;
if (ix1 < t1)
ix0 -= 1;
ix1 -= t1;
q1 += r;
}
ix0 += ix0 + ((ix1 & sign) >> 31);
ix1 += ix1;
r >>= 1;
}
/* use floating add to find out rounding direction */
if ((ix0 | ix1) != 0) {
z = one - tiny; /* trigger inexact flag */
if (z >= one) {
z = one + tiny;
if (q1 == (uint32_t) 0xffffffff) {
q1 = 0;
q += 1;
} else if (z > one) {
if (q1 == (uint32_t) KK4)
q += 1;
q1 += 2;
} else
q1 += (q1 & 1);
}
}
ix0 = (q >> 1) + 0x3fe00000;
ix1 = q1 >> 1;
if ((q & 1) == 1)
ix1 |= sign;
ix0 += (m << KK5);
INSERT_WORDS(z, ix0, ix1);
return z;
}
```
---
### 解題想法與思考
* 程式碼解讀
```clike=6
typedef union {
double value;
struct {
uint32_t lsw;
uint32_t msw;
} parts;
} ieee_double_shape_type;
```
* 建立union，為了能夠分別操作前後32bit，如下所示。
|bit|63|62 ~~~ 52|51 ~~~~~~ 32|31 ~~~~~~~~~~ 0|
|--|:---:|:---:|:---:|:---:|
|value|s|exp|frac|frac|
|parts|msw|msw|msw|lsw|
| |31|30 ~~~ 20|19 ~~~~~~ 0|31 ~~~~~~~~~~ 0|
```clike=14
/* Set a double from two 32 bit ints. */
#define INSERT_WORDS(d, ix0, ix1) \
do { \
ieee_double_shape_type iw_u = { \
.parts.msw = ix0, \
.parts.lsw = ix1, \
}; \
(d) = iw_u.value; \
} while (0)
/* Get two 32 bit ints from a double. */
#define EXTRACT_WORDS(ix0, ix1, d) \
do { \
ieee_double_shape_type ew_u; \
ew_u.value = (d); \
(ix0) = ew_u.parts.msw; \
(ix1) = ew_u.parts.lsw; \
} while (0)
```
* 這兩段 define 主要功能是使 union 中的 double 和 parts 互相轉換。
* #define在做的事，基本上就是代換，將程式其他地方出現 INSERT_WORDS(d, ix0, ix1) ，代換成下面一段程式碼，不做任何修改。通常在使用會以 `INSERT_WORDS(d, ix0, ix1);` 的樣子出現，因此在代換的程式碼中，最後不會在加 `;`。
* 因為直接代換的關係，可能會產生 dangling else ，所以需要加上 `do{...}while(0)`。
* dangling else : 在使用 if-then(-else) 時，產生語意不清的情況，如:
`if A then if B then C else D` ，會有兩種解讀:
`if A then (if B then C)else D` 或
`if A then (if B then C else D)`
* [Linux Kernel 巨集 do{...}while(0) 的撰寫](http://www.runpc.com.tw/content/content.aspx?id=108451)中，舉了一個很好的例子，將其簡化來看:
```c
#define func(A, B) \
{ \
if(A) \
B = 1; \
else \
B = 0; \
} \
int funcA(int A,int B){
if(B)
func(A, B);
else
printf("Error B");
return B;
}
```
將 func(A,B) 用define的內容代換:
```c
int funcA(int A,int B){
if(B)
{
if(A)
B = 1;
else
B = 0;
};
else
printf("Error B");
return B;
}
```
可以發現一個很尷尬的地方，else 前面有 `;` ，表示if在else之前就做完了， else 是多的，語意與原本要表達的意思完全不同，這就是dangling else的情況。
若加上 `do{...}while(0)` 即可避免這樣的情況。
```c
int funcA(int A,int B){
if(B)
do{
if(A)
B = 1;
else
B = 0;
}while(0);
else
printf("Error B");
return B;
}
```
```clike=44
/* take care of INF and NaN */
if ((ix0 & KK1) == KK2) {
/* sqrt(NaN) = NaN, sqrt(+INF) = +INF, sqrt(-INF) = sNaN */
return x * x + x;
}
```
* 考慮 INF 及 NaN 的情況，要判斷是否為這兩者，只要看exp是否全為 1 ，對 ix0 來說， bit 20~30 的位置為 exp ，因此如果是 INF 或 NaN ， ix0 必為
`x 11111111111 xxxxxxxxxxxxxxxxxxx` ，需要在意的只有 exp 那一段，因此 ++KK1 = `0 11111111111 0000000000000000000` = 0x7FF00000++
* 且， exp 必須全為 1 ，因此， ++KK2 = 0x7FF00000++ 。
* 負數開根號為 sNaN，因此，回傳值不能直接回傳 x 。要用-INF得到sNaN，要先了解INF的運算:
* +INF + +INF = +INF
* -INF + -INF = -INF
* +INF + -INF = +INF - +INF = -INF - -INF = sNaN
* +INF * +INF = +INF
* -INF * -INF = +INF
* +INF * -INF = -INF
基本上，要得到sNaN，就必須使 +INF 和 -INF 相加，而 -INF 得到 +INF 的方法只能乘以自己，因此，回傳值為 `x * x + x` 才符合要求。
* 那麼，如果想要不用乘法做回傳值，且完成sqrt(NaN) = NaN, sqrt(+INF) = +INF, sqrt(-INF) = sNaN，有什麼辦法?
* 先判斷負號，讓 -INF 在進到這個if之前，就因為負號被驅逐了，這樣就可以直接回傳 x 。
```clike=49
/* take care of zero */
if (ix0 <= 0) {
if (((ix0 & (~sign)) | ix1) == 0)
return x; /* sqrt(+-0) = +-0 */
if (ix0 < 0)
return (x - x) / (x - x); /* sqrt(-ve) = sNaN */
}
```
* 判斷為 0 或負數，+-0 的開根號為 +-0 ，負數的開根號為 sNaN 。可由 ix0 判斷，這也是為什麼 ix0 的 type 要為 int32_t 而非 uint32_t 。
* x 為 0 ，除了 sign bit 無所謂以外，其他 bit 都為 0 。拆成 msw 及 lsw 兩部份來看: `ix0 & (~sign)` 即是判斷 msw 中除了 sign bit 以外，其他 bit 是否全為 0 ， `ix1` 就是lsw全部是否為0。
* x 為負數，須回傳sNaN。
```clike=56
/* normalize x */
m = (ix0 >> 20);
if (m == 0) { /* subnormal x */
while (ix0 == 0) {
m -= 21;
ix0 |= (ix1 >> 11);
ix1 <<= 21;
}
for (i = 0; (ix0 & 0x00100000) == 0; i++)
ix0 <<= 1;
m -= i - 1;
ix0 |= (ix1 >> (32 - i));//x << i
ix1 <<= i;
}
```
* `m = (ix0 >> 20);`，表示 m 為 sign + exp 。不過，透過前面的if負數判斷， sign 必定為 0 。
* 在非正規化(exp = 0)的情況下，先將它變成正規化，透過 shift 來操作。先找出 frac 中，第一個 1 的位置j，再將其移到 exp 的最小位，並修改 m 。也就是說，原本 x = 2^1-bias^ * f ，現在 x = 2^e-bias^ * (f * 2^j^) ，因此，實際上 x = 2^m^ * 2^-bias^ * (f * 2^j^) ， m = -j+1，1代表e
* 假設原本 f 為 0001010...0，視為 0.0001010...0， shift 後， f 為 (1) 010...0 ，視為 1.010...0。不需要再考慮是否加一，這大概是非正規化數 exp 和 frac 設計巧妙的地方。
```clike=70
m -= KK3; /* unbias exponent */
ix0 = (ix0 & 0x000fffff) | 0x00100000;
```
* 之後，又進一步將 2^-bias^ 也包進 2^m^ 之中，同時，修改 ix0 exp 的部分為 0x001。此時， x = x = 2^m^ * f' (f' = f * 2^j^ ， m = e - j - bias)
因此，++KK3 = bias = 2^11-1^-1 = 1023++。
```clik=72
if (m & 1) { /* odd m, double x to make it even */
ix0 += ix0 + ((ix1 & sign) >> 31);
ix1 += ix1;
}
m >>= 1; /* m = [m/2] */
```
* 這段要做 exp 部份的開根號，基本上就是次方除以二。但要考慮到 m 可能是奇數的情況。
* m 如果是奇數的話，要將 exp+frac 的部分乘以二，以表示會被忽略的部份。
```clike=78
/* generate sqrt(x) bit by bit */
ix0 += ix0 + ((ix1 & sign) >> 31);
ix1 += ix1;
q = q1 = s0 = s1 = 0; /* [q,q1] = sqrt(x) */
r = 0x00200000; /* 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 + ((ix1 & sign) >> 31);
ix1 += ix1;
r >>= 1;
}
r = sign;
while (r != 0) {
t1 = s1 + r;
t = s0;
if ((t < ix0) || ((t == ix0) && (t1 <= ix1))) {
s1 = t1 + r;
if (((t1 & sign) == sign) && (s1 & sign) == 0)
s0 += 1;
ix0 -= t;
if (ix1 < t1)
ix0 -= 1;
ix1 -= t1;
q1 += r;
}
ix0 += ix0 + ((ix1 & sign) >> 31);
ix1 += ix1;
r >>= 1;
}
```
* 參考[二進制根式法](https://home.gamer.com.tw/creationDetail.php?sn=2032528)、[wikipedia](https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Binary_numeral_system_(base_2))、[√Square Roots](http://www.azillionmonkeys.com/qed/sqroot.html)，了解直式開根號的方法，基本上就對應了這段程式碼。兩個while迴圈，基本上就是對ix0及ix1做開根號的動作。
```clike=115
/* use floating add to find out rounding direction */
if ((ix0 | ix1) != 0) {
z = one - tiny; /* trigger inexact flag */
if (z >= one) {
z = one + tiny;
if (q1 == (uint32_t) 0xffffffff) {
q1 = 0;
q += 1;
} else if (z > one) {
if (q1 == (uint32_t) KK4)
q += 1;
q1 += 2;
} else
q1 += (q1 & 1);
}
}
```
* 這一段在做捨入， IEEE 754 預設捨入方法是向偶數捨入，但我看不懂 z 是在做什麼。
* 在需要進位的情況下，因為[q,q1]實際為一組數字，如果 q1 全為 1 的話，進位之後，會變成 1 0000...00。
* q1 + 2 之後，q 需要 +1 ，這就表示 ++KK4 = 1 000...00-2 = 111...10 = 0xfffffffe++ 。
* `q1 += (q1 & 1);` 就比較單純，表示向偶數捨入。
* 直接跑過程式，看能不能幫助理解。結果，都如我所預測，`z = one - tiny` 由於 tiny 很小，會被忽略，所以 z == one ，進入第二層 if ，在 q1沒有剛好為 `0xffffffff` 時，會直接進入第二個 else ， 因為 `z = one + tiny;` 的 tiny 也會被忽略，所以 z == one 。
* tiny = 1.0e-300，double可表達的範圍約1.7e-308~1.7e-308
:::info
- `z = one - tiny; /* trigger inexact flag */`
`if (z >= one)`的目的是什麼？在什麼情況下會出現 z > 0 ?
- 在什麼情況下，tiny 不會被忽略？
:::
```clike=131
ix0 = (q >> 1) + 0x3fe00000;
ix1 = q1 >> 1;
if ((q & 1) == 1)
ix1 |= sign;
ix0 += (m << KK5);
INSERT_WORDS(z, ix0, ix1);
return z;
```
* exp 要將 1-bias 加回去，因此加上 0x3fe00000 ，此時為 2^0^。
* m 為開根號後的 exp ，要加回 ix0 ，exp 在20~30的位置，因此，要做 shift ，且 ++KK5 = 20++ 。
* 最後再將 ix0 和 ix1 合成 z 就大功告成了~~~
### 延伸題目: 改為 float (單精度) 版本，注意應該用更短的程式碼來實作
* 想法與思考
* 要注意 ix0 和 ix1 的運算，float不需要拆成兩部分。
* 不需要拆成兩部份，是否還需要union的轉換？
float 無法做 bitwise 的操作，需要轉換成int才行。
* 先做負號判斷，再做 +INF and NaN ，這樣可以直接回傳 x 。 -INF 的部分，在負號判斷時，已回傳 sNaN 。
* bias = 2^8-1^ -1 = 127
* tiny 改為1e-30，因為 float 可表達的範圍約3.4e-38~e.4e-38
* 程式碼
```clike=
#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 & 0x7ff00000) == 0x7ff00000) {
/* 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;
}
```
* 執行結果
```
sqrt 1e+78
1.#INF00
sqrt -1e+78
-1.#IND00
sqrt 9
3.000000
sqrt 1.44
1.200000
sqrt 17453723534
132112.546875
```
## [因為自動飲料機而延畢的那一年](http://opass.logdown.com/posts/1273243-the-story-of-auto-beverage-machine-1)及學習本課程 3 週之後的感想與啟發
- 現實與想像的落差。
現實其實有更多的層面需要考量，我們時常將事情簡單化，無論是面對自以為龐大的專案，還是認為很簡單的幾個流程，實際上都是需要經過很多的研究時間、成本。我沒想過，一隻十塊錢不到的原子筆、坐上去舒服的椅子，這些生活用品是怎麼做出來的。
- 學用落差。
我曾經以為，土木系學用落差嚴重只是因為我們不可能真的找一個地方給房子、蓋橋，只能紙上談兵，那牽涉到太多人力、成本。而資訊系這方面的問題，比較沒那麼嚴重，畢竟寫程式每個人的電腦都可以做，而且馬上就可以看到結果。實際上不是，寫寫小程式不難，就像土木系在做模型時一樣，但面對大一點的就容易出現問題了。
- 全心投入。
看完這片文章，讓我想起一年前，為了做一個小房子的模型，而做過無數的實驗跟研究。那時是要蓋一個混凝土的模型，原本的想像也是很簡單，設計尺寸、做模、灌水泥、拆模。但實際下去做之後，才發現現實世界有太多問題，那時為了做這件事，整天都在想、在試。做完之後，成果也不怎麼樣，最後也是收到角落。而這之中，最寶貴的，我想，不是結果多厲害，而是在那段時間全心投入的自己。在實習過後，體會更深的，工作之餘，想要再將心力投入到其他事情上，是沒有那麼容易的，尤其是需要耗費腦力的事。希望自己能一直保有那份熱情。
- 關於課程
一開始覺得壓力很大，這種上課模式沒有遇過 : 每週都小考，還有很多課後學習的教材要看，作業還會被同學批評、觀摩 ~~(害羞)~~，一週花費超過16個小時在上面，不過，還是很多教材沒看完。但，這好像比較接近我想像中的大學課程，有點半自學的性質，也可以說很幸運的在最後一個學期修到這門課。
- 程式寫得好的人，數學都很厲害。
這三、四週下來，深深的體悟到什麼是真正在用的程式。其實很多事情的原理並不難，只是我沒有好好理解過，隨堂考的題目要解出來其實不難，但是在沒有細細看過程式前，很多都是靠感覺猜的，因為那些背後的理論我尚未融會貫通。但是，要完全看懂程式碼在做什麼，為什麼要那樣做就完全是另外一回事了，這大概就是理論跟實做的差距吧。
## 「[課後學習](http://wiki.csie.ncku.edu.tw/sysprog/schedule)」教材和 [CS:APP 3/e](http://csapp.cs.cmu.edu/)心得和提問
Sign in with Wallet
Connect another wallet