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

## 程式運作原理
### 測驗一：
程式碼運作原理如下：

$$
\begin{aligned}{{\dfrac{a}{b}}-{\dfrac{a}{b+1}}}&={\dfrac{a}{b\cdot(b+1)}}
\\&={\dfrac{\dfrac{a}{b+1}}{b}}
\end{aligned}
$$
換句話說，要計算 ${\dfrac{a}{b}}$ 時，可以換成計算 ${\dfrac{a}{b+1}} +{\dfrac{\dfrac{a}{b+1}}{b}}$ 。
* 需要討論的事項：
1. 比較該除法與一般除法的效能（實驗）
### 測驗二：
為了集中討論程式碼的核心，我將題目的程式碼縮減成如下：

#include <stdint.h>

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)

    float ieee754_sqrt(float x)
    {
        float z;
        uint32_t r;
        int32_t ix0, s0,q,m,t,i;

        EXTRACT_WORDS(ix0, x);
        m = ix0 >> 23;
        m -= 127;
        ix0 = (ix0 & 0x007fffff) | 0x00800000;
        if (m&1) {
            ix0 += ix0;
        }
        m >>= 1;
        q = s0  = 0;
        r = 0x01000000;
        ix0 += ix0;
        while (r != 0) {
            t = s0 + r;
            if (t <= ix0) {
                s0 = t + r;
                ix0 -= t;
                q += r;
            }
            ix0 += ix0;
            r >>= 1;
        }
        ix0 = (q >> 1) + 0x3f000000;
        ix0 += (m << 23);

        INSERT_WORDS(z, ix0);
        return z;

    }

根據 IEEE-754 ，浮點數格式轉換成科學記號的關係如下：$$
(-1)^{sign}\space\cdot\space(1.fraction)_2\space\cdot\space2^{exponent \space - \space bias}
$$
在 `sign` bit 為 0 (這裡我們只考慮正數)時，平方根可以用
$$\sqrt{(1.fraction)_2}\space\cdot\space2^{\dfrac{exponent \space - \space bias}{2}}$$
表示。於是，我們可以將計算拆成 `fraction` 和 `exponent` 兩部份。
以下是取得 `fraction` 和 `exponent` 的步驟：

* 第一步：取得指數

        m = ix0 >> 23;
        m -= 127;

這裡是根據 `IEEE-754` 的規範，將指數的部份取出，並減去 bias = 127 。
另外，指數有可能為奇數，之後除以 2 時要把多出的 2 乘回 `1.fraction`，這個想法對應到程式碼

        if (m&1) {
            ix0 += ix0;
        }

* 第二步：取得 `fraction`

ix0 = (ix0 & 0x007fffff) | 0x00800000;

這裡取得 `fraction` 的部份只有 `(ix0 & 0x007fffff)` 這裡，真正的目的是令 `ix0` 值為$(1.{\overbrace{fraction}^{23 個 bit}})_2$ ，binary point 在從 LSB 數來第 23 個 bit 和第 24 個 bit 之間。
到此取得了 fraction 、 exponent。平方根的指數部份能簡單求出，重點是 $\sqrt{(1.fraction)_2}$ 的求法。參考 [YLowy 的筆記](https://hackmd.io/@YLowy/HJkXHWCIv) 後，我發現我需要了解開 2 次方根的運算原理。
我查到了[開 n 次方根的直式計算與原理](http://www.hwsh.tc.edu.tw/ischool/public/resource_view/open.php?file=330d84facabdf8dff79800f4efc3fa33.pdf) 這篇。
以求 `根號2` 為例：
1. 根號2
2 的 $(1.fraction)_2$ 左移兩位。總共26 bit。此時 binary point 定在從 LSB 數來第 24 和第 25 個 bit 之間。

10|00|00|00|00|00|00|00|00|00|00|00|00

為了討論方便，在做開平方根的直式計算時，暫且將

10|00|00|00|00|00|00|00|00|00|00|00|00

視為整數來計算。
以下是計算流程：

  q                        1|  0| 1|
                        --------------------------------------------------
 ix0              √ 10|00|00|00|00|00|00|00|00|00|00|00|00
t = 1                    1
s0 = 10           
                        --------------------------------------------------
t = 101                100
s0 = 100              000
                        -------------------------------------------------
 t = 1001              10000
 s0 = 1010              1001

`q` 欄位代表平方根 ， `ix0` 欄位代表代表要被取平方根的數，`t` 欄位代表將要從 `ix0` 扣除的數， `s0` 欄位代表 `q` 左移一位的結果。
這種直式計算，用到的原理其實就是國中學的乘法公式：
若
$$ x = {(a + \beta)}^{2} = a^{2} + (2a+\beta)\beta$$
則
$$ x -a^{2} = (2a+\beta)\beta$$
若
$$\beta = b + \gamma$$
則
$$x - (a+b)^{2} = (2a + 2b + \gamma)\gamma$$
試著把上述兩個式子套到直式計算的流程：
$$a^{2} = (1{\overbrace{(00...00)}^{13 個 bit}})^{2} <= ix0 < (10{\overbrace{(00...00)}^{13 個 bit}})^{2}$$ 所以
$$ix0 = (1{\overbrace{(\beta)}^{13 個 bit}})^{2}=(a + \beta)^{2}$$
$$(a + b)^{2} = (1{\overbrace{(00...00)}^{13 個 bit}})^{2} <= ix0 < (1{\overbrace{(10...00)}^{13 個 bit}})^{2}$$
$a^{2} + 2ab+b^{2} = a^{2} + (2a+b)b$
所以
$$ix0 = (1{\overbrace{(\beta)}^{13 個 bit}})^{2}=(a + b+ \gamma)^{2} =(10 {\overbrace{(\gamma)}^{12 個 bit}})^{2}$$
由於
$$a^{2} + 2ab+b^{2} = a^{2} + (2a+b)b$$
比較 $(a+b)^{2}$是否大於 ix0 時，只要比較$(2a + b)b$ 是否大於 $ix0 - a^{2}$即可。
我們求平方根的求法是由左往右 bit by bit 的求，因此 b 永遠是當前 q 的 LSB，非 0 即 1。換句話說，只須比較 $(2a + 1)$ 是否大於 $ix0 - a^{2}$，就可以知道 b 是否為 1 。這就是

        while (r != 0) {
            t = s0 + r;
            if (t <= ix0) {
                s0 = t + r;
                ix0 -= t;
                q += r;
            }
            ix0 += ix0;
            r >>= 1;
        }

這段程式碼作的工作。另外的 `ix0 += ix0` 是為了讓 `s0` 可以不用 shift。
順帶一提，如果要更符合數學式的話，可以寫成

        while (r != 0) {
            t = s0 + r;
            if (t <= ix0) {
                ix0 -= t;
                q += r;
                s0 = q << 1;
            }
            ix0 += ix0;
            r >>= 1;
        }

### 測驗三：
從 x 開始， k 個連續正整數合為 N 的話，可以表示成如下關係式：
$kx = N \space - \space \dfrac{(k-1)k}{2}$
於是我們要檢查特定 k 是否能使上述關係式成立時，我們只須確認 $N \space - \space \dfrac{(k-1)k}{2}$ 是否為 k 的倍數即可。
接下來關於 k 的範圍，我們知道 $N \space - \space \dfrac{(k-1)k}{2} > 0$ ，故將其設為迴圈中止條件即可。