# 2025-03-18 問答簡記

## [程式設計不等同於軟體開發](https://www.ithome.com.tw/voice/89695)
* 「程式設計像是個人的作戰技巧，而軟體開發像是團隊行軍作戰，需要的不只是個人的戰技，諸如整體的戰略、陣勢、分工、武器……等等，也都相當的重要。軟體開發只講程式設計，就像兩軍交戰，我軍空有個人戰技，卻不談如何設定戰略、也不談該如何擺陣一樣。」
* 「在程式碼寫完之後，還需要一段測試及修改的時間，而這段時間通常不少於撰寫程式所花的時間，甚至倍數於撰寫程式所花的時間...錯估了軟體實際需要完成的時程，或是在時間截止時，只能交付品質不夠穩定的軟體。」
* 「我們應該要把軟體開發當做是一個獨立的學問來看待，而不是把它和程式設計給混在一起，才能夠把軟體開發做的更好。」

## 佳句偶得
* 「資訊人的本色就是做什麼，就要像什麼」 —— 洪良茂，成功大學資訊系大學部第一屆畢業生。
* 「大部分的人一輩子洞察力不彰，原因之一是怕講錯被笑。想了一點點就不敢繼續也沒記錄或分享，時間都花在讀書查資料看別人怎麼想。看完就真的沒有自己的洞察了」 —— 蔡志浩
* [能自己從程式裡找到問題，然後想辦法解決，才叫作會寫程式。僅僅實作給定的指示，不算會寫程式](https://x.com/yungyuc/status/1899457408453443923)
* [如果寫程式變得輕鬆，可以想像它的供給會變得充沛。愈能輕鬆寫的程式，就愈難獲得利潤](https://x.com/yungyuc/status/1901608433498959965)

[海底潛水遇難 沒有氧氣 他閉氣整整 38 分鐘](https://youtu.be/tSZtSpv7IjE)
> 人在極端環境下，被激發的潛能...就像某種生物的保險、一種進化的[冗餘](https://en.wikipedia.org/wiki/Redundancy)，讓我們未來面臨整體氣候和環境變化時，依然能夠延續我們的文明

[科技產業的好時光已經結束](https://www.facebook.com/PivotSoftwareEngineer/posts/pfbid0eX6zWD4t68yJPsJfCeV1qaYXzQA3YqyZ7MwUVKBQGfguSPawyxeF4mz7fBuUmV81l)

- 「我本來覺得我長很醜，直到我看到 BLACKPINK 才發現，歐! 原來大家也很看重才華。」

## 快訊
* [GCC 15 預計在今年第二季發布，Rust 程式語言的支援獲得顯著改進](https://www.facebook.com/groups/900640625606073/?multi_permalinks=946578174345651)
> [Polonius](https://github.com/rust-lang/polonius) borrow checker，其命名來自莎士比亞作品《哈姆雷特》的名言 "[Neither a borrower nor a lender be](https://literarydevices.net/neither-a-borrower-nor-a-lender-be/)"，發音: [pəˈloʊniəs](https://www.collinsdictionary.com/dictionary/english/polonius)
* [《程式設計原來不只有寫 CODE！》](https://www.facebook.com/groups/system.software2025/posts/946440737692728/): 恭喜學生成為暢銷電腦書籍作家
* [SoftBank 斥資 65 億美元收購 Ampere Computing](https://www.facebook.com/groups/system.software2025/posts/946388827697919/)
* [GIMP 3.0 正式發布](https://www.facebook.com/groups/system.software2025/posts/946327397704062/)
* [GTC 裡頭的先進 Arm 多核技術](https://www.facebook.com/groups/system.software2025/posts/945777261092409/) $\to$ [無所不在的 Arm](https://hackmd.io/@sysprog/arm-intro)
* [Andy 老師，CROWD 還給你](https://www.facebook.com/groups/system.software2025/posts/944955474507921/)
* [RISC-V 系統模擬器藉由 VirtIO GPU，將 Guest OS 的影像顯示於 Host 端](https://www.facebook.com/groups/system.software2025/posts/943783021291833/)
* [陳立武正式接任 Intel 新任執行長](https://www.facebook.com/groups/system.software2025/posts/941084441561691/)
* [SSE2NEON 專案列於 Apple 公司 Games Pathway 專頁](https://www.facebook.com/groups/system.software2025/posts/939306688406133/): 作為移轉既有 x86 程式碼到 Arm 為基礎的 Apple Silicon 的 SSE 指令集轉換工具

## Standardizing vs. Normalization
中學統計課程提過的「資料標準化」，對應到大學課程的 Z-Score，又稱 Standard Score（標準分數），這類方法稱為 standardizing, normalizing, normalization。其中 [normalization](https://en.wikipedia.org/wiki/Normalization_(statistics)) 是統計學的術語，指將不同比例或尺度的資料映射至統一尺度，便於進行比較。除了 Z-Score，另一常見的 normalization 方法是 Min-max feature scaling。

在資料分析領域，如統計學、巨量資料分析、機器學習及深度學習中，Normalization 指在分析前對資料進行預處理的步驟，又稱「資料前置處理」。然而，許多人 normalization 譯作「正規化」，或許是強調該方法意在「統一各種資料規格」，但也帶來許多混淆，因為「正規化」一詞在不同領域有不同寓意，例如在機器學習中，「正規化」可能是針對成本函式 (cost function) 的某種調整手法，而在資料庫領域，「正規化」則指對資料表進行結構化處理的過程；正規語言與自動機理論中，不但有「正規語言」的概念，還有名為「正規化」的演算法步驟。

知名人工智慧學者何愷明與 Yann LeCun 聯手，在 [Transformers without Normalization](https://jiachenzhu.github.io/DyT/) 研究中，提出僅以 9 行 Python 程式碼，成功移除 Transformer 架構中被視為不可或缺的 LayerNorm，他們提出名為「動態 Tanh」(DyT) 的新方法，徹底移除傳統 Transformer 中的正規化層，不僅簡化模型設計，亦有效減少計算資源消耗。多項任務驗證表明，移除正規化層後模型性能不減反升，顛覆 normalization 層必要性的傳統認知。
> 
> tanh 的英文名稱是 hyperbolic tangent，IPA 通常標記為 [tæn] 或 [tɑːn]，也有人發音為 [tæntʃ]，但較常讀作 tan

在 Llama 模型測試中，DyT 大幅減少推理與訓練時間，整體計算量降低約 7～8%。對於大規模模型而言，這樣的減少幅度相當可觀，尤其是在深度學習訓練成本日益提高的背景下，具備明顯實用價值。DyT 之所以能有效移除 Transformer 中的 normalization 層並提升效能，關鍵在於其解決深度學習常見的梯度消失與分布不穩定問題。在傳統 Transformer 架構中，LayerNorm 負責穩定輸入分佈與加速收斂，但本質上是一種「被動補救」，強制將模型輸出維持在標準化範圍，卻未解決輸出不穩定的根源。DyT 的創新在於使用 Dynamic Tanh（動態雙曲正切函數）取代正規化層，直接在激勵函數層面解決數值分佈問題。Tanh 本就具備內建 normalization 功能，輸出範圍固定在 $[-1, 1]$，而 DyT 加入可學習參數 $\alpha$，實作 DyT(x) = tanh(αx)，使其能根據不同層特性，調整輸出壓縮程度，穩定數值，無需 LayerNorm 修正。

除了 tanh 函式具備的 S 形壓縮特性，DyT 中的 α 還展現出一項巧妙行為 —— 它趨近於輸入張量標準差的倒數（即 α ≈ 1/std）。研究發現，在訓練過程中 α 可自動調整至接近該值，某種程度上模擬 LayerNorm 的標準化效果。這也解釋為何 DyT 即便移除 LayerNorm，依然能穩定輸出並提升效能。LayerNorm 之所以廣泛應用，部分原因是為了解決多層網路中梯度傳遞易發生離散或過度壓縮的問題。然而，這種「強制標準化」可能導致資訊（尤其是雜訊）的流失。早在 2017 年，Edward Raff 就證明過，過度套用正規化層未必總是最佳選擇。以多層 CNN 為例，ReLU 後套 normalization 層可能會削弱梯度收斂能力，影響模型學習效果。現代許多架構早已開始移除部分 LayerNorm，以減少冗餘計算與潛在副作用。更有趣的是，LeCun 早已指出，自然語言處理（NLP）問題不一定非得依賴 Transformer 架構來解決。從 DyT 的推出可見，Meta 正在布局更具延展性的模型設計策略。本研究已成功入選 CVPR 2025，並開放原始碼。DyT 在計算效率上表現尤為突出，以 Llama 7B 為例，DyT 明顯縮短推論與訓練所需時間，BF16 與 FP32 精度均呈現同樣趨勢。侷限性方面，DyT 在非 Transformer 架構（如 ResNet 中 BatchNorm）表現不佳。
> [類神經網路的 ReLU 及其常數時間實作](https://hackmd.io/@sysprog/constant-time-relu)

## 親身實踐
:::success
「紙上得來終覺淺，絕知此事要躬行」，出自《冬夜讀書示子聿》(1199 年)，作者陸游
:::

間接經驗是從書本中汲取養分，學習前人的知識的途徑，而直接經驗是親身體認到的知識，是獲取知識更重要的途徑。唯有藉由「躬行」，把書本知識變成實際知識，方可深入理解其中的道理。

從 Google 搜尋的結果來看，我們可能會看到有人提及 Linux 核心的 context switch 成本約 1 microseconds，但也有不同的數字，到底該相信哪個說法呢？其實，就該設計實驗，並且著手在特定的硬體及軟體組態下進行測量。

當然，你也許會質疑：「我們重複他人已經做過的實驗，到底有何意義？」回答此問題前，先來看諾貝爾物理學獎得主李政道博士，在美國受費米指導時期 (1946 年到 1950 年) 發生的一則軼事：費米問李政道太陽中心溫度為何，後者一開始說約一千萬度左右，費米追問是否核算過，李政道以二個關聯方程來回答，大致呼應一開始的數據，不過費米表示「你不能依靠别人的計算结果，你必須自己核算」，並建議透過計算尺查證，之後師生就一同打造木製計算尺，從而計算出符合推論的數據。

在費米指導李政道之前，費米在第二次世界大戰期間參與曼哈頓計劃。費米領導他的團隊設計並建造[芝加哥 1 號堆](https://en.wikipedia.org/wiki/Chicago_Pile-1) (位於芝加哥大學體育場 Stagg Field 的下方)。這個反應爐 1942 年 12 月 2 日進行臨界試驗，完成首次人工自持續連鎖反應，奠定日後位於洛斯阿拉莫斯國家實驗室、利用熱核反應的「超級核彈」的成功。1945 年 7 月 16 日，費米參與三位一體核試，並利用自己的方法估算爆炸當量。之後的故事大家都知道，原子彈永久地改變人類的歷史。

已是物理泰斗的費米，面對一位從中國上海留學的學生，竟不忘以身教提醒李政道，面對統計物理的相變以及凝聚體物理學的極化子的研究之前，追求真善美境界並親身體驗的態度，才是一切的根本。
> [李政道師從費米時期的一段故事](https://www.facebook.com/jay.hu.376/posts/pfbid08TTUPMXb6jTRsPhZiaTiu4z4LjtXqtZLxEjTzexzMWekMHjfjsJmgK1chKmUAQbZl)

大部份的學生很快就可「背誦」作業系統課程給予的「心法」，回答一長串，可是，其中多少人「看過」Process 和 Thread 呢？多少人知道二者的實作機制？又，為何當初電腦科學家和工程師會提出這二個概念呢？

書本說 thread 建立比 process 快，但你知道快多少？是不是每次都會快？然後這兩者的 context switch 成本為何？又，在 SMP 中，是否會得到一致的行為呢？

我們希望透過一系列的實驗，藉由統計來「看到」process 與 thread。物理學家如何「看到」微觀的世界呢？當然不是透過顯微鏡，因為整個尺度已太小。統計物理學 (statistical physics) 指的是根據物質微觀夠以及微觀粒子相互作用的認知，藉由統計的方法，對大量粒子組成的物理特性和巨觀規律，做出微觀解釋的理論物理分支。如今，當我們想要「看到」context switch, interrupt latency, jitter, ... 等微觀事件，無一不透過統計學！
> [2015 年報告: ARM-Linux](https://wiki.csie.ncku.edu.tw/embedded/arm-linux)
> [財經首長也該上的物理課](https://pansci.asia/archives/81114)

## idoleat
> [quiz1+2](https://hackmd.io/@idoleat/linux2024-homework2)

評：NVIDIA 考前衝刺

排序是一維資料附加序列位置資訊、擴增為二維，再投影至另一維之展現。

## jouae
> Pull request [#173: Improve string randomness](https://github.com/sysprog21/lab0-c/pull/173)

在 [lab0-c](https://github.com/sysprog21/lab0-c) 的 `qtest.c` 提供對於隨機字串插入佇列的功能，藉由 `RAND` 選項啟用偽隨機字串產生器。使用方式如 `ih RAND 5` ，執行 `q_insert_head` 對佇列開端插入隨機 $5$ 個字串。

在 `qtest` 中隨機字串產生的函式為 `fill_rand_string` ，對全域變數 `charset` 由 26 個小寫英文字母的組成陣列，隨機挑選介於 `MIN_RANDSTR_LEN` 和 `MAX_RANDSTR_LEN` 個字元，最後組成字串。

static void fill_rand_string(char *buf, size_t buf_size)
{
    size_t len = 0;
    while (len < MIN_RANDSTR_LEN)
        len = rand() % buf_size;

    randombytes((uint8_t *) buf, len);
    for (size_t n = 0; n < len; n++)
        buf[n] = charset[buf[n] % (sizeof(charset) - 1)];
    buf[len] = '\0';
}

其中偽隨機數的產生藉由 [randombytes](https://github.com/dsprenkels/randombytes) 實作。

而 `wuyihung` 發現在此採用 `uint8_t` 會有問題：
> However, since the number of variations is 256, which is not exact division of 26

`uint8_t` 在 [stdint.h(0p) — Linux manual page](https://man7.org/linux/man-pages/man0/stdint.h.0p.html) 提及：
> The typedef name `intN_t` designates a signed integer type with width N, no padding bits, and a two's-complement representation.

藉由此敘述得知非符號整數，以十進位表示法範圍為： $[0,255]$ 。`buf[n] % (sizeof(charset) - 1)` 從數學角度解釋，此處是對於 $x\equiv i \pmod{26}$ 選定索引的過程。

藉由[除法原理](https://en.wikipedia.org/wiki/Euclidean_division)可以得到
$$
255 = 26 \times 9 + 21.
$$

亦即，0 到 25 數字的倍數在 0 至 255 中，除了 22, 23, 24, 25 倍數之外皆出現 10 次。

**假設 `randombytes` 給定的隨機數是均勻且彼此獨立的**，在此之下 0 到 21 數字倍數出現的頻率為 10， 22, 23, 24, 25 倍數為 9 。假設 0 到 25 對應到小寫英文 a 到 z ，其離散機率分佈表示為
$$
P(X=i) =
\begin{cases}
\dfrac{10}{256}, \text{ for } i=0,1,\dots,21 ,\\
\dfrac{9}{256}, \text{ for } i=22,23,24,25.
\end{cases}
$$

也就是說，隨機產生的字元頻率不均勻。

`wuyihung` 使用 Python 設計一個[測試腳本](https://github.com/wuyihung/lab0-c/blob/master/scripts/test_prng.py)驗證上述理論上的情況。

他藉由重複執行 `ih RAND 30` $150,000/30$ 次後，將串列內所有的字串連接（ concatenate ）成一個字串，並使用 [ent](https://www.fourmilab.ch/random/) 工具分析字串的 entropy 。發現理論跟實際有些許誤差（約 $0.0009$ ）。

-    randombytes((uint8_t *) buf, len);
+    uint64_t *randstr_buf_64 = malloc(len * sizeof(uint64_t));
+    randombytes((uint8_t *) randstr_buf_64, len * sizeof(uint64_t));
    for (size_t n = 0; n < len; n++)
-        buf[n] = charset[buf[n] % (sizeof(charset) - 1)];
+        buf[n] = charset[randstr_buf_64[n] % (sizeof(charset) - 1)];

+    free(randstr_buf_64);

他改用 `uint64_t` 來改善隨機字串的產生過程。其想法很簡單，就是藉由把上述 8 bits 改成用 64 bits 表示，在上述分佈的分母由 $2^8$ 改為 $2^{64}$ 從除法原理得知：
$$
2^{64} = 709490156681136600 \times 26 + 16
$$

其離散機率分佈表示為
$$
P(X=i) =
\begin{cases}
\dfrac{709490156681136601}{2^{64}}, \text{ for } i=0,1,\dots,16 ,\\
\dfrac{709490156681136600}{2^{64}}, \text{ for } i=17,18,\dots,24,25.
\end{cases}
$$

參照他的實驗方式，此處將使用 `uint8_t` 和 `uint64_t` 的結果用 Python 套件 matplotlib 將數據以長條圖（bar plot）呈現。

採用 `uint8_t`：
總共 1050988 個字元。從 a 到 z 的頻率：
`[41080, 41161, 40861, 41197, 41150, 40723, 40969, 41191, 41006, 41099, 41278, 41312, 41146, 40890, 41113, 40846, 36915, 36894, 41141, 41231, 41035, 41269, 40993, 40905, 36696, 36887]`

採用 `uint64_t`：
總共 1049822 個字元。從 a 到 z 的頻率：
`[40629, 40651, 40560, 40750, 40237, 40244, 39793, 40428, 40261, 40372, 40540, 40571, 40310, 40795, 40222, 39952, 40452, 40058, 40700, 40048, 40713, 40398, 40332, 40431, 40047, 40328]`

`sizeof()` 在 C99 6.5.3.4 The sizeof operator 描述中也是回傳 `unsigned integer type`
> The value of the result is implementation-defined, and its type (an unsigned integer type) is size_t, defined in <stddef.h> (and other headers).

由於 $charset$ 大小固定，以 `26U` `取代 sizeof(charset) - 1`

採用 `uint8_t` 並且使用 `26U` 取代 `sizeof(charset) - 1`取模：

-        buf[n] = charset[buf[n] % (sizeof(charset) - 1)];
+        buf[n] = charset[buf[n] % 26U];

總共 1050391 個字元。從 a 到 z 的頻率：
`[40987, 41617, 40995, 41252, 41100, 40988, 40903, 41212, 40949, 40610, 40903, 41089, 41100, 40931, 41208, 40949, 40958, 40686, 41270, 41415, 40921, 40946, 37085, 36880, 36761, 36676]`

## EricccTaiwan
> [第五周測驗](https://hackmd.io/@sysprog/linux2025-quiz5)
> 03/24 21:30 線上討論
> 04/03 20:00 直播討論 (補筆記)
> 04/07 15:00 線上討論
> [Quiz5-refactor](https://gist.github.com/EricccTaiwan/acf4ead098eed57a90016d82eadea3a6) (還在修改)
> 待完成 : 改程式，針對小數點用其他表示(e.g., Q8.23)
> [google fast_tanh()](https://github.com/google/audio-to-tactile/blob/d3f449fdfd8cfe4a845d0ae244fce2a0bca34a15/src/dsp/fast_fun.h#L151)

首先，我們先針對 [Q notation](https://en.wikipedia.org/wiki/Q_(number_format)#) 的表示方式進行定義，並以此作為下方段落的基礎說明。在維基百科中提到， Q notation 主要有兩種常見的表示方式，其用途與意義雖然相同，但寫法略有不同，分別為 Texas Instruments (以下簡稱 TI ) 版本和 ARM 版本。

根據 TI 的定義， `Q` 後面會接著一組的 `m.n` ，其中 `m` 表示用於整數部分的位元數，而 `n` 則是小數部分的位元數。值得注意的是 sign bit 是不被計算在 `m` 裡面的，因此 `Qm.n` 所使用的總位元數為 $1+m+n$。因此在 [quiz5](https://hackmd.io/@sysprog/linux2025-quiz5#%E6%B8%AC%E9%A9%97-1) 中，所定義的縮放因子 `FIX16_ONE = 16^4`，任意數字經過定點化後，其對應的 TI 版本 Q notation 表示為 `Q15.16`。相對地，在 ARM 的表示法中，`m` 則包含了 sign bit，因此同樣的縮放因子 `FIX16_ONE = 16^4`，若使用 ARM 版本表達，應表示為 `Q16.16`。
**Texas Instruments version**
* Qm.n : 表示該數有 m 個整數位、n 個小數位，是 $1+m+n$ 位二補數。
> The Q notation, as defined by Texas Instruments,[1] consists of the letter Q followed by a pair of numbers m.n, where m is the number of bits used for the integer part of the value, and n is the number of fraction bits.
> The first bit always gives the sign of the value(1 = negative, 0 = non-negative), and it is not counted in the m parameter. Thus, the total number w of bits used is 1 + m + n.

**ARM version**
* Qm.n : 表示該數有 m-1 個整數位、n 個小數位，是 $m+n$ 位二補數。
> In this variant, the m number includes the sign bit. For example, a 16-bit signed integer would be denoted Q15.0 in the TI variant, but Q16.0 in the ARM variant.[2][3]

下方採用 ARM 版本進行課程討論，因為 ARM 營收比較高，所以我聽他的。

### Q1:
怎麼解釋輸出的小數點以下只有三位是精確的？

A1:
這其實是誤差累積的結果，而非固定限制。雖然從結果看起來小數點後只有三位有效數，但實際上使用的 `Q16.16` 定點數理想上可以表達到小數點後第 5 位左右的精度。

`Q16.16` 可以表示的範圍為 $-2^{15}$ 到 $2^{15} - 2^{-16}$。
* 先討論上界 $2^{15} - 2^{-16}$，為了將整數轉換為 `Q16.16` 格式，會需要乘上 `FIX16_ONE`，也就是左移 16 位元，如果傳入數值為 $2^{16}$，代表 bit 15 為 1 ，其他為 0 ，經過左移 16 位後，1-bit 會被推到 bit 31 上，而在 `int32_t` 中 bit 31 是用來表示符號的 (sign bit) ，導致 $2^{16}$ 經過轉換後會變成 $-0$ ， 造成錯誤，代表 $2^{16}$ 已經超出 `Q16.16` 可表示的範圍。
* 所以對於下界來說也是一樣，經過左移後，若超出第 31 個位元，就會 overflow。可以注意到，在 `fix16_exp()` 函式中，先檢查是否為負數，若是，先將其轉為正數進行運算，最後再以 $e^{-x} = 1 / e^x$ 計算結果，避免左移 16 位元後，還需要對 signed bit 進行附值的操作。

從原始碼中的這一行可以觀察到：

#define FIX16_ONE 0x00010000 // scaling factor

這個巨集定義了 `FIX16_ONE = 2^16 = 65536` ， 代表我們使用的是 `Q16.16` 格式的定點數表示法。這種格式的最小可表示單位是 $2^{-16} \approx 1.53 \times 10^{-5}$，也就是理論上最多可以精確到小數點後第 5 位。

然而經過實驗，發現 `fix16_tanh()` 的結果通常只精確到小數點後三位，這是由以下幾個原因造成的：
1. 定點數本身的限制與範圍問題 (已在上面解釋)
2. 相近數值相減與除法造成的精度放大 (`fix16_tanh` 中 `fix16_div` 的分子分母項)
3. 泰勒展開無法準確的逼近 (並思考是否有更快速收斂的方式)

這邊針對第 2 點和第 3 點進行討論，首先針對第 2 點， `fix16_tanh` 的函式中，會回傳 `fix16_div` ， 而這個函式的分子分母數值太接近，進而導致再小的誤差，也會被放大，至於誤差會放大多少，目前不太清楚。

/* tanh(x) is defined as (exp(x) - exp(-x)) / (exp(x) + exp(-x)). */
fix16_t fix16_tanh(fix16_t in)
{
    fix16_t e_x = fix16_exp(in);
    fix16_t m_e_x = fix16_exp(-in);
    return fix16_div(e_x - m_e_x, e_x + m_e_x);
}

那接著就要討論為何在分子分母項會有誤差產生， 可以看到 `fix16_exp` 的函式是對經過縮放後的數字取自然指數 ($e$) ， 前面的 4 個 `if` 是在對 edge case 進行處理，

/* 傳入 0 , 回傳 exp(0) = FIX16_ONE */
if (in == 0)
    return FIX16_ONE;

/* 傳入 FIX16_ONE , 回傳 exp(FIX16_ONE) = FIX16_ONE * exp(1) */
if (in == FIX16_ONE)
    return 178145;

/* 理想上 Q16.16 能最大表示的數字為 2^15 - 2^-16 ,
 * 這邊的 681391 只代表一個相對大的數， 傳入相對大的數，
 * 回傳 16 位元能表示的最大值 2147483647 (0b0111111111111111) */
if (in >= 681391)
    return 0x7FFFFFFF;

/* 理想上 Q16.16 能最小表示的數字為 -2^15 ,
 * 這邊的  -772243 只代表一個相對小的數， 傳入相對小的數，
 * 回傳 0 */ 
if (in <= -772243)
    return 0; // why not min value

首先，先對於泰勒展開進行介紹，為了逼近以 $e$ 為底數的指數函數，因此對 $e^x$ 做泰勒展開，可以得到，

$$
e^x = \sum_{n=0}^{\infty} \frac{x^n}{n!} = 1 + x + \frac{x^2}{2!} + \frac{x^3}{3!} + \cdots + \frac{x^n}{n!} + \cdots \quad \forall x \quad \text{（對所有 } x \text{ 都成立）}
$$

接著看到下方的 for 迴圈，`i=2` 開始，其實就是在算 $x\times\frac{x}{2!}$ 這一項，然後再加回去 result 中，持續下去，但會有如此誤差傳遞的原因有兩點:
1. `Q16.16` 只能表示到約 $1.53 \times 10^{-5}$，更小的數會被 truncate 成 0
2. 泰勒展開無法準確逼近真實值，收斂太慢

fix16_t result = in + FIX16_ONE; // (1 + x) * FIX16_ONE
fix16_t term = in;
for (uint_fast8_t i = 2; i < 30; i++) {
    term = fix16_mul(term, fix16_div(in, int_to_fix16(i)));
    result += term;
    /* Break early if the term is sufficiently small */
    if ((term < 500) && ((i > 15) || (term < 20)))
        break;
}

### Q2:
承上，`fix16_tanh()` 如何更精確的表示小數點? (可以表示到小數點第四位嘛? 還是沒有辦法更準確)

A2:
由於 `Q16.16` 格式先天上的精度限制，對於泰勒展開式中的高次項，若仍使用 `Q16.16` 表示，將可能因精度不足而導致誤差累積。因此，針對這些數值較小的項改用 Q1.31 格式來表示。由於 Q1.31 將 31 位元用於小數部分，其最小可表示單位為 $2^{-31}\approx 4.6 \times 10^{-10}$，能更精確地表達小數點後極小的數值。

此外，也可進一步思考是否存在比泰勒展開更適合的近似方法，以逼近 $e^x$ 的值。若能以更高精度近似 $e^x$，就能漸少在 `fix16_div` 中因前段誤差放大所導致的誤差傳遞問題。

參照 [efficient tanh computation using Lambert's continued fraction](https://varietyofsound.wordpress.com/2011/02/14/efficient-tanh-computation-using-lamberts-continued-fraction/) 和 [Rapid approximation of tanh(x) [jenkas]](https://math.stackexchange.com/a/3485944/1399226)

考慮以下實作:
> [fast_tanh 實作疑惑](https://hackmd.io/@ericisgood/fast_tanh_problem)

float fast_tanh(float v)
{
    /* Constants for polynomial approximation */
    const float c1 = 0.03138777f;
    const float c2 = 0.276281267f;
    const float log2_e = 1.442695022f; /* log2(e) */

    /* Scale input by log2(e) for exponentiation base 2 */
    v *= log2_e;

    /* Separate integer and fractional parts */
    int int_part = (int) v;
    float frac_part = v - int_part;

    /* Compute polynomial approximation terms */
    float frac_part_sq = frac_part * frac_part;

    float exp_frac_even_term = log2_e + c2 * frac_part_sq;
    float exp_frac_odd_term = frac_part + (c1 * frac_part_sq * frac_part);

    float exp_approx = exp_frac_odd_term + exp_frac_even_term;

    /* Efficiently compute scaling by 2^(int_part) */
    union {
        float f;    // holds the actual floating-point value
        uint32_t i; // used for manipulating the IEEE 754 bit representation
    } exp_scaled = {.f = exp_approx};
    exp_scaled.i += (uint32_t) (int_part << 23);

    float exp_neg_approx = exp_frac_even_term - exp_frac_odd_term;

    /* Combine terms to approximate tanh(v) */
    return (exp_scaled.f - exp_neg_approx) / (exp_scaled.f + exp_neg_approx);
}

$\tanh(x)$ 以指數函數定義如下：
$$
\tanh(x) = \frac{e^{x} - e^{-x}}{e^{x} + e^{-x}}
= \frac{e^{2x} - 1}{e^{2x} + 1}
$$

近似法本質上是在高效率地估算指數函數 ($e^x$)。

首先將指數函數的底數從 $e$ 轉成以 2 為底：
$$
e^x = 2^{x \log_2 e}
$$

輸入值 `v` 藉由常數 `log2_e` (約為 $1.442695$) 進行：

v *= log2_e;

隨後將 `v` 拆分成整數部分 (`int_part`) 與小數部分 (`frac_part`)，便於計算指數函數：
$$
2^{\text{int_part} + \text{frac_part}}
= 2^{\text{int_part}} \times 2^{\text{frac_part}}
$$

針對小數指數部分 ($2^{\text{frac_part}}$)，以多項式進行近似：
$$
2^x \approx 1 + (\ln 2)x + \frac{(\ln 2)^2}{2}x^2 + \frac{(\ln 2)^3}{6}x^3
$$

程式給定的常數 `c1` 和 `c2` 以 [Minimax approximation algorithm](https://en.wikipedia.org/wiki/Minimax_approximation_algorithm) 近似法所得，能用最少的運算獲得足夠高的精確，專門在範圍 $[0, 1)$ 的近似：

IEEE 754 浮點數表示法的形式如下：
$$
(-1)^{\text{sign}} \times 2^{(\text{exponent} - \text{bias})} \times \text{mantissa}
$$

透過直接調整指數位元，可迅速地達成 $2^{\text{int_part}}$ 運算，免去昂貴的浮點數乘法：

exp_scaled.i += (uint32_t)(int_part << 23);

此處左移 23 位元是因為 IEEE 754 單精度浮點數由 1 位元符號、8 位元指數與 23 位元尾數所組成。

此近似法可計算 $e^x$ 與 $e^{-x}$。

#### **改進 `fast_tanh()` 對於 $e^{-x}$ 的處理**
Jenkas 所提出的方法僅適用於經過轉換後 (`v*=log2_e`) 且不含整數部分的數值。這是因為在原始程式中，對於 $e^{-x}$ 的計算僅處理了 $2^{\text{frac_part}}$，並未考慮整數部分對應的縮放因子 $2^{\text{int_part}}$。因此，在 $(0, 1]$ 的範圍內，實際上僅有 $(0, \frac{1}{\log_2 e}]$ 的區間是有效的。

針對這個限制，我對程式中 $e^{-x}$ 的計算邏輯進行了修改，使其能正確處理整體的指數縮放，並拓展了適用範圍。

+int int_part_neg = (int)(-v);
float exp_neg_approx = exp_frac_even_term - exp_frac_odd_term;
+union {
+    float f;
+    uint32_t i;
+} exp_neg = {.f = exp_neg_approx};
+exp_neg.i += (uint32_t)(int_part_neg << 23);

/* Combine terms to approximate tanh(v) */
-return (exp_scaled.f + exp_neg_approx) / (exp_scaled.f - exp_neg_approx);
+return (exp_scaled.f - exp_neg.i) / (exp_scaled.f + exp_neg.i);

> [fast_tanh 修正版](https://gist.github.com/EricccTaiwan/c61e3273d7b7de18398780c1ca081d94)

$ ./jenkas # 原始 
fast_tanh(0.70) = 0.339411
$ ./modify # 修改後
fast_tanh(0.70) = 0.604368

與實際值 $\tanh(0.7) = 0.604367777$ 相比，因為已經超出 jenkas 原程式碼所能處理的有效範圍，因此誤差約 $43\%$，經過處理後的誤差小到可以忽略。

#### 實作 Minimax approximation (求 `c1` 和 `c2`)
以下是用 MATLAB 的 [CVX](https://ask.cvxr.com/) 工具對 $2^x$ 進行 Minimax Approximation，目標函數為
$$
\min(\max_{x\in[0,1)} \left| 2^x - \Bigl(1+\ln2\,x+c_2x^2+c_1x^3\Bigr) \right| )
$$

clc;clear;
cvx_begin quiet
    variables c1 c2 t
    x = linspace(0,0.99999,100)';  % 選取100個點，確保不含1
    f = 2.^x;                     % 目標函數是 2^x
    ln2 = log(2);                 % ln(2) 的精確值

    minimize(t)
    subject to
        abs( f - (1 + ln2 * x + c2 * (x.^2) + c1 * (x.^3)) ) <= t;
cvx_end

% 輸出結果
fprintf('最佳的 c1 = %.8f\n', c1);
fprintf('最佳的 c2 = %.8f\n', c2);
fprintf('最大誤差 t = %.8f\n', t);

輸出結果:

最佳的 c1 = 0.07343435
最佳的 c2 = 0.23317052
最大誤差 t = 0.00024788

接著我對於 [jenkas](https://math.stackexchange.com/a/3485944/1399226) 所提出的 `c1` 和 `c2` 計算其最大誤差，

clc;clear;

% 使用初始設定的 c1 和 c2 常數值
c1 = 0.03138777;
c2 = 0.276281267;

x = linspace(0, 0.99999, 100)';  % 選取 100 個點，確保不含 1
f = 2.^x;                        % 目標函數是 2^x
ln2 = log(2);                    % ln(2) 的精確值

% 計算使用初始常數值時的誤差
error = abs(f - (1 + ln2 * x + c2 * (x.^2) + c1 * (x.^3)));

% 找到最大誤差
t = max(error);

% 輸出結果
fprintf('使用的 c1 = %.8f\n', c1);
fprintf('使用的 c2 = %.8f\n', c2);
fprintf('最大誤差 t = %.8f\n', t);

輸出結果 :

使用的 c1 = 0.03138777
使用的 c2 = 0.27628127
最大誤差 t = 0.00684144

接著，用這兩組的 `c1` 和 `c2` 各自近似 $e^x$ ，以 $x = 0.5$ 為例
$$
error = |2^{0.5} - (1+ln2 \times 0.5 + c2 \times 0.5^2 + c1 \times 0.5^3)|
$$

| 參數來源 | $c_1$ | $c_2$ | error | `fast_tanh(0.5)` | 與真實值誤差 |
|----------|-----------------|-----------------|-----------------------------|-----------------|-----------------------------------------|
| CVX | 0.07343435 | 0.23317052 | $+1.6805\times10^{-4}$ | 0.478836 | 約 +0.01672 |
| jenkas | 0.03138777 | 0.27628127 | -0.0054 | 0.462117 | 幾乎為 0 |
> 真實值 $\tanh(0.5) = 0.46211715726$

~~這一現象顯示，僅針對 $2^x$ 的局部近似最小化並不能直接保證整個 `fast_tanh()` 演算法達到全局最佳，因此這個問題可能還需要進一步研究以尋找更合適的全局最佳係數。~~

這代表著，其實我們是要對 `tanh(x)` 去做泰勒展開的逼近，而不是去逼近 $e^x$ 。
> 補 : 貝茲曲線
### Q3:
為何 `typedef int32_t fix16_t` 改成 `typedef int16_t fix16_t` 輸出會錯誤?

A3:
在 `Q16.16` 固定點格式中，`fix16_t` 使用 32 位元整數來表示數值，其中高 16 位表示整數部分，低 16 位表示小數部分。巨集 `FIX16_ONE` 被定義為 `0x00010000`，代表實數值 1（即 $2^{16} = 65536$ 的定點表示）。

若將 `fix16_t` 改為 `int16_t`，即只使用 16 位元整數來儲存定點數，將無法容納 `FIX16_ONE` 所代表的數值。因為 `int16_t` 的範圍僅為 $-2^{15}$ 到 $2^{15} - 1$（即 -32768 到 32767），而 `0x00010000` 等於十進位的 65536，遠超出此範圍，導致 overflow，最終值會錯誤甚至變成 0。

因此，`Q16.16` 格式必須使用 32 位元整數（`int32_t`）才能正確表示與運算，否則會因為容量不足而產生錯誤結果。

## rota1001
$\displaystyle \lim_ {n\to\infty} (1+\frac{1}{n})^{n}$ 為何會收斂在 $e$ ?

證明 $e^x = \displaystyle \lim_ {n\to\infty} (1+\frac{1}{n})^{nx}$

首先從定義 $\ln x = \int_1^x \frac{1}{t}dt$ 出發，以下會在一些定理已經證明的基礎上開始敘述，所以從頭開始講一些結論：
1. $\exists!e \in \mathbb{R}\ [ \ln x = \log_ex,\ \forall x\in(0, \infty)]$
我們知道藉由對數律可以等價的定義出一個對數函數，所以只要證明 $\ln(x)$ 這個函數滿足對數律就好，而且有個實數 $e$ 使得 $\ln(e)=1$。

2. 定義 $e$ 就是在第 1 點找到的那個常數。然後因為 $\ln(x)$ 在 $(0, \infty)$ 上是 one to one 且 onto 的，所以它存在一個從 $\mathbb{R}$ 映射到 $(0, \infty)$ 的反函數 $E(x)$，接下來可以證明 $E(x)=e^x$。

準備都做完了，接下來可以進行證明：
$\text{Let}\ t=\frac{1}{n}$

$\ln(\displaystyle \lim_{n\to\infty}(1+\frac{1}{n})^n)$

$=\displaystyle \lim_{n\to\infty}\ln((1+\frac{1}{n})^n)$ （因為 $\ln(x)$ 是連續函數，從微積分第一基本定理就可以知道）

$=\displaystyle \lim_{t\to 0} \ln((1+t)^\frac{1}{t})$

$=\displaystyle \lim_{t\to 0}\frac{\ln(1+t)-\ln(1)}{t}$

$=\ln'(1)=\frac{1}{1}=1$

因為 $\ln(x)$ 在 $(0, \infty)$ 中是 one to one 的，所以 $\displaystyle \lim_{n\to\infty}(1+\frac{1}{n})^n=e$，而明顯的 $e^x = \displaystyle \lim_ {n\to\infty} (1+\frac{1}{n})^{nx}$

接下來是下一個證明， $(e^x)'=e^x$ （這裡先默認連鎖律已經證明過了）：
$(\ln(e^x))'=\ln'(e^x)(e^x)'=e^{-x}(e^x)'$

$=x'=1$

$\Rightarrow (e^x)'=e^x$

然後這裡順便來簡單的積一下前面提到的 $\int\frac{1}{e^x-1}dx$：
$$
\text{Let }u=e^x-1,\ du=e^xdx \Rightarrow dx=\frac{1}{u+1}du
$$

$$
\int\frac{1}{e^x-1}dx=\int\frac{1}{u(u+1)}du=\int (\frac{1}{u}-\frac{1}{u+1})du=\ln\lvert\frac{u}{u+1}\rvert+C=\ln\lvert 1-e^{-x}\rvert+C
$$

## devarajabc
> [從 CPU cache coherence 談 Linux spinlock 可擴展能力議題](https://hackmd.io/@sysprog/linux-spinlock-scalability)

為何常態分布 ($f(x|\mu,\sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}} e^{-\frac{(x-\mu)^2}{2\sigma^2}}$) 及泊松分佈 ($P(k) = \frac{\lambda^k e^{-\lambda}}{k!}$) 皆會用到自然底數 $e$ ?

當我們考慮大量的微小事件連續累積 (如粒子碰撞、隨機運動)，在極限情況下會自然產生指數函數。該極限過程在統計物理頗為常見，最終形成的分布正是以
$e$ 為基礎。在物理學中，熱擴散方程描述熱量如何隨時間在空間中傳播，其基本解往往是高斯分布。高斯分布的形式
$$
f(x) = \frac{1}{\sqrt{2\pi \sigma^2}} \exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)
$$
其中 $\exp$（即 $e$ 的指數函數）直接反映這種由隨機熱運動導致的能量擴散過程。亦即，在很多熱傳導和擴散過程中，隨機性和局部平衡最終導致一個以 $e$ 為底的分布形式。高斯分布在物理上也常視作許多獨立隨機變量經過無數次累積後達到的一個平衡態，這與中心極限定理相吻合。該平衡態正是由大量微小隨機變化共同作用的結果，而這些隨機變化在微分方程中通常以指數形式出現。

泊松分布描述在單位時間或空間內事件發生的次數，其機率密度函數為
$$
P(k) = \frac{\lambda^k e^{-\lambda}}{k!}
$$
本處 $e^{-\lambda}$ 的出現與自然界中「衰減」的現象相似。舉例來說，放射性衰變中，單個原子的衰變機率正是以指數衰減方式描述。這個指數項不僅提供描述隨時間 (或空間) 隨機事件衰減的數學工具，同時也是 normalize 整個分布（使得所有可能事件的總概率為 1）的必要因子。

物理中許多隨機事件 (如粒子撞擊、碰撞頻次) 在時間尺度上是離散的，但當這些事件的數量趨近無限時，我們可用泊松分布來描述這種極限行為。這與 Binomial 分布在 $n \to \infty$ 且 $np = \lambda$ 固定時收斂為泊松分布的過程類似，其中的極限運算也自然導出 $e^{-\lambda}$ 的形式。

從物理規律的角度看，$e$ 的出現反映自然界中普遍存在的指數型增長或衰減現象，這些現象無論在連續系統，抑或是離散事件中，都可藉由 $e$ 予以描述。

[A Realistic Evaluation of Memory Hardware Errors and Software System Susceptibility](https://www.usenix.org/legacy/event/atc10/tech/full_papers/Li.pdf)
> Based on the collected traces, we explore the implications of different hardware ECC protection schemes so as to identify the most common error causes and approximate error rates exposed to the software level. ...
> Statistical Error Rate Bounds
> Due to the small number of device-level errors in our
trace, the observed error rate may be quite different from the intrinsic error rate of our monitored system. To account for such inaccuracy, we use the concept of [p-value](https://en.wikipedia.org/wiki/P-value) bounds to provide a range of possible intrinsic error rates with statistical confidence. ... when the memory chips are considered identical, we can avoid this requirement. This is because in any time interval, their probability of having an error is the same, say q. Let N be the total number of memory chips operating, then the actual number of errors happening in this period, X, will be a random variable which conforms to binomial distribution: ... When N is very large (we simulated thousands of chips), we can approximate by assuming N approaches infinity. In this case the binomial distribution will turn into Poisson distribution.
$\to$ [前台大教授也搞錯 p-value 的意義](https://disp.cc/b/HatePolitics/gBj6)

[A Reliability Benchmarking Method for Linux](https://qrs24.techconf.org/download/webpub/pdfs/QRS-C2024-43b2F0XafenffERHWle5q5/656500a177/656500a177.pdf)
> The benchmarking results obtained can be used to compare the differences in reliability of different Linux distributions. The occurrence process of Linux kernel faults is modeled as Poisson process, and the timing of fault injection is determined according to the rule of exponential distribution of fault occurrence interval.

[Entropy of PRNGs and the Accuracy of Monte-Carlo Simulations in a Publicly Distributed Computing Environment](https://ieeexplore.ieee.org/document/9020334)
> To ensure the negation of the implementation of any languages proprietary Pseudo Random Number Generator (PRNG), [urandom](https://en.wikipedia.org/wiki//dev/random) was chosen as the PRNG. For most tests a single random bit suffices, however, the continuous time decay model must have a bit range higher than 1 otherwise the model will fail to fit a Poisson distribution due to too little entropy.

## BennyWang1007

Dudect 的運作原理是什麼？
- Dudect 會蒐集一段程式執行所經過的 CPU time，利用 Welch's t-test 計算 t-value，用來驗證程式是否是常數時間。

> $t=\frac{\overline{X}_1-\overline{X}_2}{\sqrt{s_{\overline{X}_1}^2+s_{\overline{X}_2}^2}},\ where\ {\displaystyle {\overline {X}}_{i}}\ and\ {\displaystyle s_{{\bar {X}}_{i}}}\ are\ the\ {\displaystyle i^{\text{th}}}\ sample\ mean\ and\ its\ standard\ error$

原本實作的問題
- 原先實作將所有測量值都納入統計，然而像是**中斷(interrupt)、分頁錯誤(page fault)、I/O 時間、CPU 排程**等原因，會導致異常值的出現，進而錯誤判定程式不是常數值間。

改進
- 首先，在測量初始化時，新增了 `percentiles[NUM_PERCENTILES]` 來儲存在測量時，應該要捨去的百分位。隨著測量進行時，我們應當捨棄越來越少的百分位，因此使用公式 $threshold(x)=1-{0.5}^{\dfrac{10\ \cdot\ x}{NUM\_PERCENTILES}}$ 來篩選是否要加入檢定。

- 此處 $NUM\_PERCENTILES$ 取 $100$，可得以下曲線。
- 最後在更新檢定的函式 `update_statistics()` 中，新增以下程式碼來篩選測試：若執行時間在 `percentile[i]` 百分位內，則納入檢定。（此處 `percentile[i]` 為該百分位的值）

    // t-test on cropped execution times, for several cropping thresholds.
    for (size_t j = 0; j < NUM_PERCENTILES; j++) {
        if (difference < percentiles[j]) {
            t_push(t, difference, classes[i]);
        }
    }

$\to$ pull request [Improve dudect test with cropped time analysis #283](https://github.com/sysprog21/lab0-c/pull/283)

有無母數統計（parametric/nonparametric）
- 有母數統計會假設樣本的分布模型，如常態分佈，用該分布的特性像是平均數、變異數等，能夠更容易的檢定出結果，如 t-test 檢測是否有差異。
- 無母數統計則相反，不假設樣本的分布模型，而使用 rank 來進行檢驗。

檢驗方式的選擇
- 首先我們考慮不同的 t-test：Welch's t-test 以及 Student's t-test，其中因為在常數時間的測試中，兩個測試群體（程式碼中的`class[]`）不一定擁有相同的變異數，因此採用 Welch's t-test 來檢驗。（Student's 假定兩樣本具有相同的標準差）

為什麼不使用標準差來篩選數據？
- 因為測試函式的時間機率分布**不一定為常態分佈**，使用標準差來篩選會有數據永遠篩選不掉、以及數據永遠都被篩選掉的情況發生。例如，假設函式的時間機率分布為均勻分布，考慮取 $thresholds = mean \pm 3 \times \sigma(X)$，大部分的情況下所有點都不會被篩選掉。

為什麼使用 CPU-time 而不用現實時間？
- CPU time 可以更好的反應出程式執行的總時間 $\sum\ instructions * CPI$，如果使用現實時間會因為**CPU頻率改變**導致結果出現波動。