# 2018q1 Homework2
contributed by <`workat60474`>

## 第 1 週測驗題 第一題
題目如下：

- 考慮以下程式碼:

#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 為多少？
`(a)` 13
`(b)` 11
`(c)` 7
`(b)` 5
`(e)` 3
`(f)` 2
* 首先，以條列式的方式，列出幾個三的倍數的二進位表示法，嘗試從其中看出規律：

| Decimal | Binary representation |
| -------- | -------- |
| 3 | 11 |
| 6 | 110 |
| 9 | 1001 |
| 12 | 1100 |
| 15 | 1111 |
| 18 | 10010 |
| 21 | 10101 |

我原先試圖以 $3$ 的二進位表示法 $(11)_2$ 來嘗試找出規律，起初我猜測其規律在於，對於所有整除 $3$ 的倍數， 其二進位表示法中的 $1$的個數和$0$的個數相等，不過很快地我發現這在21就不適用，21之二進位表示法 $(10101)_2$ 中$1$的個數為3個，而$0$的個數為2個（在這裡指的$0$是指在隨著上述程式中，右移的值還不為0之前，也就是第一個bit為1之後的0的個數），違法了我的猜測。

我重新trace了一下程式碼，發現我沒仔細研究過這行

 /* Recursive call till you get 0/1 */
    return(isMultN(abs(odd_c - even_c)));

由於整個funtion結束的條件是當 $|odd_c-even_c|$ 為 $0/1$ 的時候，當 $|odd_c-even_c|>1$的時候，funtion還未達到終止條件，而是會再次檢視上一個recursive call 傳進來的 $|odd_c-even_c|$ 是否為 $3$ 的倍數，由此不斷遞迴呼叫下去，直到達到終止條件為止。

假設猜想：
$\forall a\in N , a|3\Longleftrightarrow$ $設a=((a_k ... a_2 a_1 )_2 且 若k\in odd :(|(a_1 + a_3 +...+a_k)-(a_2 + a_4 + ... +a_{k-1})|)\mid3$
$若k\in even :(|(a_1 + a_3 +...+a_{k-1})-(a_2 + a_4 + ... +a_{k})|)\mid3$

首先觀察：
* $21_{10}$$=(10101)_2$
* even_c:0
* odd_c=3
則$|odd_c-even_c|$=$3\mid 3 \Rightarrow$ $21\mid 3$
* $15_{10}$$=(1111)_2$
* even_c:2
* odd_c=2
則$|odd_c-even_c|$=$0\mid 3 \Rightarrow$ $15\mid 3$

猜想：
考慮二進位數 $(a_2 a_1)_{2}$ ，不失一般性假設 $a_2=1$，則
$(a_2 a_1)_{2}=(1 a_1)_{2}=(10)_2*1+(01)_2*a_1$
$\Rightarrow$ $2*1+a_1$
$\Rightarrow$ $(3-1)*1+a_1$
$\Rightarrow$ $(3-1)*a_2+a_1$
$\Rightarrow$ $3*a_2+(a_1-a_2)$
$\Rightarrow$ $若 (a_1-a_2)\mid3$
$\Rightarrow$ $則 (a_2 a_1)_2 \mid 3$

考慮二進位數 $(a_3 a_2 a_1)_{2}$ ，不失一般性假設 $a_3=1$，則
$(a_3 a_2 a_1)_{2}=(1a_2 a_1)_{2}=(100)_2*1+(010)_2*a_2+(001)_2*a_1$
$\Rightarrow$ $4*1+2*a_2+a_1$
$\Rightarrow$ $(3+1)*1+(3-1)*a_2+a_1$
$\Rightarrow$ $(3+1)*a_3+(3-1)*a_2+a_1$
$\Rightarrow$ $3*(a_3+a_2)+(a_3+a_1)-(a_2)$
$\Rightarrow$ $若 (a_3+a_1-a_2)\mid3$
$\Rightarrow$ $則 (a_3 a_2 a_1)_2 \mid 3$

考慮二進位數 $(a_4 a_3 a_2 a_1)_{2}$ ，不失一般性假設 $a_4=1$，則
$(a_4 a_3 a_2 a_1)_{2}=(1a_3 a_2 a_1)_{2}=(1000)_2*1+(0100)_2*a_3+(0010)_2*a_2+(0001)_2*a_1$
$\Rightarrow$ $8*1+4*a_3+2*a_2+a_1$
$\Rightarrow$ $(9-1)*1+(3+1)*a_3+(3-1)*a_2+a_1$
$\Rightarrow$ $(9-1)*a_4+(3+1)*a_3+(3-1)*a_2+a_1$
$\Rightarrow$ $3*(3*a_4+a_3+a_2)+(a_3+a_1)-(a_4+a_2)$
$\Rightarrow$ $若(a_3+a_1-a_2-a_4)\mid3$
$\Rightarrow$ $則 (a_4 a_3 a_2 a_1)_2 \mid 3$

由上不難看出，對一個正整數 a 而言， a 是否為 3 的倍數，關鍵在於a之二進位表示法中，在奇數位(odd_bit)的1和偶數位(even_bit)的1之個數的差值，是否為3的倍數，和上方的假設猜想一致。
只要重複上述的作法，對32位元的正整數做同樣的操作，一樣能夠得到一樣的結果。

上網搜尋了一下，看到一個以前在離散數學和計算理論常常使用的DFA:
我參考了這篇：https://www.geeksforgeeks.org/check-divisibility-binary-stream/

把要判斷的數字 N 之二進位表示法$(a_{32}a_{31}...a_{1})_2$當作input：
由高位到低位依序輸入，若是input結束前，停在recept state q0 代表 N 除以 3 餘數為 0 -> $N\mid 3$
反之，若是停在 q1 代表 代表 N 除以 3 餘數為 1 -> $N\nmid 3$
同理，若是停在 q2 代表 代表 N 除以 3 餘數為 2 -> $N\nmid 3$

分析DFA：
不難得到這個DFA所接受(accepted)的字串為： {$11$}$^{*}$$\cup$ { {$10$}$^{*}$ $\cup$ {1}$^{*}$ $\cup$ {01}$^{*}$} $\cup$ {$0$}$^{*}$

但是還得注意一點，在上述的DFA中，input被定義為是由數字的高位到低位的bit依序輸入，這並不直接符合右移的行為，上述程式中右移的行為是低位到高位依序輸入，這代表我們有兩種調整方式：

1.將整個bit-pattern作reverse(這是最直觀的作法，但效率較差，除了得事先對數字做處理，當數字不大的時候，還會需要不斷處理連續的0-bit)
2.翻轉accepted string（這種做法效率較高，理由和上述相反）

翻轉後的 accepted string: {$0$}$^{*}$$\cup$ { {$10$}$^{*}$ $\cup$ {1}$^{*}$ $\cup$ {01}$^{*}$} $\cup$ {$11$}$^{*}$
（注意：翻轉後的accepted string和翻轉前相同，但這只是巧合）

> 提示:
> HackMD 支援透過 Graphviz 和 Mermaid 製圖，請多利用
> [name=jserv]

考慮N=5的狀況：
解法1
====
由於$5=2^{2}+1$$\Rightarrow$$2^2=4=5-1$ ，這代表如果一次處理兩個bits，對5來說是最有效率的，由此我們盡可能把一個數N表示成
$N=4M+k=(5-1)M+k=5M+(k-M)$ ，把(k-M)的數值當成參數，傳入進行下一次遞迴呼叫。
例： $12_{10}$=$5*2+2$

考慮二進位數 $(a_4 a_3 a_2 a_1)_{2}$
$\Rightarrow$ $(a_4 a_3 a_2)_2*(100)_2 + (a_2 a_1)_2$
$\Rightarrow$ $(a_4 a_3 a_2)_2*(5-1) + (a_2 a_1)_2$
$\Rightarrow$ $5*(a_4 a_3 a_2)_2 -(a_4 a_3 a_2)_2 + (a_2 a_1)_2$
$\Rightarrow$ $5*(a_4 a_3 a_2)_2 -[(a_4 a_3 a_2)_2 - (a_2 a_1)_2]$

int isMult5(int  n)
{
    if (!(n ^ 5) || !(n ^ 0)) return 1;
    int k = n & 3;
    n = n >> 2;
    if (n - k & 0x8000) return 0;
    return isMult5(n - k);
}

解法2
====
事實上，由於5的倍數以十進位表示時，只要檢查其個位數即可，個位數必定為 5 或 10。
我們把數字 n轉成 char 型態後檢查最後一個 byte（char)是否為'5'或是'0'即可。

int isMult5_2(int n)
{
    char str[32];
    str = itoa(n, str, 32);
    if (str[0] == '5' || str[0] == '0') return 1;
    return 0;
}

解法3
=====
透過仿造 除數=3 時的算法來對 除數=5 時做分析：
考慮二進位數$(a_5 a_4 a_3 a_2 a_1)_2$
$\Rightarrow$ $(16a_5+8a_4+4a_3+2a_2+a_1) mod 5$
$\Rightarrow$ $a_5-2a_4-a_3+2a_2+a_1$
$\Rightarrow$ $(a_5 - a_3 + a_1 ) + 2(-a_4+a_2)$
也就是説$(a_5 a_4 a_3 a_2 a_1)_2$ 是否為5的倍數，取決於 $(a_5 - a_3 + a_1 ) + 2(-a_4+a_2)$ 是否為5的倍數
這種做法提供了我們透過遞迴呼叫來求解的契機。

考慮二進位數$(a_6 a_5 a_4 a_3 a_2 a_1)_2$
$\Rightarrow$ $(32a_6 16a_5+8a_4+4a_3+2a_2+a_1) mod 5$
$\Rightarrow$ $2a_6+a_5-2a_4-a_3+2a_2+a_1$
$\Rightarrow$ $(a_5 - a_3 + a_1 ) + 2(a_6-a_4+a_2)$

將其推廣到32bit的無號數 ： $(a_{32} a_{31}...a_{2} a_1)_2$

$\Rightarrow$ $(-a_{31}+a_{29}-a_{27}+....+a_1)$ $+2(-a_{32}+a_{30}-a_{28}...-a_{4}+a_{2})$

$\Rightarrow$ $\sum_{i=0}^{15}(-1)^{i} a_{2i+1}$ $+ 2\sum_{i=1}^{16}(-1)^{i+1} a_{2i}$

透過上面的分析，嘗試寫出和本週作業類似的遞迴求解funtion。

int isMult5_3(unsigned  n) {
 int odd_diff = 0, even_diff = 0; //  variables to count odd and even SET bits 

 int subtract_bit =   0;

 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
    {
     if(subtract_bit) odd_diff--;
     else odd_diff++ ;
    }
    n >>= 1;
     if ( (n &   1) ) // odd bit is SET, increment odd_C
    {
     if(subtract_bit) even_diff--;
     else even_diff++ ;
    }
    n = n >> 1;

    subtract_bit^=1 ; //toggle subtrack bit
}
 // Recursive call till you get 0/1 
 return(isMult5_3(abs(odd_diff + 2*even_diff)));
}

事實上，上述的分析，同樣可以透過分析除以5的DFA求出來（但是需要花一點時間作accepted string reduction），DFA可以和硬體唯一對應，也就可以和程式語言對應。

## 第 2 週測驗題 第一題
題目如下：
請完成下方程式碼，依循 IEEE 754 單精度規範，輸出 2 x 2x 的浮點數表達法，而且考慮兩種狀況：

- 當 x x 過小的時候，回傳 0.0 0.0
- 當 x x 過大的時候，回傳 +∞ +∞

注意：這裡假設 `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);
}

將 Y0 ~ Y7 填入相對應的數值，使得該函式依input: $x$ output得到 $x^2$ 的浮點數值。

首先考慮 input x 太小 (too small )的情況，在什麼情況下，我們將 x 稱作“太小”呢？
根據[IEEE754](https://zh.wikipedia.org/zh-hant/IEEE_754)定義：


當某個浮點數其exponent field為 0 且 fraction field 全為 0 時 ，被視為$\pm0$
其exponent field為 0 但 fraction field 不為 0 時 ，被視為$\pm Denormalized-number$
而最小的可表示的 $positive-Denormalized-number$為 $(0.00000000000000000000001)_2*2^{1-127}$
$\Rightarrow$ $(1.0)_2 *2^{-23}*2^{-126}$
$\Rightarrow$ $2^{-149}$
因此任何小於$2^{-149}$的在single precision都被視為too small而無法表示

/* too small */
    if (x < 2 - pow(2, Y0) - 23) {
        exp = 0;
        frac = 0;
    }

因此若$2^x < 2^{-149}$則無法表示
$\Rightarrow$ $x<-149=2-2^{Y_0}-23$
$\Rightarrow$ $Y_0=7$

再來考慮以下這段

/* denormalize */
    } else if (x < Y1 - pow(2, Y2)) {
        exp = 0; 
        frac = 1 << (unsigned) (x - (2 - pow(2, Y3) - Y4));
    }

如上面所述，exponent field為 0 但 fraction field 不為 0 時 ，被視為$\pm Denormalized-number$
而$\pm Denormalized-number$能夠表示的範圍為
$(0.111...11)_2*2^{1-127} 到 (0.000...01)_2*2^{1-127}$
$\Rightarrow$ $(2-2^{-23})*2^{-127} 到 2^{-149}$ 之間
$\Rightarrow$ $2^{-126}-2^{-149}$ 到 $2^{-149}$ 之間
由此可見，首先當$n\leq2^{-126}-2^{-149}$時會被視為$\pm Denormalized-number$
$\Rightarrow$ $2^x=n\leq2^{-126}-2^{-149}$
$\Rightarrow$ $x<-126$
$\Rightarrow$ $x<2-2^{7}=-126$
$\Rightarrow$ $Y_1=2$ 且 $Y_2=7$

之後，由於$n$為非正規化數，要注意：

非正規化數計算公式為： $(-1)^{sign}*(0.fraction)*2^{-126}$
而exp被設定為0（exp=0的狀況是特別保留給非正規化數和$\pm0$），這時候我們盡量左移(x+149)來還原fraction
$\Rightarrow$ $2-2^{Y3}-Y4 = -149$
$\Rightarrow$ $Y_3=7$ $Y_4=23$

接下來來考慮正規化數的情況：

/* normalized */
    else if (x < pow(2, Y5)) {
        exp = pow(2, Y6) - 1 + x;
        frac = 0;

    } 

最大的正規化數字： $(1.11..1)_2 * 2^{127} < 2^{128} = 2^{2^7}$
$\Rightarrow$ $Y_5= 7$
而在計算exponent的時候，因為IEEE754是以無號數來看待exp，得將exp扣去bias(bias在單精度中為127)
$\Rightarrow$ $exp= x - 127 = x - 2^{7} +1$
$\Rightarrow$ $Y_6=7$

最後考慮too large的情況：
在IEEE754中，當exp=255以及 farc=0時被視為 $\pm\infty$

 /* too large */
    } else {
        exp = Y7;
        frac = 0;
    }

由上所述，顯然$Y_7＝255=0xff$

**延伸題: 給定一個單精度的浮點數值，輸出較大且最接近的$2^{x}$值，需要充分考慮到邊界**
此時，考慮x為float-type
首先，考慮幾點
1. $2^{+\infty}$ return $+\infty$
2. $2^{-\infty}$ return $+0$
3. $2^{\pm0}$ return $1$
4. $2^{NaN}$ return $NaN$
5. $2^{x}>((2-2^{23})*2^{127})$ return $overflow$
6. $2^{x}<(2^{-149})$ return $out$ $of$ $range$

float nth_root(float a,int n)
{
    const int K = 6;
    flaot x[K] = {1};
    for(int k=0;k<K-1;k++)
    {
        x[k+1] = (1.0/n) * ((n-1) * x[k]+ a/ pow(x[k],n-1));
    }
    return x[K-1];
}
float exp2_fp_hw2(float x) {
     if(x==0) return 1;
    else if(x>127) return FLT_MAX ;
    else if(x<-149) return FLT_MIN ;
    else if(x<0) return 1/pow(2,-1*x);
    else if(x>0 && x<1)  return nth_root(2,1/x);
    else if(!(int)x%2) 
    {
        float half_pow = pow(2,x/2);
        return half_pow * half_pow;
    }
    else 
    {
        return 2 * pow(2,x-1);
    }
    
}