# 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;
}
```
利用二分法逼近 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}$ 的情況已在第九行時就有判斷過了，所以成立。
Sign in with Wallet
Connect another wallet