---
# System prepended metadata

title: 2020q3 Homework5 (quiz5)

---

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

## 測驗 1
```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 / 2, od ? (slots + 1) >> 1 : slots >> 1);
    if (od)
        result += divop(result, slots);
    return result;
}
```
這題除法依照除數 slots 的奇偶來分成兩個 result。

slots 為偶數時 : (這裡將 orig 設成 m, slots 設成 n)
$\dfrac{m / 2}{n / 2}= \dfrac{m}{n}$

slots 為奇數時 : 
$\dfrac{m/2}{(n+1)/2} + \dfrac{\dfrac{m/2}{(n+1)/2}}{n} = \dfrac{n \times m/2  + m/2}{\dfrac{n \times (n + 1)}{2}} = \dfrac{(n + 1) \times m/2}{\dfrac{n \times(n+1)}{2}} = \dfrac{m}{n}$

接著一直做遞迴直到終止條件 `(slots == 1 || orig == 0)`。
## 測驗 2
:::spoiler 題目
```c=
#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;
}
```
:::
這題做的是浮點數的開根號，EXTRACT_WORDS 將 float 改用 int 表示(位元不變動)，INSERT_WORDS 則相反。這題先將 float 改成 int 來方便做運算，算好值之後再轉換成 float。
```c=36
    /* 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 */
    }
```
這段在處理 0 和負數的問題，如果除了 sign bit 以外都是 0，代表此數是 0 (雖然 ix0 是 int，但實際上表示的是 float，所以有 +-0)，直接 return 原數。而如果是負數的話，則 return NaN。
```c=43
    /* take care of +INF and NaN */
    if ((ix0 & 0x7f800000) == 0x7f800000) {
        /* sqrt(NaN) = NaN, sqrt(+INF) = +INF,*/
        return x;
    }
```
這段對 infinity 和 NaN 做處理，由於這兩數開根號後依然不變，所以也是直接返回原數。
```c=48
    /* 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] */
```
`m = (ix0 >> 23);` 取出 exponent 部分。
由於浮點數有 normalize 和 denormalize 模式，50 ~ 54 行將 denormalize 轉換成 normalize，
`ix0 <<= 1;` 將 denormalize 的數乘以二(由於這裡一定是正數，所以這邊不用將 sign bit 補回)，
`m -= i - 1;` 則是將乘過的次數補回來，之所以會 +1 是因為 denormalize 時是 $\times 2^{-126}$，但轉換成 normalize 後為 $\times 2^{exponent-127}$，以上部分可以參考[浮點數運算和定點數操作](https://hackmd.io/@NwaynhhKTK6joWHmoUxc7Q/H1SVhbETQ)

`m -= 127;` 將原本的 exponent 改成更直觀的寫法，`ix0 = (ix0 & 0x007fffff) | 0x00800000;` 把 fraction 部分取出並加上 1，到這邊為止，被開根號的數表示為 $ix0 \times 2^m，1 \leq ix0 < 2$

開根號後，我們可以輕易地處理 $2^m$ 的部分。
如果 m 是偶數，則最後的結果是 $\sqrt{ix0} \times 2^{\lfloor m/2 \rfloor}$
如果 m 是奇數，則最後的結果是 $\sqrt{2 \times ix0} \times 2^{\lfloor m/2\rfloor}$
也就是 57 ~ 60 行的內容。
```c=62
/* 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;     //t = s_i + 2^(-(i+1))
        if (t <= ix0) { //s_i + 2^(-(i+1)) <= y_i
            s0 = t + r; //s_(i+1) = s_i + 2^(-i)
            ix0 -= t;   //y_(i+1) = y_i - s_i - 2^(-(i+1)) (* 2 later)
            q += r;     //q_(i+1) = q_i + 2^(-(i+1))
        }
        ix0 += ix0;     //y_(i+1) = y_i * 2
        r >>= 1;
    }
```
這裡部分參考 [RinHizakura](https://hackmd.io/@RinHizakura/SJnXiSewD#%E6%B8%AC%E9%A9%97-2) 的解釋。
第 62 行我想了蠻久的，我覺得這裡應該不是對數值本身乘二，而是將小數點往左移一格，所以目前的小數點在(右邊數來)第 24, 25 位之間，因此 
r = 0x01000000;

$q_i$ 表示為 $\sqrt y$ 到小數點後第 i 位的近似值，$q_0$     = 1。
而依照 ${(q_{i} + 2^{-(i + 1)})}^2 \leq y$ 判斷式來決定 $q_{i+1}$ 的值，若上式成立時，$q_{i+1} = q_i + 2^{-(i+1)}$ (下一位補 1)，反之 $q_{i+1} = q_i$ (下一位補 0)。

把上式做展開後，可得到
$$
q_i^2 + q_i \times2^{-i} + 2^{-2(i+1)}\leq y
\\
q_i \times2^{-i} + 2^{-2(i+1)}\leq y -q_i^2
\\
q_i \times 2 + 2^{-(i+1)}\leq (y-q_i^2)\times 2^{i+1}
$$
令以下兩個變數
$$
s_i = q_i\times 2
\\
y_i = (y-q_i^2)\times 2^{i+1}
$$
於是原本的不等式變成
$$
s_i + 2^{-(i+1)}\leq y_i\space
$$
若以上式子成立時
$$
s_{i+1} = s_i + 2^{-i}
\\
y_{i+1} = (y-q_{i+1}^2)\times 2^{i+2} = (y - q_i^2 - q_i \times2^{-i} - 2^{-2(i+1)})\times2^{i+2} = (y_i-s_i-2^{-(i+1)})\times 2
$$
不成立時
$$
s_{i+1} = s_i
\\
y_{i+1} = y_i \times 2
$$
```c=78
    /* 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);
        }
    }
```
這裡在判斷系統是哪種進位方式，一開始先把 `z = one - tiny`，由於 tiny 遠小於 one，所以如果此時 `z >= one` 成立，代表系統可能為無條件進位或四捨五入，再透過 `z = one + tiny` 以及 `z > one` 來分辨無條件進位以及四捨五入。

由於後面會將 `q >> 1`，所以這裡需要考慮進位問題，如果系統是無條件進位的話，就直接把第 2 位，也就是最後結果的第 1 位加 1，
而如果是四捨五入的話，就會做 `q += (q & 1)` ，如此一來，如果第 1 位是 1 的話會進位，反之則不進位。

```c=90
    ix0 = (q >> 1) + QQ1;
    ix0 += (m << QQ2);
```
經過上面的了解後，最後這部分就比較容易回答了。由於前面把小數點往前移了一格，所以這裡做 `q >> 1` 把小數點移到 float 格式中的位置，又因為原本把 exponent 部分減 127，所以這裡要加回來，因為 `q >> 1` 的 MSB 必定是第 24 位 ($1\leq \sqrt y < 2$)，所以這裡的 QQ1 = 0x3f000000。再把 $2^m$ 乘回來，所以 QQ2 = 23
### 延伸問題

#### 練習 LeetCode 69. Sqrt(x)，提交不需要 FPU 指令的實作
```c=
int mySqrt(int x) {
    if (x == 0 || x == 1)
        return x;
    int a = x / 2;
    int b = 1 << (((32 - __builtin_clz(x)) >> 1) - 1);
    while(b <= a) {
        int c = (a + b) / 2;
        if (c == x / c)
            return c;
        if (c < x / c)
            b = c + 1;
        if (c > x / c)
            a = c - 1;
    }
    return a;
}
```
![](https://i.imgur.com/dOKlPwk.png)

利用二分法逼近 x 的平方根，在 x >= 2 條件下，x / 2 >= $\lfloor \sqrt x\rfloor$，
網路上查到的方法中，b 多半設為 1，這裡用了這種算法可以逼近得更快。
## 測驗 3
給定一個正整數 N，問 N 能寫成多少種連續正整數之和，假設正整數從 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}
$$
```c=
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;
} 
```
第七行的 i 就是上面的 k，從 2 開始試，因為 k == 1 一定成立 (x == N)。

$\dfrac{(k-1)k}{2} = 1 + 2 + ... +(k - 1)$，所以函式中的 $X = N - \dfrac{(k-1)k}{2} > 0$ (初始值為N，之後在第八行加 (1 - i))(為了避免跟前面的正整數起始 x 混淆，這裡用大寫表示)，而第九行在判斷 $N - \dfrac{(k-1)k}{2}$ 是否等於 $kx$，如果是的話代表 k 等於該值的時候有解，所以 `ret++`

在第七行有一個 `i < X` 的終止條件，因為 $k <= N - \dfrac{(k-1)k}{2}$，且 $k = N - \dfrac{(k-1)k}{2}$ 的情況已在第九行時就有判斷過了，所以成立。