# 2020q3 Homework5 (quiz5)
contributed by < `Rsysz` >.
###### tags: `sysprog`

>[測驗題](https://hackmd.io/@sysprog/2020-quiz5)

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

將浮點數 `orig` 除以大於 `0` 的整數 `slots`
若 `slots` 為偶數
$$
\frac{orig/2}{slots/2}=\frac{orig}{slots}
$$

若 `slots`為奇數，透過 `result += divop(result, slots)` 補償
$$
\frac{orig/2}{(slots+1)/2}+\frac{\frac{orig/2}{(slots+1)/2}}{slots}=\frac{orig/2}{(slots+1)/2}\times(1+\frac{1}{slots})=\frac{orig/2}{(slots+1)/2}\times\frac{(slots+1)}{slots}=\frac{orig}{slots}
$$

對奇數 `slots` 加 `1` 後右移與補償，巧妙地加速計算速度，D1 = `(c) 2`, D2 = `(d) 1`
## 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;
}

首先看到 `do while(0)` 封裝的 `marco`，這種用法是為了確保替換後程式正確編譯與減少內部 `Bug` 的發生，再來看到各 `marco` 的功能。
* ` INSERT_WORDS ` 將 `int ix0` 轉為 `float d`
* ` EXTRACT_WORDS ` 將 `float d` 轉為 `int ix0`
* 因無論 `unsigned` 與否，內部數值皆為相同的，因此皆以 `int` 表示

接著看到整體程式碼，首先將 `float x` 轉為 `int ix0`
* `ix0 <= 0` 也就是 `float x` 為負數時有兩個結果
* `(ix0 & (~sign)) == 0` 表示除 `sign` 以外皆是 `0`，代表 `float x`為 `+-0` 開根號還是 `0`，回傳 `x`
* 先前條件不滿足且 `ix0 < 0`，表示 `float x` 為負浮點數，無法開根號，因此回傳 `NaN`

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

接著處理 `INF` 與 `NaN`，先前條件不滿足且 `(ix0 & 0x7f800000) == 0x7f800000` 表示 `float x` 的 `Exponent` 部分皆為 `1`
* `+INF` 表示為 $1_{sign}11111111_{exponent}..._{fraction}$ `fraction` 皆為 `0`
* `NaN` 表示為 $0_{sign}11111111_{exponent}..._{fraction}$ `fraction` 部分皆可，此外先前 `(ix0 < 0)` 已排除 `sign` 為 `1` 的結果

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


* 對浮點數做正規化，首先取出 `ix0` 的前九個bit 也就是 `float x` 的 `sign` 與 `exponent` 部分
* 若 `m == 0` 則對 `ix0` 做左移，直到 `exponent` 部分的 $2^0$ 為 `1`
* 因 `m` 表示 `exponent` 部分，根據上圖對浮點數的正規化，`m = m + 1 - i`

* 接著將 `m - 127` 得到對應 `exponent` 的數字，並以 `ix0 = (ix0 & 0x007fffff) | 0x00800000` 擷取 `ix0` 的後 `23` bit並將第 `24` bit 設為 `1`，得到$1.fraction$
* 最終得到$2^{m} \times ix0，1 \leq ix0 < 2$ 接著為了將$\sqrt{2^{m}} = 2^{\frac{m}{2}}$，若 `m` 為奇數，將多的 `2` 乘至 `ix0` 部分，再將 `m` 除以 `2` 得到$2^{\frac{m}{2}} \times (2 \times 1.fraction)，1.fraction \leq ix0 \leq (2 \times 1.fraction)$
* 又 $1 \leq 1.fraction < 2，所以1 \leq ix0 < 4$，整理如下$$2^{\frac{m}{2}} \times ix0，1 \leq ix0 < 4，m \in even$$

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

* 接著因 $1 \leq ix0 < 4$ 表示 $1 \leq \sqrt{ix0} < 2$，使得可以從 `1` 開始逼近 $\sqrt{ix0}$ 的數值。
* 此外前面將 `fraction` 部分擴大至 `24` bit(新增 `1` 的部分)，這裡又將 `ix0 += ix0`，因此範圍變成$2 \leq ix0 < 8$，所以為了保證 `ix0` 最小值 `1.fraction` 能正確逼近，這裡要從 `25` 位開始做判斷， QQ0 = `(a) 0x01000000`

:::warning
想請問這裡一下，這裡 `ix0 += ix0` 與 `ix0 << 1` 在各種方面上有差異嗎?
:::

* 令 $q_i = \sqrt{ix0}$ 取小數點後 $i$ 位的精度$q_0 = 1$， $s_i = 2q_i$且$ix0_i = 2^{i+1} (ix0 - q_i^2)$
* 為了透過$q_i$計算$q_{i+1}$，若$(q_i+2^{-(i+1)})^2 \leq ix0$成立則$q_{i+1} = q{i} + 2^{-(i+1)}(補1) 否則 q_{i+1} = q_i(補0)$
* 接著透過代換得到$s_i+2^{-(i+1)} \leq ix0_i$
* 若成立，$s_{i+1} = s_{i} + 2^{-i}，ix0_{i+1} = ix0_i - s_i - 2^{-(i+1)}$
* 否則，$s_{i+1} = s_{i}，ix0_{i+1} = ix0_i$

* 這裡以圖解說明逼近過程，首先假設 `ix0 = 1.10...`

digraph FP{
	node[shape=record]
	
	BF16 [label="<f0>0|<f1>0|<f2>0|<f3>0|<f4>0|<f5>0|<f6>0|<f7>0|<f8>1|..."]
    BF16:f0->Sign
	BF16:f8->24
}

`ix0 += ix0`

digraph FP{
	node[shape=record]
	
	BF16 [label="<f0>0|<f1>0|<f2>0|<f3>0|<f4>0|<f5>0|<f6>0|<f7>1|<f8>1|..."]
    BF16:f0->Sign
	BF16:f8->24
}

`q = s0 = 0`
`r = 0x01000000`
`t = s0 + r` 也就是 $s_i+2^{-(i+1)} \leq ix0_i$ 首項 `s` 為 `0` 是因為還沒進入小數位
`t <= ix0` `t = 0x01000000` < `ix0 = 0x018XXXXX`
`s0 = t + r` 也就是 $s_{i+1} = s_{i} + 2^{-i}$ `s0 = 0x02000000`
`ix0 -= t` 也就是 $ix0_{i+1} = ix0_i - s_i - 2^{-(i+1)}$ `ix0 = 0x008XXXXX`

digraph FP{
	node[shape=record]
	
	BF16 [label="<f0>0|<f1>0|<f2>0|<f3>0|<f4>0|<f5>0|<f6>0|<f7>0|<f8>1|..."]
    BF16:f0->Sign
	BF16:f8->24
}

`q += r` `q = 0x01000000`
`ix0 += ix0` `ix0 = 0x01XXXXXX`

digraph FP{
	node[shape=record]
	
	BF16 [label="<f0>0|<f1>0|<f2>0|<f3>0|<f4>0|<f5>0|<f6>0|<f7>1|<f8>0|..."]
    BF16:f0->Sign
	BF16:f8->24
}

`r >>= 1` `r = 0x00800000`

* 透過 `ix0 += ix0` 對 `ix0` 左移，透過 `r >> 1` 對 `r` 右移，不斷逼近 $\sqrt{ix0}$ 並將逼近後每位數值存於 `q`

* 剩下的數值將會根據不同系統上的精度做 `rounding`，四捨五入

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

* 最後將計算結果還原成浮點數表示
* `q` 表示的是小數點落在 24~25 bit間的編碼，因此右移一次回到正常範圍
* QQ1 = `(d) 0x3f000000`，照理來講應該加上 `0x3f800000` ，但因 `q = 1.fraction` 其最高位的 `1` 右移時補上了缺口
* `m` 表示 `exponent` 部分，QQ2 = `(g) 23` 左移回正確位置

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

    INSERT_WORDS(z, ix0);
    return z;

### 延伸問題

[LeetCode 69. Sqrt(x)](https://leetcode.com/problems/sqrtx/)

typedef union {
    float value;
    uint32_t v_int;
} ieee_float_shape_type;

#define EXTRACT_WORDS(ix0, d)        \
    do {                             \
        ieee_float_shape_type ew_u;  \
        ew_u.value = (d);            \
        (ix0) = ew_u.v_int;          \
    } while (0)

int mySqrt(int x)
{
    int32_t sign = 0x80000000;
    int32_t ix0, m, 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 */
    }   
    /* 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] */
    
    unsigned int lo = 1 << m, hi = 1 << (m + 1);
    unsigned int mid = (lo + hi) >> 1;
    while(1){
        if(x > (mid + 1) * (mid + 1)){
            lo = mid;
            mid = (lo + hi) >> 1;
        }
        else if (x < mid * mid){
            hi = mid;
            mid = (lo + hi) >> 1;
        }
        else
            return mid;
    }
}

因本題無須深入求小數點位的精度要求，因此借用本題方式轉換為$2^{\frac{m}{2}} \times 1.fraction$接著利用$2^{\frac{m}{2}}$做 `Binary search` 以此快速縮小範圍，最終回傳 `mid`。

## 3.

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

給定正整數N，求N能寫成多少種連續正整數的和，假設從 `x` 開始，且個數共有 `k` 個
根據等差數列的求和公式
$$
kx = N - \frac{(k-1)k}{2}
$$
* 求有幾種組合`k` 與 `x` 能滿足 `N - {0 ~ k}的和 = kx`
* 首先看到 `x = N` 對應著公式中的 `N`，代表 `ZZZ` 對應著 `{0~k}的和`，因此 `i` 對應 `k`
* `for` 迴圈則是用來計算不同終點的等差數列，所以ZZZ = `(d) 1 - i`
* 舉例來說，若 `N = 3`
* 第一個迴圈`N - {0~1}的和 = 2` 也就是 `x += 1 - i`，因此檢查 `x % i` 也就是$x = \frac{N - \frac{(k-1)k}{2}}{k}$就能知道是否有可行的解
* 第二個迴圈`N - {0~2}的和 = 0` 因 `x` 於第一個迴圈已被更新，因此算出結果仍然等於 `x += 1 - i`，且 `x % i == 0` 因此得到 `x = 0` 的整數解一個，以此類推，直到 `k = x`。