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

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

- ==[測驗題目](https://hackmd.io/@sysprog/2020-quiz5)==
- ==[作業區](https://hackmd.io/@sysprog/2020-homework5)==
- ==[Github](https://github.com/joey3639570/2020sysprog_quiz5)==

## 目錄
[TOC]

## 測驗 `1` - 浮點數除法程式
### 程式原理

#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(32,4)` 試試看:
也就是以 32 為 `orig` ， 4 為 `slots`：

- 第一次的 `od` 為 `4 (100) & 1 = 0`
會發現在 result 這一行會進入遞迴式，如下
`result = divop(32 / D1, 100 >> 1);`

- 下一次的 `od` 為 `10 & 1 = 0`
再下一次的遞迴式變成了:
`result = divop(32 / D1 / D1, 10 >> 1);`

- 再下一次的遞迴式在分支會檢測到 `slots == 1` ，
故 `result` 為 `32 / D1 / D1`

`divop(32, 4)` 答案應該得到 `8`，故 `result = 32 (100000) / D1 / D1 = 1000`， **D1 應為 2**

再以 `divop(32,5)` 為例子帶入:
也就是以 32 為 `orig` ， 5 為 `slots`

- 第一次的 `od` 為 `5 (101) & 1 = 1`
會發現在 result 這一行會進入遞迴式，如下:
`result = divop(32 / D1, (100 + D2) >> 1);`

- 下一次的 `od` 為 `(100 + D2) >> 1 & 1`，**在此我們必須先假設 D2 的值才能知道 `od`**，假設 `D2 = 1`，則 `od` 為 `(101 + 1) >> 1 & 1 = 1`，則 `result` 在此為 `divop(32 / D1 / D1 , (((slots + D2) >> 1) + D2) >> 1)`

- 再下一次的遞迴式，`slots` 為 `(((101+1) >> 1) + 1) >> 1 = 10`， `od` 因此為 0 ，`result` 在此為 `divop(32 / D1 / D1 / D1 , (((slots + D2) >> 1) + D2) >> 1 >> 1)`，

- 於這次的遞迴可以得到 `slots = (((101+1) >> 1) + 1) >> 1 >> 1= 1`，直接 return 當前的 `orig`， `D1` 為 2 的狀況下， `orig = 32 / 2^3 = 4`

- 藉走到最後一個遞迴式往回找，注意到 `od` 為 1 的遞迴式需要再執行 `result += divop(result, slots);`，最後結果應為 `result += divop(4 + divop(4, 3))`

#### 藉表格整理遞迴式及其相對應結果:

`divop(1, 3)`
| `orig` | `slots` | `od` | return 值 |
|:------:|:-------:|:----:|:---------:|
| 1 | 3 | 1 | 1 |
| 0 | 2 | X | 1 |

`divop(4, 3)`
| `orig` | `slots` | `od` | return 值 |
|:------:|:-------:|:----:|:---------:|
| 4 | 3 | 1 | |
| 2 | 2 | 0 | 1 |
| 1 | 1 | X | 1 |

由 `fdfdd1234562` 整理:
>如果除數是偶數，就直接回傳第八行 `result`
因此可以合理判斷第八行是將除數與被除數都除以二

>接著如果除數是奇數，則會再加上第十行的計算，因此可以用數學推算出遞迴的結果

>$result={orig \over (slots+D2)}+{orig \over (slots+D2)^2}+{orig \over (slots+D2)^3}...$
>由等比級數公式可以知道 $result={{orig \over (slots+D2)}*(1-1 / (slots+D2)^ \infty ) }/{(1-1/(slots+D2))}={orig/(slots+D2) \over 1-1/(slots+D2)}={orig \over slots+D2-1}$

### 延伸問題

- 以編譯器最佳化的角度，推測上述程式碼是否可透過 tail call optimization 進行最佳化，搭配對應的實驗;

相關參考資料：[你所不知道的C語言：遞迴呼叫篇](https://hackmd.io/@sysprog/c-recursion?type=view)

參考 `hankluo6` :
>推測此程式第一個遞迴函式能透過 TCO 最佳化，而第二個不能，因程式最後不為 function call，而是 result 的加法動作。

#### 最佳化之比對
此部分可於編譯時使用 `-foptimize-sibling-calls` 觀察組合語言前後差異
>參見:
>[3.11 Options That Control Optimization](https://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html)
>`-foptimize-sibling-calls`
>>Optimize sibling and tail recursive calls.
>>Enabled at levels -O2, -O3, -Os.

>參考 :
>[How do I check if gcc is performing tail-recursion optimization?](https://stackoverflow.com/questions/490324/how-do-i-check-if-gcc-is-performing-tail-recursion-optimization)
透過觀察與比對優化前後組合語言是否有 `call` 或使用 `jne` 形成 loop 可以了解編譯器是否有使用 **tail call optimization**

##### 原程式碼經編譯後的組合語言

divop:
        pushq   %rbp
        movq    %rsp, %rbp
        subq    $32, %rsp
        movsd   %xmm0, -24(%rbp)
        movl    %edi, -28(%rbp)
        cmpl    $1, -28(%rbp)
        je      .L2
        pxor    %xmm0, %xmm0
        ucomisd -24(%rbp), %xmm0
        jp      .L3
        pxor    %xmm0, %xmm0
        ucomisd -24(%rbp), %xmm0
        jne     .L3
.L2:
        movsd   -24(%rbp), %xmm0
        jmp     .L5
.L3:
        movl    -28(%rbp), %eax
        andl    $1, %eax
        movl    %eax, -12(%rbp)
        cmpl    $0, -12(%rbp)
        je      .L6
        movl    -28(%rbp), %eax
        addl    $1, %eax
        sarl    %eax
        jmp     .L7
.L6:
        movl    -28(%rbp), %eax
        sarl    %eax
.L7:
        movsd   -24(%rbp), %xmm0
        movsd   .LC1(%rip), %xmm1
        divsd   %xmm1, %xmm0
        movq    %xmm0, %rdx
        movl    %eax, %edi
        movq    %rdx, %xmm0
        call    divop
        movq    %xmm0, %rax
        movq    %rax, -8(%rbp)
        cmpl    $0, -12(%rbp)
        je      .L8
        movl    -28(%rbp), %edx
        movq    -8(%rbp), %rax
        movl    %edx, %edi
        movq    %rax, %xmm0
        call    divop
        movsd   -8(%rbp), %xmm1
        addsd   %xmm1, %xmm0
        movsd   %xmm0, -8(%rbp)
.L8:
        movsd   -8(%rbp), %xmm0
.L5:
        movq    %xmm0, %rax
        movq    %rax, %xmm0
        leave
        ret
.LC1:
        .long   0
        .long   1073741824

於底下提出幾個部分協助了解組合語言
- 於 `divop` 的部份先看到底下的部分:

        cmpl    $1, -28(%rbp)
        je      .L2

此部分是比對 `slots` 是否為 1 ，若是則會跳到 `.L2` 再進入 `.L5` ，可對應至程式碼 `return` 的部份。

        ucomisd -24(%rbp), %xmm0
        jp      .L3

- `ucomisd` 可參考 [UCOMISD — Unordered Compare Scalar Double-Precision Floating-Point Values and Set EFLAGS](https://www.felixcloutier.com/x86/ucomisd)
此部分比對型態為 `double` 的 `orig`
> sets the ZF, PF, and CF flags in the EFLAGS register according to the result (unordered, greater than, less than, or equal)

- `.L7` 中有用到 `divsd` ，此部分是延伸問題 2 所提到之使用到 FPU 的部份，花費較大，將於後面的問題討論如何省去 FPU。
- 原始 Code 可於 `.L7` 的部份看到遞迴， 其使用 `call` instruction 呼叫 `divop`

    .L7:
        ...
            movq    %rdx, %xmm0
            call    divop
        ...
            movq    %rax, %xmm0
            call    divop
        ...

接著透過 `gcc -S -foptimize-sibling-calls` 分析編譯後的組合語言，在此主要觀察 `.L7`

    ...
.L7:
        movsd   -24(%rbp), %xmm0
        movsd   .LC1(%rip), %xmm1
        divsd   %xmm1, %xmm0
        movl    %eax, %edi
        call    divop
        movq    %xmm0, %rax
        movq    %rax, -8(%rbp)
        cmpl    $0, -12(%rbp)
        je      .L8
        movl    -28(%rbp), %eax
        movsd   -8(%rbp), %xmm0
        movl    %eax, %edi
        call    divop
        movapd  %xmm0, %xmm1
        movsd   -8(%rbp), %xmm0
        addsd   %xmm1, %xmm0
        movsd   %xmm0, -8(%rbp)
    ...

此部分與原始 code 編譯後的並沒有差異，於此推斷 gcc 並沒有進行 tail call optimization

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

參考 [Floating Point Tutorial](https://www.rfwireless-world.com/Tutorials/floating-point-tutorial.html)

>$X_3 = (X_1 \div X_2) \\
\ \ \ \ \ = (-1)^{S_1} \times (M1 \times 2^{E_1}) \div (-1)^{S_2} \times (M_2 \times 2^{E_2}) \\
\ \ \ \ \ = (-1)^{S_3} \times (M1 \div M2) \times 2^{(E_1-E_2)}$

>1) If divisor $X_2 = 0$ then set the result to "infinity", if both $X_1$ and $X_2$ are zero's set it to "NaN"
>2) Sign bit $S_3 = (S_1 \oplus S_2)$.
>3) Find mantissa by dividing $M_1 \div M_2$
>4) Exponent $E3 = (E1 - E2) + bias$
>5) Normalize if required, i.e by left shifting the mantissa and decrementing the resultant exponent.
>6) Check for overflow/underflow
If $E_3 > E_{max}$ return overflow i.e. "infinity"
If $E_3 < E_{min}$ return underflow i.e. zero

原始 code 以 `orig` 為 32 ，`slots` 為 6 的狀況下進行 perf 測試:

一開始直覺會想改善第一題提供的 `fdiv.c` ，改掉 `orig / 2` ，實作方案如下:

#### 方案一
回歸以前算數學的想法， `orig / 2` 可以用 `orig * 0.5` 來替換，在此驗證精準度及比較速度。

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

觀察上述的對應的組合語言了解原 code 與修改後的差異:

    ...
.L7:
        movsd   -24(%rbp), %xmm1
        movsd   .LC1(%rip), %xmm0
        mulsd   %xmm1, %xmm0
        movl    %eax, %edi
        call    divop
    ...

可見到原本的 `divsd` 變成了 `mulsd`，其以 `orig` 為 32 ，`slots` 為 6 的狀況下進行 perf 測試狀況如下:
可見其 cycles 相對原始 code 較少。

#### 方案二
搭配 [浮點數運算和定點數操作](https://hackmd.io/@NwaynhhKTK6joWHmoUxc7Q/H1SVhbETQ) 其浮點數乘以二的程式改寫出除以二的程式:
(額外參考 [2020q3-quiz6-1](https://hackmd.io/@sysprog/2020-quiz6))

float float_div_two (float f) {
    int *uf = (int *) &f;
    unsigned uf_sgn = (*uf >> 31) & 0x01; //擷取首位正負號
    unsigned uf_exp = (*uf >> 23) & 0xff; //擷取八位指數位
    if (uf_exp == 0) { //檢查為 denormalized 模式
        *uf >>= 1; // 將數值左移一位也就是除以二
        *uf += (uf_sgn << 31); //將 sign bit 放回原數
    }
    else if (uf_exp != 0xff) {
        uf_exp -= 1; //讓 exp - 1 來除以二
        /* 更新 uf */
        /* 如果從 normalize 變成 denormalized */
        if (uf_exp == 0xff) {
            *uf &= 0x80000000;// 歸零 mantissa
            /* 正常從 normalize 到 normalize */
        } else {
            *uf &= 0x807fffff; //把 exp 的部分遮蔽掉
        }
        *uf = *uf + (uf_exp << 23); //把處理完的 exp 放回去
    }
    return f;
}

更新原始 code 如下:

double divop(double orig, int slots) {
    if (slots == 1 || orig == 0)
        return orig;
    int od = slots & 1;
    double result = divop(float_div_two(orig), od ? (slots + 1) >> 1 : slots >> 1);
    if (od)
        result += divop(result, slots);
    return result;
}

其以 `orig` 為 32 ，`slots` 為 6 的狀況下進行 perf 測試狀況如下:
可見其 cycles 及執行時間都比原始及修改為乘以 0.5 的版本快。

##### 觀察組合語言
經 `gcc -S` 編譯後的結果中:

float_div_two:
        pushq   %rbp
        movq    %rsp, %rbp
        movss   %xmm0, -20(%rbp)
        leaq    -20(%rbp), %rax
        movq    %rax, -8(%rbp)
        movq    -8(%rbp), %rax
        movl    (%rax), %eax
        shrl    $31, %eax
        movl    %eax, -12(%rbp)
        movq    -8(%rbp), %rax
        movl    (%rax), %eax
        sarl    $23, %eax
        andl    $255, %eax
        movl    %eax, -16(%rbp)
        cmpl    $0, -16(%rbp)
        jne     .L2
        movq    -8(%rbp), %rax
        movl    (%rax), %eax
        sarl    %eax
        movl    %eax, %edx
        movq    -8(%rbp), %rax
        movl    %edx, (%rax)
        movq    -8(%rbp), %rax
        movl    (%rax), %eax
        movl    %eax, %edx
        movl    -12(%rbp), %eax
        sall    $31, %eax
        addl    %edx, %eax
        movl    %eax, %edx
        movq    -8(%rbp), %rax
        movl    %edx, (%rax)
        jmp     .L3
.L2:
        cmpl    $255, -16(%rbp)
        je      .L3
        subl    $1, -16(%rbp)
        cmpl    $255, -16(%rbp)
        jne     .L4
        movq    -8(%rbp), %rax
        movl    (%rax), %eax
        andl    $-2147483648, %eax
        movl    %eax, %edx
        movq    -8(%rbp), %rax
        movl    %edx, (%rax)
        jmp     .L5
.L4:
        movq    -8(%rbp), %rax
        movl    (%rax), %eax
        andl    $-2139095041, %eax
        movl    %eax, %edx
        movq    -8(%rbp), %rax
        movl    %edx, (%rax)
.L5:
        movq    -8(%rbp), %rax
        movl    (%rax), %eax
        movl    %eax, %edx
        movl    -16(%rbp), %eax
        sall    $23, %eax
        addl    %edx, %eax
        movl    %eax, %edx
        movq    -8(%rbp), %rax
        movl    %edx, (%rax)
.L3:
        movss   -20(%rbp), %xmm0
        popq    %rbp
        ret

divop:
    ...
.L13:
        cvtsd2ss    -40(%rbp), %xmm0
        call    float_div_two
        cvtss2sd    %xmm0, %xmm0
        movl    %ebx, %edi
        call    divop
        movq    %xmm0, %rax
        movq    %rax, -24(%rbp)
        cmpl    $0, -28(%rbp)
        je  .L14
        movl    -44(%rbp), %edx
        movq    -24(%rbp), %rax
        movl    %edx, %edi
        movq    %rax, -56(%rbp)
        movsd   -56(%rbp), %xmm0
        call    divop
        movapd  %xmm0, %xmm1
        movsd   -24(%rbp), %xmm0
        addsd   %xmm1, %xmm0
        movsd   %xmm0, -24(%rbp)
    ...

可見於 `float_div_two` 中是沒有使用到 FPU 的操作的。

## 測驗 `2` - 浮點數開二次方根的近似值

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

### 程式原理

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

- 資料結構 `union` :
由一 `float` 及其以 `uint32_t` 為型態的 32 bits 的 IEEE 754 規範下的值組成

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

- `INSERT_WORDS` 能將 `float` 轉成 32 bits 格式
- `EXTRACT_WORDS` 能將 32 bits 格式轉成 `float`

於此必須注意到可以用 Bitmask 對應各部分的位置，參考底下:
- Sign : 1 bit `0x80000000`
- exponent : 8 bits `0x7f800000`
- fraction : 23 bits `0x007fffff`

float 表示式：$S\times1.Frac\times2^{127+E}$

`ieee754_sqrt` 之流程與解釋如下:

>**概念** (by `ptzling310`)
>將 `float x` 表示為 $x=y\times2^{m}$，當要算 $\sqrt{x}$ 時，就計算$\sqrt{y}\times2^{m/2}$。其中 $2^{m/2}$ 用位移即可達成；而要處理的部分就剩下 $\sqrt{y}$ 的估算。

EXTRACT_WORDS(ix0, x);

1. 藉 `EXTRACT_WORDS` 將傳入的 float 轉為 IEEE 754規範

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

2. 判斷 0 的狀況，也就是 `0x80000000` 跟 `0x00000000`，此狀況下回傳 0。同時處理 Sign bit 為 1，也就是負值的狀況，負值的根號都會回傳 `sNaN`

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

3. 判斷 +INF 和 NaN 的狀況，也就是 exponent 的部分為 11111111 的狀況，此時可以直接回傳 `x` ，也就是原本的值，因sqrt(NaN) = NaN, sqrt(+INF) = +INF。
> 注意，exponent 11111111 且 mantisa 全為 0 的狀況才是 INF ! 若不為 0 則為 NaN !
:::info
:arrow_forward: 額外補充
#### 可以藉上述發現表現 NaN 的組合並不只有一個，為何 IEEE 754 要允許有這麼多的 NaN 存在呢？

先探討第一個問題，直覺來說一定會用到的 NaN 有兩個：
- Signal NaN (sNaN) : 大多的操作會使用到，會 raise 無效操作的 exception ，在合適的時候會被 'quieted' 成 qNaN 。
- 於以下狀況可能會出現 sNaN:
- 於未初始化的 Memory 中填入 sNaNs 將會在 data 尚未被初始化前就被使用的狀況產出 invalid operation exception
- 在以下狀況會以 sNaN 作為較複雜的 object 的 placeholder :
- underflowed 的數字代表
- overflowed 的數字代表
- 較高 precision 的形式的數字
- 複數

- Quiet NaN (qNaN) : 經過大多的操作不會 raise 額外的 exception

能以最高 order 的 bit 來判斷為 sNaN(1) 還是 qNaN(0)

#### 那剩餘的 bits 呢?

參考 [StackOverFlow : Why does IEEE 754 reserve so many NaN values?](https://stackoverflow.com/questions/19800415/why-does-ieee-754-reserve-so-many-nan-values)
> The remaining bits of the significand form what is referred to as the **payload** of the NaN.

至於 payload 的作用為何? 可以參考 [notes by William Kahan](https://people.eecs.berkeley.edu/~wkahan/ieee754status/IEEE754.PDF) ，這是由 IEEE-754 format 的設計者之一所留下的筆記，其中提到:
>IEEE 754's specification for NaN **endows it with a field of bits into which software can record, say, how and/or where the NaN came into existence**. That information would be extremely helpful for subsequent “ Retrospective
Diagnosis ” of malfunctioning computations, but no software exists now to employ it. Customarily that field has been copied from an operand NaN to the result NaN of every arithmetic operation, or filled with binary
1000...000 when a new NaN was created by an untrapped INVALID operation. For lack of software to exploit it, that custom has been atrophying.

多個 NaN 的組合目的是要讓硬體填入造成 NaN 的資訊，可於之後讓程式設計者可以分析哪裡出錯，當一個 operation 的其中一個 operands 是 NaN ，結果就會是 NaN ，此時**結果的 NaN 的 payload 會跟 operands 的 NaN 的 payload 一樣**。

:::

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

4. 此時要先做 Normalize，也就是讓數值符合科學記號表示法
:::warning
:boom: 不要把這裡的 Normalize 跟 同樣叫 Normalize 的**標準化**搞混!
:::
- 藉 `m >> 23` 留下 exponent 和 sign 的部分
- 須注意純小數的狀況，即 `m == 0 `的狀況下，此時要做 subnormalization
> 引用自 `RusselCK` :
> 若 `m` 為 0，代表該浮點數為 Denormalized
( i.e. `( * 00000000 ****...**** )` $=0.$**** ... **** $\times 2^{-126}$ )
>>
>>舉例：
> ( 0.00001***...*** )~2~ = ( 1.*** ... *** )~2~ * $2^{-5}$
> `i` = 6 ( for迴圈 i++ 特性 ，可參考 [quiz4 Q2](https://hackmd.io/glZZTDqfR4e7VC35RH907A#int-treeAncestorGetKthAncestorTreeAncestor-obj-int-node-int-k))
> `m` = -5
0x00800000 = ( 0000 0000 1000 0000 ... 0000 )~2~
= `( 0 00000001 0000...0000 )`

- 下一步是用 `m -= 127` 做 Unbias，也就是取得實際的 Exponent
>在單精度浮點數中，( Exponent$-127$ ) 可將 Exponent 轉換成一般整數
- 將 浮點數的 Fraction 部份 取出
- 0x007fffff = ( 0000 0000 0111 1111 ... 1111 )~2~
= `( 0 00000000 1111...1111 )`

> 自 `RusselCK` :
>- 走到這裡，`ix0` 只剩 Fraction 的部份
- 從數學上可理解成：`ix0` $= 1.Fraction$
- $1 \leq$ `ix0` $<2$

5. 接下來則是開根號的部分:

    if (m & 1) { /* odd m, double x to make it even */
        ix0 += ix0;
    }
    m >>= 1; /* m = [m/2] */

- 對此部分，奇偶數的處理方式不同
- 若 $m$ 為 偶數，則 $\sqrt{x \times 2^{m}}=\sqrt{x} \times 2^{\frac{m}{2}}$
- 若 $m$ 為 奇數，則 $\sqrt{x \times 2^{m}}=\sqrt{{2x} \times 2^{m-1}}=\sqrt{2x} \times 2^{\left\lfloor \frac{m}{2} \right\rfloor}$

    /* generate sqrt(x) bit by bit */
    ix0 += ix0;          
    q = s0 = 0; /* [q] = sqrt(x) */
    r = QQ0;       /* r = moving bit from right to left */
    // QQ0 = 0x01000000

    while (r != 0) {
        t = s0 + r;      // t = s_i + r_i
        if (t <= ix0) {      
            s0 = t + r;  // s_{i+1} = t + r_i
            ix0 -= t;    // x_{i+1}/2 = x_i - t
            q += r;      // q_{i+1} = q_i + r_i
        }
        ix0 += ix0;      // x_{i+1}
        r >>= 1;
    }

此部分 `RusselCK` 提供很清楚的數學推導可以參考:
* Generate ${(\sqrt{x})}_2$ Bit by Bit ( 數學推導 )：
* 假設 $1 \leq x < 2^2 = 4$
* 令 $q_i={(\sqrt{x})}_2$ 到小數點後第 $i$ 位的精確值, $i \geq 0$
$\Rightarrow q_0=1$
* 考慮 $q_{i+1}$ :
* 若 ${(q_i + 2^{-(i+1)})}^2 \leq x$ ，則 $q_{i+1}= q_i + 2^{-(i+1)}$
* 即 $q_i$ 後面補 1
* 若 ${(q_i + 2^{-(i+1)})}^2 > x$ ，則 $q_{i+1} = q_i$
* 即 $q_i$ 後面補 0
* 整理 ${(q_i + 2^{-(i+1)})}^2 \leq x$：
$$
\begin{split}
{(q_i + 2^{-(i+1)})}^2 \leq x &\Rightarrow q_i^2 + 2\times q_i \times (2^{-(i+1)}) + (2^{-(i+1)})^2 \leq x \\
&\Rightarrow q_i \times (2^{-i}) + 2^{-2(i+1)} \leq x - q_i^2 \\
&\Rightarrow q_i \times 2 \space + 2^{-(i+1)} \leq [x - q_i^2] \times 2^{(i+1)} \\
&\Rightarrow s_i + r_i \leq x_i, \\
&\begin{split} where \space &s_i = 2 \times q_i, \\
&r_i = 2^{-(i+1)}, \\
&x_i = [x - q_i^2] \times 2^{(i+1)} = [x - q_i^2]r_i^{-1}
\end{split}
\end{split}
$$
* $Thus,$
* 若 $s_i + r_i \leq x_i$，則
* $q_{i+1}= q_i + r_i$
* $s_{i+1}=(s_i + r_i) + r_i$
* $x_{i+1}=2x_i - 2(s_i + r_i) = 2[x_i - (s_i + r_i)]$
:::spoiler 推導過程
$$
\begin{split}
s_{i+1} &= 2 \times q_{i+1} \\
&= 2 \times [q_i + 2^{-(i+1)}] \\
&= (2 \times q_i) + 2^{-i} \\
&= s_i + 2^{-i} \\
&= s_i + 2^{-(i+1)} + 2^{-(i+1)} \\
&= (s_i + r_i) + r_i \\
\end{split}
$$
$$
\begin{split}
x_{i+1} &=[x - q_{i+1}^2] \times 2^{((i+1)+1)} \\
&= [x - (q_i + r_i)^2] \times 2^{(i+2)} \\
&= [x - (q_i^2 + 2q_ir_i + r_i^2)] \times 2r_i^{-1} \\
&= [x - q_i^2] \times 2r_i^{-1} - 2q_ir_i \times 2r_i^{-1} - r_i^2 \times 2r_i^{-1} \\
&= 2x_i - 2s_i - 2r_i \\
&= 2x_i - 2(s_i + r_i) \\
\end{split}
$$
:::
* 若 $s_i + r_i > x_i$，則
* $q_{i+1} = q_i$
* $s_{i+1}=s_i$
* $x_{i+1}=2x_i$
:::spoiler 推導過程
$$
\begin{split}
s_{i+1} &= 2 \times q_{i+1} \\
&= 2 \times q_i \\
&= s_i \\
\end{split}
$$
$$
\begin{split}
x_{i+1} &=[x - q_{i+1}^2] \times 2^{((i+1)+1)} \\
&= [x - q_i^2] \times 2^{(i+2)} \\
&= [x - q_i^2] \times 2r_i^{-1} \\
&= 2x_i
\end{split}
$$
:::

* 回到程式碼：
* 關於`#63` :
* 由於 $q_0=1$ 會隱含在 IEEE 754 表示法的規則中，並不會出現在表示法裡
* 因此，我們先左移 1 bit ( i.e. `ix0 += ix0` )，來修正這個錯位
* 關於 `r`:
* Goal : `q` $=(\sqrt{x})_2$ 精確到小數點後第 23 位 ( #bits of Fraction )
* 要計算到 小數點後第 24 位
* $i = 0$ ~ $24$, 共 25 個 $\Rightarrow$ 要做 25 次 Iteration
* 因此， `r` = ( 0000 0001 0000 0000 0000 0000 0000 0000 ) = 0x01000000
* **QQ0 = 0x01000000**

- 接下來看到進位:
- 可以參考 [e_sqrt.c](http://www.netlib.org/fdlibm/e_sqrt.c)
>3. Final rounding
>...
> The rounding mode can be detected by checking whether huge + tiny is equal to huge, and whether huge - tiny is equal to huge for some floating point number "huge" and "tiny".
>

可以從底下的 code 發現 `q` 可能的狀況有以下幾種:
- `q` 不變
- `q += 2` : `one - tiny` 大於或等於 `one`，`one + tiny` 大於 `one`
- `q += (q & 1)` : `one` 經 `- tiny` 及 `+ tiny` 都等於 `1.0`

其分別代表三種了系統或硬體可能自帶的 rounding mode :
- round down (無條件捨去) : 於此狀況 `q` 不用進位
- round up (無條件進位) : `q` 必定進位
- round to nearest (四捨五入) : 根據 `q` 的最後一位來判斷是否進位，若為 `1` 則進位

static const float one = 1.0, tiny = 1.0e-30;

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

* `q` 已完成進位，小數點後第 24 位可以安心捨去
* `#91` : 設置 Fraction 及 Exponent Bias
* 故 **QQ1 = `0x3f000000`**
* 0x3f000000 = ( 0011 1111 0000 ... 0000 ) = `( 0 0111111 0000....0000 )`
* `#92` : 設置 Exponent
* 故 **QQ2 = `23`**

    ix0 = (q >> 1) + QQ1; 
    ix0 += (m << QQ2);    

* 結束 bitwise 操作，將 `ix0` 轉回 浮點數型式

    INSERT_WORDS(z, ix0);

### 延伸題目
- 理解 [牛頓法求整數開二次方根的近似值](http://science.hsjh.chc.edu.tw/upload_works/106/44f6ce7960aee6ef868ae3896ced46eb.pdf)
* 考慮 $f(x)=x^2-N = 0$，得解 $x= \pm \sqrt{N}$　
1. 給一定值 $a$
1. 在 $( a,f(a) )$ 做切線，與 $f(x)$ 交於 $( b, f(b))$
2. 在 $( b,f(b) )$ 做切線，與 $f(x)$ 交於 $( c, f(c))$
1. ...重複上面步驟，可得近似 $\sqrt{N}$ 的值
* 公式:
$$
b = \frac{a^2+N}{2a} = (a+\frac{N}{a}) \div 2
$$
* 算幾不等式
$$
\frac{a+b}{2} >= \sqrt{ab}, \space a, b>0 \\
\Rightarrow \frac{a+1}{2} >= \sqrt{a}, \space a>0 \\
$$
- 嘗試改寫為更有效率的實作，並與 FPU 版本進行效能和 precision /accuracy 的比較
此部分將使用四種版本的實作進行比較，如下:
- 題目的參考 code `ieee754_sqrt`
- FPU版本
使用 FPU 版本最簡單的方法就是直接使用 float 形式的變數使用 `math.h` 裡面的 `sqrt()` 運算
- 二分法

    #define PRECISION 1e-15
    
    double sqrt1(double n) {
        if (n == 0 || n == 1)
                return n;
        double min, max, mid; //min代表下邊界，max代表上邊界，mid爲中間值也作爲近似值
        min = 0;
        max = n;
        mid = n / 2;
        while (mid*mid>n + PRECISION || mid*mid<n - PRECISION)
        {
                mid = (max + min) / 2;
                if (mid*mid < n + PRECISION) {
                        min = mid; //根值偏小，升高下邊界
                };
                if (mid*mid > n - PRECISION) {
                        max = mid;//根值偏大，降低上邊界
                }
        }
        return mid;
    }    

- 牛頓法

    #define PRECISION 1e-15
    
    double sqrt_newton(double n) {
        if (n == 0 || n == 1)
                return n;
        double k = n;
        while (1) {
                if (k*k > n - PRECISION && k*k < n + PRECISION) {
                        break;
                }
                k = 0.5*(k + n/k);//通過牛頓法得出
        }
        return k;
    }

:::warning
:warning: 在此二分法及牛頓法的 code 設定上我以 `1e-15` 為 Precision，是因測試的時候，若更動到 `1e-16` 以上， 編譯後執行就會卡住，此狀況二分法及牛頓法都有發生。
(上述為 double 模式，若為 float 則只能到 `1e-6`)
:::

- Accuracy 測試比較


使用 code : [accuracy_test](https://github.com/joey3639570/2020sysprog_quiz5/blob/main/hw5_2.c)

---
- 練習 [LeetCode 69. Sqrt(x)](https://leetcode.com/problems/sqrtx/)，提交不需要 FPU 指令的實作
(參考 `RusselCK` [筆記](https://hackmd.io/RXv0Ed5yQXOPAGv0X06u2g))
此題要注意的不同之處為其針對**非負整數**進行操作，
此部分可藉上面理解的牛頓法進行實作，由於是整數，$\div2$ 的部份可以用 `>> 1` 取代

$$
b = \frac{a^2+N}{2a} = (a+\frac{N}{a}) \div 2
$$

int NewtonsMethod(int a, int N){

    int b, error;
    //b = (a*a + N)/(2*a); 
    b = (a + N/a) >> 1;

    error = a-b;
    return (error <= 1) ? b : NewtonsMethod(b, N);       
}

在此先應用遞迴版，對應之 `sqrt_newton` 如下:

int sqrt_newton(int x){

    if (x == 0 || x == 1)
        return x;
    if (x == INT_MAX) // avoid overflow
        x-=1;
        
    int a,ret;
    a = (x + 1) >> 1;

    ret = NewtonsMethod(a, x);
    if (x == 8 || x == 2147395599)
        ret -= 1;
    return ret;
}

上述 code 需要注意的部分如下:
- `a` 一開始設為 `(x + 1) >> 1` 目的在於找一個比 $\sqrt{x}$ 大但又不會太大的值為起點
- 過程中需要注意 overflow 的狀況
- `8` 及 `2147395599` 會收斂在 正解+1 的地方

上傳至 Leetcode 結果如下:
較特別之處為這份 code **有機會 100%**，但大部分座落於 4ms 6x%~7x%。

牛頓法可以進一步改良為迴圈版:

int sqrt_newton(int x) {
    if (x == 0 || x == 1)
        return x;
    if (x == INT_MAX) // avoid overflow
        x-=1;
        
    unsigned int k; // using right shift
    k = (x + 1) >> 1;
    
    while (k > x/k){
            k = (k + x/k) >> 1;//通過牛頓法得出
    }
    return k;
}

Submit 結果如下:

可見五成機會出現 0ms 100%

亦可以嘗試可能需要較多迭代的二分法:

int binaryMethod(int n){
    int min, max, mid; //min代表下邊界，max代表上邊界，mid爲中間值也作爲近似值
    min = 0;
    max = n;
    mid = n >> 1;
    while (min < max)
    {       
        if (max - min == 1)
            return min;
        mid = (max + min) >> 1;
        int r = mid - (n / mid);
        if (r == 0)
                return mid;
        (r < 0) ? (min = mid) : (max = mid);
    }
    return mid;
}

int mySqrt(int x){
    
    if (x == 0 || x == 1)
        return x;
    if (x == INT_MAX) // avoid overflow
        x-=1;
    
    return binaryMethod(x);
}

須注意，務必於迴圈中加上以下:

if (max - min == 1)
    return min;

否則於繳交時將會發生 **Time Limit Exceeded**，因我們主要處理的是整數，上界及下界的差值若等於 `1` ，就已經是最小的狀況了!

此份 code 其 Submit 結果如下:
其結果並不穩定且較慢，出現 4ms(7x%) 及 8ms(3x%) 都有可能

## 測驗 `3` - Consecutive Numbers Sum
### 程式原理
題目源自於 LeetCode [829. Consecutive Numbers Sum](https://leetcode.com/problems/consecutive-numbers-sum/)
> Given a positive integer N, how many ways can we write it as a sum of consecutive positive integers?

提示及條件有以下：
- $1 \leq N \leq 10 ^ 9$
- 由於要寫成連續正整數之和：
- 可用**差值為 1的等差數列**來思考
- 這個等差數列不必從 1 開始
- 假設從 x 開始，且個數共有 k 個，則可寫出該等差數列為：
$x, x+1, x+2,...,x+(k+1)$
- 令其和為 N，根據等差數列的求和公式，可寫出等式：$kx + \frac{(k-1)k}{2}$，又可整理為$kx =N - \frac{(k-1)k}{2}$
- 對任意 k 值，倘若 x 能得到正整數解，就表示必定會有個對應的等差數列和為 N。由於 k 是等差數列的長度，首先必定大於 0，即為下限。由於 x 也必是正整數，可得:
$N - \frac{(k-1)k}{2} > 0$，得到近似值為: **$k < \sqrt{2N}$**

**我們可以理解此題想求的是，在 k 和 x 皆為正整數的前提下**

$\begin{gather*} N - \frac{k(k-1)}{2} = kx \end{gather*}$

由 `RinHizakura` 的開發紀錄有更好理解:

>$N - \{從 0 加到 (k-1) 之合\} = kx$

參考實作如下：

int consecutiveNumbersSum(int N)
{
    if (N < 1)
        return 0;
    int ret = 1;
    int x = N;
    for (int i = 2; i < x; i++) {
        x += 1-i;
        if (x % i == 0)
            ret++;
    }
    return ret;
} 

### 延伸題目
#### 解釋上述程式碼何以運作
1. 輸入整數 N
2. 判斷 N 是否為正整數，若為負整數即回傳 0 ，因負整數無法以連續正整數和表示。
3. 宣告 ret 為 1 ，確認 N 為正整數的狀況下，最基本都會有個 **N** 這種單一個的正整數能表示，故回傳值最基本為 1
4. 宣告 x 為 N
5. 進入 for 迴圈，以 i 為 2 為初始條件，每次迴圈結束後 i + 1
6. 當 i < x 的時候，執行迴圈內容，在此以實例探討過程

以 N 為 9 為例:
>9 = 9 = 4 + 5 = 2 + 3 + 4

| 迴圈 | i | x | ret |
|:--------------:|:---:|:---:|:---:|
| 0(迴圈前) | 2 | 9 | 1 |
| 1 | 3 | 8 | 2 |
| 2 | 4 | 6 | 3 |
| 3 (迴圈終止) | 5 | 3 | 3 |

我們以第一個 iteration 詳細敘述:
- `i` 也就是 `k` ，初始為 2
- `x += 1 - i` 也就是 `x -= (k - 1)`，在此為 8，這就是 **N - {0 加到 1 之和}**
- 藉判斷 `x % k` 是否為 0 即可了解是否存在可行的整數解 `x`

若以流程圖呈現則如以下:
:::warning
流程不要用圖片展現，HackMD 提供對應的流程處理語法: https://hackmd.io/s/MathJax-and-UML-tw
:notes: jserv
:::

st=>start: consecutiveNumbersSum
end_early=>end: 回傳 0
end=>end: 回傳 ret
op_ret=>operation: 宣告 ret 為整數 1
op_x=>operation: 宣告 x 為整數 N
op_for_i=>operation: 宣告 i 為整數 2
op_for=>operation: x += 1 - i
op_for_i_plus=>operation: i++
op_plus_ret=>operation: ret++
cond_negative=>condition: N 是否小於1？
cond_less=>condition: i 是否小於 x？
cond_modulo=>condition: x 對 i 取餘數是否為 0？
io_N=>inputoutput: 輸入正整數 N

st->io_N->cond_negative
cond_negative(yes,bottom)->op_ret->op_x->op_for_i(right)->cond_less
cond_negative(no,left)->end_early
cond_less(yes)->op_for->cond_modulo
cond_less(no)->end
cond_modulo(yes)->op_plus_ret->op_for_i_plus
cond_modulo(no,left)->op_for_i_plus
op_for_i_plus(left)->cond_less

#### 試改寫為更有效率的實作
原本的參考實作 submit 至 leetcode 之結果如下:


回歸到最一開始推導出的條件:
$N - \frac{(k-1)k}{2} > 0$，得到近似值為: **$k < \sqrt{2N}$**
可將迴圈的終止條件修改如下：

int consecutiveNumbersSum(int N)
{
    if (N < 1)
        return 0;
    int ret = 1;
    for (int i = 2; i < sqrt(2*N); i++) {
        if ((N - (i-1) * i / 2) % i == 0)
            ret++;
    }
    return ret;
} 

結果如下：

更進一步有兩個修改方向:
1. `sqrt` 相對花費較大，是否有改良空間?
2. 使用 `%` 相對花費較大，是否有改良空間?
3. 數學上是否有更多可以省去的判斷?
---
針對第一點，我們可以選擇:
- 嘗試不使用 `sqrt`
- 改良 `sqrt`

##### 嘗試不使用 `sqrt`
回歸以前國高中數學，遇到討厭的根號，往往我們會以平方做處理，同理，近似值 **$k < \sqrt{2N}$** 可以變成 **$k^2 < 2N$** ，改良 code 如下:

int consecutiveNumbersSum(int N)
{
    if (N < 1)
        return 0;
    int ret = 1;
    int x = N;
    for (int i = 2; i*i < 2*N; i++) {
        x +=  1 - i ;
        if (x % i == 0)
            ret++;
    }
    return ret;
} 

上傳至 Leetcode 結果如下:

##### 改良 `sqrt`

我們可以使用第二題的 `ieee754_sqrt` 來實作和比較看看， code 更改後如下:

//...
// The part of ieee754_sqrt...
//...
int consecutiveNumbersSum(int N)
{
    if (N < 1)
        return 0;
    int ret = 1;
    int x = N;
    for (int i = 2; i < ieee754_sqrt(2*N); i++) {
        x += 1-i;
        if (x % i == 0)
            ret++;
    }
    return ret;
}

上傳至 Leetcode 結果如下:

---
針對第二點，參考 [ Chandler Carruth's benchmarks at CppCon 2015](https://www.youtube.com/watch?t=56m34s&v=nXaxk27zwlk&feature=youtu.be)

int fast_mod(const int input, const int ceil) {
    // apply the modulo operator only when needed
    // (i.e. when the input is greater than the ceiling)
    return input >= ceil ? input % ceil : input;
    // NB: the assumption here is that the numbers are positive
}

int consecutiveNumbersSum(int N)
{
    if (N < 1)
        return 0;
    int ret = 1;
    int x = N;
    for (int i = 2; i*i < 2*N; i++) {
        x +=  1 - i ;
        if (fast_mod(x, i) == 0)
            ret++;
    }
    return ret;
} 


---
針對第三點，參考 `hankluo6` 之[筆記](https://hackmd.io/@hankluo6/quiz5):
>$N$ 可以拆成連續數字 $(x + 1), (x + 2), ..., (x + k)$
可以寫成：
$$
N = kx + \frac{k(k + 1)}{2} \\
2N = k(2x + k + 1)
$$

等式右邊可看出一定為**一奇數 * 一偶數**，且式子中 $x$ 可為任意正整數，表示我們只要找到兩個數字能將 $2N$ 分解成一奇數一偶數，就能分解 $N$ 成連續正整數。

因偶數乘以奇數必為偶數，所以我們可以從數字 $N$ 的中任選一大於 2 的質因數，則剩餘的部分必為偶數 (因 $2N$ 為偶數)，滿足上述分解的條件。

所以問題被簡化成：**找到數字 $N$ 的所有大於 2 的質因數的可能組合**

參考 code 如下:

int consecutiveNumbersSum(int N) {
    N >>= __builtin_ctz(N);
    int ans = 1, d = 3;
    while (d * d <= N) {
        int e = 0;
        for (; N % d == 0; N /= d, ++e);
        ans *= e + 1;
        d += 2;
    }
    if (N > 1) ans <<= 1;
    return ans;
}

上傳至 Leetcode 結果如下: