---
tags: linux2022
---
# 2022q1 Homework2 (quiz2)
contributed by < [howardjchen](https://github.com/howardjchen) >
### 測驗 1
考慮以下對二個無號整數取平均值的程式碼:
```c
#include <stdint.h>
uint32_t average(uint32_t low, uint32_t high)
{
return low + (high - low) / 2;
}
```
接著我們可改寫為以下等價的實作:
```c
#include <stdint.h>
uint32_t average(uint32_t a, uint32_t b)
{
return (a >> 1) + (b >> 1) + (EXP1);
}
```
- **想法&思路**
```(a >> 1) + (b >> 1)``` 可以看成 a/2 + b/2 但不等於 (a+b)/2 因為會涉及到進位的問題
如果 a 或 b 有兩個都是奇數的話, a/2 + b/2 跟 (a+b)/2 就會相差 1
因此 EXP1 的作用是如果 a 和 b 都是奇數的時候需要加 1
```
EXP1 = a & b & 1
```
所以可以改寫成
```c
#include <stdint.h>
uint32_t average(uint32_t a, uint32_t b)
{
return (a >> 1) + (b >> 1) + (a & b & 1);
}
```
我們再次改寫為以下等價的實作:
```c
uint32_t average(uint32_t a, uint32_t b)
{
return (EXP2) + ((EXP3) >> 1);
}
```
- **想法&思路**
這個直觀上很像 ```return low + (high - low) / 2;```
因此```EXP2 = low``` ```EXP3 = high - low``` 但是並不知道誰是 low 跟 high
另一個想法是從[你所不知道的 C 語言：數值系統篇](https://hackmd.io/@sysprog/c-numerics) 有提到為了避免 overflow
- ```(x + y) / 2``` 這樣的運算，有個致命問題在於 (x + y) 可能會導致 overflow (考慮到 x 和 y 都接近 [UINT32_MAX](https://msdn.microsoft.com/en-us/library/mt764276.aspx)，亦即 32-bit 表示範圍的上限之際)
* 於是我們可以改寫為以下:
```cpp
(x & y) + ((x ^ y) >> 1)
```
* 用加法器來思考: `x & y` 是進位, `x ^ y` 是位元和, `/ 2` 是向右移一位
* 位元相加不進位的總和: `x ^ y`; 位元相加產生的進位值: `(x & y) << 1`
* `x + y = x ^ y + ( x & y ) << 1`
* 所以 (x + y) / 2
= `(x + y) >> 1`
= `(x ^ y + (x & y) << 1) >> 1`
= `(x & y) + ((x ^ y) >> 1)`
因此
```
EXP2 = a & b
EXP3 = a ^ b
```
```c
uint32_t average(uint32_t a, uint32_t b)
{
return (a & b) + ((a ^ b) >> 1);
}
```
> 2. 比較上述實作在編譯器最佳化開啟的狀況，對應的組合語言輸出，並嘗試解讀 (可研讀 CS:APP 第 3 章)
我們可以透過下列指令做編譯最佳化並輸出處和語言 ```test.s```
```shell
gcc -S -O3 test.c
```
**(high + low) / 2;**
```c=
endbr64
leal (%rsi,%rdi), %eax
shrl %eax
ret
```
[What does the LEAL assembly instruction do?](https://stackoverflow.com/questions/11212444/what-does-the-leal-assembly-instruction-do)
> LEA (load effective address) just computes the address of the operand, it does not actually dereference it. Most of the time, it's just doing a calculation like a combined multiply-and-add for, say, array indexing.
- 簡單來說 leal 可以看成相加, 所以把兩個變數相加後放在 ```%eax```
- shrl 是 shift left, 把```%eax``` 除 2
**low + (high - low) / 2;**
```c=
endbr64
movl %esi, %eax
subl %edi, %eax
shrl %eax
addl %edi, %eax
ret
```
- 前兩行就是相減
- 第四行向右移 1
- 最後加上 low
**(a >> 1) + (b >> 1) + (a & b & 1);**
```c=
endbr64
movl %edi, %eax
movl %esi, %edx
andl %esi, %edi // a & b
shrl %eax // a >> 1
shrl %edx // b >> 1
andl $1, %edi // a & b & 1
addl %edx, %eax // (a >> 1) + (a & b & 1)
addl %edi, %eax // (a >> 1) + (b >> 1) + (a & b & 1)
ret
```
**(a & b) + ((a ^ b) >> 1);**
```c=
endbr64
movl %edi, %eax
andl %esi, %edi // a & b
xorl %esi, %eax // a ^ b
shrl %eax // (a ^ b) >> 1
addl %edi, %eax // (a & b) + ((a ^ b) >> 1)
ret
```
> 研讀 Linux 核心原始程式碼 include/linux/average.h，探討其 Exponentially weighted moving average (EWMA) 實作，以及在 Linux 核心的應用
### 測驗 2
改寫 [解讀計算機編碼](https://hackmd.io/@sysprog/binary-representation) 一文的「不需要分支的設計」一節提供的程式碼 min
```c=
int32_t min(int32_t a, int32_t b) {
int32_t diff = (a - b);
return b + (diff & (diff >> 31));
}
```
針對 max 的情形我們可以改寫成
```c=
int32_t max(int32_t a, int32_t b) {
int32_t diff = (a - b);
return a - (diff & (diff >> 31));
}
```
若要更近一步改寫成題目的形式
```c=
#include <stdint.h>
uint32_t max(uint32_t a, uint32_t b)
{
return a ^ ((EXP4) & -(EXP5));
}
```
所以 EXP4 = a ^ b
- if a > b, then return a
那就表示 ```((a ^ b) & -(EXP5)) = 0```
-(EXP5) 需要等於 0
也就是 EXP5 = 0
- if a <= b, then return b
那就表示 ```((a ^ b) & -(EXP5)) = a ^ b```
-(EXP5) 需要等於 0xffffffff'h = -1'd
也就是 EXP5 = 1
所以 EXP5 = -((a - b) >> 31)
因此我們可以重新把這個 function 寫成
```c=
#include <stdint.h>
uint32_t max(uint32_t a, uint32_t b)
{
return a ^ ((a ^ b) & ((a - b) >> 31));
}
```
- 延伸問題:
2. 針對 32 位元無號/有號整數，撰寫同樣 branchless 的實作
不管是有號整數或是無號整數， (a-b) 都是 -1，所以實作一樣
:::success
4. Linux 核心也有若干 branchless / branch-free 的實作，例如 lib/sort.c:
```c
/*
* Logically, we're doing "if (i & lsbit) i -= size;", but since the
* branch is unpredictable, it's done with a bit of clever branch-free
* code instead.
*/
__attribute_const__ __always_inline
static size_t parent(size_t i, unsigned int lsbit, size_t size)
{
i -= size;
i -= size & -(i & lsbit);
return i / 2;
}
```
請在 Linux 核心原始程式碼中，找出更多類似的實作手法。請善用 git log 檢索。
:::
### 測驗 3
考慮以下 64 位元 GCD (greatest common divisor, 最大公因數) 求值函式:
```c
#include <stdint.h>
uint64_t gcd64(uint64_t u, uint64_t v)
{
if (!u || !v) return u | v;
while (v) {
uint64_t t = v;
v = u % v;
u = t;
}
return u;
}
```
改寫為以下等價實作:
```c
#include <stdint.h>
uint64_t gcd64(uint64_t u, uint64_t v)
{
if (!u || !v) return u | v;
int shift;
for (shift = 0; !((u | v) & 1); shift++) {
u /= 2, v /= 2;
}
while (!(u & 1))
u /= 2;
do {
while (!(v & 1))
v /= 2;
if (u < v) {
v -= u;
} else {
uint64_t t = u - v;
u = v;
v = t;
}
} while (COND);
return RET;
}
```
> 1. if both x, y are 0, $gcd(0,0) = 0$
> 2. $gcd(x,0)=x$ and $gcd(0,y)=y$
```c
if (!u || !v) return u | v;
```
> 3. If x and y are both even, $gcd(x,y) = 2 * gcd(\dfrac{x}{2}, \dfrac{y}{2})$
```c
for (shift = 0; !((u | v) & 1); shift++) {
u /= 2, v /= 2;
}
```
所以$gcd(u, v) = 2^n * gcd(\dfrac{u}{2}, \dfrac{v}{2})$
$n = shift$
> 4. If x is even and y is odd, $gcd(x,y) = 2 * gcd(\dfrac{x}{2}, y)$
> 5. If x is odd and y is even, then $gcd(x,y) = 2 * gcd(x, \dfrac{y}{2})$
```c
while (!(u & 1))
u /= 2;
```
```c
while (!(v & 1))
v /= 2;
```
這兩個表示 u 跟 v 不會同時是偶數
是利用輾轉相除法的原理，只是用了連續減法來找出相除的餘數結果。
輾轉相除法就是除到餘數為 0 就是答案。
最後一步 ```u = v``` 的時候, 就是餘數為 0 的時候，此時 ```uint64_t t = u - v```會讓 t 為 0, 然後被指派給 v.
所以 COND 就是餘數 v , 而 u 在這裡就是最後的那個最大公因數。
但要記得真正的答案是 $2^n * gcd(\dfrac{u}{2}, \dfrac{v}{2})$
$n = shift$
因此 RET 就是 u << shift.
> 在 x86_64 上透過 __builtin_ctz 改寫 GCD，分析對效能的提升;
__builtin_ctz 並不是 C 的標準函式庫，只有在 gcc 中可以使用
:::success
int __builtin_ctz (unsigned int x)
int __builtin_ctzl (unsigned long)
int __builtin_ctzll (unsigned long long)
返回右起第一個1之後的0的個數
Returns the number of trailing 0-bits in x, starting at the least significant bit position. If x is 0, the result is undefined.
:::
似乎可以改善
```c
int shift;
for (shift = 0; !((u | v) & 1); shift++) {
u /= 2, v /= 2;
}
```
成為
```c
int shift_u = __builtin_ctz(u);
int shift_v = __builtin_ctz(v);
int shift = (shift_u > shift_v)? shift_v : shift_u;
```
另一方面還有
```c
while (!(u & 1))
u /= 2;
```
成為
```c
u = u >> __builtin_ctz(u);
```
因此可以改寫程式碼如下：
```c
#include <stdint.h>
uint64_t gcd64(uint64_t u, uint64_t v)
{
if (!u || !v) return u | v;
int shift_u = __builtin_ctz(u);
int shift_v = __builtin_ctz(v);
int shift = (shift_u > shift_v)? shift_v : shift_u;
u >>= shift_u;
do {
v >>= _shift_v;
if (u < v) {
v -= u;
} else {
uint64_t t = u - v;
u = v;
v = t;
}
} while (v);
return (u << shift);
}
```
各執行了 100000 次得到結果，使用```__builtin_ctz```可以提升約 1.49 倍的效能
```c
size_t rand64()
{
time_t t;
srand((unsigned) time(&t));
return (rand() % 0xffffffff);
}
float execute_cost(size_t (*func)(size_t, size_t), int times)
{
clock_t start = clock();
int cnt = 0;
while (times--)
{
size_t ans = func(testbench[cnt], testbench[cnt+1]);
cnt++;
}
clock_t end = clock();
return (float)(end - start) / CLOCKS_PER_SEC;
}
```
```bash
gcd64 cost: 0.000961
gcd64_f cost: 0.000641
```
> Linux 核心中也內建 GCD (而且還不只一種實作)，例如 lib/math/gcd.c，請解釋其實作手法和探討在核心內的應用場景
Sign in with Wallet
Connect another wallet