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

## 測驗一

#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 / 2, od ? (slots + 1) >> 1 : slots >> 1);
    if (od)
        result += divop(result, slots);
    return result;
}

### 程式碼討論

假設舉例 $orig=2.4$ 、 $slots=3$ 我們預期得到的結果會是 0.8 。
根據遞迴的方式，我們發現浮點數 2.4 除 3 之後的結果會是
`divop(2.4, 3)` + `divop(0.6, 3)` +`divop(0.3, 3)` +...=0.8

>我們發現第二項的參數傳入其實是由第一項計算出來的答案繼續往下遞迴的，事實上，如果 $slots$ 中有除了2的冪次以外的因數，其除法就不能直接以二進制表示，需要遞迴運算直到 `result` 的 `fraction` 紀錄不下為止。如果以 $2.4/3$ 為例，就運算了 **1067** 次之多。

分析程式碼的部份:
先看當 $slots$ 是偶數的情形

$divop(orig/2, slots>>1) = \frac{\frac{orig}{2}}{\frac{slots}{2}}$

我們可以把它理解為被除數和除數先提出一個 2 的公因式再繼續下去運算，代表的值沒有改變，接下來看奇數的情形:

$divop(orig/2, slots+1>>1)$

$= \frac{\frac{orig}{2}}{\frac{slots+1}{2}}$

$= \frac{(orig)slots}{(slots+1)slots}$

$= \frac{orig(slots+1)}{(slots+1)slots}-\frac{orig}{(slots+1)slots}$

$= \frac{\frac{orig}{2}}{\frac{slots}{2}}-\frac{orig}{(slots+1)slots}$

$= result$
需要加上 $\frac{orig}{(slots+1)slots}$ 得到的才會是我們真正想要的解 $\frac{\frac{orig}{2}}{\frac{slots}{2}}$
$\frac{orig}{(slots+1)slots}$ 可以再表示成 $\frac{\frac{orig}{2}}{(slots+1)\frac{slots}{2}}=\frac{\frac{orig}{2}}{\frac{slots+1}{2}}/slots$ ，也就是 $\frac{result}{slots}$ ，因此改寫成 $divop(result, slots)$

Ans:
**D1 = 2
D2 = 1**

:::info
:bulb:當除數過大或是一連串奇數相乘，遞迴呼叫過多次效率反而比直接使用除法(/)降低許多，探究其原因是因為多層遞迴呼叫的關係，不曉得使用時機是什麼時候？或許是當除數比較小的時候所實現的快速浮點數除法？
:::

---
### 比照浮點數運算和定點數操作，改寫上述程式碼，提出效率更高的實作，並與使用到 FPU 的實作進行效能比較
參照 IEEE754 對於 `double` 的規範，實做較為迅速的 `div2` :

double div2(double dou) {
    u_int64_t ud = *(u_int64_t *)&dou;
    // assert(sizeof(u_int64_t), sizeof(long long int));
    printf("uint64 %llu\n", (long long int)ud);
    u_int64_t ud_sgn = (ud >> 63) & 0x01; //extract sign bit
    u_int64_t ud_exp = (ud >> 52) & 0x7ff;//extract exponent
    if (ud_exp == 0) {//denormalized
        ud = ud >> 1;
        ud &= 0x7fffffffffffffff;
        ud += (ud_sgn <<63);
    }
    else if (ud_exp != 0x7ff) {//+Inf or NaN
        ud_exp -=1;
        ud &= 0x800fffffffffffff;
        ud = ud + (ud_exp << 52);
        printf("uint64 %d\n", ud);
    }
    return *(double*)&ud;
}

這邊直接使用了 double 的 biteise 運算實做除法，按照以下步驟完成:
1. 浮點數轉型成 `uint64` 的方式儲存，方便我們直接對位元操作
2. 用 `bit mask` 分別抓出 `sign bit` 和 `exponent` 項
3. 如果 `exponent` 項全為零，代表 denormalized ，直接把 `fraction` **右移**處理然後把 `sign bit` 放回
4. 如果 `exponent` 項不全為1，代表為 normalized number ，**`exponent` -1** 後將 `sign bit` 和 `exponent` 項放回
5. 其他種情形(Inf、NaN)直接回傳

#### 效能檢測

---
## 測驗二
- 表一：IEEE754浮點數規範

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

這邊的兩個巨集分別是
1. 將32 bit int 的 `ix0` 轉成 float 的 `d`
2. 將float 的 `d` 轉成 32 bit int 的 `ix0`

float z;
    int32_t sign = 0x80000000;
    uint32_t r;
    int32_t ix0, s0, q, m, t, i;

接著是變數的宣告，並且我們注意到 `sign` 代表的意思是 32 bit int 的第一個 `sign bit` 的 mask ，而 `ix0` 代表的是我們要把浮點數轉成整數型態儲存

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

處理傳入浮點數是正負零或是負數的情況。如果是正負零，我們用 `(~sign)` 判斷除了 `sign bit` 以外的其他位元是不是都為零。如果全為零，則直接回傳0。

如果是負數，回傳 0/0 ，也就是 `NaN`。

  /* take care of +INF and NaN */
    if ((ix0 & 0x7f800000) == 0x7f800000) {
        /* sqrt(NaN) = NaN, sqrt(+INF) = +INF,*/
        return x;
    }

現在的 `ix0` 已經是一個正數，但是我們還需要處理原本的浮點數如果是(正)無限大或是 `NaN` 的情況，根據表一，擔這兩種狀況發生的時候，中間的 Exponent 項會全部都是 1，所以我們需要用一個 mask `0x7f800000` 去判斷 `ix0` 指數項以排除(正)無限大或是 `NaN` 的情況。

/* normalize x */
    m = (ix0 >> 23);
    if (m == 0) { /* subnormal x */
        for (i = 0; (ix0 & 0x00800000) == 0; i++)
            ix0 <<= 1;
        m -= i - 1;
    }

右移了 23 個位元後得到 `m` 是浮點數的指數項( exponent )，如果全部為零，排除掉0的可能性只可能會是 `denormalized` 浮點數，將其正規化，方式是將 `fraction` 的部份右移直到他的第一個位元是 1 ，而這時候我們的 `m` 也會相對應的變小。

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

減掉了 bias 的 127 之後才是真正的指數次方。
透過`ix0 = (ix0 & 0x007fffff) | 0x00800000` 擷取 `fraction` 的部份，並且設定第 24 bit 為 1 ，得到 $1.fraction$

根據[以牛頓法求整數開二次方根的近似值](http://science.hsjh.chc.edu.tw/upload_works/106/44f6ce7960aee6ef868ae3896ced46eb.pdf)中得到 $\sqrt{N}$ 在 x 軸的第一個近似值為:
$\sqrt{N}$
$=\frac{a(n-1)+\frac{N}{a^{n-1}}}{n}$
$=\frac{a+\frac{N}{a}}{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);
	

## 測驗三

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

---
### 程式碼解析
從原方程式的

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

我們將它理解為

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

也就是 $N$ 減掉從 1 到 k-1 之間連續正整數的和，假如我們能夠找到 $k=i$ ，且 $k$ 為正整數，使得等號成立，那必定存在一組連續正整數的和等於 $N$ ，也就是說， $N-\{1+2+...+(k-1)\}$ 能夠被 $x$ 整除，這個狀況下的 $k$ 便是一個可行解。

回到題目本身
**第3行:** 當 $N<1$ 的情況無法形成任何連續正整數的和，因此直接回傳0。
**第5行:** 對於任何大於 1 的正整數來說，至少都能夠表示成只有自己的連續正整數和，因此起始值為一，這邊其實有另一個意思，也就是說當 $k=1$ 時，$x=N-\{(1-1)\}$ 的時候，只要令 $x=N$ ，則頂事必定會成立。
我們用 9 當作例子：
|iteration| x | i | ret |
| --- | --- | --- | --- |
|1|9|2|1|
|2|8|3|2|
|3|6|4|3|
最後回傳的答案是 3，為什麼呢？我們令 $N=9$ 以及 $k=2$ (這邊的 $k$ 和 $i$ 是相同意思，因為 $k=1$ 我們已經討論過必定成立，所以程式碼從 $k=2$ 開始)代入看看:
- **第一次遞迴 $k=2$**

$2x=9-\{(2-1)\}$
因為可以整除所以 `ret+1`

- **第二次遞迴 $k=3$**

$3x=9-\{(2-1)+(3-1)\}$
也可以整除所以 `ret+1`

從這邊我們已經可以看出端倪，正好可以解釋程式碼的實做方式:
**第8行:** 一開始的初始值是 $N$ ，接著，每次的迴圈中比上一次等號右手邊的值還要再多減去一個 $(k-1)$ ，因此我們便不用每個迴圈都重新計算前面的整數和。
**第9行:** 接著看向等號的左手邊，如果右手邊的值，也就是程式碼中的 $x$ ，能夠整除等號左邊的 $k$ ，也就是程式碼裡面的 $i$ ，我們可以知道這是一個可行解，因此 `ret` 要加一。
**由此我們可以得到 ZZZ= -(k-1)**

- **第三次遞迴 $k=4$**

$4x=9-\{(2-1)+(3-1)+(4-1)\}$
在程式碼當中，如果當前的 $x$ 已經小於 $i$ 了，意味著這次遞迴當中等號右手邊的值甚至會比 $i$ 還要小，因此不可能被整除。

---

### 嘗試更有效率的實做
程式碼本身已經是精簡之後的結果了，以演算法本身來說想不太到比老師的作法更加快速的方法，參考了 [huang-me](https://hackmd.io/@huangme/embedded_quiz5) 同學筆記下老師的建議，分析 `%` 是否能夠被省略。
使用之前的方法，假定:

u_int64_t M =((u_int64_t)(__INT64_C(0xFFFFFFFFFFFFFFFF) / i + 1));
if ( x * M <= M-1)
    ret++;

結果出來時間並沒有比較快，應該是因為還是用到除法的關係。
試過其他方法，結果都超時。

:::info
TODO:繼續尋找更快速的方法
:::