# 2020q3 Homework5 (quiz5)
contributed by < `JimmyLiu0530` >

###### tags: `進階電腦系統理論與實作`

## 測驗1: 浮點數除法
考慮到以下浮點數除法程式: (fdiv.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;
}

假設 `divop()` 的第二個參數必為大於 0 的整數，而且不超過 `int` 型態能表達的數值上界。請補完程式碼。

### 程式說明

原理:

$\dfrac{N}{d}$ = $\dfrac{N}{d+1}*\dfrac{d+1}{d}$ = $\dfrac{N}{d+1}*(1+\frac{1}{d})$ = $\dfrac{N}{d+1}+\dfrac{N}{d(d+1)}$

e.g. $\dfrac{N}{5}=\dfrac{N}{5+1}+\dfrac{N}{5(5+1)}$

Line 8 的作用在於:

若除數 `slots` 是偶數，那麼被除數 `orig` 與除數 `slots` 可以先同時除以 2，再進行運算，因此 `D1 = 2`。
若 `slots` 為奇數，則根據上述的原理 $\dfrac{orig}{slots}$ = $\dfrac{orig}{slots+1}+\dfrac{orig}{slots(slots+1)}$ 告訴我們兩件事:
1. 首先就是 `D2 = 1`。又因為 `slots` 原本是奇數加 1 後變成偶數，所以也可以與 `orig` 同除以 2。
2. 再來就是最終結果 `result` 除了計算 $\dfrac{orig}{slots+1}$ 外還要再加上 $\dfrac{orig}{slots(slots+1)}$ ，因此 line 10 需要 `divop(result, slots)`。
### 延伸問題

#### 1. 以編譯器最佳化的角度，推測上述程式碼是否可透過 [tail call optimization](https://en.wikipedia.org/wiki/Tail_call) 進行最佳化，搭配對應的實驗;
> 搭配閱讀 [C 語言: 編譯器和最佳化原理](https://hackmd.io/@sysprog/c-compiler-optimization?type=view) 及 [C 語言: 遞迴呼叫](https://hackmd.io/@sysprog/c-recursion?type=view)

#### 2. 比照 [浮點數運算和定點數操作](https://hackmd.io/@NwaynhhKTK6joWHmoUxc7Q/H1SVhbETQ)，改寫上述程式碼，提出效率更高的實作，同樣避免使用到 FPU 指令 (可用 `gcc -S` 觀察對應的組合語言輸出)，並與使用到 FPU 的實作進行效能比較

------------------------------------------------
## 測驗2: 平方根
延續 [從 √2的運算談浮點數](https://hackmd.io/@sysprog/ryOLwrVtz?type=view)，假設 float 為 32-bit 寬度，考慮以下符合 IEEE 754 規範的平方根程式，請補上對應數值，使其得以運作：

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

Reference:
- [Floating point division and square root algorithms and implementation in the AMD-K7/sup TM/ microprocessor](https://ieeexplore.ieee.org/document/762835/)
- [√Square Roots](http://www.azillionmonkeys.com/qed/sqroot.html)
- [以牛頓法求整數開二次方根的近似值](http://science.hsjh.chc.edu.tw/upload_works/106/44f6ce7960aee6ef868ae3896ced46eb.pdf) (第 19 頁)

請補完程式碼，使上述程式碼得以運作。

### 程式說明
- **牛頓法**求整數開根號近似值的**核心概念**:
- 假設估算 $\sqrt{N}$，則令$f(x)=x^2-N$，利用 $f(x)=0$ 求得 $x=\pm\sqrt{N}$。
- 接著，取一個估算值 a，對 $(a,f(a))$ 做切線，並與 x 軸相交於 b。一樣再對 $(b,f(b))$ 做切線與 x 軸相交於 c，重複上述步驟即可逼近 $\sqrt{N}$。

### 延伸問題

#### 1. 嘗試改寫為更有效率的實作，並與 FPU 版本進行效能和 precision /accuracy 的比較

#### 2. 練習 LeetCode [69. Sqrt(x)](https://leetcode.com/problems/sqrtx/)，提交不需要 FPU 指令的實作

------------------------------------------------
## 測驗3: LeetCode [829. Consecutive Numbers Sum](https://leetcode.com/problems/consecutive-numbers-sum/)

給定一個正整數 N，問 N 能寫成多少種連續正整數之和，例如 9 可寫成
$4+5$，或者 $2+3+4$。由於要寫成連續正整數之和，則可用等差數列來思考，且差值為 1，這個等差數列不必從 1 開始，假設從 x 開始，且個數共有 k 個，則可寫出該等差數列為：
$x, x+1, x+2,..., x+(k-1)$

令其和為 N，根據等差數列的求和公式，可寫出以下等式：

$kx+\dfrac{(k-1)k}{2}=N$

整理後可得：

$kx=N-\dfrac{(k-1)k}{2}$

對任意 k 值，倘若 x 能得到正整數解，就表示必定會有個對應的等差數列和為 N。由於 k 是等差數列的長度，首先必定大於 0，即為下限。由於 x 也必是正整數，可得:

$N-\dfrac{(k-1)k}{2}>0$

從而得到近似解：

$k<\sqrt{2N}$

參考實作如下:

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

請補完程式碼，使上述程式碼得以運作。

### 程式說明
根據題目的敘述及提示可知 $k$ 為正整數且 $k<\sqrt{2N}$，有了 k 值的範圍後，在這範圍內一一去檢驗 $kx=N-\dfrac{(k-1)k}{2}$ 是否成立就可以知道 N 能不能寫成 k 個連續整數相加。

因為任何正整數一定能寫成 1 個數相加(即自己)，所以直接從 k=2 開始看。

接著用一個例子來說明如何將數學式轉換成程式碼: (N=9為例)
| k | $kx=N-\dfrac{(k-1)k}{2}$ | x |是否能寫成 k 個連續整數相加 |
| - | - | - | - |
| 2 | $2x=9-\dfrac{1*2}{2}=8$ | 4 | :O: |
| 3 | $3x=9-\dfrac{2*3}{2}=6$ | 2 | :O: |
| 4 | $4x=9-\dfrac{3*4}{2}=3$ | $\dfrac{3}{4}$ | :X: |

根據上面例子會發現
1. 每一次迭代，$N-\dfrac{(k-1)k}{2}$ 的值都會比上一次等差地減少。
- e.g. k=2 時，值為 8 會是上一輪的 9 減掉 **1**; k=3 時，值為 6 會是上一輪的 8 減掉 **2**, 以此類推，因此 `ZZZ = 1-i`。
2. $2x = 8$ 有整數解釋因為 8 被 2 整除，然而 $4x = 3$ 無整數解是因為 3 不被 4 整除，因此可以推廣到 $kx=m$ 有整數解就代表 m 要能被 k 整除，呼應 line 9 的判斷式。

- 執行效能:

### 延伸問題
#### 1. 嘗試改寫為更有效率的實作

int consecutiveNumbersSum(int N){
    int res = 1, count;
        while (N % 2 == 0) N /= 2;
        for (int i = 3; i * i <= N; i += 2) {
            count = 0;
            while (N % i == 0) {
                N /= i;
                count++;
            }
            res *= count + 1;
        }
        return N == 1 ? res : res * 2;
}

> TODO: 透過 shift 改寫
> :notes: jserv

- 此方法的原理：
由上面的推論可得 $kx=N-\dfrac{(k-1)k}{2}$ ，變換一下得到 $2N=k(2x+k-1)$，因為 $2N$ 是偶數，代表著 $k$ 或是 $(2x+k-1)$ 要是偶數，而當 k 偶數，$(2x+k-1)$ 就會是奇數; 當 k 奇數，$(2x+k-1)$ 就會是偶數。

也就是說，N 的每一個 odd factor 都會對應一組答案。因此，此問題就變成了 N 有幾個 odd factor。

- 計算 odd factor 數量:
> 若 $N=3^a*5^b*7^c*...$，則 N 共有 (a+1)(b+1)(c+1)... 個因數
> ( 因為在 line 3 會將所有2因數移除，所以才可以如此假設 N )

- 程式碼說明:
上面這個新的程式碼就是用來找 `N` 有多少個 odd factor。首先，因為我們只關心 odd factor，所以先移除 `N` 的所有 2 因數，再對所有小於 `N` 的奇數去看是否為 `N` 的因數。若是，則去計算該因數有多少次方(即 `count`)，加 1 後再乘入 `res` 中。

最後，for 迴圈結束檢查若 `N` 為 1，那麼直接回傳 `res`; 否則代表 N 中還剩一個**最大質因數** P，且其次方必為1，因此回傳 `res * (1+1)`。

> 因為 line 4 for 迴圈的中止條件為 `i * i <= N`，所以執行完 for 迴圈後 N 有可能剩下的就只有其最大質因數 P (因為 (1) 若為合數，那麼早在 for 迴圈就被除掉了。(2) 假設 N 還存在另一個質因數 Q，Q < P，則代表 for 迴圈沒有執行到 `i = Q`，即 Q * Q > N。又 ${N}\geq{P}*Q>Q*Q$，很明顯矛盾，因此 N 中不存在另一個質數)，使得 P * P > N，也因為 P * P > N，所以 P 的次方必為 1。

- 執行效能: