---
# System prepended metadata

title: 2020q3 Homework5 (quiz5)
tags: [sysprog2020]

---

# 2020q3 Homework5 (quiz5)
contributed by < `joe-U16` >
###### tags: `sysprog2020`

## 測驗 `1`
考慮到以下浮點數除法程式: (fdiv.c)
```cpp=
#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;
}
// D1 = 2
// D2 = 1
```
假設 `divop()` 的第二個參數必為大於 `0` 的整數，而且不超過 `int` 型態能表達的數值上界。請補完程式碼。

**程式碼運作原理:**
* `int od = slots & 1;`
    * 如果 slots 為奇數，則 `od = 1`
* `double result = divop(orig / D1, od ? (slots + D2) >> 1 : slots >> 1);`
    * 先看 D2 ，如果 `od = 1` 則下一次的 slots 會是 `(slots + D2) >> 1`
    * 如果 D2 是 2、3、4 的話，那每一個遞迴函式的 `slots` 都不會變為 1 ，這會讓最前面 `if (slots == 1 || orig == 0)` slot判斷沒有意義，所以 `D2 = 1`
* 再來是這個程式碼做除法的方式
    * 把被除數和除數一起除以 2 ，當除數 slots 為 1 時即代表完成，因為任何數除以 1 都是自己。
    * 可以合理推得 `D1 = 1`

## 測驗 `2`

先看到程式的兩個巨集：
```c=
/* 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)
```
上面巨集如同註解所述，把 int 轉換成 float ，且不經過換算，是直接轉換。

e.g. 0x40800000，以 int 表示即為 0x40800000 ，轉成以 float 表示 0x40800000 ，需使用 ieee754 標準，會使其所代表的數值為 4。

```c=
/* 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)

```
上面巨集如同註解所述，把 float 轉換成 int ，且不經過換算，是直接轉換。

e.g. 數值 4 以 float 表示需使用 32 bit，以浮點數表示成 0x40800000，直接轉成 int 為 0x40800000。

```c=
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 */
    }
```
使用巨集 **EXTRACT_WORDS(ix0, x)** 以 整數 `ix0` 儲存 `x` 的浮點數表示法。

* `if ((ix0 & (~sign)) == 0)` 判斷 Exponent 和 Mantissa 是否全為 `0` ，成立則代表 `ix0` 為 `+-0`，即可 return
* `if (ix0 < 0)` return sNaN

```c=
/* take care of +INF and NaN */
    if ((ix0 & 0x7f800000) == 0x7f800000) {
        /* sqrt(NaN) = NaN, sqrt(+INF) = +INF,*/
        return x;
```
* `if ((ix0 & 0x7f800000) == 0x7f800000)` 判斷 Exponent 是否全為 `1`

```c=
/* 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] */
```
* 做 `sqrt(x)` 可將 x 表示成 $x = y \cdot 2^{2k}$ ，那麼 sqrt(x) = sqrt(y) $\cdot 2^k$
    * 要確保 `k` 為整數， `m = 2k`， `m` 需為偶數。
    * y 為 `1.fraction`，為了確保 m 為偶數，可能做成 $x = (y\cdot2)\cdot2^{2k-1}$
    * 故 y 的範圍會在 `(1, 4]`
* `m = ix0 >> 23` 取到前 9 個bit，即為 Sign + Exponent。
* `if (m == 0)` 判斷是否為 Denormalize
* 若是，則 `i` 儲存 `ix0` fraction 部分的 leading zero
* `ix0 = (ix0 & 0x007fffff) | 0x00800000;` 將 `ix0` 取成 `1.fraction` 的形式，用到後面24個bit
* 浮點數公式為： $(-1)^s\cdot2^{e-127}\cdot(1.m)$ 
    * `s: sign bit`, `e = exponent`,  `m = fraction`
* `if (m & 1) ix0 += ix0` ，這邊的 `m` 為 exponent ，當 `m` 為奇數，將 `ix0` 乘 2 ，可看成將 exponent 的 2 丟給 `ix0`，再將 `m` 減 1，將 m 變成偶數，當然因為待會會做 right shift ，有沒有減 1 不影響 right shift 的結果，故可將減 1 這個動作省略


```c=
/* 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) {  // if t <= y_i
            s0 = t + r;  // s_{i+1} = s_i + 2^(-i)
            ix0 -= t;    // y_{i+1} = y_i - s_i - 2{-(i+1)}
            q += r;
        }
        ix0 += ix0;      // 
        r >>= 1;
    }
// QQ0 = 0x01000000
```

讓 $q_i$ = sqrt(y) ， $q_0$=1 ，  e.g.  $q_0.q_1q_2q_3q_4q_5...$ ， `1.fraction`

$s_i = 2\cdot{q_i}\ \ \ \  {and} \ \ \ y_i = 2^{i+1}\cdot{(y - q_i^2)}\tag{1}$
To compute $q_{i+1}$ from $q_i$,one checks whether

${q_i+2^{-(i+1)}}^2 \le y\tag{2}$

If (2) is false, then $q_{i+1} = q_i$; otherwise $q_{i+1} = q_i + 2^{-(i+1)}$
With some algebric manipulation, it is not difficult to see that (2) is equivalent to

$s_i+2^{-(i+1)} \le y_i\tag{3}$

The advantage of (3) is that $s_i$ and $y_i$ can be computed by the following recurrence formula:

$s_{i+1} = s_i\ \ ,\ \ \ \ y_{i+1} = y_i \cdot2;\tag{4}$
:::warning
原本作者是寫：$s_{i+1} = s_i\ \ ,\ \ \ \ y_{i+1} = y_i ;\tag{4}$

:::

otherwise,

$s_{i+1} = s_i+2^{-i}\ \ ,\ \ \ \ y_{i+1} = y_i-s_i-2^{-(i+1)}\tag{5}$

One may easily use induction to prove (4) and (5)

* QQ0 = 0x01000000 等於可以做25次判斷，最開始的 `ix0+=ix0` 是為了讓最後一個位元可以拿來做rounding。
* 最下面的 `ix0 += ix0;` 即為上述式子(4)

---

* `ix0 = (q >> 1) + QQ1; //QQ1 = 0x3f000000` 
    * 為計算好的 q 加上 exponent的偏移值 `127`
* `ix0 += (m << QQ2); QQ2 = 23` 
    * 把算好的 exponent 放到bit 24-31 ，即浮點數表示法 exponent 的位子



**練習 LeetCode [69. Sqrt(x)](https://leetcode.com/problems/sqrtx/submissions/)，提交不需要 FPU 指令的實作**
使用 binary search 的方式找到 sqrt(x)
```c=
int mySqrt(int x){
    int l = 0;
    int h = 65535; //2^16 - 1 = 65535
    unsigned int m = h;
    
    while(l < h) {
        m = (l + h + 1) / 2;

        if (m * m == x) {
            return m;
        } else if(m * m > x){
            h = m - 1;
        } else {
            l = m;
        }
    }
    return h;
}
```

**Submission Result:**
![](https://i.imgur.com/NTXVcX2.png)



## 測驗 `3`
[LeetCode 829. Consecutive Numbers Sum](https://leetcode.com/problems/consecutive-numbers-sum/) 給定一個正整數 N，問 N 能寫成多少種連續正整數之和，例如 9 = 9 = 4 + 5 = 2 + 3 + 4，所以需 return `3`。

參考實作如下 [1]:
```c=
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 = 1-i
```

**解釋程式原理：**
這個等差數列不必從 1 開始，假設從 x 開始，且個數共有 k 個，則可寫出該等差數列為：
`x + (x+1) + (x+2)+...+ k`
`kx + k*(k-1)/2 = N` implies
`kx = N - k*(k-1)/2`
那麼當 `k` 可整除 `kx` 即代表可以從 x 開始每次加 `1` ，共 k 個連續整數相加可得 正整數 N。

Leetcode 實作如下：
```c=
int consecutiveNumbersSum(int N)
{
    if (N < 1)
        return 0;
    int ret = 1;
    int x = N;
    for (int k = 2; k < sqrt(2 * x); k++) {
        if ((x - (k * (k - 1) / 2)) % k == 0) ret++;
    }
    return ret;
} 
```

可以觀察到上下兩個程式碼原理相同，差別在於程式碼[1] 在 for 迴圈裡不停的加上一個負值，且此負值呈等差級數下降。
* 觀察到下面程式碼，當 k 值上升時的差值。
${k(k+1) \over 2} - {k(k-1) \over 2} = k$
* 可以得知程式碼 [1] 為每次扣掉差值 k ，才做後續判斷。



**Submission Result:**
![](https://i.imgur.com/8oFlqj5.png)




## Reference
https://www.netlib.org/fdlibm/e_sqrt.c
https://leetcode.com/problems/consecutive-numbers-sum/discuss/129015/5-lines-C%2B%2B-solution-with-detailed-mathematical-explanation.


