# 2023-06-18 測驗
> 受測者: POCHUN-CHEN
## 測驗 `1`
用 Python 計算開平方根:
```python
import math
print(math.sqrt(16))
print(math.sqrt(63))
```
以下程式碼嘗試計算開平方根: (版本一)
```c
#include <math.h>
int isqrt(int N)
{
int msb = (int) log2(N);
int a = 1 << msb;
int result = 0;
while (a != 0) {
if ((result + a) * (result + a) <= N)
result += a;
a >>= 1;
}
return result;
}
```
對應的測試程式:
```c
#include <stdio.h>
int main()
{
int N = 36;
printf("%d\n", isqrt(N));
return 0;
}
```
編譯方式:
```shell
$ gcc -o isqrt isqrt.c -lm
```
改寫上述 `isqrt` 函式,使其不依賴 [log2](https://man7.org/linux/man-pages/man3/log2.3p.html) 並保持原本的函式原型 (prototype) 和精度 (precision)。
:arrow_down_small: 作答區 (完整的 `isqrt` 函式,版本二)
```c
int isqrt(int N)
{
int msb = 0;
int n = N;
while (n > 1) {
n >>= 1;
msb++;
}
int a = 1 << msb;
int result = 0;
while (a != 0) {
if ((result + a) * (result + a) <= N)
result += a;
a >>= 1;
}
return result;
}
```
Linux 核心的 [lib/math/int_sqrt.c](https://github.com/torvalds/linux/blob/master/lib/math/int_sqrt.c) 採用 [Digit-by-digit calculation](https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Digit-by-digit_calculation),也就是常見的直式開方法。程式碼改寫如下: (版本三)
```c
int isqrt(int x)
{
if (x <= 1) /* Assume x is always positive */
return x;
int y = 0;
for (int m = 1UL << ((31 - __builtin_clz(x)) & ~1UL); m; m >>= 2) {
int b = y + m;
y >>= 1;
if (x >= b)
x -= b, y += m;
}
return y;
}
```
直式開方法是把要開方的數拆成 2 的冪相加 (若是十進位就是 10 的冪),就可以拆成以下的形式:
$$
x = (000b_0b_1b_2\cdots b_{n-1}b_{n})^2
$$
上面的 $000b_0b_1b_2\cdots b_{n-1}b_{n}$ 為 bits pattern ,其中 $b_0$ 為最高位的 1,寫成數學式如下:
$$
x = (2^{n}b_0+2^{n-1}b_1+\cdots+2^1b_{n-1}+2^0b_{n})^2 \\
a_i = 2^{n-i}b_i \\
\begin{aligned}x
&=(a_0+a_1+\cdots+a_{n-1}+a_{n})^2 \\
&= a_0^{2}+2a_0a_1 + a_1^2 + 2(a_0+a_1)a_2 + a_2^2 + \cdots + a_{n-1}^2 + 2\left(\sum\limits_{i = 0}^{n-1}{a_i}\right)a_n + a_n^2 \\
&=a_0^2 + (2a_0+a_1)a_1 + [2(a_0+a_1)+a_2]a_2 + \cdots + [2\left(\sum\limits_{i = 0}^{n-1}{a_i}\right)+a_n]a_n
\end{aligned}
$$
最後的形式為 $[2\left(\sum\limits_{i = 0}^{n-1}{a_i}\right)+a_n]a_n$ 的相加,其中每項就對應到上述程式碼的每一輪迴圈
```c
if (x >= b) {
x -= b;
y += m;
}
```
要判斷並可能減掉的 `b`。
原理是將欲求平方根的 $N^2$ 拆成:
$$
N^2=(a_{n}+a_{n-1}+a_{n-2}+...+{a_0})^2, a_m=2^m\ or\ a_m=0
$$
將其乘開得到
$$
\left[
\begin{array}{ccccc}
a_0a_0 & a_1a_0 & a_2a_0 & ... & a_na_0 \\
a_0a_1 & a_1a_1 & a_2a_1 & ... & a_na_1 \\
a_0a_2 & a_1a_2 & a_2a_2 & ... & a_na_2 \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
a_0a_n & a_1a_n & a_2a_n & ... & a_na_n \\
\end{array}
\right]
$$
分成幾個部份觀察
$$
\left[
\begin{array}{c}
a_0a_0 \\
\end{array}
\right]
$$
$$
\left[
\begin{array}{cc}
& a_1a_0 \\
a_0a_1 & a_1a_1
\end{array}
\right]
$$
$$
\left[
\begin{array}{ccc}
& & a_2a_0 \\
& & a_2a_1 \\
a_0a_2 & a_1a_2 & a_2a_2
\end{array}
\right]
$$
原式可整理成
$$
\begin{split} \\
N^2 =& a_n^2+[2a_n+a_{n-1}]a_{n-1}+[2(a_n+a_{n-1})+a_{n-2}]a_{n-2}+...+[2(\displaystyle\sum_{i=1}^{n}a_i)+a_0]a_0 \\
=& a_n^2+[2P_n+a_{n-1}]a_{n-1}+[2P_{n-1}+a_{n-2}]a_{n-2}+...+[2P_1+a_0]a_0 \\
P_m =& a_n+a_{n-1}+...+a_m \\
\end{split}
$$
$P_0=a_n+a_{n-1}+...+a_0$ 即為所求平方根 $N$
顯然 $P_m = P_{m+1} + a_m$。我們的目的是從試出所有 $a_m$,$a_m=2^m$ or $a_m=0$。從 $m=n$ 一路試到 $m=0$,而求得 $a_m$ 只要試驗 $P_m^2 \leq N^2$ 是否成立,二種可能狀況:
$$
\begin{cases}
P_m = P_{m+1} + 2^m, & \text{if $P_m^2 \leq N^2$}\\
P_m = P_{m+1}, & \text{otherwise}
\end{cases}
$$
但我們不可能每輪去計算 $N^2 - P_m^2$ 去試驗 $a_m$,因為成本太高。一種想法是從上一輪的差值 $X_{m+1}$ 減去 $Y_m$ 得到 $X_{m}$
- $X_m = N^2 - P_m^2 = X_{m+1} - Y_m$
- $Y_m = P_m^2 - P_{m+1}^2 = 2P_{m+1}a_m+a_m^2$
也就是紀錄上一輪的 $P_{m+1}$ 來計算 $Y_m$
進一步調整,將 $Y_m$ 拆成 $c_m, d_m$
$$
c_m = P_{m+1}2^{m+1} \\
d_m = (2^m)^2 \\
Y_m=\left.
\begin{cases}
c_m+d_m & \text{if } a_m^2=2^m \\
0 & \text{if } a_m=0
\end{cases}
\right.
$$
藉由位元運算從 $c_m, d_m$ 推出下一輪$c_{m-1}, d_{m-1}$
- $c_{m-1}=P_m2^m=(P_{m+1}+a_m)2^m=P_{m+1}2^m+a_m2^m=\begin{cases}
c_m/2+d_m & \text{if }a_m=2^m \\
c_m/2 & \text{if }a_m=0
\end{cases}$
- $d_{m-1}=\dfrac{d_m}{4}$
結合上述方法撰寫演算法來尋求$a_n+a_{n-1}+...+a_0$,假設從 $a_n$ 開始往下試,至於如何求得 $a_n$ 後面再說明。因為是從 $a_n$ 開始,所以:
- $X_{n+1}=N$
- 一開始都沒猜,差值是 N
- $n$ 是最大的位數,使得 $a_n^2=(2^n)^2= d_n^2 = 4^n \leq N^2$ 所以用迴圈逼近 $d_n$,
```c
int32_t d = 1 << 30; // The second-to-top bit is set.
while (d > n)
d >>= 2;
```
- $c_n = 0$
- $c_m=P_{m+1}2^{m+1}$。$a_n$ 已是已知最大,所以 $P_{n+1} = 0 \to c_n=0$
因為 $d_m$ 恆大於 $0$,使用位移操作 $d_m$ 到最後必為 $0$,我們可以使用 `d!=0` 作為條件。又 $c_{-1} = P_0 = a_n+a_{n-1}+...+a_0$,算到最後的 $c_{-1}$ 即為所求 $N$,$a_n$ 也一併得知。
對應程式碼
由上述的理論可知,求出根號需要求得 $c_{-1}$ 即為所求 $N$ 。故寫出迭代的關係,即可求得解。
首先,先判斷輸入的數值,是否符合求解根號的範圍。
```c
if (x <= 1) /* Assume x is always positive */
return x;
```
由於開根號,本身在負數與 `0` 是無意義的,故不討論。且小數不能被整數型別 `int` 儲存,所以也不用討論。
`int m` 對應上述理論中的 $d_{m}$ ,
由於首項 $d_n = (2^n)^2$,可知
```c
int m = 1UL << ((31 - __builtin_clz(x)) & ~1UL)
```
> [`__builtin_clz(x)`](https://gcc.gnu.org/onlinedocs/gcc/Other-Builtins.html) 為 GCC 的內建函式。
> Returns the number of leading 0-bits in x, starting at the most significant bit position. If x is 0, the result is undefined.
> 其功能為找到 x 最高位數前面有幾個 0 , x 為 0 時,未定義。
> `1UL` 為 `(unsigned long) 1` 。
`int m` 為了求得 32 位元下,輸入參數 `x` 的最高有效位數,並將其向下取成偶數。
:::warning
詳細待補的部分:
1. 並將其向下取成偶數,這裡的寫法我還是有點不確定,為什麼奇數會有問題,目前推估是在做 Bitwise >> 2 時的問題。
2. 首項 $d_n = (2^n)^2$ ,但程式碼中,似乎只有取出最高位數($2^n$),並沒有做平方。
:::
`m >>= 2` 為 $d_{m-1}=\dfrac{d_m}{4}$ 除以 4 的計算。
由程式與理論對應可得知,
$b=Y_{m}$
$y=c_{m}$
$x=X_{m}$
最終經由迭代,得知$c_{-1}$ 即為所求 $N$ 。
### 效能測試
使用 `perf stat` 來進行測試。
為了不讓 I/O 影響到測試的時間,測試時先移除 `printf("%d\n", isqrt(N));`,並重複執行 100000 次,比較執行時間。
編譯後,執行 `perf stat`
```shell
$ perf stat -r 100000 ./isqrt1
```
版本一
```
Performance counter stats for './isqrt1' (100000 runs):
0.29 msec task-clock # 0.561 CPUs utilized ( +- 0.01% )
0 context-switches # 0.000 /sec
0 cpu-migrations # 0.000 /sec
62 page-faults # 211.765 K/sec ( +- 0.00% )
99,4673 cycles # 3.397 GHz ( +- 0.01% )
104,6590 instructions # 1.06 insn per cycle ( +- 0.00% )
19,6820 branches # 672.253 M/sec ( +- 0.00% )
5070 branch-misses # 2.58% of all branches ( +- 0.03% )
0.000525393 +- 0.000000108 seconds time elapsed ( +- 0.02% )
```
版本二
```
Performance counter stats for './isqrt2' (100000 runs):
0.24 msec task-clock # 0.500 CPUs utilized ( +- 0.01% )
0 context-switches # 0.000 /sec
0 cpu-migrations # 0.000 /sec
51 page-faults # 211.142 K/sec ( +- 0.00% )
79,3884 cycles # 3.287 GHz ( +- 0.01% )
85,2080 instructions # 1.05 insn per cycle ( +- 0.00% )
15,9847 branches # 661.772 M/sec ( +- 0.00% )
4015 branch-misses # 2.51% of all branches ( +- 0.04% )
0.000470606 +- 0.000000101 seconds time elapsed ( +- 0.02% )
```
版本三
```
Performance counter stats for './isqrt3' (100000 runs):
0.24 msec task-clock # 0.502 CPUs utilized ( +- 0.01% )
0 context-switches # 0.000 /sec
0 cpu-migrations # 0.000 /sec
51 page-faults # 211.153 K/sec ( +- 0.00% )
79,5358 cycles # 3.293 GHz ( +- 0.01% )
84,9166 instructions # 1.04 insn per cycle ( +- 0.00% )
15,9826 branches # 661.720 M/sec ( +- 0.00% )
4148 branch-misses # 2.59% of all branches ( +- 0.03% )
0.0004700040 +- 0.0000000973 seconds time elapsed ( +- 0.02% )
```
單從執行 100000 次的測試,可知版本三的指令數量 (instructions) 最少,且執行時間最短(time elapsed)。
---
## 測驗 `2`
以下程式碼計算整數除 `10`,不使用 `/` (除法) 或 `%` (modulo) 運算子。
```c
static inline int32_t div10(const int32_t x)
{
return (int32_t) (((int64_t) x * INT64_C(/* Your answer */)) >> INT64_C(35));
}
```
補完上述程式碼,使其得以運作。
:arrow_down_small: 作答區 (完整的 `div10` 函式)
```c
static inline int32_t div10(const int32_t x)
{
return (int32_t) (((int64_t) x * INT64_C(0xCCCCCCCD)) >> INT64_C(35));
}
```
敘述 `>> INT64_C(35)` 實際作用為「除 $2^{35}$」,故 `INT64_C(/* Your answer */))` 為 $2^{35}\times\dfrac{1}{10}$ ,這時候答案會是 `0xCCCCCCCC` ,但是經由測試,當**整除**發生時,商數會比預期的小 1。
誤差原因:
`0xCCCCCCCC` = 3435973836~10~,且 $2^{35} = 34359738368$,利用 $2^{35} \times \dfrac{1}{10}$ 得到 `0xCCCCCCCC` 時,已捨棄掉最低一位數,導致問題出現。
舉例來說,若 $x = 10 \times y$,商為 `y` ,則
$$
\begin{split}
x \times 3435973836 \div 34359738368 &= y \times 10 \times 3435973836 \div 34359738368 \\
&= y \times 0.099999999976717 \\
\end{split}
$$
這時候即可知 $2^{35}\times\dfrac{1}{10}$ 捨去掉的最低位數,會出現計算商數比正確商數少一的情況 (原本應為 $y\times0.1$)。
`0xCCCCCCCD` 則是在被捨去掉最低一位數字的數值 (`0xCCCCCCCC`) 的狀況,加上 `1` ,確保捨去掉的位數,不會讓原先的數值低於捨去後(此做法建立在商數為整數時,不會因為捨去最低位數,而被無條件捨去成小於正確商數值減 `1` 的商值)。
負數商:
經由改寫成 `0xCCCCCCCD` 時,帶入負數 `-20` ,答案會出現 `-3`。
負數情況,由於商數規則是捨去掉小數位,原本 `-5` 的商值依照二補數,捨去後,計算的商值會是 `-1` ($-10 = -1 \times 10$) ,但這不符合一般數學負數商數的規定,故 `+1` 符合負數商對應正數商的關係。
根據 [C11 規格書](https://www.open-std.org/jtc1/sc22/WG14/www/docs/n1570.pdf)第 6.5.5 節 Multiplicative operators 的描述:
> 5. The result of the / operator is the quotient from the division of the first operand by the second; the result of the % operator is the remainder. In both operations, if the value of the second operand is zero, the behavior is undefined.
>
> 6. When integers are divided, the result of the `/` operator is the algebraic quotient with any fractional part discarded. If the quotient `a/b` is representable, the expression `(a/b)*b + a%b` shall equal `a`; otherwise, the behavior of both a/b and a%b is undefined.
根據以上兩個規則敘述,可知:
1. 當除數為 `0` ,為未定義行為。
2. 商數運算是捨棄掉除法運算後的小數部分,如果商數運算不超出資料型別的表示範圍,則符合 `(a / b) * b + a % b = a` 也就是 「$商數 \times 除數 + 餘數 = 被除數$」。超出資料型別範圍則未定義。
因為商數運算是**捨棄掉除法運算後的小數部分**,可以從此規則上得知
3. 商的結果會朝向 `0` 取整數。
4. 餘數一定和被除數同正負號。
故以上調整是合理的,程式碼中加入負數的處理,如下。
```c
static inline int32_t div10(const int32_t x)
{
int32_t result = (int32_t) (((int64_t) x * INT64_C(0xCCCCCCCD)) >> INT64_C(35));
return (x < 0) ? (result + 1): result;
}
```
將上式改成 branchless 寫法
```c
static inline int32_t div10(const int32_t x)
{
int32_t result = (int32_t) (((int64_t) x * INT64_C(0xCCCCCCCD)) >> INT64_C(35));
return result + (x >> 31 & 1);
}
```
---
## 測驗 `3`
延續上題,針對除數為 `10`、給定被除數,輸出商與餘數,不得使用除法或 modulo 運算子。
```c
typedef struct {
int32_t quotient;
int32_t rem;
} divmod_t;
static inline void divmod10(const int32_t x, divmod_t *ret)
{
/* Your answer here */
}
```
:arrow_down_small: 作答區 (完整的 `divmod10` 函式)
先以不改寫上題題目為情況下寫此題
```c
static inline void divmod10(const int32_t x, divmod_t *ret)
{
ret->quotient = (int32_t) (((int64_t) x * INT64_C(0xCCCCCCCD)) >> INT64_C(35));
ret->rem = (int32_t) (x - ret->quotient * 10);
}
```
:::warning
TODO:
1. 提供必要的錯誤檢查
2. 參照 [第 7 週課堂問答簡記](https://hackmd.io/@sysprog/SyOhklg-3),改寫上述程式碼,探討執行成本。
:::