# 2022q1 Homework2 (quiz2)
contributed by < `Uduru0522` >

###### tags: `linux2022`
## 問題一：Average Of Two Integer

最原始的寫法，逕直將兩樹相加除以二。倘若 `a + b > UINT_MAX` ，會造成 overflow。

#include <stdint.h>
uint32_t average(uint32_t a, uint32_t b)
{
    return (a + b) / 2;
}

> 不考慮以下所有效能增進，再者平均必定介於兩者之間，
> 其實也可以用一個 `uint64_t` 暫存 `a + b` 也可以避免 ovetflow 就是了(笑)

假設已知兩數的大小，可以改寫成下方的樣式，由於平均 `high` 及 `low` 必定介於兩者之間，可以避免 overflow；且已知大小可以避免在相減時造成 underflow。

#include <stdint.h>
uint32_t average(uint32_t low, uint32_t high)
{
    return low + (high - low) / 2;
}

### EXP1
首先，從嘗試從可免除 Overflow 的實作改寫成無除法的形式：

#include <stdint.h>
uint32_t average(uint32_t a, uint32_t b)
{
    return (a >> 1) + (b >> 1) + (EXP1);
}

從數式 $\text{average}=\text{low}+\frac{\text{high} - \text{low}}{2}$，可以推導出 $\text{average}=(\text{low} - \frac{\text{low}}{2}) + \frac{\text{high}}{2}=\color{red}{\frac{\text{low}}{2} + \frac{\text{high}}{2}}$

此時 `low` 及 `high` 之間不再進行計算，可以忽視他們的大小關係，可以 `a` 和 `b` 代替。
而 `a` 及 `b` 各自除以二則可以以 right-shift 代替除法。

然而，僅以 `(a >> 1) + (b >> 1)` 計算時由於 LSB 會被無條件捨棄，有可能造成 0.5 的誤差，當 `a` 和 `b` 雙方皆為奇數時，會導致計算結果相對於正確答案少了一的結果；因此，`EXP1` 便需要用來應對這個狀況。

:::info
舉例來說，當 `a = 3, b = 9`，`(a >> 1) + (b >> 1) = 1 + 4 = 5`，然而兩數平均皆為 6。
:::

綜上所述，`EXP1` 應該有以下的計算結果：
+ 當 `a`、`b` 皆為奇數結果為 **1**。
+ 其他狀況結果為 **0**。

可以得到 `(a & 1) & (b & 1)`，化簡可得 `a & b & 1`。
因此 `define EXP1 a & b & 1`。

### EXP2、EXP3
接著，嘗試更激進的改寫：

uint32_t average(uint32_t a, uint32_t b)
{
    return (EXP2) + ((EXP3) >> 1);
}

模仿 [binary adder](https://www.electronics-tutorials.ws/combination/comb_7.html)，我們可以用 XOR 及 OR 運算替代掉加法：
$\require{cancel}\\
\begin{align}\\
\text{average}&=\frac{a+b}{2}\\
&= \frac{(a\ \mathtt{\text{XOR}}\ b)+(a\ \mathtt{\text{AND}}\ b) \ll 1}{2}\\
&= \frac{a\ \mathtt{\text{XOR}}\ b}{2}+\frac{(a\ \mathtt{\text{AND}}\ b)\cancel{\ll 1}}{\cancel{2}}\\
&= \frac{a\ \mathtt{\text{XOR}}\ b}{2}+(a\ \mathtt{\text{AND}}\ b)
\end{align}$
+ $(a\ \mathtt{\text{AND}}\ b)\times 2$ 計算**進位**，倘若在第 31 個 bit 的計算結果是 1，乘以2 (向左位移 1) 時便會造成溢位，因此利用取平均時除以二的動作消除這個位移。
+ $a\ \mathtt{\text{XOR}}\ b$ 計算**無進位加法**。

再將兩者分別除以 2 後相加，便可以得到兩數平均。接著探討這次相加有沒有可能造成溢位，以下以 $S$ 代稱 $a\ \mathtt{\text{XOR}}\ b\over 2$，$C$ 代稱 $(a\ \mathtt{\text{AND}}\ b)$：
+ ==溢位的原因為 $S$，$C$ 及進位的第 31 bit，三個之中有兩個以上是 1。==
- 例如：`b'01... +　b'11...`
+ ==$S$ 由於除以 2 等同右移一位的效果，第 31 bit 必定是 0。==
==可以確定若有溢位，必定為 C 的第 31 bit 及 兩者第 30 bit 發生進位。==

從上述基礎開始推導：
+ 若 $C$ 第 31 bit 為 1，則 `a`，`b` 第 31 bit 皆為 1。
+ 由於 `a`，`b` 第 31 bit 皆為 1，$S$ 的第 30 bit 為 0。為使第 30 bit 發生進位，$C$ 的第 30 bit 以及第 29 bit 的進位必須為 1。
+ 若 $C$ 第 30 bit 為 1，則 `a`，`b` 第 30 bit 皆為 1。
+ 由於 `a`，`b` 第 30 bit 皆為 1，$S$ 的第 29 bit 為 0。為使第 29 bit 發生進位，$C$ 的第 29 bit 以及第 28 bit 的進位必須為 1。
+ ... ... ... ...
+ 若 $C$ 第 2 bit 為 1，則 `a`，`b` 第 2 bit 皆為 1。
+ 由於 `a`，`b` 第 1 bit 皆為 1，$S$ 的第 1 bit 為 0。為使第 1 bit 發生進位，$C$ 的第 1 bit 以及第 0 bit 的進位必須為 1。
+ 若 $C$ 第 1 bit 為 1，則 `a`，`b` 第 1 bit 皆為 1。
+ 由於 `a`，`b` 第 1 bit 皆為 1，$S$ 的第 0 bit 為 0，不符合前述條件，**可確定 $C$ 的第 31 bit 為 1 的狀況下第 30 bit 不可能進位**。

得證**本次加法不可能造成溢位**。

因此，`#define EXP2 a & b`，`#define EXP3 a ^ b`。

### 比較三者編譯後的輸出
編譯以下程式碼對三個函式利用 Compiler Explorer 進行比較：
+ 編譯器使用 x86-64 gcc 11.2。
+ 優化選項使用 `-O3`
+ 為了使三個函式的參數傳入規則相同，第一個版本加入了將變數更改為 `a > b` 的流程，以求公平。

#include <stdint.h>
uint32_t average_1(uint32_t a, uint32_t b)
{
    if(a < b){
        a ^= b;
        b ^= a;
        a ^= b;
    }
    return b + (a - b) / 2;
}

uint32_t average_2(uint32_t a, uint32_t b)
{
    return (a >> 1) + (b >> 1) + (a & b & 1);
}

uint32_t average_3(uint32_t a, uint32_t b)
{
    return (a & b) + ((a ^ b) >> 1);
}

編譯輸出：

average_1:
        cmp     edi, esi
        jnb     .L2
        mov     eax, edi
        mov     edi, esi
        mov     esi, eax
.L2:
        sub     edi, esi
        shr     edi
        lea     eax, [rdi+rsi]
        ret
average_2:
        mov     eax, edi
        shr     eax
        mov     edx, esi
        shr     edx
        add     eax, edx
        and     edi, esi
        and     edi, 1
        add     eax, edi
        ret
average_3:
        mov     eax, edi
        xor     eax, esi
        shr     eax
        and     edi, esi
        add     eax, edi
        ret

#### XOR swap / Normal swap
實驗途中，發現不論在 `-O1`、`-O2`、`-O3` 參數下，將 `a`、`b` 交換的流程都會將原始程式碼中的 XOR Swap 更改使用一個暫存器及 `mov` 指令來進行， 而不會用三個 `xor` 指令。

假設 `a` `b` 兩者原先都是存在於記憶體中，`mov` 指令的方式會以以下方式執行，==總共四個指令，四個記憶體操作==：
+ 讀取 `a` 到 `$1` (1)
+ 讀取 `b` 到 `$2`
+ 儲存 `$2` 到 `a`
+ 儲存 `$1` 到 `b`

而假設我們用 `xor` 指令進行處理，則會以以下方式執行，==總共七個指令，四個記憶體操作==：
+ 讀取 `a` 到 `$1`
+ 讀取 `b` 到 `$2`
+ **XOR** `a`、 `b` 到 `a`
+ **XOR** `a`、 `b` 到 `b`
+ **XOR** `a`、 `b` 到 `a`
+ 儲存 `$1` 到 `a`
+ 儲存 `$2` 到 `b`

但是此處的資料並不是從記憶體中提取，而是直接對 register 做操作，
編譯成 `xor %edi, %esi; xor %esi, %edi; xor %edi, %esi;` 貌似並不會影響執行速度，為何還需要更改成 `mov` 指令的樣式呢？

由於實在是無法整理出一個理由，決定直接到 StackOverflow 提問：
[Why is the XOR swap optimized into a normal swap using the MOV instruction?](https://stackoverflow.com/questions/71382441/why-is-the-xor-swap-optimized-into-a-normal-swap-using-the-mov-instruction)

而歸納出幾個原因：
+ **對 Instruction-level parallelism(ILP) 不友善：**
- 倘若我們使用 `mov` 指令，在將 `a`(`%edi`) 之內容移動至 `temp`(`%eax`) 之後，兩個 `mov` 指令可以平行執行；而 XOR swap 的方式則需要等待依序執行。
+ **Intel x86 架構下有更好的選擇：**
- *x86 指令集中存在 `XCHG`，作用為交換兩個 register 之中的內容*，速度比 XOR swap 的三個指令還要快。即使沒有多餘的 register 可以使用也應該使用這個指令。
+ **MOV-elimination 的加速：**
- intel 在 Ivy Bridge 以後的架構內建，在執行 `mov` 指令時與其直接移動 register 內的內容，而是修正內部的 register file 直接使其指向正確的值。比起執行速度很快的 `xor`，`mov` 直接**不執行**，速度更快。

此外，和 Clang 的輸出結果比較：

average_1:                              # @average_1
        cmp     edi, esi
        mov     eax, esi
        cmovb   eax, edi
        cmovb   edi, esi
        sub     edi, eax
        shr     edi
        add     eax, edi
        ret

可以發現 CLang 利用 [Tail Duplication](http://elegantc.blogspot.com/2016/05/tail-duplication-with-llvm.html)，有效的讓不論 `a > b` 或是 `a < b` 都可以在相同的指令數及 latencty 下執行，而 GCC 在 `-O1`, `-O2`, `-O3` 等配置下並未這麼做；因此我們可以手動進行這項優化：

uint32_t average_1_taildup(uint32_t a, uint32_t b)
{
    if(a < b){
        return a + (b - a) / 2;
    }
    else{
        return b + (a - b) / 2;
    }
}

可以產出

average_1_taildup:
        cmp     edi, esi
        jnb     .L5
        sub     esi, edi
        shr     esi
        lea     eax, [rsi+rdi]
        ret
.L5:
        sub     edi, esi
        shr     edi
        lea     eax, [rdi+rsi]
        ret

犧牲些許的程式碼大小達成不論狀態皆為 5 個 uop 的結果。
:::info
根據 GCC 說明[(連結)](https://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html#:~:text=O2%2C%20-O3%2C%20-Os.-,-ftracer,-Perform%20tail%20duplication)，Tail duplication 優化需要啟用 `-ftracer`，而它並不被包含在 `-O1`、`-O2`、`-O3`、`-Og`、`-Os` 的任意一個優化選項裡。添加之後則可以產出類似 `average_1_taildup()` 的編譯結果。
:::

### 研讀 linux/include/linux/average.h

#define DECLARE_EWMA(name, _precision, _weight_rcp)			\
	struct ewma_##name {						\
		unsigned long internal;					\
	};								\
	static inline void ewma_##name##_init(struct ewma_##name *e)	\
	{								\
		BUILD_BUG_ON(!__builtin_constant_p(_precision));	\
		BUILD_BUG_ON(!__builtin_constant_p(_weight_rcp));	\
		/*							\
		 * Even if you want to feed it just 0/1 you should have	\
		 * some bits for the non-fractional part...		\
		 */							\
		BUILD_BUG_ON((_precision) > 30);			\
		BUILD_BUG_ON_NOT_POWER_OF_2(_weight_rcp);		\
		e->internal = 0;					\
	}								\
	static inline unsigned long					\
	ewma_##name##_read(struct ewma_##name *e)			\
	{								\
		BUILD_BUG_ON(!__builtin_constant_p(_precision));	\
		BUILD_BUG_ON(!__builtin_constant_p(_weight_rcp));	\
		BUILD_BUG_ON((_precision) > 30);			\
		BUILD_BUG_ON_NOT_POWER_OF_2(_weight_rcp);		\
		return e->internal >> (_precision);			\
	}								\
	static inline void ewma_##name##_add(struct ewma_##name *e,	\
					     unsigned long val)		\
	{								\
		unsigned long internal = READ_ONCE(e->internal);	\
		unsigned long weight_rcp = ilog2(_weight_rcp);		\
		unsigned long precision = _precision;			\
									\
		BUILD_BUG_ON(!__builtin_constant_p(_precision));	\
		BUILD_BUG_ON(!__builtin_constant_p(_weight_rcp));	\
		BUILD_BUG_ON((_precision) > 30);			\
		BUILD_BUG_ON_NOT_POWER_OF_2(_weight_rcp);		\
									\
		WRITE_ONCE(e->internal, internal ?			\
			(((internal << weight_rcp) - internal) +	\
				(val << precision)) >> weight_rcp :	\
			(val << precision));				\
	}

#### Exponentially weighted moving average (EWMA)

## 問題二
## 問題三
## 問題四
## 問題五
## 問題六
## 問題七：FizzBuzz
### 程式說明

#### 在不使用求餘運算下檢測數字可否被 3 或 5 整除
觀察程式碼，
可以注意到**以下的不等式用來得知一 32-bit 整數是否能被 d($d\in\{3,5\}$) 整除**：
> $S_{\mathtt{i}, 64}=\mathtt{\text{0xiiiiiiiiiiiiiiii}}, \mathtt{i}\in{[\text{0}_{16},\text{F}_{16}]}$

$n\times (S_{F,64}\div 3+1)\leq S_{F,64}\div 3$
經過化簡，可得
$n\times S_{5,64}+n\leq S_{5,64}$

接著，考慮三個不同 n 的值：整除、除三餘一、除三於二：
| $n\equiv$ | 0 (mod 3) | 1 (mod 3) | 2 (mod 3) |
| :-------: | :-------: | :-------: | :-------: |
| $n\times S_{F,64}+n$ | $\xrightarrow{\text{overflow}} n$ | $\xrightarrow{\text{overflow}}S_{5, 64}+n$ | $\xrightarrow{\text{overflow}}S_{A,64}+n$ |

利用 overflow 來取代除法運算中重複減去的過程，節省時間，
其中只有 $n\equiv 0(\text{mod 3})$ 時， 等式才成立。

:::info
此等式可以應用於檢測所有 $S_{F,64}$ 之因數 $k$ 且 $(S_{F,64}\div k) > L, L=$ **檢測範圍 bit 長度** 的
整數的因數上，舉例來說，
若要檢測是否可被 15 整除，$M_{15}$ 即可設成 $S_{1,64}+1$。
:::

**<div style="color:red">優點：$M_k$ 可提前計算，大幅減少除法運算的次數(僅一次)</div>**

#### 字串輸出優化
程式中字串輸出以『同個字串，不同起點和終點』來輸出，如以下表格：
| **字串內容** | F | i | z | z | B | u | z | z | %u | \0 |
| :-: | :-: | :-: | :-: | :-: | :-: | :-: | :-: | :-: | :-: | :-: |
| **起始點 (included)** | $0(mod 3)$ $0(mod 15)$| | | | $0 (mod 5)$ | | | | Others | |
| **終點 (excluded)** | | | | | $0(mod 3)$ | | | | $0 (mod 5)$ $0(mod 15)$ | Others |

#### 目標
由 `div5` 以及 `div3` 兩個值計算出**起點**以及**字串長**，`div5` $,$`div3` $\in\{0,1\}$
將上方表格轉換成數值：
| `div3` | `div5` | **起點** | `length` |
| :------: | :------: | :------: | :------: |
| 0 | 0 | 8 | $\geq$ 1 |
| 0 | 1 | 4 | 4 |
| 1 | 0 | 0 | 4 |
| 1 | 1 | 0 | 8 |

:::info
**起點方程式**
透過觀察，可以發現滿足以下方條件的方程式可以為解：

0. $\text{start_point}(0,0)=8$
1. **只要 `!div3` 便為0**
2. `div5` 會造成除以 2 的效果

可得 `start_point(div3, div5) = !div3 & (8 >> div5)`。(我的算法)
與題目中解答的相異之處為條件 1 的處理，題目使用的方式為
『若 `div3`，向右移 4-bit』，
無論運算值為何 (8 或 4) 皆會被消除。

#### `length` 方程式
透過觀察，可以發現滿足以下方條件的方程式可以為解：

0. $\text{length}(0,0)\geq 1$
1. `div3` 和 `div5` 任一為 1，結果為 4
2. `div3` 和 `div5` 同時成立會造成額外乘以 2 的效果

可得 `length = 2 << div3 << div5`。
:::