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

###### tags: `sysprog2020`

---
## 測驗 `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` 型態能表達的數值上界。請補完程式碼。

D1 = (c\) 2
D2 = (d) 1
### 解題想法

:::success
延伸問題:
1. 解釋上述程式碼運作原理;
:::
參考[浮點數除法程式](https://hackmd.io/@billsun/B1gii716N)

已知 $\dfrac{被除數}{除數}$ 等價於 $\dfrac{\dfrac{被除數}{n}}{\dfrac{除數}{n}}$ 。

第5行：遞迴終止條件為`除數==1`或`被除數==0`
第7行：od 作用為判斷 slot (除數) 是否為奇數
此題第 8 行的行為是:
若 slot 為奇數，則將 `被除數/2` 與 `(除數+1)/2` 做為遞迴參數傳入，且須修正誤差 (第9~10行)；
若 slot 為偶數，則將 `被除數/2` 與 `除數/2` 做為遞迴參數傳入。
因此 ==D1 = (c\) 2== ， ==D2 = (d) 1==

誤差修正:
舉例: `1 / 3`
會先化成 1/4 = 0.25
>再算 0.25/3
>再化成 0.25/4 = 0.0625
>>再算 0.625/3
>>再化成0.625/4 = 0.015625
>>>再算 0.15625/3
>>>再化成 0.15625/4 = 0.00390625
>>>>再算 0.00390625/3
>>>>再化成 0.00390625/4 = 0.0009765625
>>>>>再算 0.0009765625/3
>>>>>...依此類推

最終結果為 0.25 + 0.0625 + 0.015625 + 0.00390625 + 0.0009765625
=0.333078125
約為 $0.\overline{3}$

:::success
延伸問題:
2. 以編譯器最佳化的角度，推測上述程式碼是否可透過 tail call optimization 進行最佳化，搭配對應的實驗;
- 搭配閱讀 C 語言: 編譯器和最佳化原理 及 C 語言: 遞迴呼叫
3. 比照 浮點數運算和定點數操作，改寫上述程式碼，提出效率更高的實作，同樣避免使用到 FPU 指令 (可用 gcc -S 觀察對應的組合語言輸出)，並與使用到 FPU 的實作進行效能比較

:::

//TODO

---
## 測驗 `2`
延續 從 $\sqrt{2}$ 的運算談浮點數，假設 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
- √Square Roots
請補完程式碼，使上述程式碼得以運作。

QQ0 = (a) 0x01000000
QQ1 = (d) 0x3f000000
QQ2 = (g) 23

### 解題想法
單精度float32格式如下圖:
來源:[你所不知道的 C 語言: 浮點數運算](https://hackmd.io/@sysprog/c-floating-point)

34行：將 int32 轉成 float 格式，結果放至 ix0
36~47行：例外處理 0, +INF, NaN
49行：m = sign + exponent 部分
50~54行：處理小於 1 大於 0 的數 (`m==0`的情況)，需進行正規化
>把 `0.000123...` 轉成 `1.23...` 左移 i 次，所以m的正確值應為 `m-i-1`
>>-1原因:令數=0.1，i=0，正確應為m-1

//TODO

91行 要把 m 存的 exponent 存回原數的正確位置，因此應該要放在 31~23 bit 上，所以這邊選擇左移 23 bit，所以答案為 ==(g) 23==

:::success
延伸問題:
1. 解釋上述程式碼何以運作
- 搭配研讀 以牛頓法求整數開二次方根的近似值 (第 19 頁)
2. 嘗試改寫為更有效率的實作，並與 FPU 版本進行效能和 precision /accuracy 的比較
3. 練習 LeetCode 69. Sqrt(x)，提交不需要 FPU 指令的實作
:::

//TODO

---
## 測驗 `3`
LeetCode 829. 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+\frac{(k−1)k}{2}=N$

整理後可得：
$kx=N−\frac{(k−1)k}{2}$

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

從而得到近似解：
$k<\sqrt{2N}$

確認 k 的範圍後，即可走訪數值。

參考實作如下:

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

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

作答區

ZZZ = (d) 1 - i

### 解題想法

根據題目定義
要找出輸入的數能有幾種連續正整數和的寫法。
要能判斷是否能寫成連續正整數和，只要滿足 $(N-\frac{(k−1)k}{2})\ mod\ k=0$ 即可。
題目中`x += ZZZ`就是在做 $(N-\frac{(k−1)k}{2})$ 的部分。
`k`在程式中的變數名稱為`i`。
$\frac{(k−1)k}{2}$ 代表從`1`加到`k-1`，同義於程式中的`1`加到`i-1`。
$(N-\frac{(k−1)k}{2})$ 代表 N-1, N-2, ... 往下遞減。
`i`從2開始。

因此選 ==ZZZ = (d) 1 - i==。
:::success
延伸問題:
1. 解釋上述程式碼何以運作
:::
第3~4行：若 N 小於 1 則直接回傳 0，因為 1 沒辦法寫成其他組連續整數和。
>個人想法:這邊應該改成 `N<2` 才對。 2 好像沒辦法表示成連續正整數和，最小應該是 3 = 1+2。
>
`ret`用來計算能寫成幾組正整數解
`x`為輸入值`N`
第7~11行：從 N-1 往下加，若加總能整除 i 則表示能表示成一組連續正整數解 $(x)+(x+1)+(x+2)+...+(x+(i−1))$ ，`ret`加 1
最後回傳ret的值

:::success
延伸問題:
2. 嘗試改寫為更有效率的實作

:::

//TODO