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

[作業要求](https://hackmd.io/s/BJ0v7eBFz)

## 第一週測驗題 測驗`1`

考慮以下程式碼:
其作用為檢查輸入整數是否為 N 的倍數，那麼 N 為多少？

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

一開始看到這個題目可能會不知道從哪裡下手，那就嘗試從最後面一層一層推上去
先看到`if (n == 0)`代表最後一次遞迴的 n 需要為 0 才是題目要求的 **N 的倍數**

再往上一層看到`return(isMultN(abs(odd_c - even_c)));`，發現 odd_c 和 even_c 必須數值相同

再往上一點，看到 while 迴圈裡面，嘗試解釋一下迴圈的作用
* 先判斷 LSB 是否為 1，是的話 odd_c + 1
* n 向右位移一個 bit
* 再判斷 LSB 是否為 1，是的話 even_c + 1
* n 向右位移一個 bit

可以發現 while 迴圈的用途是數奇數位和偶數位 bit 為 1 的數量，如果要讓 odd_c - even_c 為 0，以一個 8 bit 數字來說，可能的數字有
* $11111111_2 = 255_{10}$
* $00001111_2 = 15_{10}$
* $01100110_2 = 102_{10}$
* ...還有很多就不一一列出了

從題目的要求可以反問**什麼的倍數奇數位元 - 偶數位元**會是 0 ？
發現是 3 ，因為三的二進位表示是 $11_2$
做一些簡單的乘法就會發現
* $3\times2$ 是向左位移一個 bit
* $3\times3$ 是向左位移兩個 bit
* $3\times10$ 是向左位移三個 bit，再加上 $3\times2$

都滿足奇數位元 - 偶數位元等於 0，因此答案是 `N=3`

### 延伸問題：將 N 改為其他質數再改寫為對應的程式碼，一樣使用遞迴

參考[這篇文章](http://mypaper.pchome.com.tw/iustlovefish/search/%E4%BA%8C%E9%80%B2%E4%BD%8D%E5%88%A4%E6%96%B73%E7%9A%84%E5%80%8D%E6%95%B8)裡對 3 的倍數的分析，依樣畫葫蘆找出 5 的判斷方式

#### 判斷 5 的倍數
* $2^0\%5 = 1$
* $2^1\%5 = 2$
* $2^2\%5 = 4$
* $2^3\%5 = 3$
* $2^4\%5 = 1$
* $2^5\%5 = 2$
* .....之後一直以餘數為 1 2 4 3 循環下去

觀察一下循環的餘數，如果想要讓數字被 5 整除，必須讓餘數也被 5 整除
會產生這六種組別：(1, 4) (2, 3) (1, 2, 2) (1, 1, 3) (1, 1, 1, 2) (1, 1, 1, 1, 1)
步驟
* 4 個 bit 一組，計算有幾個
* 看第一個 bit 和 第二個 bit 的兩倍有沒有相同
* 有 -> 5 的倍數
* 沒有 -> 不是 5 的倍數

> 發現試著找規律總是會混入其他數字（不只 5 的倍數）

靈機一動！（其實之前靈機一動很多次了，但是都失敗告終）
繼續上面的思路，每 4 個 bit 一組，依序取 5 的餘數分別是 1, 2, 4, 3。

在前面我嘗試要把所有組合的可能列出來，再一個一個判斷，完全忘了可以使用遞迴！
我們將這些 bit 按照餘數賦予權重，然後遞迴檢查是否是 5 的倍數，這樣應該是可行的，接下來我用程式實做一次看看。

int is5(int x) {
    int ret = 0;

    if (x == 0 || x == 5) return 1;
    else if (x < 8) return 0;

    while(x) {
        if (x & 1)  ret += 1;
        x >>= 1;
        if (x & 1)  ret += 2;
        x >>= 1;
        if (x & 1)  ret += 4;
        x >>= 1;
        if (x & 1)  ret += 3;
        x >>= 1;
    }   
    return is5(ret);                                                            
}

實做後發現有思慮不周全的地方，有些數字按照上述的方法生成的結果會等於自己，這會導致無窮迴圈。這些數字是 0, 1, 2, 3, 4, 5, 6, 7，額外篩選掉這些特別的數字就可以成功判斷是否為 5 的倍數了。
:::info
但是我還沒想到為什麼只有 0～7 會產生這種特別的狀況
:::
:::success
找到了！因為循環到 4 之後就會因為被取餘數而打亂權重，詳細解釋後面有寫。
:::

終於想到怎麼解決 5 的倍數讓我信心爆棚，我決定試試看其他數字，但是在這之前我要先確認是不是對於所有的數字 A ，讓 $2^n\ \% \ A$ 產生的數字一定會循環。

#### 簡易的證明：

首先，我們可以確定 $2^0\ \%\ A = 1, \forall A \in N$
接下來的餘數我們可以推導為 $2^n$，但是餘數不能大於除數，所以正確的餘數形式應該會是 $2^n\ \%\ A$

然後可以發現當我們產生了一個之前已經產生過的數字，這個數列必定會循環（但是不保證循環的數字和數量）

這個並不難理解，假設 $A = 11$，產生的前11個數字會是
`1 2 4 8 5 10 9 7 3 6 1`
注意到最後一個數字是 1 ，在前面已經產生過了 1 ，因此後面必定會按照前面的數列循環下去。
因為這樣，我們可以很確定對於所有正整數 A ，$2^n\ \%\ A$ 必定會循環，而且循環數列大小會小於等於 A-1

#### 判斷 13 的倍數
依樣化葫蘆，我們先找出 13 的循環數列
`1 2 4 8 3 6 12 11 9 5 10 7`
然後再找到會造成無窮迴圈的數字 0~15，就可以完成判斷是否為 13 的倍數了。

### 思考：這麼做的好處
除法運算在四則運算中最慢，如果可以只用加減或是位元運算子來代替`%`，在速度上會快很多。<s>缺點是這種方式沒有規律</s>，必須要根據不同的數字找到對應的規則。
:::warning
做到最後面回來打自己嘴巴，也許 % 運算子真的比較慢，但是我們繞了那麼一大圈只為了求得某個數字是否是某個質數的倍數，我覺得不太值得。
:::

後來想了一下這其中應該是有規律的，只是有沒有找到而已。**這個方法應該能夠被推廣到其他的質數**，只是需要更嚴謹的證明。

:::info
還沒找到造成無窮迴圈的數字的規律
:::
:::success
找到了！會造成循環的數字都有一個共同的特性：乘上權重後仍然是自己。
從這個角度觀察，就可以知道唯有權重還沒有被取過餘數的數字，乘上權重會是自己。

以數字 17 為例子：
循環數列是 `1 2 4 8 16 15 13 9`，發現在 16（含）之前都沒有被取過餘數（取過之後數字沒變），因此數字 0~31(32-1) 會是造成無窮迴圈的數字。
:::
> 終於想到了好開心 :notes:

都做到這裡了，來模仿[`blake11235`](https://hackmd.io/s/Syaf4cqYM#) 同學和 mod 比較時間差異好了（不然寫出來沒有比較快不知道可以幹麻 QQ)。
:::info
現在才發現這樣有點本末倒置的感覺...我用 % 運算子找出一些特定的數字來協助我不要用 % 運算子來計算 Orz，不過都做了就做下去吧！(其實可以用減法和加法來代替)
:::

:::danger
還是有機會變快的，考量到數值範圍超過 64-bit (如密碼學常見的 256-bit) 時，上述的手段就會跟一般 mod 運算有差距。來實驗吧！
:notes: jserv
:::

放一下程式碼當參考(部份程式碼直接使用 `blake11235` 同學的)

#include <stdio.h>
#include <time.h>
#define MAX_LIST_NUM 100000

static double diff_in_second(struct timespec t1, struct timespec t2) {
    struct timespec diff;
    if (t2.tv_nsec-t1.tv_nsec < 0) {
        diff.tv_sec = t2.tv_sec - t1.tv_sec - 1;
        diff.tv_nsec = t2.tv_nsec - t1.tv_nsec + 1000000000;
    } else {
        diff.tv_sec = t2.tv_sec - t1.tv_sec;
        diff.tv_nsec = t2.tv_nsec - t1.tv_nsec;
    }
    return (diff.tv_sec + diff.tv_nsec / 1000000000.0);
}

/* build cycle list */
void build_list(int p, int *list, int *cnt) {
    int db = 1;
    do {
        list[(*cnt)++] = db;
        db += db;
        if (db > p) db -= p;
    } while(db != 1);
}

int isMultN(int num, int p, int *list, int cnt) {

    while (1) {
        if (num == 0 || num == p) return 1;
        else if (num < p+p) return 0;

        int ret = 0;

        while (num) {
            for (int i=0; i<cnt; ++i) {
                if (num & 1) ret += list[i];                                                  
                num >>= 1;
            }
        }
        num = ret;
    }
}

int main() {
    int p;

    puts("please enter a prime number (less than 100000):");
    while(scanf("%d", &p)){
        int list[MAX_LIST_NUM];
        int list_cnt = 0;
        int x = 10000; // test number: x % p

        build_list(p, list, &list_cnt);
        puts("---list---");
        for (int i=0; i<list_cnt; ++i) printf("%d ", list[i]);
        puts("\n----------\n");

        struct timespec start1, end1, start2, end2;
        clock_gettime(CLOCK_REALTIME, &start1);
        printf("%s\n", isMultN(x, p, list, list_cnt)?"yes":"no");
        clock_gettime(CLOCK_REALTIME, &end1);

        clock_gettime(CLOCK_REALTIME, &start2);
        printf("%s\n", (x%p)?"no":"yes");
        clock_gettime(CLOCK_REALTIME, &end2);

        printf("time using mod: %.10f sec\n", diff_in_second(start2, end2));
        printf("time    no mod: %.10f sec\n", diff_in_second(start1, end1));
    }
    return 0;
}

結果：

time using mod: 0.0000032460 sec
time    no mod: 0.0000049340 sec

繞了一大圈最後還是比用 % 運算子還要慢，自打嘴巴 QQ
> 不要急著下結論，思考大數運算 [name=jserv]

## 第二週測驗題 測驗`3`

考慮到某些實數的二進位表示形式如同 $0.yyyyy...$ 這樣的無限循環小數，其中 $y$ 是個 $k$ 位的二進位序列，例如 $\frac{1}{3}$ 的二進位表示為 $0.01010101...$ (y = `01`)，而 $\frac{1}{5}$ 的二進位表示為 $0.001100110011...$ (y = `0011`)，考慮到以下 y 值，求出對應的十進位分數值。

* y = `010011` => $\frac{19}{X1}$
* y = `101` => $\frac{5}{X2}$
* y = `0110` => $\frac{2}{X3}$

一開始看到題目有點摸不著頭緒，國高中有看過類似的題目，但是也只有按照書上寫的規則依樣畫葫蘆回答問題。想要解答這個問題的話，先嘗試從二進位轉成十進位下手。

$\frac{1}{3}$ 的二進位表示為 $0.01010101...$ (y = `01`)，換算成十進位表示的方式會是

$0\times2^{-1}+1\times2^{-2}+0\times2^{-3}+1\times2^{-4}+0\times2^{-5}+1\times2^{-6}...$之後按照規則循環下去

接下來利用分配律化成比較好看的形式：

$\ \ \ \ (0\times2^{-1}+1\times2^{-2})(1+2^{-2}+2^{-4}+2^{-6}.....)$

$=(0\times2^{-1}+1\times2^{-2})\sum_{i=0}^\infty(2^{-2i})$

$=(1\times2^{-2})(\frac{1}{1-2^{-2}})$

$=\frac{1}{4}\times\frac{3}{4}$

$=\frac{1}{3}$

可以手寫轉換了，但是還要找出一個公式比較好判斷答案。
從上面的轉換過程可以發現答案會是 $x \times \frac{1}{1-2^{-n}}$ ，其中 n 是循環的位數。
在 x 的地方，我們只要向右位移 n 位就可以得到整數的值，記得最後的答案要除回去。

以 $0.01010101...$ (y = `01`) 為例子，$01_2 = 1_{10}$ ，$\frac{1}{1-2^{2}} = \frac{4}{3}$ ，$\frac{1 \times \frac{4}{3}}{4} = \frac{1}{3}$ ，得到答案是 $\frac{1}{3}$ 。

其他數字按照上面寫的方式就可以算出答案了！

## 第二週測驗題 測驗`1`

請完成下方程式碼，依循 IEEE 754 單精度規範，輸出 $2^x$ 的浮點數表達法，而且考慮兩種狀況：
* 當 $x$ 過小的時候，回傳 $0.0$
* 當 $x$ 過大的時候，回傳 $+\infty$

注意：這裡假設 `u2f` 函式返回的浮點數值與其無號數輸入有著相同的位數，也就是至少 32-bit

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

`exp2_fp` 會依序判斷太小，非正規數，正規數，太大的數字，並適當的回傳 $2^x$ 的值（如果表示得出來）

非正規數指的是那些介於 1 和 -1 的數字，他們的 exponent 會是 127
正規數指的則是所有不包含 非正規數，NaN，和 INF 的所有表現得出來的數字

我們先列出在這些情況下的數值，方便之後判斷
* sign bit 為 1 的正規數下限：`1 11111110 00000....`
(`exponent` 為 `00000000` 是非正規數)
* sign bit 為 1 的正規數上限：`1 00000001 11111....`
* 非正規數下限：`1 00000000 111111111....`
* 非正規數上限：`0 00000000 111111111....`
* sign bit 為 0 的正規數下限：`0 00000001 00000....`
* sign bit 為 0 的正規數上限：`0 11111110 11111....`
(`exponent` 為 `11111111` 是 inf 或 nan )

### 解題：
* 第一個判斷需要篩選出太小的數字。
要篩選出太小的數字可以從 sign bit 為 0 的非正規數下限下手。
> 因為要找到最小能夠表示的正數

轉成二進位表示後，發現下限為 $2^{-23} \times 2^{-126} = 2^{-149}$
比下限還要更小的數就無法用單精度浮點數表示了，所以 `2 - pow(2, Y0) - 23 = -149`
:::info
https://en.wikipedia.org/wiki/IEEE_754-1985#Examples
根據 wikipedia 在 denomalized 時 Actual exponent 為 -126
也因此 $2^{-23} \times 2^{-126} = 2^{-149}$
:::
> 謝謝 vulxj0j8j8 同學，這回答到我下面的問題了！我發現因為不習慣 `2 - pow(2, Y0) - 23` 的表示方式讓我被困在某種思維中離不開。按照平常的想法我會把這樣的數字看成 `-23 + -(pow(2, Y0) - 2)`，這樣的意思就是兩個負數相加！
>
>另外，在[成大資工 Wiki](http://wiki.csie.ncku.edu.tw/) 的首頁也可以看到，[wiki](https://en.wikipedia.org/wiki/Wiki) 不等於 wikipedia。我覺得應該要說明白是根據 wikipedia 而不是 wiki ，因為 wiki 涵蓋的範圍太廣了。

* 第二個判斷需要篩選非正規數
因為題目不會出現負數的情況，所以可以直接用 $<$ 把非正規數篩選出來而不用限定一個區間。
再參考上面列出的上下限，非正規數上限的二進位表示會是 $(2-2^{-23}) \times 2^{0-126} = 2^{-125}-2^{-149} > 2^{-126}$，得到`Y1-pow(2,Y2) = -126`

* 把 $2^x$ 轉成非正規數形式表示
因為非正規數的 exponent 可以確定是 `00000000`，所以我們只需要討論 mantissa 的部份就好。

首先，mantissa 有 23 bit 可以儲存資料。注意到 $1101_2$ 轉換成 10 進位數字會是 $2^{-1}+2^{-2}+2^{-4}$，和在 $x>0$ 的情況剛好相反。意思是如果想要用 bitwise shift ，就需要移動 23+n 位而不是 n 位。

現在為到非正規數的表示： $2^x = 2^{-126} \times 2^n$ ，整理後得到 $n = x+126$
再把 n 代入剛剛得到的結果，如果我們想要用 bitwise shift 來達到目的，需要向左位移 $23+x+126 = x+149$ 個 bit。
在和原本的式子互相比對後就可以得到答案了。
:::info
感覺這個方法不是最好的方法，雖然有算出需要向左位移 x+149 個 bit，但真正的答案還是有點湊出來的感覺
從答案看到數字 23 就可以聯想到 mantissa 有 23 位，不知道是否有什麼關聯？
:::
> 已回答在上面

* 第三個判斷需要篩選正規數
因為剩下來的數字只有正規數和太大兩種而已，所以可以很容易的知道條件一定要是 x 小於 exponent 的極限。
搭配上面列出的 sign bit 為 0 的正規數上限，發現最大值會是 $(2-2^{23}) \times (2^{254-127}) < 2^{128}$

* 把 $2^x$ 轉成正規數形式表示
因為都是在二進位之間的轉換，所以 mantissa 的部份一定會是 0 。
在 exponent 的部份只要注意先把 bias 加上去就可以了
（bias 在單精度是 127，在雙精度是 1023）

因為在前面我們都把 exponent 和 mantissa 分開求值，最後要合併時只要把 exponent 推到正確的位置 (<<) 再把兩個合併 (|) 就可以得到結果了。
另外，因為不會有負數的情況發生，所以不用考慮 sign bit 的問題。

> 這一題其實並不需要什麼技巧或是什麼想法才寫的出來，只要照著 IEEE 754 的標準走就可以了。
> 但是就是這樣讓我在這裡花了好多時間，寫了才發現自己似乎根本沒有熟悉這個標準過，但是之前遇到類似問題（浮點數表示法）時總是裝作沒看到。
> 遇到沒搞懂的東西真的要馬上去試著理解他而不是放過就算了 QQ
:::success
學海無涯，回頭(看書)是岸
:notes: jserv
:::