# 2018q1 Homework3 (assessment)
contributed by <`PunchShadow`>

## 測驗題重作

### 第三週第一題

[題目](https://hackmd.io/s/SknkEfVFf)：主要程式目的為進行$\sqrt{2}$ 的 double(兩個32 bits) float point 運算

#### 程式理解及解題想法

typedef union {
    double value;
    struct {
        uint32_t lsw;
        uint32_t msw;
    } parts;
} ieee_double_shape_type;

首先將一個 double floating point 之值存於 value 中，並將其 IEEE754 表現法拆成 32 bits 的兩部份（lsw 和 msw)，並採用 Little Endian 的方式

> 
> msw 和 lsw 分別代表 Most Significiant Word 和 Least Significiant Word
> 以此圖為例，若資料為 0x0A0B0C0DFFFFFFFF 則若採 Little Endian 方式儲存，則 0A 會存於 memory address 最小的地方，而其 msw 和 lsw 則分別是 0x0A0B0C0D 和 0xFFFFFFFF

由上圖可以看出 lsw 為 double-precision 的前半段（相當於single-precision 的長度），其中包含了 sign bit、exponential 和一部分的 fraction，而 msw 則全部都是 fraction 的部份。

/* 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)

此段程式碼是定義從32 bit ints 轉成 double-precision floating point 的方式，延續上面的 type define，將 .parts.msw 和 .parts.lsw 為 [designated initializor](https://gcc.gnu.org/onlinedocs/gcc/Designated-Inits.html) 宣告在 structure type iw_u 中 msw 和 lsw 的值為 ix0 和 ix1(==這兩個搞了超久才知道原來是一個 int32_t(ix0)和 uint32_t(ix1) 變數而已:sob: 原本還以為是16進位的某個表示法勒==），並將 iw_u.value assign 到 d 中，

:::info
不知道是否可以將這個 function 看成 **輸入 ix0 和 ix1 兩個 int32_t 得到 一個 double floating point d ?** 用 **define 和 function 有甚麼差別嗎** ?
:::

這邊為了瞭解 Macro 的使用及為何這個 function 可以將兩個 int32 轉換成 double ，因此找到了這個程式再 Android 內的[原始碼](https://chromium.googlesource.com/android_tools/+/af1c5a4cd6329ccdcf8c2bc93d9eea02f9d74869/ndk/sources/android/support/src/msun/math_private.h)，看到很多 Macro 的形式都長得像這樣，讓我有了一個想法
:::info
有沒有辦法利用兩個連續記憶體位置的 int32 去表示一個 double-precision floating point ?
於是就有了以下的實驗
:::

#### 2個連續 int32 轉成 double-precision floating point

#include <stdio.h>
#include <stdint.h> 

/* 使用 sqrt function 內 union 的 data structure */
typedef union
{
  double value;
  struct
  {
    /* CPU 為 i5-4200h 為 Little Endian 的儲存方式 */
    /* 宣告的 lsw 記憶體位置就會小於 msw   */ 
    uint32_t lsw;
    uint32_t msw;
  } parts;
} ieee_double_shape_type;


int main(){
    ieee_double_shape_type test;
    
    /* 0x4034000000000000 = 20.0 */
    test.parts.lsw = 0x00000000;
    test.parts.msw = 0x40340000;
	
    /* Want to read the continuous address 
     * of lsw and msw in double floating type. */
    double *d = &test.parts.lsw;
	
    printf("lsw's address : %d\n", &test.parts.lsw);
    printf("msw's address : %d\n", &test.parts.msw);
    printf("d is %lf\n", d);
	
    return 0;
} 

跟預想的一樣，lsw 的記憶體位置在 msw 前面 4 個 bit
但 d 卻沒有像預想中一樣印出，看來想用 pointer 存取兩個連續的 int32 是行不通的

lsw's address : 343809464
msw's address : 343809468
d is 0.000000

後來在 Stackoverflow 上面看到有關 [bit sequential combination](https://stackoverflow.com/questions/11824762/convert-two-32-bit-floats-to-one-64-bit-number-and-vice-versa) 的問題，於是就將 d 直接改用 union 內宣告的 double value

double d = test.value

想不到真的可以成功，而且 `test.value` 的初始 memory address 竟然是跟 test.parts.lsw 一樣，這讓我有點驚訝，於是就去求助 C99 規格書啦！

>(6.7.2.1-15) Within a structure object, the non-bit-field members and the units in which bit-fields reside have addresses that increase in the order in which they are declared. A pointer to a structure object, suitably converted, points to its initial member (or if that member is a bit-field, then to the unit in which it resides), and vice versa. There may be unnamed padding within a structure object, but not at its beginning.

這裡很明確的指出會有一個 pointer 指向 object 的起始，且後來宣告的會 padding 在 object 的後方，不過還是覺得奇怪，為什麼 double 不會在一開始就 allocate 了呢？
後來查了一下，發現原來問題出在 **Union** 上面阿，原來 Union 和 Struct 雖然功能很像，不過差別就差在 **Memony Allocation** 上面阿， Union 一次只准一個 member access ，若不同 member 同時 access，則會造成 **data garbage** 的現象，請參考 [C Programming Unions](https://www.programiz.com/c-programming/c-unions)

這樣就很蠻合理的，由於一次只能 access 一個 member，所以 address 會一樣，原本應該要被視為 garbage 的兩個 int32 ，卻可以當成轉換成 double-precision floating point 的兩部份，真是高招！

---
拉回來～拉回來～

#define EXTRACT_WORDS(ix0, ix1, d)

這個 function 則是與 INSERT_WORDS 相反，將一個 double floating point `d` 轉成兩個 int32_t `ix0` 和 `ix1`

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;
    }
.......

* x : 想進行平方根的數
* z : 平方根後的結果
* sign : 0x80000000 寫成二進位的話即是取出 32 bits 中的第一位
* **ix0 : 代表 x 的 Most Significant Word 部份**
* **ix1 : 代表 x 的 Least Significant Word 部份**
* if ((ix0 & ==KK1==) == ==KK2==) : 這邊是要判斷 INF 和 NaN 的部份，首先我們先從 INF 和 NaN 的定義下手
上圖為 single-precision 在 IEEE 754 中定義的 INF 和 NaN 皆是 exponent 部分全為一，而 double-precision 亦是如此，因此要 找出部分(1~11個 bits)需要 ix0 和 ==0x7ff00000 (KK1)== 做 AND，且得出的結果應該也是 ==0x7ff00000 (KK2)==
* return x * x + x : sqrt(NaN) = NaN, sqrt(+INF) = +INF, sqrt(-INF) = sNaN
這邊還蠻好玩的，先稍微簡介一下在 IEEE 754 中所定義的 +-INF 和 NaN
* INF :
* 數學運算上出現無限時 e.g : 1/0 = INF
* 當運算大於所定義的最大值時(single-precision : $\ {3.4*10^{34}}$ , double-precision : $\ 1.797 * 10^{304}$) e.g : $\ 10^{1000}$ = INF
* 有正負 INF
* NaN :
* 未定義的數學運算 e.g : 0/0 = NaN
* 分為 `quiet NaN` 和 `signal NaN` : 前者會在程式執行完成後再進行檢查，後者則會產生 exception 需要進入 exception handle 進行例外排除
* INF 和 NaN 運算 :
* INF + 10 → INF
* INF/0 → INF
* 10 - INF → -INF
* 100/0 = INF → True
* ==INF - INF → NaN==
* INF/INF → NaN
* ==INF * -INF = -INF==
* ==INF - INF → NaN==
* 各種例外情況 :
* x 為 NaN : NaN * NaN = NaN, NaN + NaN = NaN
* x 為 +INF : (+INF) * (+INF) = (+INF), (+INF) + (+INF) = (+INF)
* x 為 -INF : (-INF) * (-INF) = (+INF), (+INF) + (-INF) = NaN (引發 exception)
* **所以 x * x + x 就可以用上面黃底的運算得到所定義的內容，真是太強大惹** :+1: :+1:
>參考 : [INF, NAN, and NULL - Exception values](http://wiki.analytica.com/index.php?title=INF,_NAN,_and_NULL_-_Exception_values)

/* 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 */
}

這邊是要處理零的情況 (不是靈的轉移 :laughing: )
想要是 0，exponent 的部分就必須要是 0 才可以，一開始有點疑惑， ix0 是 int32_t 的形式，所以若 sign bit 是 1 的話不論 exponent 部分為多少，條件皆會成立，進入迴圈判斷後才決定是 0 或是負數，雖然邏輯上面是沒有錯，不過這樣是不是有點多此一舉？
:::info
加上 **if (ix0 <= 0)** 的目的是希望能減少 branch 嗎? 不然在程式邏輯上，這條件和下面兩個條件好像有點重疊
:::

:::danger
程式碼是給你修改實驗用，不是拿來背誦，有想法很好，就動手改！
順便研究該如何有效測試程式
:notes: jserv
:::

/* 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] */

一開始看到這邊的時候有被 >> 困惑住，讓我開始想這個 shift right operator 是 arithmetic 或是 logical 的，後來查了 C11 規格書後，在(6.5.7)中很明確的告訴我們 shift right 是 arithemtic
>The result of E1 >> E2 is E1 right-shifted E2 bit positions. If E1 has an unsigned type or if E1 has a signed type and a nonnegative value, the value of the result is the integral part of the quotient of E1 / 2E2. If E1 has a signed type and a negative value, the resulting value is implementation-defined.
><C11 Standard> (6.5.7)

這下我們就知道 **m 是 ix0 扣除掉 fraction 的值，也就是剩下 sign bit 和 exponent 的值**，而此時 m 若等於 0 代表其 exponent 和 sign bit 皆為零，因此其值應該要以 Subnormal(Demornal) 的形式呈現，這邊原本有個小小疑問，若為**負的 subnormal**，則 m 應為 0xfffff800 != 0 則不會進入判斷 subnormal 中，但後來發現，**負數的情形竟然在上面已經討論過了**，以下已經不會再有負數的情況出現，所以==條件的編排順序也是很重要的阿!!!==

一開始想到讀到 Submornal 正確值，應該要主要的部份是 fraction 內的值，所以應該是想辦法找出 ix0[12:31] 和 ix1[0:31] 然後其值應為 0.xxxx * $\ 2^{-1024}$(xxxx為 fraction 部分)

while (ix0 == 0) {
    m -= 21;
    ix0 |= (ix1 >> 11);
    ix1 <<= 21;
}

上面這段還蠻玄的，想了很久終於想通倒底在做甚麼
* while (ix0 == 0) : 已經知道 0~11 bit 皆未 0 ，那用此條件的用意在哪？其實很簡單，就是要找出 fraction 前半段(12~31 bit) 是否為 0
* m -= 21 : 若 fraction 前 21 個 bits 皆為 0 ，那麼若要正規化前 21 個 bits 皆為零的二進位，就要乘以 $\ 2^{-21}$


* ix0 |= (ix1 >> 11) : 進行這段動作的原因是，將 ix0 的 fraction 部分(12~31 bit)套上 ix1 的 0~19 個 bits (如上圖)
* ix1 <<= 21 : 因為已經 * $\ 2^{-21}$ ，所以位數上升 21 位
* while 用意 : 可能 ix1 之 0~19 bits 皆為零，此時就要再進行一次 bit operation，而離開迴圈時，可以將 ix0[12:31] 填入非 0 值，並得到 $\ *2^{m}$

for (i = 0; (ix0 & 0x00100000) == 0; i++)
            ix0 <<= 1;
        m -= i - 1;
        ix0 |= (ix1 >> (32 - i));
        ix1 <<= i;

* for 迴圈的重點在於 0x00100000，是將 ix0 向左 shift 直到**第12個 bit 為 1**為止，目的在將 fraction 的第一個 bit 對齊為 1
* m -= i-1 : shift 了多少個 bit 就要將 m 減去多少，以維持正確的數值
* ix0 |= (ix1 >> (32 - i)) : 由於 ix0 shift 完後最後面幾個 bits 會都是 0 ，所以比需將 fraction 部分從 ix1 中遞補上來
* ix1 <<= i : ix1 遞補了幾個 bits 就要 shift right 幾個 bits
:::success
處理到這裡完，已經將 Subnormal number 化成以下形式 :

真是可喜可賀阿 :100: :100: :100:
不過沒想到ㄚㄚㄚㄚ~~~
這樣處理完之後竟然可以**和 normailized number 一起處理**了呢，真是太神奇了 :+1: :+1:
:::
> 第12個應為編號11才對喔
> [name=HexRabbit][color=orange]
>> 你是說圖的部份嗎？double 的 exponent 部份不是從編號 1~11 嗎？還是 fraction 部份錯了？
>> [name=PunchShadow][color=green]
>>> double 的 exponent 的確是 1~11，
>>> 但是這個操作正如你所說的，
>>> 會將 ix0 向左 shift 直到**第12個 bit 為 1**，
>>> 所以 1 的位置編號應是 11
>>> [name=HexRabbit][color=purple]
>>>> 沒錯，編號 11 的位置是 1，感謝 H大，已更新圖
>>>> [name=PunchShadow][color=green]

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] */

當 subnormalized number 和 normalized number 都化成同一個形式之後，我們就可以一起來處理啦
* m -= ==KK3== : 這邊可以從上面處理 subnormalized number 可以得知 m 就是所謂的 exponential bias， 公式為
$\ {2^{ exponent 位數 - 1} - 1}$
所以single-precision 是 $\ {2^{8- 1} - 1=128}$, double-precision 則是 $\ {2^{ 11 - 1} - 1 = 1023}$
得到 ==KK3 = 1023==
>參考 : [floating point representation](http://fourier.eng.hmc.edu/e85_old/lectures/arithmetic/node12.html)

* ix0 = (ix0 & 0x000fffff) | 0x00100000 : 將 ix0 中除了第 11 個 bit 和 fraction 部分外，都變成 0 ，因為已經取得 bias ，所以將 ix0 都變成 **0x001xxxxx** 的形式
* if (m & 1) : 考慮 m 為奇數的情況，假設 m = -5 則在進行++第7行++程式 m >>= 1 時， [m/2] = -3，所以要將 fraction 內值 ***2**，其數值才會正確
* ix0 += ix0 + ==((ix1 & sign) >> 31)== : 這邊會這樣寫是要考慮當 ix1 兩倍時，是否會進位的問題，所以才要 (ix1 & sign) 去取得 ix1 的第 1 個 bit，如果是 1 則會進位，因此要在 ix0 加上 1，而這種寫法可以**省去一個判斷式，應該就是有品味的程式了吧！！**

:::success
至此終於全部解釋完要開根號的前置作業 (主菜都還沒上就快不行了 :sob:) ，然後一噴就是一個下午，這就是修行阿 :cry:
:::

/* 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 */

* r = 0x0020000 : 設定第10個 bit 為 1，以方便後面 check bit by bit
* q = q1 = s0 = s1 = 0 : 為牛頓法的範圍

:::info
這邊有個疑問，為何 ++2~3++ 還要將 ix0、ix1 還要乘以 2 ? 不是已經處理完奇數的部份了，應該是不用這兩行才對
:::

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;
    }

這邊就是核心部份了，然後.......尷尬想不出來.....只能先往後面發展....

20180326 : 看了 [HexRabbit ](https://hackmd.io/s/BkMaaAmcM#)共筆後，有種豁然開朗的感覺，現在終於可以理解為何老師一再強調**合作的重要性！！**

* 首先要先去理解 2 進位的找平方根方法
這邊用 [Binary numeral system (base 2)](https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Binary_numeral_system_(base_2)) 上面的範例解釋這段程式碼
1. **r** : 可以看成一個指標，從0x00200000，每次 shift right 一個 bit 直到 0，每經過一個 bit 就檢查對應的 ix0 夠不夠減去 t
2. **t** : 如圖所示，是代表左邊進行短除法的數字，若商數(q)出現 1 ，則 t 之最左邊的 bit 需要 *2 加回去 t 當成下一次的 t
3. **s0** : 為進行完一次除法後，要準備下次的暫存值

r = sign; /* 看成從第0個 bit 往左邊掃*/
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;
}

* 這邊可以看成結合 ix0 和 ix1 的平方根，有幾種情況:
* **ix0 為完全平方** : 在上段程式結束後就會歸零，所以 ==if ((t < ix0) || ((t == ix0) && (t1 <= ix1)))== 就會不會成立，直到 ix1 前面進位上來才會執行 if 內的操作
* **ix0 不為完全平方** : ix0 內的剩餘值會和 ix1 做結合，當 ix0 一直 shift left 與 t1 對齊時就會進入++第5行條件式++
* **t1** : 對 ix1 相減的除數，相當於 **t** 對 ix0 的功能，提供更精細的除法用
* **s1** : 相當於上一段程式中 **s0** 的角色，負責當作邊界，然後加回去 **t1**
* **s0** : 承接上一段程式的值，繼續對 ix0 進行平方根，因為 ix0 會持續 shift left ，所以還是要繼續減去此值以維持正確
* **段落邏輯分析** :

		if (((t1 & sign) == sign) && (s1 & sign) == 0) 
			s0 += 1;

* (t1 & sign) == sign : t1 的第 0 個 bit 為 1，代表 t1 以無法再更大，若要進位只能去更改 t
* (s1 & sign) == 0 : s1 的第 0 個 bit 不為 1，因為 s1 = t1 + r ，若 t1 = $\ {1xxxxxx.....}_2$ ，加上 r 卻讓第一個 bit 為 0，代表說要進位，因此要去將 s0 += 1，其值才會正確

		ix0 -= t;

* 如同上方說的， ix0 會一直左移，因此每次進行減法時，都要考慮 t 進去

		if (ix1 < t1)
			ix0 -= 1; 

* if (ix1 < t1) 情況 : 當 t < ix0 正要進行減法時，但 ix1 卻又小於 t1 無法直接進行相減，所以就要從 ix0 借一位後方可減去

		q1 += r;

* **q1** : 相當於 ix1 的**商數**，想法和 q 有異曲同工之妙

		ix0 += ix0 \+ ((ix1 & sign) >\> 31); 
		ix1 += ix1; r >>= 1;

* 將 ix0 和 ix1 shift left 1 bit (和上面做法一樣)
:::info
以上這段完全沒有概念，只覺得有點像快速除法的感覺，++8-9++應該 ix0 和 ix1 *2 的樣子，覺得 *2 很詭異，完全不知道這段的用意，帶數字 trace 一遍還是不曉得，總覺得是在實現 [Binary numeral system (base 2)
](https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Binary_numeral_system_(base_2))不過完全沒有想法
:::
:::success
太好了，再花了一個下午和看了 `HexRabbit` 的分析後，很多東西都慢慢可以稍微推導出來了，開心 :100:
:::

於是崩潰的跳過了無數行程式，終於來到我比較看得懂的地方，這裡是要結束的部分，想說從後面反推搞不好會有意想不到的進展

ix0 = (q >> 1) + 0x3fe00000;
ix1 = q1 >> 1;
if ((q & 1) == 1)
    ix1 |= sign;
ix0 += (m << KK5);
INSERT_WORDS(z, ix0, ix1);
return z;

從這段可以看出來，最後 output 的是 z ，而 z 是由 ix0 和 ix1 結合而來的，而 ix0 又是從 (q>>1) + 0x3fe00000 而來，所以我推測，==前面的 q 一定是與平方根數有絕對的關係
* 這邊可以先猜測一下 **KK5**的值，前面得到的 m 是 bias 開根號後的值，也就是說在最後面要合成 z 的時候要將其加回去，因為要加到 exponent 的部分，所以推測是要 **shift left ==20==** 才有辦法對齊 ix0 之 exponent 部分

---

## 因自動飲料機而延畢的那年心得 (未完成)

在看完了全部23篇的網誌後，心中其實是百感交集的，有很多熱血感動的情緒，但有更多的是對未來的惶恐和不安。

同樣身為大四生，同樣都是飲料的愛好者，同樣希望自己有朝一日能夠改變世界，哪怕只是一小部份也好，也希望為這個世界貢獻一己之力，但令我值得欽佩的是，作者願意勇敢踏出創業的那一小步。大家都在嚷攘著大學教育無法與實做接軌，但又有多少人，是真正去想辦法把想做的東西做出來？又或者說，**有多少人真正有想做的事或者想實現的目標**，這讓我回想起在看 Jserv 的 [你所不知道的C語言：指標篇(上)](https://www.youtube.com/watch?v=G7vERppua9o&t=6742s)中老師所提到的，為什麼有人會對C語言的本質有相當多的誤解和不熟悉，那是因為你**只是把它當作一個語言在學習，而沒有實際想利用c做的事情**。

一聽到這句話，當下其實很心有戚戚焉，不只是c語言，想要精通某事物或達到某個目的，所需要的絕對不僅只是一個**想學好的念頭**，而是一個很強烈的==動機==，而這個動機，必須是與它的本質相符合的。以本身為例，經歷過了研究所考試，雖然旁人都覺得我看起來很辛苦很努力，不過其實我念得很開心，當我很明確的知道我為何要讀這些，讀這些東西對我以後做的事會有什麼幫助後，一切都變得相當充實、相當理所當然，甚至會自發性的去深掘一些背後的理論和延伸，這些東西在我看來，其實一點都不辛苦，因為我知道，這都是為了我以後想做的事情在打的基礎，但也有很多人只是為了考上研究所而讀書，常常就會顯得力不從心（好像有點講多了，套句老師的口頭禪：`拉回來～拉回來` :laughing: :laughing: ）

回到飲料機，在現實世界中

To be continued .....

---

## 課後學習及 CS:APP 3/e 心得、紀錄