# 歐拉數 $e$：描述連續變化的基石
> 資料整理: [jserv](https://wiki.csie.ncku.edu.tw/User/jserv), weihsinyeh


## 點題
你一定聽過圓周率 $\pi$，但你知道在數學和科學領域中，還有個同等重要的常數，歐拉數 $e$ 嗎？之所以如此重要，是因為 $e$ 是用來描述連續變化的常數。試想銀行利息不斷複利累積、熱咖啡逐漸變涼、乃至於病毒在人群中擴散等等現象背後的數學模型，都有 $e$ 的身影。

本文回顧 $e$ 的起源，一探數學家如何定義它、發現它的奇妙性質 (例如 $e^x$ 的導數仍是自身)，及它何以稱為「自然對數」的底。讀者也將理解 $e$ 在碳 14 測定法、微積分、金融、機率統計、物理領域所奠定的根基，並探討其在嵌入式系統中的數值方法，以及在機器學習領域 (sigmoid、softmax、交叉熵) 的關鍵作用。

## 用來解釋連續變化的常數

:::success
"The most powerful force in the Universe is compound interest." (宇宙中最強大的力量是複利)
> 常被歸因於愛因斯坦，確切出處不詳
:::

歐拉數 ([Euler's number](https://en.wikipedia.org/wiki/E_(mathematical_constant)), $e \approx 2.71828$，注意 Euler 讀作 [ˈɔɪlɚ]) 廣泛出現於數學、物理、生物、化學和經濟等領域，特別是與指數成長、衰減以及自然對數相關的問題。除了稱為歐拉數，$e$ 也稱為自然常數、自然底數 ([自然對數](https://en.wikipedia.org/wiki/Natural_logarithm)的底數)、Napier 常數，這些命名都反映出其作用。不過，歐拉數不能稱呼為「[歐拉常數](https://en.wikipedia.org/wiki/Euler%27s_constant)」(Euler's constant，又名 Euler-Mascheroni constant)，後者約為 $0.57721$，定義為調和級數與自然對數之差在 $n \to \infty$ 時的極限，即 $\gamma = \lim_{n\to\infty}\bigl(\sum_{k=1}^n \frac{1}{k} - \ln n\bigr)$。論及物理與工程系統時，[Euler number](https://en.wikipedia.org/wiki/Euler_number_(physics)) (縮寫為 $Eu$) 與數學中 Euler's number 完全不同，前者出現在流體力學與應用物理中，後者則源自數學分析。

人類自古以來就對未來充滿疑問：下一刻會發生什麼？人們發現有些變化的規律，當數學家從單純的計數問題 (即自然數系統的建立) 進入探討連續變化的領域時，他們需要一個新的工具來描述這種自然的、持續的變化，這便是歐拉數出現的契機。

$e$ 最初是在研究複利的過程中發現：假設你在銀行存入 $1$ 元 (下圖藍色圓圈)，且恰逢嚴重通貨膨脹，導致銀行存款利率飆升至驚人的 $100\%$。一般情況下，銀行一年才給付一次利息。當滿一年後，你能獲得 $1$ 元利息 (下圖藍綠色圓圈)

存款餘額成為 $1 \times (1 + 100\%) = 2$

若銀行改為每半年給付一次利息，並將付給你的利息自動轉為本金，下半年便能和本金一起生息 (下圖紅色圓圈):

存款餘額成為 $1 \times \Bigl(1 + \frac{100\%}{2}\Bigr)^2 = 2.25$

以此類推，若銀行一年給付 $n$ 次利息 (下圖為 $n=3$ 的情形，紅色、紫色圓圈):

存款餘額成為 $1 \times \Bigl(1 + \frac{100\%}{n}\Bigr)^n$。
> 圖片來源: [An Intuitive Guide To Exponential Functions & e](https://betterexplained.com/articles/an-intuitive-guide-to-exponential-functions-e/)

你或許會貪婪地推論：倘若 $n$ 趨近無窮大，存款金額是否也會趨近無窮大呢？試算如下：
| $n$ | $(1 + \frac{1}{n})^n$ |
|-----------|--------------|
| 1 | 2 |
| 2 | 2.25 |
| 3 | 2.37 |
| 5 | 2.488 |
| 10 | 2.5937 |
| 100 | 2.7048 |
| 1,000 | 2.7169 |
| 10,000 | 2.71814 |
| 100,000 | 2.718268 |
| 1,000,000 | 2.7182804 |
| ... | ... |


一年約有 $31,536,000$ 秒 (以 365 天計)，假設每秒收到的利息皆立即再投入本金：
$$
\biggl(1 + \frac{1}{31{,}536{,}000}\biggr)^{31{,}536{,}000} \approx 2.71828
$$

累積下來的利滾利餘額約為 $2.7182817785$ 元，並非原本貪婪設想的無窮大。換言之，上述計算得到的數值相當接近歐拉數 $e$。


這種趨近的極限即定義為 $e$，數學上可表示為：
$$
e = \lim_{x \to +\infty} \left(1 + \frac{1}{x}\right)^x
$$

金融領域廣泛使用歐拉數來計算連續複利，公式如下：
$$
FV = PV \times e^{rt}
$$

其中：
- $FV$ 為未來價值
- $PV$ 為現值 (本金)
- $r$ 為年利率
- $t$ 為時間 (年)

例如，若你有 1,000 美元，以年利率 2% 的連續複利投資 3 年，最終金額為 $1,000 \times e^{0.02 \times 3} \approx 1,061.84$ 美元。該數值略高於離散方式 (例如每月複利一次) 的計算結果，隨著本金增大、利率增高或時間延長，差距將更加明顯。

當銀行年利率為 $5\%$，我們想知道 100 元存款何時可翻倍，只要將問題轉化為以下方程式：
$$
100 \,\times e^{\,0.05 \,t} \;=\; 200
$$

並加以求解，可得
$$
t \;=\; \frac{\ln(2)}{0.05} \;=\; \frac{0.693}{0.05} \;\approx\; 13.86
$$

單位是「年」，亦可近似計算 $\frac{72}{5} \;\approx\; 14.4$，這裡的 $72$ 源自 $100 \times \ln(2) \approx 69.3$，選用 72 的考量是其因數較多，便於心算。這是金融領域的「[72 法則](https://zh.wikipedia.org/zh-tw/72%E6%B3%95%E5%89%87)」，亦即估算資產翻倍所需的時間。


公元 1683 年，瑞士數學家 Jacob Bernoulli 研究複利頻率對本金增長的影響，顯示上述極限值介於 2 與 3 之間，但未能精確計算其值 (參見 Bernoulli《[Ars Conjectandi](https://en.wikipedia.org/wiki/Ars_Conjectandi)》附錄)。約五十年後，Leonhard Euler 將此數精確計算至 23 位小數 (見《Introduction》第 6 章，§122)。字母 $e$ 最早見於其約 1727/28 年的未發表手稿，公元 1731 年致 [Goldbach 的書信](https://en.wikipedia.org/wiki/Euler%27s_number#History)是目前已知最早使用此符號的公開文獻，而公元 1736 年出版的《[Mechanica](https://en.wikipedia.org/wiki/Mechanica)》則是 $e$ 首次出現在印刷品中的紀錄。Euler 於公元 1737 年以連分數方法證明 $e$ 為無理數 (該證明於 1744 年發表，見 [E71: De fractionibus continuis dissertatio](http://eulerarchive.maa.org/pages/E071.html))，並在公元 1748 年的著作《[Introductio in Analysin Infinitorum](https://en.wikipedia.org/wiki/Introductio_in_analysin_infinitorum)》中系統性地詳述 $e$ 的性質與級數表示，後人遂以 Euler 之名將該常數命名為「歐拉數」。

## 自然界與歐拉數的關聯
人類對自然界的理解，最初依賴宗教與神話。然而，自文藝復興時代起，對時間的精確測量與數學分析逐漸讓人類意識到，自然界的變化並非隨機無序，而是遵循一定的規律，且可用數學表達。當伽利略透過斜面實驗與水鐘計時，系統性地研究落體運動並發現等加速度規律時，人類正式踏上以數學和科學解析自然的道路。

在自然界中，許多現象的變化速率與其目前狀態成正比，例如：
- 熱水的冷卻速度與熱水與環境溫度的溫差成正比 ([牛頓冷卻定律](https://episte.math.ntu.edu.tw/applications/ap_cooling/))
- 生物族群的成長速度與目前族群大小成正比 ([指數成長模型](https://en.wikipedia.org/wiki/Exponential_growth))
- 放射性元素的衰變速率與剩餘原子數目成正比 ([放射性衰變](https://en.wikipedia.org/wiki/Radioactive_decay))

這類現象可用微分方程描述：
$$
\frac{dy}{dx} = ky
$$

其通解為指數函數：
$$
y = Ce^{kx}
$$

其中的 $e$ 就是歐拉數，描述這類自然變化模式的常數，其數學定義源於極限與無窮級數，表示如下：
$$
e = \sum_{n=0}^{\infty} \frac{1}{n!}
$$

此級數與前述極限定義的等價性，可由二項式定理推導。對 $\bigl(1+\frac{1}{n}\bigr)^n$ 展開：
$$
\left(1+\frac{1}{n}\right)^n = \sum_{k=0}^{n} \binom{n}{k} \frac{1}{n^k}
= \sum_{k=0}^{n} \frac{1}{k!} \cdot \frac{n(n-1)(n-2)\cdots(n-k+1)}{n^k}
$$

後方因子可改寫為：
$$
\frac{n(n-1)\cdots(n-k+1)}{n^k} = 1 \cdot \left(1 - \frac{1}{n}\right) \cdot \left(1 - \frac{2}{n}\right) \cdots \left(1 - \frac{k-1}{n}\right)
$$

對固定的 $k$，當 $n \to \infty$ 時，每個因子趨近 $1$，故整體趨近 $\frac{1}{k!}$。級數的尾部 ($k > K$ 的部分) 可用 $\sum_{k>K} 1/k!$ 控制，因此：
$$
\lim_{n \to \infty} \left(1+\frac{1}{n}\right)^n = \sum_{k=0}^{\infty} \frac{1}{k!}
$$

對實數或複數 $z$ 進行推廣，可定義：
$$
\exp(z) := \lim_{n \to \infty} \left(1+\frac{z}{n}\right)^n = \sum_{k=0}^{\infty} \frac{z^k}{k!} = e^z
$$

以同樣的二項式展開，$\bigl(1+\frac{z}{n}\bigr)^n = \sum_{k=0}^{n} \binom{n}{k} \frac{z^k}{n^k}$，當 $n \to \infty$ 時各項中的因子 $\prod_{j=0}^{k-1}(1-j/n)$ 皆趨近 $1$，即得 $\sum z^k/k!$。此推廣使 $e$ 從一個數值常數提升為指數函數的「單位」，正如 $\pi$ 是圓的度量單位、$i$ 是旋轉的度量單位。任意底數 $b > 0$ 的實數指數函數皆可歸約至 $e$：
$$
b^x = e^{x \ln b}
$$

這個展開式橫跨離散數學中的「階乘」與連續數學中的「無窮和」，使 $e$ 成為離散與連續世界的橋樑。當物理學家研究熱水冷卻的速率、生物學家分析細菌繁殖的過程時，他們最終發現，這些變化的數學模型往往依賴以 $e$ 為底的指數函數。


在統計力學中，[波茲曼因子](https://en.wikipedia.org/wiki/Boltzmann_distribution) $e^{-E/kT}$ 描述熱平衡狀態下粒子處於能量 $E$ 的相對機率，其中 $k$ 為[波茲曼常數](https://en.wikipedia.org/wiki/Boltzmann_constant)，$T$ 為絕對溫度。建立在此因子之上的[馬克士威-波茲曼分佈](https://en.wikipedia.org/wiki/Maxwell%E2%80%93Boltzmann_distribution)描述氣體分子速度的機率分佈，反映高能分子出現的可能性隨速度平方指數遞減的特性。隨機運動方面，[布朗運動](https://en.wikipedia.org/wiki/Brownian_motion)是微小顆粒在流體中因分子碰撞而產生的隨機運動，其擴散機率密度呈高斯分佈 $\propto e^{-x^2/(4Dt)}$ (其中 $D$ 為擴散係數)，同樣以 $e$ 為主體。

在核物理學中，[放射性衰變](https://en.wikipedia.org/wiki/Radioactive_decay) (radioactive decay) 描述不穩定的原子核自發性地轉變成其他種類原子核的自然現象。這種轉變遵循著明確的數學規律，也就是指數衰減定律 (exponential decay law)，其數學形式如下：
$$
N(t) = N_0 e^{-\lambda t}
$$

其中：
* $N(t)$：代表在時間 $t$ 時，樣品中還剩下未衰變的原子核數量
* $N_0$：代表在初始時刻 ($t=0$) 的原子核數量
* $\lambda$：稱為衰變常數 (decay constant)，是個正數，代表該種原子核發生衰變的速率或可能性。$\lambda$ 的值越大，表示原子核衰變得越快

> 放射性原子核數量隨時間呈指數遞減

原子核的數量會隨著時間以指數方式減少，亦即每經過一段固定的時間，剩下的原子核數量佔原本數量的比例都是相同的。當一個原子核透過放出輻射 (例如 $\alpha$ 粒子、$\beta$ 粒子或 $\gamma$ 射線) 釋放能量並轉變為另一個原子核時，這個過程就稱為衰變。

幾乎所有原子序大於 82 的重元素 (如鈾、釷) 的原子核都會發生放射性衰變，而某些較輕元素的不穩定同位素也具備同樣的衰變行為。其中，碳-14 ($^{14}\text{C}$) 是一個重要的例子：宇宙射線產生的高能中子撞擊大氣中的氮原子，依循反應式 $^{14}_{7}\mathrm{N} + ^{1}_{0}\mathrm{n} \rightarrow ^{14}_{6}\mathrm{C} + ^{1}_{1}\mathrm{H}$ 持續產生放射性的 $^{14}\text{C}$。這些 $^{14}\text{C}$ 原子迅速氧化成 $^{14}\text{CO}_2$，並與大氣中主要的 $^{12}\text{CO}_2$ 混合。植物透過光合作用吸收這些二氧化碳，再經由食物鏈傳遞給動物，使生物體內的 $^{14}\text{C}$ 與 $^{12}\text{C}$ 比例大致與大氣保持平衡。當生物死亡後，碳交換停止，體內的 $^{14}\text{C}$ 因 $\beta$ 衰變 ($^{14}\text{C} \rightarrow ^{14}\text{N} + e^{-} + \bar{\nu}_e$) 而逐漸減少，半衰期約為 5730 年。

> 圖片出處: [Free Webinar: Choosing Optimal Bone Samples for Analysis](https://www.radiocarbon.com/optimal-bone-samples-radiocarbon-dating/)

[放射性碳定年法](https://en.wikipedia.org/wiki/Radiocarbon_dating) (radiocarbon dating) 正是利用上述衰變規律來測定古老有機材料的年代。測量樣品中殘留的 $^{14}\text{C}$ 比例，再代入指數衰減公式 $N(t) = N_0 e^{-\lambda t}$，即可反推生物死亡至今的時間，此方法可有效測定約 5 萬年內的樣本年代。自公元 1949 年 Arnold 與 Libby 證實其可行性後，$^{14}\text{C}$ 定年法對考古學產生革命性影響，為遺址和文化演進提供精確的時間框架，在地質學中也廣泛應用於測定晚第四紀以來的地質事件年代。

半衰期是描述衰變快慢的常用概念，符號為 $T_{1/2}$，定義是樣品中一半數量的原子核發生衰變所需要的時間。舉例來說，假設某種放射性同位素的半衰期是 10 年。這就代表：
* 經過 10 年後，原子核數量會剩下原來的 $\frac{1}{2}$
* 再經過 10 年 (總共 20 年) ，數量會剩下 $\frac{1}{2} \times \frac{1}{2} = \frac{1}{4}$
* 再經過 10 年 (總共 30 年) ，數量會剩下 $\frac{1}{4} \times \frac{1}{2} = \frac{1}{8}$

半衰期 $T_{1/2}$ 與衰變常數 $\lambda$ 之間關係如下，說明指數衰減的行為：
$$
N(T_{1/2}) = N_0 e^{-\lambda T_{1/2}} = \frac{1}{2} N_0
$$

從上式可推導出：
$$
e^{-\lambda T_{1/2}} = \frac{1}{2} \implies -\lambda T_{1/2} = \ln\left(\frac{1}{2}\right) = -\ln(2)
$$

$$
\lambda = \frac{\ln(2)}{T_{1/2}} \approx \frac{0.693}{T_{1/2}}
$$

亦即，衰變越快 ($\lambda$ 越大)，則半衰期 $T_{1/2}$ 就越短。放射性物質的活性 (activity) 指單位時間內發生衰變的次數 (即衰變速率)，通常用符號 $A(t)$ 表示。活性 $A(t)$ 正比於當時存在的原子核數量 $N(t)$，即 $A(t) = \lambda N(t)$。因此，活性也會隨時間以完全相同的指數速度下降：$A(t) = A_0 e^{-\lambda t}$，其中 $A_0 = \lambda N_0$ 是初始活性。

我們知道放射性衰變遵循 $N(t) = N_0 e^{-\lambda t}$ 這個規律，但這個公式是怎麼來的呢？能否從更基本的假設出發來推導？為何非得要用到 $e$ 不可呢？

試著從微觀角度來思考。考慮每個不穩定的原子核，它發生衰變是個隨機事件。但在一個很短的時間間隔 $dt$ 內，任何一個尚未衰變的原子核，都有個發生衰變的微小機率，該機率應該正比於時間間隔 $dt$ 的長度，我們可描述為：單一原子核在 $dt$ 時間內衰變的機率 $p = \lambda dt$，此處的 $\lambda$ 就是前面提到的衰變常數，它反映該種原子核內在的不穩定性。

現在，考慮一大堆 (例如 $N$ 個) 相同的放射性原子核。在 $dt$ 這個微小的時間間隔內，每個原子核獨立衰變，衰變總數服從二項分佈。根據大數法則，當 $N$ 足夠大時，衰變數量的期望值主導平均行為，亦即:
$$
E[dN] = - N \times p = - N \lambda \, dt
$$

其中負號表示 $N$ 隨時間減少。在平均行為下 (即以期望值取代隨機變數)，當 $dt$ 趨近於零時，上式成為確定性的微分方程 $\frac{dN}{dt} = -\lambda N$，亦即原子核數量的變化速率 $\frac{dN}{dt}$ 與目前的數量 $N$ 成正比，比例常數為 $-\lambda$。

接著求解這個微分方程。利用變數分離法，將所有跟 $N$ 有關的項移到等號左邊，跟 $t$ 有關的項移到等號右邊: $\frac{dN}{N} = -\lambda dt$，然後二邊同時進行積分: $\int \frac{dN}{N} = \int -\lambda dt$

積分結果為 $\ln N = -\lambda t + C$，這裡的 $C$ 是積分常數。$\ln N$ 是以 $e$ 為底的自然對數。為了得到 $N(t)$ 的顯式表達式，我們對等號兩邊同時取指數 (以 $e$ 為底):
$e^{\ln N} = e^{-\lambda t + C}$
$N = e^{-\lambda t} e^C$

由於 $e^C$ 也是個常數，我們可用新的符號 $A$ 來代表它 ($A = e^C$) 。所以得到通解：
$N(t) = A e^{-\lambda t}$

這個常數 $A$ 的值由初始條件決定。假設在 $t=0$ 時，原子核的數量為 $N(0) = N_0$。將 $t=0$ 代入通解：
$$
N(0) = A e^{-\lambda \times 0} = A e^0 = A \times 1 = A
$$

於是常數 $A$ 就是初始原子核數量 $N_0$。最終解就是 $N(t) = N_0 e^{-\lambda t}$

在上述推導與求解過程中，不僅解釋放射性衰變為何符合指數變化規律，也揭開 $e$ 在此現象中扮演的角色 (源自於 $\int \frac{dN}{N} = \ln N$) 。

指數函數 $e^{-\lambda t}$ (或更一般的 $e^{kt}$) 與微分運算之間，尚存在非常重要的特性。觀察微分方程 $\frac{dN}{dt} = -\lambda N$ 的解 $N(t) = N_0 e^{-\lambda t}$。對解進行微分：
$$
\frac{d}{dt} N(t) = \frac{d}{dt} (N_0 e^{-\lambda t}) = N_0 \frac{d}{dt} (e^{-\lambda t}) = N_0 (-\lambda e^{-\lambda t}) = -\lambda (N_0 e^{-\lambda t}) = -\lambda N(t)
$$

這驗證上方的解，同時也突顯關鍵性質：
$\frac{d}{dt} (e^{-\lambda t}) = -\lambda (e^{-\lambda t})$

亦即，對指數函數 $e^{-\lambda t}$ 進行微分運算，得到的結果是原函數自身再乘上一個常數 $-\lambda$。

倘若我們將微分運算子 $\frac{d}{dt}$ 看作是作用在函數上的「操作」，那麼指數函數 $e^{-\lambda t}$ 就好似該操作的「特殊函數」：它在被操作後，其「形態」(指數函數的形式) 保持不變，僅乘上一個數值 $-\lambda$。

這與線性代數中的特徵值 (eigenvalue) 和特徵向量 (eigenvector) 的概念極為相似。對於一個矩陣 $A$ 和一個非零向量 $\mathbf{v}$，若它們滿足關係式：
$A\mathbf{v} = \lambda \mathbf{v}$

則 $\lambda$ 就被稱為矩陣 $A$ 的特徵值，而 $\mathbf{v}$ 則是被稱為對應於 $\lambda$ 的特徵向量。矩陣 $A$ 作用在它的特徵向量 $\mathbf{v}$ 上，效果等同於將 $\mathbf{v}$ 乘以純量 $\lambda$。

接著我們可類比：
* 將微分運算子 $\frac{d}{dt}$ 看作是線性代數中的「矩陣 $A$」
* 將指數函數 $e^{-\lambda t}$ 看作是線性代數中的「向量 $\mathbf{v}$」

於是，關係式 $\frac{d}{dt} (e^{-\lambda t}) = -\lambda (e^{-\lambda t})$ 就完全對應於 $A\mathbf{v} = \lambda \mathbf{v}$ 的形式。

在這裡：
* 指數函數 $e^{-\lambda t}$ 扮演特徵函數 (eigenfunction) 的角色
* 常數 $-\lambda$ 扮演特徵值 (eigenvalue) 的角色

這種微分運算子與其特徵函數、特徵值之間的關係，不只是個有趣的類比，還揭露指數函數之所以特殊和重要的原因：在物理學和工程學的許多領域中，基本定律往往以微分方程的形式出現，而這些方程的解，常與特定「運算子」(例如 $\frac{d}{dt}$ 或更複雜的運算子) 的特徵函數和特徵值緊密相關。

特徵函數通常代表系統的一種「固有」行為模式或狀態。例如，在振動系統中，它們可能代表特定的自然振動頻率。在量子力學中，$e$ 的角色更為根本：
* 定態薛丁格方程式 $H\psi = E\psi$ 是個典型的特徵值問題，能量運算子 $H$ 作用在穩態波函數 $\psi$ (特徵函數) 上，得到能量值 $E$ (特徵值) 乘以該波函數
* [時間演化算符](https://en.wikipedia.org/wiki/Time_evolution#In_quantum_mechanics) $U(t) = e^{-iHt/\hbar}$ 描述量子態隨時間的演化，其中 $H$ 為系統的 Hamiltonian，$\hbar$ 為約化 Planck 常數。此算符之所以取指數形式，根源於薛丁格方程 $i\hbar \frac{\partial}{\partial t}\psi = H\psi$ 與前述 $\frac{dy}{dt} = ky$ 的結構一致：解為 $\psi(t) = e^{-iHt/\hbar}\psi(0)$
* 自由粒子的[平面波](https://en.wikipedia.org/wiki/Plane_wave)解 $\psi = e^{i(kx - \omega t)}$ 同時涉及歐拉公式 $e^{i\theta} = \cos\theta + i\sin\theta$，將波動的振幅與相位統一於單一指數函數中

這些例子揭示 $e$ 從古典的連續變化 (複利、衰變、冷卻) 延伸至量子態的機率波演化，其主體始終是微分方程 $\frac{dy}{dt} = ky$ 的指數解結構。此類比在本文後段「橫跨連續現實與離散計算的橋樑」一節中將再次浮現；離散系統的移位運算子 $z$ 與連續微分運算子 $D$ 之間，正是透過 $z = e^{hD}$ 建立特徵值的對映關係。

$\frac{dN}{dt} = -\lambda N$ 這個微分方程以及相應的特徵函數結構，並不侷限於放射性衰變，而是貫穿前文列舉的複利、冷卻、族群成長等現象。$e$ 並非僅是個數學常數，而是連續變化規律的根本法則。

> "How can it be that mathematics ... is so admirably adapted to the objects of reality?" — 物理學家 [Wigner](https://en.wikipedia.org/wiki/The_Unreasonable_Effectiveness_of_Mathematics_in_the_Natural_Sciences), 1960

延伸閱讀: [e^(iπ) in 3.14 minutes, using dynamics](https://youtu.be/v0YEaeIClKY)

## 自然底數的由來
前文藉由極限定義 $e$，但科學史上研究自然對數的起點遠早於此。對數的出現甚至早於指數的形式化定義。

> [John Napier](https://en.wikipedia.org/wiki/John_Napier) [ˈneɪpiər], 1550–1617

17 世紀初，John Napier 率先發表對數表 (詳見後文)，而直到公元 1637 年，法國數學家笛卡兒才引入現代指數記號系統。真正釐清二者之間的關係，則已是百餘年後的公元 1748 年，由歐拉首次明確指出：以現代觀點看，對數是指數的反函數。

對數之所以先於指數而被發明，根本原因在於其應用價值：17 世紀的歐洲，隨著天文學與航海技術的發展，人們面對大量數值計算，亟需更有效率的處理方式，對數因應而生。在對數問世以前，天文學與航海領域高度依賴三角函數表來進行繁複的計算，這些表格詳列各角度對應的正餘弦值，精確至小數點後多位，於是數學家利用以下三角恆等式：
$$
\cos(x) \cos(y) = \frac{\cos(x+y) + \cos(x-y)}{2}
$$
將二角的乘積轉化為和差形式，也就是著名的「積化和差公式」(prosthaphaeresis)。這種手法雖能間接完成乘法，但本身並非對數，而是對數問世前的替代方案，操作步驟繁瑣：需查閱三角表二次、執行二次反查 (從 $\cos$ 值回推角度)、再搭配數次加減與平均操作。

對數的誕生汲取相同的主體想法，將「乘法化為加法」，卻以更通用、更易操作的形式實現。在以 $10$ 為底的常用對數 (記為 $\log_{10}$，有別於後文以 $e$ 為底的自然對數 $\ln$) 下，$x \cdot y = 10^{\log_{10} x + \log_{10} y}$，流程極為簡潔：
- 查對數表二次，得 $\log_{10} x$ 與 $\log_{10} y$
- 加總二值
- 查一次反對數表，即得乘積

以計算 $567.89 \times 3141.59$ 為例：
$$
\begin{aligned}
\log_{10}(x) &= \log_{10}(567.89 \times 3141.59) \\
&= \log_{10}(567.89) + \log_{10}(3141.59) \\
&= \log_{10}(10^2 \times 5.6789) + \log_{10}(10^3 \times 3.14159) \\
&= 5 + \log_{10}(5.6789) + \log_{10}(3.14159)
\end{aligned}
$$

查表得 $\log_{10}(5.6789) \approx 0.7543$、$\log_{10}(3.14159) \approx 0.4971$，合計 $\log_{10}(x) \approx 6.2514$。再查反對數表得 $10^{0.2514} \approx 1.784$，因此 $x \approx 10^6 \times 1.784 = 1{,}784{,}000$，與直接乘法結果 $1{,}784{,}077.5451$ 的相對誤差僅約 $0.004\%$。整個過程只需查表二次、加法一次、再反查一次，即可完成大數乘法，而對於兩個任意大數的乘積 $x_1 x_2$，令 $x_1 = a_1 \times 10^{n_1}$、$x_2 = a_2 \times 10^{n_2}$ ($1 \le a_1, a_2 < 10$)，有：
$$
\log_{10}(x_1 x_2) = n_1 + n_2 + \log_{10}(a_1) + \log_{10}(a_2)
$$

將 $a_1$、$a_2$ 的值四捨五入到最近的四位小數，查對數表得 $\log_{10}(a_1)$、$\log_{10}(a_2)$，依上式加總即得 $\log_{10}(x_1 x_2)$。將結果拆分為整數部分 $N$ 和 $[0, 1)$ 的小數部分 $f$，再以反對數表求 $10^f$，最終 $x_1 x_2 \approx 10^N \times 10^f$。

一張對數表足以大幅加速乘除與開根號等運算，其後演化出的[對數尺](https://en.wikipedia.org/wiki/Slide_rule) (slide rule)，在尺面上以對數刻度標示數字，使乘法化為刻度距離的相加 (依據 $\log_b(p \times q) = \log_b p + \log_b q$)，工程師只要滑動尺身、對齊刻度，即可直接讀取乘積，不必查表或筆算，堪稱攜帶型類比計算器，猶如現今理工科學生獲贈的科學計算器。
然而，當時的對數計算主要用[常用對數](https://en.wikipedia.org/wiki/Common_logarithm) (底數為 $10$)，這已足以解決日常計算需求，不會特別涉及 $e$。

Jost Bürgi 曾擔任天文學家克卜勒的助手，負責大量繁瑣的天文計算工作，受前人 Michael Stifel 研究等差與等比數列關係的啟發，Bürgi 發展對數概念，並於公元 1610 年左右完成對數表，但直到公元 1620 年才正式發表於《Arithmetische und Geometrische Progress Tabulen》(等差與等比數列表) 一書。

遠在蘇格蘭的 Napier 也展開類似的工作，他出身貴族，是愛丁堡附近 Merchiston 城堡的第八代地主，未曾有過正式的職業，但對天文與數學充滿熱情，完全獨立於 Bürgi 的成果。古希臘哲學家提出「自然」(φύσις) 的概念，認為萬物自然而然地變化，哲學家赫拉克利特引入 λόγος (英語: Logos) 一詞，表達自然變化背後所隱含的秩序與規律性。λόγος 原本含義廣泛，可指言語、演講、敘事、論理，或是原則與法則。在哲學語境中，它逐漸轉化為代表萬物背後的尺度、分寸與比例，即強調萬物之間的數量關係。Napier 結合 Logos (比例) 與 Arithmos (數) 二詞並創造出 Logarithm 一詞，意指「比例之數」，也就是後來的「對數」。

Napier 最初構想的對數，與我們今日以指數函數為基礎的對數概念迥然不同，他不是從代數公式出發，而是以幾何的方式定義對數：透過二個以不同速度移動的質點，其行進距離之間的對應關係，來刻畫數量變化的比例性，如下圖：

圖中展示 Napier 定義對數的運動學模型：想像二個質點沿線段運動，上方 Moving $b$ Particle，等速運動。下方的速度 Moving $\beta$ Particle 是等比例遞減。接下來想要找一個數學函數去描述這兩個質點位置的關係。先寫成簡單的表格。
| 時間 $t$ | 位置 $x(t)$ | 速度 $v(t)$ | $y = \log x$ | $\Delta y$ |
| --- | --- | --- | --- | --- |
| 0 | $x$ | (等比遞減) | $\log(x)$ | (等速) |
| 1 | $x - \frac{x}{2} = \frac{x}{2}$ | $\frac{x}{2}$ | $\log(\frac{x}{2})$ | $\log(\frac{x}{2}) - \log(x)$ |
| 2 | $x - \frac{x}{2} - \frac{x}{4} = \frac{x}{4}$ | $\frac{x}{4}$ | $\log(\frac{x}{4})$ | $\log(\frac{x}{4}) - \log(\frac{x}{2})$ |
| 3 | $x - \frac{x}{2} - \frac{x}{4} - \frac{x}{8} = \frac{x}{8}$ | $\frac{x}{8}$ | $\log(\frac{x}{8})$ | $\log(\frac{x}{8}) - \log(\frac{x}{4})$ |

如此可發現若質點以「等比例遞減的速度」移動，另一質點以「等速度」移動，兩者之間的距離關係就是 $y = \log(x)$ 。

Napier 的原始定義是幾何式的 (參見[〈指數與對數發展史簡介〉](https://www.sec.ntnu.edu.tw/uploads/asset/data/66cf6cde0e0b30672f4bad39/1983-056-10_58-73_.pdf))：取線段 $\overrightarrow{AB}$ (長度 $10^7$) 與射線 $\overrightarrow{CD}$，質點 $P$ 從 $A$ 出發沿 $\overrightarrow{AB}$ 移動，速度等於剩餘距離 $\overrightarrow{PB}$ (等比遞減)；質點 $Q$ 從 $C$ 出發沿 $\overrightarrow{CD}$ 等速移動。二者同時出發，Napier 定義 $\overrightarrow{CQ}$ 為 $\overrightarrow{PB}$ 的對數。


分析此定義：

若 P 與 Q 兩質點的初速度均為 $10^7$，則可發現
- 質點 $P$ 在位置 $A$ 時的速度為 $\overrightarrow{AB}$ 的長度
- 質點 $P$ 到位置 $P$ 時的速度為 $\overrightarrow{PB}$ 的長度

接下來把質點 $P$ 從位置 $A$ 移到位置 $B$ 的過程分成 $n$ 段，每一段經過的時間為 $t, (0 < t < 1)$。
第 $i$ 段的起始位置為 $P_{i-1}$，第 $i$ 段的末位置為 $P_{i}$。
因為 $t$ 很小，所以質點 $P$ 在 $\overrightarrow{P_{i-1}P_{i}}$ 移動的過程等速，且該速度與質點 $P$ 在位置 $P_{i-1}$ 時一樣，也就是 $\overrightarrow{P_{i-1}B}$ 。因此寫成下面的式子。
(距離) $\overrightarrow{P_{i-1}P_{i}} = t \times \overrightarrow{P_{i-1}B}$ (時間 x 速度)。
用兩個向量相減的計算得出 : $\overrightarrow{P_{i}B} = \overrightarrow{P_{i-1}B} - \overrightarrow{P_{i-1}P_i}$，再將上面 $\overrightarrow{P_{i-1}P_{i}}$ 代入後得下面的式子。
$\overrightarrow{P_{i}B} = \overrightarrow{P_{i-1}B} - t \times \overrightarrow{P_{i-1}B}=(1-t)\overrightarrow{P_{i-1}B}$ 。突然發現每一段相同時間經過的距離會等比例遞減，也就是速度會等比例遞減，跟上面推導的結果一致。也就是說若一個質點的速度是等比例遞減，另一個質點的速度為等速，即可用上述表格描述。這樣兩個速度的變化就是 $\overrightarrow{CQ} = \log(\overrightarrow{PB})$
最後再把示意圖中位置 $A$、$B$、$C$、$D$ 與動點 $P$、$Q$ 都標註在圖上，即可得知首圖的由來。
雖然現在透過對數原理可知一個質點目前速度為 $x$，前一刻速度為 $\frac{x}{1-t}$。代入 $y = \log(x)$ 求另一質點的距離 $y$，$y' = \log(\frac{x}{1-t})$ 即可得知該質點在固定時間內會多行進 $-\log(1-t)$ (即 $\log\!\bigl(\frac{1}{1-t}\bigr)$，為正的常數)，呈現等速運動。但當時是為了運用數學描述這兩者間的關係，才發展出此對數表達式。

上圖將前述二個質點的運動軌跡合併呈現。上方等速質點的軌跡對應 $y = \log(x)$，下方等比遞減質點的軌跡對應 $x = \sin(\theta)$，比例常數 $R = \text{Sinus Totus} = 10{,}000{,}000$。上方標記點 $A$–$I$ 與下方的 $\alpha$–$\omega$ 構成一對一映射，以垂直線段連結，展現等差運動與等比運動之間的幾何對應。該模型將對數視為描述運動變化的比例關係，並與當時天文計算所需的三角函數緊密結合。

Napier 歷時近二十年編製出史上第一份公開發表的對數表，於公元 1614 年發表其鉅作《Mirifici Logarithmorum Canonis Descriptio》(對數奇妙表的說明)，該書甫一出版，便在歐洲科學界引起廣泛關注，倫敦數學家 Henry Briggs 在公元 1615 年前往愛丁堡拜訪 Napier，他們一致認為，對數的使用不應僅限於三角函數表的計算需求，更可延伸至一般十進位數的運算。

Briggs 改以 10 為底的對數，這便是「常用對數」。Briggs 後來成為劍橋大學聖約翰學院院士與牛津大學薩維爾幾何學教授，他以現代化的記號與計算方法，構造出更實用的常用對數表。公元 1624 年，克卜勒親自編輯並補充出版 Napier 的對數表，使其更適合用於天文與數學領域。Napier 與 Briggs 的合作成果，成功將對數從單一的天文計算工具擴充為通用的數值分析方法。正因如此，數學家拉普拉斯讚嘆：「對數的發明得以簡化計算，相當於延長天文學家的壽命。」

以下範例具體展示對數表在天文學中的威力。假設地球以圓形軌道繞太陽運動，已知地球與太陽的平均距離 $R = 1.496 \times 10^{11}$ m、地球公轉週期 $T = 3.156 \times 10^{7}$ s、萬有引力常數 $G = 6.672 \times 10^{-11}$ m$^3$/(s$^2 \cdot$ kg)，根據萬有引力定律與牛頓運動定律可知太陽質量：
$$
M = \frac{4\pi^2 R^3}{GT^2}
$$

不用計算器，僅以對數表做加減法即可求出。分別查表並計算：
$$
\begin{aligned}
\log(4\pi^2) &= 2\log(2\pi) = 2 \times 0.7982 = 1.5964 \\
\log(R^3) &= 3\log(R) = 3 \times (11 + 0.1750) = 33.5250 \\
\log(GT^2) &= \log(G) + 2\log(T) = (-11 + 0.8242) + 2 \times (7 + 0.4991) = 4.8224
\end{aligned}
$$

$$
\log(M) = 1.5964 + 33.5250 - 4.8224 = 30.2990
$$

因此 $M = 10^{30.2990} \approx 1.991 \times 10^{30}$ kg，與目前公認值 $(1.98855 \pm 0.00025) \times 10^{30}$ kg 的相對誤差約 $0.1\%$。僅靠一張對數表與四則運算，天文學家就能完成如此龐大的數值計算。

數學家在編撰對數表的過程中，面對繁瑣的計算步驟，意識到若採用以下形式作為底數：
$$
\Bigl(1 + \frac{1}{n}\Bigr)^n, \quad n \gg 0
$$
可顯著簡化計算。以現代觀點回溯，Napier 採用的底數可近似理解為 $\Bigl(1 - \frac{1}{10^7}\Bigr)^{10^7} \approx 1/e$，而 Bürgi 則選擇 $\Bigl(1 + \frac{1}{10^4}\Bigr)^{10^4} \approx e$ (參見 Cajori《[A History of Mathematics](https://archive.org/details/historyofmathema00caborich)》關於 Napier 與 Bürgi 的章節)。須注意，當時尚未明確建立「底數」的概念，Napier 的系統更為複雜，並非直接使用現代意義的底數，上述近似是後人以指數框架重新詮釋的結果。

或許讀者會好奇，當初 Napier 為何要選擇如此特別的對數底數？事實上，這是為了提高對數表在計算上的實用性。若底數選擇得太大，則對數表的資料密度 (即計算精度) 將變得稀疏，導致許多數字在表中找不到適合的對數值，內插誤差便難以控制，嚴重影響對數表的效能。

舉例來說，若以 $2$ 作為對數的底數，當對數值從 $1$ 增加到 $10$ 時，對應的原始數值將從 $2^1=2$ 急遽增長到 $2^{10}=1024$。在此情況下，若要查找如 $x=798$ 這樣的數值，對數表內可能僅有附近的 $512$ (即 $2^9$) 與 $1024$ (即 $2^{10}$) 這二個參考點，內插時誤差就會非常大。

因此，我們想尋找一個理想的底數，其性質是當指數以固定幅度增加時，原始數值不會呈倍數暴增，而是以穩定、近似線性的方式成長。換言之，新的問題是：能否將乘法導致的等比變化，轉換為加法產生的等差變化？

我們假設存在某個底數 $a$，使得其指數函數滿足以下對應關係：
$$
a^{b + \Delta b} = y + \Delta y
$$

儘管 $a^b$ 屬於等比級數，但我們希望其對應的值 $y$ 則呈現等差級數成長，亦即：
$$
\begin{aligned}
a^{b-1} &\longrightarrow y - \Delta y \\
a^b &\longrightarrow y \\
a^{b+1} &\longrightarrow y + \Delta y \\
a^{b+2} &\longrightarrow y + 2\Delta y
\end{aligned}
$$

此時，左側的變化是每次乘以固定倍數 $a$，而右側則是每次加上固定增量 $\Delta y$，也就是將原本的等比變化對應為等差變化。這種對應不僅有助於簡化計算，更提供一種理解指數成長的加法視角，最終形成對數的概念。

然而，絕大多數底數並不具備這種特性。若 $a$ 過大，$a^b$ 的變化將過於劇烈，無法對應出平滑的等差變化。

數學家如何克服上述問題，從而建立以 10 為底的對數表呢？

根據換底公式：
$$
\log_a x = \frac{\log_b x}{\log_b a}
$$
只要選定一個適合的底數，並針對該底構造出一份覆蓋範圍廣泛的對數表 (例如從 $\log(0.0001)$ 至 $\log(9.9999)$)，就可用換底公式，輕易推導出其他底數的對數值。

初看之下，$10$ 似乎是個合理選擇。若取 $a = 10$，雖然 $\log_a x$ 的間距可設為極小值如 $0.0001$，但其對應的真數 (即對數運算中的輸入值，如 $\log_a{x}$ 的 $x$) 值就變得難以計算：需對 $10, 100, 1000,\dots$ 等數開 $10000$ 次方根，這在當時不可行。

為簡化真數的生成，我們可反向思考，取某個數的 $10000$ 次方為底，例如 $a = 10^{10000}$。此時雖可輕易列出真數 $1, 10, 100, 1000\dots$，但由於真數的變化過於劇烈，在高位時增幅失控。

接著嘗試以較小的底，如 $a = 2^{10000}$，則表中真數如下：

| $\log_a x$ | $x$ |
|-------------|--------|
| 0.0000 | 1 |
| 0.0001 | 2 |
| 0.0002 | 4 |
| 0.0003 | 8 |
| 0.0004 | 16 |

雖然表中的數值間距縮小，但仍增加過快，不適合廣泛採用。

直到我們選用更接近 $1$ 的底數，例如 $a = 1.5^{10000}$，情況才有所改善：

| $\log_a x$ | $x$ |
|-------------|-------------|
| 0.0000 | 1.0000 |
| 0.0001 | 1.5000 |
| 0.0002 | 2.2500 |
| 0.0003 | 3.3750 |
| 0.0004 | 5.0625 |

我們從中發現規律：上述各表的底數皆形如 $a = b^N$ (其中 $N = 10000$)，表中相鄰真數的比值恰為 $b$。當 $b$ 越接近 $1$，相鄰真數的間距就越均勻，計算上也更合理。若我們選擇 $b = 1.0001$ (即 $a = 1.0001^{10000}$)，結果如下：

| $\log_a x$ | $x$ |
|-------------|--------------------------|
| 0.0000 | 1.0000 |
| 0.0001 | 1.0001 |
| 0.0002 | 1.00020001 |
| 0.0003 | 1.00030003 |
| ... | ... |
| 0.9999 | 2.71787413941... |
| 1.0000 | 2.71814592682... |

這樣的表格同時具備二項優點：對數和真數皆呈均勻遞增，計算上只需重複乘以 $1.0001$。

Bürgi 正是採用該方法。由於每筆新真數 $(1.0001)^{N}$ 可由前一筆 $(1.0001)^{N-1}$ 乘以 $1.0001$ 得到，填滿 $10000$ 筆表格僅需約 $10000$ 次乘法。然而每次乘法涉及多位數的手工運算，且需維持高精度以避免誤差累積，工作量仍然極為可觀。這樣的設計讓早期的對數表建構者得以確保高精度，進而利用以下換底公式推導出常用對數：
$$
\log_{10}(x) = \frac{\log_a x}{\log_a 10}
$$

其中 $\log_a(10)$ 的值 $L$ 滿足 $a^L = 10$。這個 $L$ 可以透過類似二進位展開的方式 $L = \sum_{k} c_k 2^k$ (這裡 $c_k$ 為 $0$ 或 $1$) 來逼近，而計算 $a^{2^k}$ 只需要對 $a$ 反覆進行平方運算，計算 $a^{1/2^k}$ 則需要反覆開平方根。透過比較 $a^L$ 與 10 的大小來逐步確定係數 $c_k$，可將 $L$ 計算到所需精度。

換言之，構造高精度對數表的訣竅，在於選擇形如 $b^N$ 的底數，其中 $b$ 趨近 $1$ 而 $N$ 足夠大 (如 $10000$ 或以上)，如此方能兼顧計算便利與精度要求。

接著用現代語言分析。設底數 $a = 1 \pm c$，其中 $c$ 為極小正數 (例如 $c = 10^{-7}$)。由二項式展開 $(1 \pm c)^k \approx 1 \pm kc$ (忽略 $c^2$ 以上高階項)，等比數列 $(1 \pm c)^0, (1 \pm c)^1, (1 \pm c)^2, \ldots$ 的近似值恰好構成等差數列 $1, 1 \pm c, 1 \pm 2c, \ldots$，正是對數「將乘法化為加法」的關鍵想法。

為找出此映射的極限形式，設 $x = (1 \pm c)^b$ 並令 $b = y/c$，使 $y$ 的單位變化對應 $b$ 的 $1/c$ 變化。分別代入二種底數：
- $a = 1 + c$: $x = (1+c)^{y/c} = \bigl[(1+c)^{1/c}\bigr]^y \to e^y$，反函數 $y = \ln x$
- $a = 1 - c$: $x = (1-c)^{y/c} = \bigl[(1-c)^{1/c}\bigr]^y \to e^{-y}$，反函數 $y = -\ln x$

(此處利用 $\lim_{c \to 0^+}(1+c)^{1/c} = e$ 及 $\lim_{c \to 0^+}(1-c)^{1/c} = e^{-1}$。) 以 $c = 10^{-7}$ 驗證，近似函數幾乎與 $\ln x$ 重合：

> 圖一: $x = 1.0000001^{10^7 y}$ 與 $y = \ln{x}$ 的比較

> 圖二: $x = 0.9999999^{10^7 y}$ 與 $y = -\ln{x}$ 的比較

直觀理解：設 $y = f(x)$ 為上述映射，以 $a = 1 + c$ 為例，取差商可得 $f'(x) \approx (1+c)/x \to 1/x$，積分並取 $f(1) = 0$ 即得 $f(x) = \ln x$；$a = 1 - c$ 分支的差商為 $-(1-c)/x \to -1/x$，對應 $f(x) = -\ln x$。(此為啟發式論證，嚴格推導需處理極限交換條件。)

上述推導顯示，當底數趨近 $1$ 時，「將幾何級數近似為算術級數」的出發點自然引出 $\ln x$ 作為反函數。Napier 在微積分尚未成形的時代，獨力以二十年的反覆運算與觀察，在建構對數表的過程中觸及 $e$ 的極限形式，這也是 $e$ 有時稱為「Napier 常數」的由來 (儘管不甚精確)。他過世後出版的《De arte logistica》(直譯為「計算之術」) 展現其關注如何透過數學工具實際提升計算效率，與現代數學強調結構與理論的取向相當不同。此後，介於 Napier 與牛頓、萊布尼茲之間的數學家，逐步將 $e$ 納入連續變化與極限計算的架構中，使其從早期的計算工具演變為現代分析、機率與數值方法的常數。
> 延伸閱讀: [對數與約翰．納皮爾 (John Napier)](https://episte.math.ntu.edu.tw/articles/mm/mm_03_4_07/)

### 從雙曲線面積到自然底數
Napier 與 Bürgi 的工作奠定對數表的實用基礎，但 $e$ 作為常數的理論地位，是在後續一個世紀中由多位數學家逐步確立的。這條路線從雙曲線 $y = 1/x$ 的面積問題開始，經惠更斯、Mercator、Johann Bernoulli，最終由 Euler 收束。

費馬 (Fermat) 在公元 1636 年之前就已推導出
$$
\int_{1}^{a} x^n \, dx \;=\; \frac{a^{n+1}-1}{n+1}, \quad n \neq -1
$$

於是人們自然聯想到：
$$
\int_{1}^{a} \frac{1}{x} \, dx = \;? \tag{1}
$$

在公元 1649 年以前，二位耶穌會會士 Grégoire de Saint-Vincent 與 Alphonse Antonio de Sarasa 研究雙曲線 $xy = 1$ 的求積法，藉由計算雙曲扇形的面積，他們推導出後來稱為[雙曲對數](https://en.wikipedia.org/wiki/Hyperbolic_functions) (hyperbolic logarithm) 的函數，該函數的性質與現今的自然對數相同。這二位耶穌會會士發現，雙曲線 $y=1/x$ 下方的面積 $A(a) = \int_1^a \frac{1}{t} dt$ 具有對數的關鍵性質 $A(ab) = A(a) + A(b)$，也就是說，$y = \tfrac{1}{x}$ 曲線下方的面積和 $y$ 的對數之間呈正比。

惠更斯 (Christiaan Huygens) 是 17 世紀科學革命的重要人物，橫跨數學、物理、天文、工程與哲學等多個領域，堪稱一代巨擘。受到 Grégoire de Saint-Vincent 的影響，惠更斯在探索等軸雙曲線的過程，串起 Napier 發明的對數與 Grégoire de Saint-Vincent 發現的幾何結構之間的關聯。
> 作為物理學家，惠更斯建立向心力定律，提出動量守恆原理與光學中的惠更斯原理；作為天文學家，他發現土衛六 (Titan)、獵戶座大星雲、火星極冠與土星環的結構；作為工程師，他改良望遠鏡鏡片設計，發明精密擺鐘與火藥引擎 (內燃機的雛形)；作為數學家，他發展漸屈線理論，提出單擺周期公式，並被視為機率論的早期奠基者。

在上述對等軸雙曲線 $xy = 1$，惠更斯注意到幾何性質：若取曲線下某區間的面積，只要滿足兩區間的比值關係一致，其對應面積亦相等。若區間 $[p, q]$ 與 $[r, s]$ 滿足 $\frac{p}{q} = \frac{r}{s}$，則二者在曲線 $y = \frac{1}{x}$ 下方所圍成的面積相同。

等軸雙曲線 $y = \frac{1}{x}$ 的圖例中，斜線標注的區域表示從 $x = 1$ 到 $x = a$ 所圍成的面積，即為 $\ln(a)$，若 $a < 1$，則面積為負值，表示 $\ln(a) < 0$。曲線 $y = \frac{1}{x}$ 下從 $x = 1$ 到 $x = a$ 的面積正好等於 $\ln(a)$，也就是自然對數的定義：
- 函數 $f(x) = \frac{1}{x}$
- 面積定義：$\ln(a) = \int_1^a \frac{1}{x} \, dx$
- 若定義 $A(x) = \int_1^x \frac{1}{t} \, dt$，則滿足對數的加法性質 $A(xy) = A(x) + A(y)$

這類滿足乘積轉加法規律的函數，即為對數函數。惠更斯從此幾何積分觀點出發，辨識出以 $e$ 為底的自然對數與常用十進位對數之間的本質差異。他也是最早以高精度計算 $e$ 的數學家之一，將其精確至小數點後第 17 位。

關於自然對數的最早記錄，出現在 Nicholas Mercator 於公元 1668 年出版的《Logarithmotechnia》一書，推動自然對數的形式化和應用。在許多實驗與自然現象中，被觀察對象的變化率 (增減速率) 往往與其當下的數量成正比，例如前述的複利和投資獲利，而所謂「自然對數」中的「自然」(natural) 一詞，並非因例子取於自然界，而是數學家在建構這些模型時，發現以 $e$ 為底的對數，能讓公式最簡潔、最直觀地表達出這種成長或衰減的特性。因此，這個「自然」是指在數學表示中最自然、無須多餘調整的選擇，讀者不用糾結於「自然」一詞和大自然的直接對應，一如[有理數](https://en.wikipedia.org/wiki/Rational_number) (rational number) 著重於「比例」，而非「有理」與否。
> 延伸閱讀: [資訊科技詞彙翻譯](https://hackmd.io/@sysprog/it-vocabulary)

Johann Bernoulli (Jacob Bernoulli 之弟) 時代，積分問題進一步推廣為下式 (2)：
$$
\int \frac{dx}{a x^2 + b x + c}, \tag{2}
$$

顯然可透過配方與換元法，將其化為與式 (1) 類似的形式。關鍵在於如何解方程 $a x^2 + b x + c = 0$，其解可能會是複數。

Johann Bernoulli 研究 $\lim_{n \to \infty} \Bigl(1 + \frac{1}{n}\Bigr)^n$ 的性質，認為這是個重要常數，他的學生 Euler 也深受啟發。

另一方面，在 Johann Bernoulli 求解下式 (3) 這個特例時：
$$
\int \frac{dx}{b^2 + x^2} \tag{3}
$$

發現透過巧妙的代換 (如涉及複數 $i$ 和變數 $t$ 的分式代換)，可將此積分與對數函數聯繫起來，進而在三角函數 ($\arctan$) 與複數對數之間建立關聯，也就是：
$$
\frac{-\,dt}{\sqrt{-1}\,\cdot 2\,b\,t} \tag{4}
$$

式 (3) 的積分結果是 $\arctan$，而式 (4) 的結果是含虛數的對數，進而在三角函數與複數對數之間建立關聯。需要注意的是，Johann Bernoulli 對複數的理解仍停留在 Gerolamo Cardano 時期的水準，計算過程中常會避免複數的出現。

公元 1740 年，Euler 發現
$$
\cos z = \frac{e^{\,\mathrm{i}z} + e^{-\mathrm{i}z}}{2}
$$
這個指數形式滿足與 $\cos z$ 相同的基本性質 (或微分方程)，因此確立彼此之間的等價關係。接著，Euler 在 1743 年指出三角函數與複指數的關係：
$$
\cos x = \frac{e^{\,\mathrm{i}x} + e^{-\mathrm{i}x}}{2}, \quad
\sin x = \frac{e^{\,\mathrm{i}x} - e^{-\mathrm{i}x}}{2\,\mathrm{i}}
$$

並在公元 1748 年最終提出著名的[歐拉公式](https://en.wikipedia.org/wiki/Euler%27s_formula):
$$
e^{\,\mathrm{i}x} = \cos x + \mathrm{i}\,\sin x
$$

從 Johann Bernoulli 與 Euler 的研究中，我們不難發現，冪級數展開對許多函數而言都至關重要。此技術在 Brook Taylor (1715) 與 Colin Maclaurin (1742) 的著作中獲得系統化，成為後續分析學的基本工具。

延伸閱讀:
* [自然對數漫談](https://www.math.sinica.edu.tw/media/pdf/d231/23111.pdf)
* [指數與對數發展史簡介 (上)](https://www.sec.ntnu.edu.tw/uploads/asset/data/66cf6cde0e0b30672f4bad39/1983-056-10_58-73_.pdf)
* [指數與對數發展史簡介 (下)](https://www.sec.ntnu.edu.tw/uploads/asset/data/66cf6a940e0b305152a537b2/1983-057-08_32-50_.pdf)

從現代觀點看 $e$ 的定義：$e^x$ 是導數等於自身的指數函數，即 $f(x) = e^x \Rightarrow f'(x) = f(x)$。在探究指數函數 (例如 $2^x$) 的導數時，數學家發現一個極為特殊的底數：所對應的指數函數在 $x = 0$ 處的切線斜率恰為 $1$，導數在任意點皆與函數值相等，滿足以下微分方程：
$$
f'(x) = f(x), \quad f(0) = 1
$$

為理解這個底數的特殊性，我們比較常見底數如 $2$ 與 $3$ 所構成的指數函數：
- 對於 $y = 2^x$，其在 $x = 0$ 處的導數約為 $0.7$
- 對於 $y = 3^x$，其在 $x = 0$ 處的導數約為 $1.1$

二者的變化率皆與其函數值不一致，顯示無法滿足 $f'(x) = f(x)$ 的條件。然而，存在一個介於 $2$ 與 $3$ 之間的底數 $e$，使得：
$$
\frac{d}{dx} e^x = e^x
$$

> 圖片取自 [Stewart Calculus Textbooks and Online Course Materials](https://www.stewartcalculus.com/)

亦即，函數 $y = e^x$ 是唯一導數恆等於自身的指數函數。特別是在 $x = 0$ 時，其切線斜率為 $1$，與函數值一致，展現出極致的對稱與自然性。

該結論可從指數函數的一般導數公式中推得。考慮 $f(x) = a^x$，其導數為：
$$
f'(x) = \lim_{\Delta x \to 0} \frac{a^{x + \Delta x} - a^x}{\Delta x}
= a^x \cdot \lim_{\Delta x \to 0} \frac{a^{\Delta x} - 1}{\Delta x}
$$

這顯示 $a^x$ 的導數總是等於其自身乘上一個常數，該常數與 $a$ 有關，與 $x$ 無關。為了使導數與函數值完全一致，我們自然會希望此常數為 $1$，亦即：
$$
\lim_{\Delta x \to 0} \frac{a^{\Delta x} - 1}{\Delta x} = 1
$$

唯一滿足此條件的底數就是 $e$，也可透過另一個等價的極限來定義：
$$
e = \lim_{\Delta x \to 0} (1 + \Delta x)^{1 / \Delta x} \approx 2.71828
$$

因此，$e^x$ 的導數為：
$$
\frac{d}{dx} e^x = 1 \cdot e^x
$$

該特性使 $e^x$ 在數學分析與微積分不可或缺。接著，對任意底數 $a > 0$，可藉由換底公式將其轉寫為以 $e$ 為底的形式：
$$
a^x = e^{x \ln a}
$$

應用鏈式法則可得其導數為：
$$
\frac{d}{dx} a^x = \frac{d}{dx} \bigl(e^{x \ln a} \bigr) = \ln a \cdot e^{x \ln a} = \ln a \cdot a^x
$$

因此，任意指數函數 $a^x$ 的導數皆為其自身乘上一個比例因子 $\ln a$。部分常見底數的對應值如下：
$$
\begin{aligned}
\frac{d}{dx} 2^x &= \ln(2) \cdot 2^x \approx 0.6931 \cdot 2^x \\
\frac{d}{dx} 3^x &= \ln(3) \cdot 3^x \approx 1.0986 \cdot 3^x \\
\frac{d}{dx} 5^x &= \ln(5) \cdot 5^x \approx 1.6094 \cdot 5^x \\
\frac{d}{dx} 10^x &= \ln(10) \cdot 10^x \approx 2.3025 \cdot 10^x
\end{aligned}
$$

換言之，所有指數函數皆可透過自然對數 $\ln a$ 表現其導數比例，而 $e$ 的特殊地位正來自 $\ln e = 1$，使得 $e^x$ 是唯一導數與積分形式最簡潔的底數。

延伸閱讀:
* [Earliest Uses of Symbols for Constants](https://jeff560.tripod.com/constants.html)
* [Logarithms: The Early History of a Familiar Function](https://old.maa.org/press/periodicals/convergence/logarithms-the-early-history-of-a-familiar-function-introduction)

## 證明極限存在
前述已從歷史脈絡揭示 $e$ 的來源與其作為自然底數的地位。然而，$e = \lim_{n\to\infty}(1+1/n)^n$ 這個極限是否真的收斂，仍需嚴格的分析學論證。證明的策略在於利用「單調有界數列收斂原理」：先證明數列 $a_n = (1+1/n)^n$ 單調遞增，再證明其存在上界。

### 步驟一：證明數列單調遞增
首先，我們利用廣義二項式定理或更簡單的代數不等式。對於任意 $b > a > 0$，考慮以下恆等式：
$$
b^n - a^n = (b - a)(b^{n-1} + b^{n-2}a + \cdots + a^{n-1})
$$

由於右側括號內共有 $n$ 項，且每項均滿足 $b^{n-k}a^{k-1} < b^{n-1}$ (因為 $a < b$)，故有：
$$
b^n - a^n < (b - a) \cdot n \cdot b^{n-1}
$$

整理可得重要不等式：
$$
a^n > b^{n-1} \bigl[ b - n(b - a) \bigr] \tag{5}
$$

現在，為了構造數列 $a_n$ 的遞迴關係，我們令：
$$
a = 1 + \frac{1}{n}, \quad b = 1 + \frac{1}{n - 1} \quad (\text{對 } n > 1)
$$

代入不等式 (5)，計算括號內各項：
$$
b - a = \left(1 + \frac{1}{n-1}\right) - \left(1 + \frac{1}{n}\right) = \frac{1}{n(n-1)}
$$
$$
b - n(b - a) = \left(1 + \frac{1}{n-1}\right) - n \cdot \frac{1}{n(n-1)} = 1 + \frac{1}{n-1} - \frac{1}{n-1} = 1
$$

代入 (5) 式後得到：
$$
\left(1 + \frac{1}{n}\right)^n > \left(1 + \frac{1}{n-1}\right)^{n-1} \cdot 1 = \left(1 + \frac{1}{n-1}\right)^{n-1}
$$

這證明數列 $x_n = (1+1/n)^n$ 嚴格單調遞增。

### 步驟二：證明數列有上界
接著，我們證明該數列不會發散至無窮大。再次利用不等式 (5)，但選取不同的參數。令：
$$
a = 1, \quad b = 1 + \frac{1}{2(n-1)} \quad (\text{對 } n > 1)
$$

代入 (5) 式。先計算括號內的值：
$$
b - a = \frac{1}{2(n-1)}, \quad n(b - a) = \frac{n}{2(n-1)}
$$
$$
b - n(b - a) = 1 + \frac{1}{2(n-1)} - \frac{n}{2(n-1)} = 1 + \frac{1 - n}{2(n-1)} = 1 - \frac{1}{2} = \frac{1}{2}
$$

因此不等式 (5) 給出：
$$
1 > \left(1 + \frac{1}{2(n-1)}\right)^{n-1} \cdot \frac{1}{2}
$$

整理得：
$$
\left(1 + \frac{1}{2(n-1)}\right)^{n-1} < 2
$$

令 $k = n - 1$ (其中 $k \geq 1$)，上式改寫為 $\left(1 + \frac{1}{2k}\right)^{k} < 2$。兩邊平方：
$$
\left(1 + \frac{1}{2k}\right)^{2k} < 4
$$

注意到左側恰為 $x_{2k} = \left(1 + \frac{1}{2k}\right)^{2k}$，因此偶數項子序列 $x_2, x_4, x_6, \ldots$ 皆小於 $4$。又因數列 $\{x_m\}$ 單調遞增，對於任何 $m$，總能找到偶數 $2k > m$，使得 $x_m < x_{2k} < 4$，故數列被 $4$ 所界定。

綜合上述結果，由單調有界原理可知極限 $\lim_{n \to \infty} (1 + 1/n)^n$ 必然存在。我們將此極限定義為歐拉數 $e$。歐拉隨後證明此極限與無窮級數 $\sum 1/n!$ 等價，後者收斂速度極快，為 $e$ 的精確計算奠定基礎。

## 計算歐拉數
只要有極限，任何運算都能用加減乘除處理，就算如下方這樣的普通計算器 (calculator，並非運算能力更強的電腦 [computer])：

其中，自然對數特別容易計算，接下來我們將介紹一種方法，主要靠計算器「開平方」功能。

自然對數 $\ln(1+x)$ 的泰勒展開如下：
$$
\ln(1 + x) = x - \frac{x^2}{2} + \frac{x^3}{3} - \frac{x^4}{4} + \cdots \quad (|x| < 1)
$$

對於普通計算器來說，實在難以計算上方的無窮級數。我們換個思路：回顧函數 $y = a^x$ (其中 $a > 0$ 且 $a \neq 1$) 的導數的推導過程：
$$
\frac{d}{dx} a^x = \lim_{\Delta x \to 0} \frac{a^{x + \Delta x} - a^x}{\Delta x} = \lim_{\Delta x \to 0} a^x \frac{a^{\Delta x} - 1}{\Delta x} = a^x \lim_{\Delta x \to 0}\frac{a^{\Delta x} - 1}{\Delta x} \tag{6}
$$

可見，最後面那一項極限和 $x$ 無關，而只跟 $a$ 有關，也就是個常數。該常數是多少呢？只要令 $x = 0$ 就能得到：
$$
\left. \frac{d}{dx} a^x \right|_{x=0} = \lim_{\Delta x \to 0} \frac{a^{\Delta x} - 1}{\Delta x}
$$

根據指數函數的圖形：

由式 (6) 可知，$a^x$ 的導數等於 $a^x$ 乘以一個僅取決於 $a$ 的常數 $\lim_{\Delta x \to 0} \frac{a^{\Delta x} - 1}{\Delta x}$。若令該常數恰為 $1$，則導數處處等於函數本身，滿足此條件的唯一底數即定義為 $e$。依據定義，我們知道：
$$
\frac{d}{dx} a^x = \frac{d}{dx} \bigl(e^{(\ln a) x}\bigr) = e^{(\ln{a}) x} \ln{a} = a^x \ln{a}
$$

和式 (6) 比較，可得：
$$
\ln{a} = \lim_{\Delta x \to 0} \frac{a^{\Delta x} - 1}{\Delta x} \tag{7}
$$

式 (7) 的左邊是對數，右邊是乘冪，離我們設定的目標仍很遠，接下來的關鍵想法是把指數式轉換成根式。
$$
\begin{aligned}
\ln{a} &= \lim_{x \to 0} \frac{a^x - 1}{x} \\
&= \lim_{k \to +\infty} \frac{a^{2^{-k}} - 1}{2^{-k}} \\
&= \lim_{k \to +\infty} \Bigl(a^{\frac{1}{2^k}} - 1\Bigr) \, 2^k \\
&= \lim_{k \to +\infty} \Bigl(\sqrt[2^k]{a} - 1\Bigr) \, 2^k \\
&= \lim_{k \to +\infty} \Bigl(\underbrace{\sqrt{\cdots \sqrt{a}}}_{k \text{times}} - 1\Bigr)\, 2^k
\quad (k \in \mathbb{N}).
\end{aligned}
$$

也就是把 $a^x$ 替換成一連串的「反覆開平方」，最後再用 $2^k$ 這個因子修正。對於只有開平方功能的普通計算器，這種方法極易上手，只要反覆開平方就能逼近結果。

例如，取 $k = 8$ 來估算 $\ln(2)$，我們可得：

計算器顯示 $0.694086413$，而事實上 $\ln(2) \approx 0.693147$，誤差約為 $0.000939413$。

若改算 $\ln(3)$，則約為 $1.100972987$，實際值為 $1.09861$，誤差約 $0.002363$。

可發現，當 $x$ 越接近 $1$ 時，這種方法計算 $\ln{x}$ 的誤差越小。若把 $y = \ln{x}$ 與 $y = (\sqrt[2^k]{x} - 1) \times 2^k$ 繪製在同一坐標系上，能明顯看到二者密切貼合。

至此，只要我們懂得利用反覆開平方，再結合極限的概念，就能在普通計算器近似計算歐拉數和自然對數，精度隨開平方次數 $k$ 的增大而提升。

### 在 Minecraft 中近似計算歐拉數
計算歐拉數 $e$ 的方法多樣，除了傳統的數學分析，我們甚至可在風靡全球的遊戲 [Minecraft](https://www.minecraft.net/) 中，透過模擬機率實驗來近似 $e$ 的值。正如 Molly Lynch 與 Michael Weselcouch 在論文〈[Approximating Mathematical Constants using Minecraft](https://arxiv.org/pdf/2411.18464)〉中所展示，我們可利用 Minecraft 遊戲內的機制來模擬隨機排列，並藉此估算 $e$。

該方法的數學基礎源自組合數學中的[錯位排列](https://en.wikipedia.org/wiki/Derangement) (derangement)：若一個排列中，所有元素皆未出現在其原本的位置上，則稱為錯位排列。例如，排列 `312` 是 `123` 的一個錯位排列，因為其中 1, 2, 3 都未位於原始位置；而 `321` 則不是，因為數字 2 仍然處於原來的第二個位置。其關鍵的特性是，當元素數量 $n$ 趨近於無限大時，一個隨機產生的 $n$ 元素排列是錯位排列的機率，會趨近於 $1/e$，亦即：
$$
\lim_{n \to \infty} \frac{!n}{n!} = \frac{1}{e}
$$
其中 $!n$ 代表 $n$ 個元素的錯位排列總數，$n!$ 代表 $n$ 個元素的全部排列總數。

根據這個極限關係，我們可反過來利用機率實驗近似 $e$：只要產生大量 ($N$ 次) 的隨機排列，並計算其中錯位排列出現的次數 ($D$ 次)，則它們的比率 $D/N$ 會近似於 $1/e$。因此：
$$
e \approx \frac{N}{D} = \frac{\text{隨機排列的總試驗次數}}{\text{觀察到的錯位排列次數}}
$$

我們該如何在 Minecraft 中模擬上述呢？可用遊戲中的「紅石」(Redstone) 系統來搭建自動化裝置：
1. 產生隨機排列：使用「[發射器](https://minecraft.wiki/w/Dropper)」(Dropper)。發射器內可放置最多 9 種不同的、可堆疊數量為 1 的物品 (例如 9 種不同顏色的羊毛，或 9 種不同的唱片)，每種物品代表集合 $[9] = \{1, 2, \dots, 9\}$ 中的一個數字。當發射器被觸發時，它會隨機射出一件物品。連續觸發 9 次，記錄下物品射出的順序，就構成一次對 $[9]$ 元素的隨機排列。
2. 判斷是否為錯位排列：這需要設計[紅石電路](https://minecraft.wiki/w/Redstone_circuits)。我們可設定一個系統，記錄第 $k$ 次 ($k$ 從 1 到 9) 噴出的物品是代表數字 $j$ 的物品。如果對於所有的 $k$ (從 1 到 9)，都有 $j \neq k$ (即第 $k$ 次噴出的物品所代表的數字不等於 $k$)，於是該次產生的排列就是個錯位排列。這可運用比較器 (comparator) 和邏輯閘 (logic gates) 來實作。
3. 計數與統計：利用「[漏斗](https://minecraft.wiki/w/Tutorial:Hopper)」(Hopper) 和儲存裝置 (如「[箱子](https://minecraft.wiki/w/Chest)」) 來計數。設計二組計數器，一組記錄總共產生多少次排列 ($N$)，另一組則只在偵測到錯位排列時才計數 ($D$)。
4. 自動化迴圈：將上述步驟組合成可自動重複運行的系統，讓它可長時間、大量地產生排列並進行統計，以獲得足夠大的樣本數，從而提高近似的準確度。


在論文提及的實驗中，產生 647 次排列，其中有 238 次是錯位排列。根據上述公式，得到的 $e$ 的近似值為：
$$
e \approx \frac{647}{238} \approx 2.71849
$$
誤差僅約 0.0076%。不僅展示 $e$ 與錯位排列機率之間的關聯，也說明 Minecraft 作為計算模擬與實驗平台的潛力，將抽象的數學概念轉化為遊戲中高度視覺化的趣味過程。
展示影片: [Approximation of e in Minecraft](https://youtu.be/YH-wx5YQftY) 及[播放清單](https://www.youtube.com/playlist?list=PLscpLh9rN1Rf-dqLO4r3GAwvm1O_xL7D1)。

### 以 binary splitting 計算 $e$ 精確至百萬位數
為高精度計算自然對數底 $e$，我們採用下列無窮級數：
$$
e = \sum_{k=0}^{\infty} \frac{1}{k!}
$$

該級數具有快速收斂性，若期望小數點後達 $d$ 位正確，我們需確保截尾誤差 $\sum_{k>n} 1/k! < 10^{-d}$。由於尾和可用等比級數估計 $\sum_{k>n} 1/k! < \frac{1}{n!} \cdot \frac{1}{n}$ (對 $n \gg 1$)，以 $1/n! < 10^{-d}$ 作為充分條件即可。因此該找出最小整數 $n$，使得：
$$
\frac{1}{n!} < 10^{-d} \quad\Rightarrow\quad \log_{10}(n!) > d
$$

為估算 $n!$ 的對數，可用 Stirling 近似公式：
$$
\log_{10}(n!) \approx \frac{1}{\ln 10} \left(n \ln n - n + \frac{1}{2}\ln(2\pi n)\right)
$$

藉由二分搜尋，可有效找到滿足上述條件的最小 $n$，將此作為級數展開的項數上限。以 $10^6$ 位精度為例，計算所需項數約為 $205{,}023$。我們採取 [GMP](https://gmplib.org/) 多精度整數 (`mpz_t`) 與浮點數 (`mpf_t`) 為基礎進行實作，主要流程如下：
* 使用遞迴函式 `calc(mpz_t a, mpz_t b, mpz_t p, mpz_t q)` 以 [binary splitting](https://en.wikipedia.org/wiki/Binary_splitting) 方式計算級數部分和 (其中 `p`, `q` 為輸出參數)。此函式維護以下不變量：
- 不變量：對 $a > 0$ 的區間 $[a, b)$，輸出滿足 $q = b!/a!$，$p/q = \sum_{k=a+1}^{b} a!/k!$
- 基底情況 ($a > 0$, $b = a+1$)：回傳 $p=1$, $q=b$，此時 $q = b!/a! = b$，$p/q = 1/b$
- 基底情況 ($a=0$, $b=1$)：回傳 $p=2$, $q=1$，將 $1/0! + 1/1! = 2$ 合併為一個葉節點
- 驗證：以 `calc(0, 3)` 為例。基底：$[0,1)$ → $p=2, q=1$；$[1,2)$ → $p=1, q=2$；$[2,3)$ → $p=1, q=3$。合併 $[1,3)$：$p=1 \cdot 3 + 1 = 4$, $q=2 \cdot 3=6$。合併 $[0,3)$：$p=2 \cdot 6 + 4 = 16$, $q=1 \cdot 6=6$，故 $p/q = 8/3$，正好等於 $\sum_{k=0}^{3} 1/k! = 1 + 1 + \frac{1}{2} + \frac{1}{6} = \frac{8}{3}$
* 遞迴分解：將區間 $[a, b)$ 對半分割為 $[a, m)$ 與 $[m, b)$，其中 $m = \lfloor(a+b)/2\rfloor$
* 遞迴計算：
* 呼叫 `calc(a, m, p_L, q_L)` 計算左半 (程式碼中為 `p`, `q1`)
* 呼叫 `calc(m, b, p_R, q_R)` 計算右半 (程式碼中為 `p2`, `q`)
* 合併結果：由於 $q_L = m!/a!$ 與 $q_R = b!/m!$，合併後 $q_{\text{new}} = q_L \cdot q_R = b!/a!$。分子方面，左半的部分和需乘以 $q_R$ 將分母統一至 $b!/a!$，再加上右半分子：$p_{\text{new}} = p_L \cdot q_R + p_R$ (對應程式碼 `mpz_mul(p, p, q); mpz_add(p, p, p2); mpz_mul(q, q, q1)`)
* 精度保持：所有中間步驟皆以 GMP 大整數 (`mpz_t`) 執行精確的有理數運算，僅在最終計算完成後，才用 `mpf_div()` 將總和的 $p/q$ 轉換為高精度浮點數 (`mpf_t`)，從而減少浮點捨入誤差的累積

GMP 浮點型別需指定位元精度，公式為：
$$
\text{bits} = \left\lceil d \cdot \log_2(10) \right\rceil \approx \left\lceil d \cdot 3.3219 \right\rceil
$$

因此，百萬位數需要約 3.3 百萬位元，即約 0.4 MB 記憶體。對一億位數則需約 41.5 MB。參考程式碼如下:

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <inttypes.h>
#include <gmp.h>
#include <math.h>
#include <time.h>
#include <stdint.h>

#define FILENAME "e.out"
#define PI 3.141592653589793238462643L

/* progress tracking */
typedef struct {
    uint64_t count, total;
    int8_t last_percentage;
    clock_t start_time;
} progress_t;

/* Global progress tracker */
static progress_t prog = {0, 0, -1, 0};

/* Compute log10(n!) using Stirling's approximation */
static long double log_factorial(uint64_t n)
{
    if (n <= 1) return 0;
    long double ln_n = logl(n);
    return (logl(2.0L * PI) / 2 + ln_n / 2 + n * (ln_n - 1)) / logl(10.0L);
}

/* Binary search for number of terms */
static uint64_t terms(uint64_t digits)
{
    uint64_t lower = 0, upper = 1;
    long double f;

    while ((f = log_factorial(upper)) <= digits) {
        lower = upper;
        upper <<= 1;
    }

    while (upper - lower > 1) {
        uint64_t n = lower + (upper - lower) / 2;
        f = log_factorial(n);
        if (f > digits) {
            upper = n;
        } else {
            lower = n;
        }
    }
    return upper;
}

static void write_e(mpf_t e, uint64_t digits)
{
    FILE *fp = fopen(FILENAME, "w");
    if (!fp) {
        fprintf(stderr, "Failed to open %s\n", FILENAME);
        return;
    }

    /* Write header */
    fprintf(fp, "Calculation of e to %" PRIu64 " digits\n\n", digits);
    fprintf(fp, "e = ");

    mp_exp_t exp_out;
    char *str = mpf_get_str(NULL, &exp_out, 10, digits + 1, e);
    if (!str) {
        fclose(fp);
        return;
    }

    fprintf(fp, "%c.", str[0]);
    for (uint64_t i = 1; i <= digits; i++) {
        fputc(str[i], fp);
        if (i % 10 == 0 && i != digits)
            fputc(' ', fp);
        if (i % 50 == 0) {
            fprintf(fp, "   [%6" PRIu64 "]\n", i);
            if (i < digits)
                fprintf(fp, "      ");
        }
    }
    if (digits % 50 != 0)
        fputc('\n', fp);

    free(str);
    fclose(fp);
}

static void show_progress_bar(void)
{
    if (prog.total == 0) return;

    int percentage = (int)(100.0 * prog.count / prog.total);
    if (percentage <= prog.last_percentage) return;

    prog.last_percentage = percentage;
    int width = 50;
    int filled = (int)(width * prog.count / prog.total);

    clock_t now = clock();
    double elapsed = (double)(now - prog.start_time) / CLOCKS_PER_SEC;
    double eta = (prog.total - prog.count) * (elapsed / (prog.count ? prog.count : 1));

    printf("\r[");
    for (int i = 0; i < width; i++) {
        printf(i < filled ? "█" : " ");
    }
    printf("] %d%% (ETA: %.0fs)", percentage, eta);
    fflush(stdout);

    if (prog.count == prog.total) printf("\n");
}

/* Binary splitting: p/q = sum_{k in [a,b)} 1/k! */
static void calc(mpz_t a, mpz_t b, mpz_t p, mpz_t q)
{
    assert(mpz_cmp(a, b) < 0);

    mpz_t m, p2, q1;
    mpz_inits(m, p2, q1, NULL);

    mpz_sub(m, b, a);
    if (mpz_cmp_ui(m, 1) == 0) {
        if (mpz_cmp_ui(a, 0) == 0) {
            mpz_set_ui(p, 2);  /* 1/0! + 1/1! */
            mpz_set_ui(q, 1);
        } else {
            mpz_set_ui(p, 1);
            mpz_set(q, b);
        }
    } else {
        mpz_add(m, a, b);
        mpz_fdiv_q_ui(m, m, 2);

        calc(a, m, p, q1);
        calc(m, b, p2, q);

        mpz_mul(p, p, q);
        mpz_add(p, p, p2);
        mpz_mul(q, q, q1);
    }

    prog.count++;
    show_progress_bar();
    mpz_clears(m, p2, q1, NULL);
}

/* Main calculation function for e */
static void calculate_e(uint64_t digits) {
    uint64_t d = digits + 1;
    uint64_t num_terms = terms(d);
    uint64_t precision = ceill(d * (logl(10.0L) / logl(2.0L))) + 1;

    prog.total = num_terms * 2 - 1;
    prog.start_time = clock();

    printf("Calculating e to %" PRIu64 " digits\n", digits);
    printf("Terms: %" PRIu64 ", Precision: %" PRIu64 " bits\n", num_terms, precision);

    mpf_set_default_prec(precision);
    printf("Actual precision: %lu bits\n", mpf_get_default_prec());

    mpz_t a, p, q, n;
    mpf_t pf, qf;
    mpz_inits(a, p, q, n, NULL);
    mpf_inits(pf, qf, NULL);

    mpz_set_ui(n, num_terms);
    calc(a, n, p, q);

    mpf_set_z(pf, p);
    mpf_set_z(qf, q);
    mpf_div(pf, pf, qf);

    printf("Writing result to %s...\n", FILENAME);
    write_e(pf, digits);

    mpz_clears(a, p, q, n, NULL);
    mpf_clears(pf, qf, NULL);
}

int main(int argc, const char *const argv[])
{
    uint64_t digits = 1000000ULL;
    const uint64_t precision_max = (sizeof(long) == 4) ? 
        1292913975ULL : 5553023288523357111ULL;

    if (argc == 2) {
        digits = strtoull(argv[1], NULL, 10);
        if (digits == 0 || digits > precision_max) {
            fprintf(stderr, "Digits must be 1 to %" PRIu64 "\n", precision_max);
            return 1;
        }
    } else if (argc > 2) {
        fprintf(stderr, "Usage: %s [#digits]\n", argv[0]);
        return 1;
    }

    calculate_e(digits);
    return 0;
}

上述程式碼可高效率完成小數點百萬位的精確計算，對記憶體需求有限，且透過進度條顯示，最終輸出名為 `e.out` 的檔案:

e = 2.7182818284 5904523536 0287471352 6624977572 4709369995    [    50]
      9574966967 6277240766 3035354759 4571382178 5251664274    [   100]
      2746639193 2003059921 8174135966 2904357290 0334295260    [   150]
      5956307381 3232862794 3490763233 8298807531 9525101901    [   200]
...

讀者可指定更多有效位數進行運算。此外，由於 GMP 已針對大整數乘法實作 Karatsuba、Toom-Cook 與 FFT 最佳化，上方程式無須自行處理底層運算效率的議題。

延伸閱讀: [從尤拉數 e 到 Stirling 常數](https://episte.math.ntu.edu.tw/articles/mm/mm_20_1_03/)

## 快速計算指數函數
前節探討如何精確計算常數 $e$ 的數值 (一個定值，可離線運算任意位數)；本節則聚焦於在即時系統中快速近似整個函數 $e^x$，二者的設計取捨截然不同。

標準函式庫中的 [exp](https://man7.org/linux/man-pages/man3/exp.3.html) 通常能達到接近 1 [ULP](https://en.wikipedia.org/wiki/Unit_in_the_last_place) 的精確度，然而許多應用情境，如即時系統、訊號處理、機器學習或遊戲引擎等，通常更注重運算時間，且能接受略低的精度，於是就會運用標準函式庫以外的指數函數實作。

〈[留意浮點數運算的陷阱](https://hackmd.io/@sysprog/floating-point)〉以 $e$ 的運算為例，浮點數的捨入誤差不一定需要累積才能導致問題，關鍵在於識別誤差產生的關鍵環節。以下嘗試對效能與精確度進行權衡 (tradeoff)，適合於可容忍小誤差、但對速度敏感的情境，並從快速但低精度的版本開始，逐步改進到更精確的實作，探討其中的數學背景及技巧。

為快速近似計算指數函數 $e^x$，我們用以下數學恒等式：
$$
e^x = 2^{x \cdot \log_2(e)}
$$

透過此恒等式，我們將指數函數從自然底數 $e$ 轉換到以 2 為底的計算，更適合電腦的二進位架構。

### 基礎實作
將輸入 $x$ 分成整數部分 $i$ 與小數部分 $f$：
$$
2^{i+f} = 2^i \cdot 2^f
$$

其中
- 整數部分使用 IEEE-754 浮點格式的指數欄位直接進行位元操作
- 小數部分 $2^f$ 以低次多項式近似

接著，我們將 $t = x \cdot \log_2 e$ 拆成整數與小數二部分：
$$
t = i + f, \quad i = \lfloor t \rfloor,\, f = t - i
$$

於是：
$$
e^x = 2^t = 2^{i + f} = 2^i \cdot 2^f
$$

步驟 1：以低次多項式逼近 $2^f$

運用[極小極大逼近多項式](https://en.wikipedia.org/wiki/Minimax_approximation_algorithm) (minimax polynomial) 來估算 $2^f$。對於 $f \in [0,1)$ 的範圍，該手法可有效近似：
$$
2^f \approx P(f) = (a \cdot f + b) \cdot f + c
$$

這裡的係數經過調整，可將最大相對誤差控制在 $1.73 \times 10^{-3}$:
* $a = 0.3371894346$
* $b = 0.657636276$
* $c = 1.00172476$

步驟 2：透過位元操作計算 $2^i$

在 IEEE-754 單精度浮點數格式中，浮點數的記憶體配置如下：
> 符號位元 (1 個位元), 指數 (8 個位元), 尾數 (23 個位元)

一個 float 可視為：
$$
f = (-1)^s \cdot 1.m \cdot 2^{e - 127}
$$

我們要產生一個浮點數，滿足以下：
- 符號位元 (`s`) 為 0
- 指數部分 (`e`) 設為 $i + 127$
- 尾數部分 (`m`) 來自 $2^f$ 的估算值，而指數部分則透過位元運算設定

於是可透過將整數 $i$ 左移 23 個位元並加到原始 float 的位元表示中，達成對 $2^i$ 的計算，而不用乘法或呼叫 [pow](https://man7.org/linux/man-pages/man3/pow.3.html) 函式。

運算流程如下：
1. 計算 $t = x \cdot \log_2 e$
2. 拆成 $t = i + f$
3. 以多項式 $P(f)$ 估算 $2^f$
4. 藉由浮點數位元操作實作 $2^i$，得到 $2^t$

因此可得 $\exp(x) \approx P(f) \cdot 2^i$，並藉由位元運算以結合二者。

對應的 C 程式碼：

#include <math.h>
#include <stdint.h>
#include <string.h>

/* Fast approximate exp(x), x ∈ [-87, 88], max relative error ≤ 1.73e-3 */
float fast_exp(float x)
{
    const float LOG2_E = 1.442695041f;
    const float A = 0.3371894346f;
    const float B = 0.657636276f;
    const float C = 1.00172476f;

    float t = x * LOG2_E;
    float fi = floorf(t);
    float f = t - fi;
    int32_t i = (int32_t)fi;

    /* Polynomial approximation for 2^f */
    float exp2_frac = (A * f + B) * f + C;

    uint32_t bits;
    memcpy(&bits, &exp2_frac, sizeof(bits));
    bits += ((uint32_t)i << 23);

    float result;
    memcpy(&result, &bits, sizeof(result));
    return result;
}

不難發現，程式碼不依賴成本高昂的數學函式庫呼叫，如 [expf](https://man7.org/linux/man-pages/man3/exp.3.html) 和 [powf](https://man7.org/linux/man-pages/man3/pow.3.html)，唯一用到的函式呼叫是 [floorf](https://man7.org/linux/man-pages/man3/floor.3.html) 以達成浮點數對整數的轉換，隨後可換成處理器對應的指令，在 x86 的 SSE/AVX 或 Arm NEON 都存在這樣的指令。主要的運算由乘法、加法與位移運算構成。

### 進階實作
基礎實作針對單精度浮點數，進階的實作則針對倍精度 (`double`)，適合需要更高精度及較大範圍的科學或一般運算情境，可運用以下常見技巧：
* 使用多項式擬合 (polynomial fitting)
* 範圍縮減 (range reduction)
* 查表法 (lookup table)

以下逐項探討。

#### 多項式擬合

為了近似特定區域 (segment) 的數學函式，最常見的技巧是使用多項式來擬合 (polynomial fitting) 該函數的局部區間。我們的目標是找出適當的係數，以建構出近似於 $\exp(x)$ 的多項式函數。以二次多項式為例：
$$
f(x) = C_2 \cdot x^2 + C_1 \cdot x + C_0
$$

不過，由於多項式函數的形狀與 $\exp(x)$ 不同，我們僅能期望多項式在特定的小區間內能有效近似原函數。為了求得多項式的係數，我們可透過以下二種方式。

方法一：[SciPy](https://scipy.org/) 函式庫

利用 Python 的科學計算函式庫 [SciPy](https://scipy.org/) 中所提供的 [Levenberg–Marquardt algorithm](https://en.wikipedia.org/wiki/Levenberg%E2%80%93Marquardt_algorithm) 方法，透過非線性最小平方擬合 (nonlinear least-squares fitting) 來求得多項式係數。
> 利用三次多項式進行曲線擬合，參見 [exp.py](https://github.com/nadavrot/fast_log/blob/main/src/exp.py)

方法二：[Sollya](https://www.sollya.org/) 極小極大多項式

使用專門針對浮點數函式近似而設計的數值工具 [Sollya](https://www.sollya.org/)。Sollya 可產生[極小極大多項式](https://en.wikipedia.org/wiki/Minimax_approximation_algorithm)，此多項式的特性為最大化降低近似函數的「最大誤差」，理論上可提供最佳的近似效果。

以下用 Sollya 找出近似 $\exp(x)$ 的係數：

display = decimal; Q = fpminimax(exp(x), 6, [|D...|], [-0.000001, 1.001]); Q;
display = decimal; Q = fpminimax(exp(x), 5, [|D...|], [-0.0039, 0.0039]); Q;

上述 Sollya 腳本分別產生 6 次與 5 次的極小極大多項式，適用於不同區間寬度。實務中需根據目標精度與運算成本選擇合適的次數；下方 C 程式碼為追求速度，採用 3 次多項式 (由 SciPy 擬合取得)，其在區間 $[0, 1)$ 的最大誤差約為 $0.002$，離開該區間後誤差會迅速增加。

#### 範圍縮減
多項式近似只能保證在有限區間內的準確度，因此我們需要額外處理更廣的數值範圍。我們透過範圍縮減的技巧 (range reduction) 來處理這個問題。

利用指數函數的特性 $e^{a+b} = e^a \cdot e^b$，我們可將任意數值 $x$ 拆分為整數部分與小數部分。例如，將數字 3.1415 拆分為整數部分 $3$ 及小數部分 $0.1415$，亦即 $e^{3.1415} = e^{3} \cdot e^{0.1415}$。透過此種拆分方式，我們可分別使用二種方法進行處理：
* 整數部分 (如 3) ：透過預先建立的查詢表 (lookup table) 直接取得精確值
* 小數部分 (如 $0.1415$) ：利用上述的多項式擬合方式計算其近似值。

#### 建立查詢表
對於 IEEE-754 倍精度 (double precision) 的浮點數，$e^{709} \approx 8.22 \times 10^{307}$ 仍在可表示範圍內，但 $e^{710}$ 已超出 `DBL_MAX` 而溢位。因此，我們可建立一個範圍為 $[-709, 709]$ 的表格，以提供快速精確的 $\exp(i)$ 計算結果：

#define TABLE_MIN -709
#define TABLE_MAX 709
#define TABLE_SIZE (TABLE_MAX - TABLE_MIN + 1)

static double EXP_TABLE[TABLE_SIZE];
/* precomputed exp(i) */
__attribute__((constructor))
static void init_exp_table(void)
{
    for (int i = 0; i < TABLE_SIZE; ++i)
        EXP_TABLE[i] = exp((double)(i + TABLE_MIN));
}

注意 `__attribute__((constructor))` 是 GCC/Clang 的 C 語言擴充，指定 `init_exp_table` 函式在程式碼載入時期，早於 `main` 函式予以執行，若你不用 GCC 或 Clang 相容編譯器，應當在使用 `fast_exp` 函式前，先呼叫 `init_exp_table` 函式。

#### 處理負數問題的技巧
負數可藉由以下恒等式處理：
$$
e^{-x} = \frac{1}{e^x}
$$

然而，實務上若透過條件分支 (branch) 處理，可能造成 CPU 分支預測失誤 (branch misprediction)，導致效能下降。因此，更好的方式是將查表法範圍直接延伸至負數，如上述範例，從而避免分支指令的使用。

#### 誤差評估
若要評估多項式近似函數的表現，我們可將近似值與實際值相減，計算近似函數的絕對誤差 (absolute error) 或相對誤差 (relative error)，並將結果繪製成誤差圖 (error plot)，以直觀評估近似函數在整個區間內的表現。

以下 C 程式目標是在實務中提供足夠的精確度與效能： (沿用上方 `init_exp_table` 函式和 `TABLE_` 開頭的巨集)

#include <math.h>
#include <stdint.h>

/* Improved fast approximate exp(x) */
double fast_exp(double x)
{
    double i_part = floor(x);
    double f_part = x - i_part;

    /* Polynomial approximation for exp(f), f ∈ [0,1) */
    const double c0 = 0.28033708;
    const double c1 = 0.425302;
    const double c2 = 1.01273643;
    const double c3 = 1.00020947;

    /* Horner's method for polynomial evaluation */
    double poly = c3 + f_part * (c2 + f_part * (c1 + f_part * c0));

    if (i_part < TABLE_MIN) return 0.0;
    if (i_part > TABLE_MAX) return HUGE_VAL;
    int i = (int) i_part;

    return poly * EXP_TABLE[i - TABLE_MIN];
}

透過 [Horner 表示法](https://en.wikipedia.org/wiki/Horner%27s_method)，多項式計算效率大幅提升：
$$
a x^3 + b x^2 + c x + d = d + x(c + x(b + a x))
$$

此方法運算量少，且易於轉換成中央處理器的 FMA (fused multiply-add) 指令，Intel 與 AMD 處理器均支援 [FMA](https://en.wikipedia.org/wiki/FMA_instruction_set) 指令集。

延伸閱讀: [A Fast, Compact Approximation of the Exponential Function](https://nic.schraudolph.org/pubs/Schraudolph99.pdf)

### 針對無 FPU 嵌入式系統的數值方法
前述 `fast_exp` 實作依賴 IEEE-754 浮點數格式與位元操作，適用於具備 FPU (浮點運算單元) 的處理器 (如 Arm Cortex-M4F/M7)。然而，在更低功耗或低成本的嵌入式系統 (如 Arm Cortex-M0/M3) 中，通常缺乏 FPU，直接使用浮點運算會導致巨大的軟體模擬開銷。此時，工程師會採用定點數 (fixed-point) 或 [CORDIC](https://en.wikipedia.org/wiki/CORDIC) 演算法。

#### 定點數與查表內插
在定點數系統中，數值以整數表示，並預設一個固定的比例因子 (scaling factor，如 $2^{16}$)。計算 $e^x$ 時，常見手法是結合查表與線性內插。
$$
e^x \approx e^{x_k} + (x - x_k) \cdot e^{x_k}
$$
其中 $x_k$ 是查表中最接近 $x$ 的基準點。由於 $e^{x_k}$ 的導數仍是其自身，線性項的係數直接取自表中的函數值，節省額外的斜率儲存空間。

#### CORDIC 演算法
[CORDIC](https://en.wikipedia.org/wiki/CORDIC) (Coordinate Rotation Digital Computer) 是僅使用「位移」與「加減法」即可計算超越函數的強大工具。對於雙曲函數 (進而推導出 $\exp$ 與 $\ln$)，其遞迴式如下：
$$
\begin{aligned}
x_{i+1} &= x_i + \sigma_i y_i 2^{-i} \\
y_{i+1} &= y_i + \sigma_i x_i 2^{-i} \\
w_{i+1} &= w_i - \sigma_i \tanh^{-1}(2^{-i})
\end{aligned}
$$
其中 $\sigma_i$ 根據目標值 $w$ 的正負號決定，且迭代從 $i = 1$ 起算，因為 $\tanh^{-1}(2^0) = \tanh^{-1}(1)$ 發散。實務上雙曲模式的 CORDIC 還需在特定迭代 (如 $i = 4, 13, 40, \ldots$) 重複執行以確保收斂，並乘以相應的增益因子。由於 $2^{-i}$ 對應位元位移，而 $\tanh^{-1}(2^{-i})$ 可預先存入微小的查表，該方法在 FPGA 與早期微處理器中極其重要。

### 數值穩定性：expm1 與 log1p
前述的快速近似法著重於速度，但在數值分析中，另一類問題更為棘手：當 $x$ 極接近 $0$ 時，直接計算 $e^x - 1$ 會遭遇[災難性抵消](https://en.wikipedia.org/wiki/Catastrophic_cancellation) (catastrophic cancellation)。

原因在於 IEEE 754 浮點數的精度有限。以雙精度為例，$e^{10^{-15}} \approx 1.000000000000001$，而浮點數將其表示為 $1.0 + \epsilon$，其中 $\epsilon$ 僅剩約 1 位有效數字。若此時計算 $e^x - 1$，兩個幾乎相等的數值相減，大量有效位元被抵消，結果的相對誤差可達數個數量級。

標準函式庫提供 [expm1](https://man7.org/linux/man-pages/man3/expm1.3.html) 與 [log1p](https://man7.org/linux/man-pages/man3/log1p.3.html) 來解決此問題：

: `expm1(x)` 直接計算 $e^x - 1$，在 $|x| \ll 1$ 時利用泰勒展開 $e^x - 1 = x + \frac{x^2}{2} + \frac{x^3}{6} + \cdots$ 避免先算 $e^x$ 再減 $1$ 的抵消
: `log1p(x)` 直接計算 $\ln(1+x)$，在 $|x| \ll 1$ 時利用 $\ln(1+x) = x - \frac{x^2}{2} + \frac{x^3}{3} - \cdots$ 避免先算 $1+x$ 再取對數的精度損失

以下範例展示差異：

#include <math.h>
#include <stdio.h>

int main(void)
{
    double x = 1e-15;
    printf("exp(x)-1  = %.17e\n", exp(x) - 1);   /* 約 1.11e-15，誤差 ~11% */
    printf("expm1(x)  = %.17e\n", expm1(x));       /* 約 1.00e-15，正確 */
    printf("log(1+x)  = %.17e\n", log(1 + x));     /* 精度損失 */
    printf("log1p(x)  = %.17e\n", log1p(x));        /* 正確 */
    return 0;
}

這兩個函式在工程實務中出現的典型情境包括：
- 金融計算：連續複利 $FV = PV \cdot e^{rt}$ 中，短期利率 $r \cdot t \ll 1$ 時計算淨利息 $FV - PV = PV \cdot (e^{rt} - 1)$ 需要 `expm1`
- 機器學習：計算極小機率的交叉熵 $-\ln(q_i)$，當 $q_i \approx 1$ 時改寫為 $-\text{log1p}(q_i - 1)$ 以保持精度
- 統計學：增長率近似 $\ln(Z_t/Z_{t-1}) \approx X_t$ (詳見後文「增長率的近似與經濟學詮釋」)，在實作中即為 $log1p((Z_t - Z_{t-1}) / Z_{t-1})$
- 數值積分：微分方程求解器中的指數積分因子 $(e^{h\lambda} - 1)/\lambda$ 在 $h\lambda \to 0$ 時需要 `expm1`

從本文的主題來看，`expm1` 與 `log1p` 的存在正反映 $e^x$ 與 $\ln(1+x)$ 在 $x = 0$ 附近的泰勒展開結構：$e^x - 1 \approx x$ 且 $\ln(1+x) \approx x$，二者互為反函數的線性近似在此交會，而浮點運算的有限精度迫使我們必須直接利用級數而非組合運算。

## 歐拉數在機器學習中的角色
前節討論如何高效計算 $e^x$，而機器學習正是這類計算最密集的應用領域之一。現代神經網路的訓練與推論過程中，幾乎每個環節都仰賴以 $e$ 為主的函數。$e^x$ 的導數恆等於自身這項性質，使得反向傳播 (backpropagation) 中的梯度計算格外簡潔，不需額外的常數修正，連鎖律直接作用於函數本身。以下列舉三個關鍵函數。

### Sigmoid 函數
[Sigmoid 函數](https://en.wikipedia.org/wiki/Sigmoid_function)將任意實數映射至 $(0, 1)$ 區間，常用於[邏輯迴歸](https://en.wikipedia.org/wiki/Logistic_regression)與二元分類：
$$
\sigma(x) = \frac{1}{1 + e^{-x}}
$$

其導數可用函數自身表達：$\sigma'(x) = \sigma(x)(1 - \sigma(x))$，不需回頭查詢 $x$ 的原始值，反向傳播只要保留前向傳播的輸出即可完成梯度計算。這種「導數以自身表達」的特性，正是源自 $\frac{d}{dx} e^x = e^x$。

### Softmax 函數
[Softmax](https://en.wikipedia.org/wiki/Softmax_function) 將模型輸出的原始分數 (logits) 轉換為多類別的機率分佈：
$$
\text{softmax}(x_i) = \frac{e^{x_i}}{\sum_{j} e^{x_j}}
$$

大型語言模型 (LLM) 產生的每個 token 機率都經過此函數。使用 $e^{x_i}$ 而非其他正值函數 (如 $2^{x_i}$) 的原因在於：以 $e$ 為底時，$\ln$ 與 $\exp$ 互為反函數，搭配交叉熵損失函數後，梯度化簡為 $\hat{y}_i - y_i$ (預測值減目標值)，沒有額外的 $\ln a$ 係數。

### Log-Sum-Exp 技巧：數值穩定性的守護者
在實作 Softmax 或計算機率分佈的總和時，常遇到一組數值 $\{x_1, \dots, x_n\}$，我們需要計算：
$$
\text{LSE}(x) = \ln \left( \sum_{i=1}^n e^{x_i} \right)
$$
若 $x_i$ 的值很大 (例如 $x_i = 1000$)，直接計算 $e^{1000}$ 會導致浮點數溢位 (overflow)。相反地，若 $x_i$ 很小且為負，則會發生下溢 (underflow)。

為解決此問題，實務中會使用 **Log-Sum-Exp trick**：
$$
\text{LSE}(x) = a + \ln \left( \sum_{i=1}^n e^{x_i - a} \right), \quad a = \max(x_1, \dots, x_n)
$$
藉由減去最大值 $a$，指數項內的最大值恆為 $e^0 = 1$，保證計算過程中不會發生溢位。這個小技巧是所有現代機器學習框架 (如 PyTorch, TensorFlow) 在處理交叉熵與 Softmax 時的標準做法，再次體現 $e$ 在工程實作中對數值邊界的敏感性。

### 交叉熵損失函數
[交叉熵](https://en.wikipedia.org/wiki/Cross-entropy) (cross-entropy) 是訓練分類器與語言模型最常用的損失函數，定義為：
$$
H(p, q) = -\sum_{i} p_i \ln q_i
$$

其中 $p$ 為目標分佈，$q$ 為模型預測的分佈 (通常由 softmax 產生)。此處的 $\ln$ 即以 $e$ 為底的自然對數。選用 $\ln$ 使得交叉熵與 KL 散度、資訊熵之間保持一致的單位，且對 softmax 輸出求導後，梯度形式簡潔，適合大規模梯度下降最佳化 (關於對數為何是唯一合理選擇、以及 $e$ 為底的結構性優勢，詳見後文「資訊理論中的 $\ln$ 與自然單位」一節)。

從 1683 年 Bernoulli 對複利問題的探索，到當代每個神經網路的每次梯度更新，$e$ 貫穿始終。Bernoulli 當年思考的是利率，無從預見自己正在為梯度下降奠基。純粹數學的危險之處在於：它看似無用，直到忽然無處不在。

## 深受 $e$ 啟發的平衡三進位系統
現代電腦普遍採用二進位系統，但這並非唯一的選擇。事實上，三進位計算機曾在歷史上留下珍貴足跡，並於近年重新受到關注。

1950 年代末至 1960 年代，莫斯科國立大學的研究人員設計出早期的三進位電子計算機「[Сетунь](https://en.wikipedia.org/wiki/Setun)」(英語發音: Setun) 與其後繼型號「Сетунь 70」。「Сетунь」這個名字來自莫斯科附近的一條河。「Сетунь」電子計算機的設計由數學家 Sergei Sobolev (С·Л·Соболев) 等人於 1956 年左右發起，目的是提供高價效比的計算資源給大專院校與科研機構。研究團隊在資源有限的情況下，設計出基於鐵氧體磁芯 (ferrite cores) 與半導體二極體 (diodes) 的平衡三進位計算機，其邏輯電路利用正電壓 (代表 `+1`) 、零電壓 (代表 `0`) 與負電壓 (代表 `-1`) 表示 3 個狀態，這種設計一度被認為具有潛在的速度、功耗和可靠性優勢。

這台採用時序邏輯的機器配備快速乘法器，使用小型鐵氧體環記憶體 (ferrite ring memory) 作為暫存器/累加器，搭配磁鼓 (magnetic drum) 作為主要記憶體，提供 24 道指令，其中包含預留指令，原型機於 1958 年完成，經過測試後於 1960 年代初期小規模生產。儘管「Сетунь」據稱表現穩定、應用廣泛，甚至獲得國外關注，但由於未能契合當時蘇聯主流的二進位計算機發展規劃，該專案最終未能大規模發展。總共生產約 50 台，分佈於蘇聯各地，用於工程計算、工業控制與教學等。

1970 年，原團隊的部分成員推出改進的「Сетунь 70」機種，並建立更系統化的三進位運算結構，包括：
- 定義 "Tryte" (三進位的位元組)，由 6 個 trit (三進位元) 構成，約等於 $\log_2(3^6) \approx 9.5$ 個二進位位元 (bit) 的資訊量 (此為 Сетунь 70 的設計；原始 Сетунь 採用 18-trit 字組)
- 支援三進位指令集與可變運算長度 (例如 1 至 6 tryte)

然而，由於缺乏持續的支持，「Сетунь 70」也僅生產極少量，三進位計算機的發展就此沉寂。「[Сетунь](https://en.wikipedia.org/wiki/Setun)」最關鍵的特點正是採用平衡三進位的對稱編碼 $({ -1, 0, +1 })$，能自然且統一地表示正負數，無需額外的符號位元或二補數轉換。

電腦科學家 Donald E. Knuth 在《[The Art of Computer Programming](https://en.wikipedia.org/wiki/The_Art_of_Computer_Programming)》第 2 卷說：
> "Perhaps the prettiest number system of all is the balanced ternary notation"
> (也許所有數值系統中最佳雅者，就是平衡三進位表示法)

此處 "ternary" 指「三元的」或「以三為基底的」(base-3)。標準三進位使用 $({0, 1, 2})$，而 [balanced ternary](https://en.wikipedia.org/wiki/Balanced_ternary) (平衡三進位) 則使用對稱的 $\{-1, 0, +1\}$ 作為基本單元 (trit)，常簡記為 $({ T, 0, 1 })$ 或 $({ -, 0, + })$。
> the ternary values as being "balanced" around the mid-point of 0. The same rules apply to ternary as to any other numeral system: The right-most symbol, R, has its own value and each successive symbol has its value multiplied by the base, B, raised to the power of its distance, D from R.
> (三進位數值圍繞中點 $0$ 平衡分佈。其規則與其他數值系統相同：最右邊的符號 R 代表其自身的值，往左每個符號的值需乘以基數 B 的 D 次方，D 為該符號與 R 的距離)

### 為何選擇三進位？
一項用來度量進位系統效率的指標是「基數經濟性」([radix economy](https://en.wikipedia.org/wiki/Optimal_radix_choice))，它試圖量化表示給定範圍數值所需的「成本」，通常定義為基數 $b$ 與所需位數 $k$ 的乘積: $E(b, N) = b \times k$，此處 $k$ 表示 $N$ 個不同數值所需的最少位數，即 $k = \lceil \log_b N \rceil$。該成本函數 $E(b, N)$ 試圖在「單一位數的複雜度」 (基數 $b$ 越大，單一位數能表示的狀態越多，但也可能越複雜) 與「表示的總長度」 (位數 $k$ 越少越好) 之間取得平衡。

倘若我們暫時忽略位數 $k$ 必須是整數的限制，使用 $k \approx \log_b N$，並利用對數換底公式 $\log_b N = \frac{\ln N}{\ln b}$，於是表示成本的近似理論值可寫為：
$$
C(b, N) \approx b \cdot \log_b N = b \cdot \frac{\ln N}{\ln b} = \ln N \cdot \frac{b}{\ln b}
$$

對於固定的數值範圍 $N$ (此時 $\ln N$ 為常數) ，要最小化成本 $C(b, N)$，就等同於最小化函數 $f(b) = \frac{b}{\ln b}$。

我們可藉由微分找出 $f(b)$ 的最小值：
$$
f'(b) = \frac{d}{db} \left( \frac{b}{\ln b} \right) = \frac{\ln b - 1}{(\ln b)^2}
$$

令導數 $f'(b) = 0$，得到 $\ln b - 1 = 0$，解出 $b = e \approx 2.718...$。由於 $f(b)$ 在定義域 $b > 1$ 上，於 $(1, e)$ 區間遞減、於 $(e, \infty)$ 區間遞增 (可由 $f'(b)$ 的分子 $\ln b - 1$ 在 $b < e$ 時為負、$b > e$ 時為正確認)，故 $b = e$ 為全域最小值。換言之，從理論上看，以 $e$ 為基底效率最高。

下圖基於 $E(b, N) = b \lceil \log_b N \rceil$ 計算 (另一種成本定義，考慮位數必須是整數)，比較不同基底表示特定範圍 $N$ 的成本：

該圖顯示，雖然 $e$ 是理論最佳值，但對於實際的整數基底，$3$ 比 $2$ 或 $4$ 更接近 $e$ 所代表的最佳效率點。接著比較整數基底的理論效率，也就是基底 $b_1$ 和 $b_2$ 對應的成本 (忽略常數 $\ln{N}$) ：
$$
\frac{C(b_1)}{C(b_2)} \approx \frac{b_1 / \ln(b_1)}{b_2 / \ln(b_2)} = \frac{b_1 \ln(b_2)}{b_2 \ln(b_1)}
$$
計算 $f(b) = b / \ln(b)$ 的值：
* $f(2) = 2 / \ln(2) \approx 2.885$
* $f(3) = 3 / \ln(3) \approx 2.731$

由於 $f(3) < f(2)$，這表明基底 $3$ 在理論上比基底 $2$ 更經濟。因此，三進位是整數基底中最具儲存效率的選擇，這正是三進位系統 (特別是平衡三進位) 引人入勝的理論基礎。

### 平衡三進位的運算特性
平衡三進位使用 $({-1, 0, +1})$ 作為 trit (三進位元) 的值，其數值表示圍繞 0 對稱：
$$
n = \sum_i a_i \cdot 3^i, \quad a_i \in \{-1, 0, +1\}
$$

平衡三進位採用與二進位相似的位置表示法概念，但其關鍵優勢在於負數表示極為簡便：只需將正數表示式中所有非零的 trit 反轉 (亦即 `+1` 變 `-1`，`-1` 變 +1) 即可得到其負數。
* 二進位 (使用二補數系統的 4 位元)
$(6)_{10} = (0110)_2$
$(-6)_{10} = (\text{反轉 } 0110 \rightarrow 1001; \text{ 加 } 1 \rightarrow 1010)_2$
* 平衡三進位 (簡稱: `bal3`)
$(6)_{10} = (1T0)_{bal3} = 1 \cdot 3^2 + (-1) \cdot 3^1 + 0 \cdot 3^0 = 9 - 3 + 0 = 6$
$$
\begin{flalign*}
(-6)_{10} &= (\text{反轉 } 1T0 \rightarrow T10)_{bal3} && \\
&= (-1) \cdot 3^2 + 1 \cdot 3^1 + 0 \cdot 3^0 && \\
&= -9 + 3 + 0 && \\
&= -6 &&
\end{flalign*}
$$

由上可知，在平衡三進位中獲取一個數的負數，只需簡單地將所有非零 trit 進行符號反轉，無需像二進位那樣進行位元反轉再加 $1$ (二補數系統) 的操作，運算本質上更為快速簡單。考慮以下平衡三進位範例 (使用 $114 = 81+27+9-3$：
$$
\begin{flalign*}
(114)_{10} &= (111T0)_{bal3} \\
&= 1 \cdot 3^4 + 1 \cdot 3^3 + 1 \cdot 3^2 + (-1) \cdot 3^1 + 0 \cdot 3^0 \\
&= 81 + 27 + 9 - 3 + 0 = 114 &&
\end{flalign*}
$$

其負數 $-114$ 的表示法為：
$$
\begin{flalign*}
(-114)_{10} &= (\text{反轉 } 111T0 \rightarrow TTT10)_{bal3} \\
&= (-1) \cdot 3^4 + (-1) \cdot 3^3 + (-1) \cdot 3^2 + 1 \cdot 3^1 + 0 \cdot 3^0 \\
&= -81 - 27 - 9 + 3 + 0 = -114 &&
\end{flalign*}
$$

這種對稱性代表正負數使用統一的規則表示，無需像在二進位系統中區分 `signed` 和 `unsigned` 類型。

平衡三進位不僅能一致地表示整數，也可用於表示分數 (近似浮點數) 。以下是十進位 `0.2` 的一個平衡三進位近似表示：
$$
\begin{flalign*}
(0.1TT1)_{bal3} &= 1 \times 3^{-1} - 1 \times 3^{-2} - 1 \times 3^{-3} + 1 \times 3^{-4} \\
&= \frac{1}{3} - \frac{1}{9} - \frac{1}{27} + \frac{1}{81} = \frac{(27 - 9 - 3 + 1)}{81} \\
&= \frac{16}{81} \approx 0.1975 &&
\end{flalign*}
$$

如同十進位的 $\frac{1}{3}$ 之於二進位，$0.2$ 在三進位中是無限循環小數 $(0.\overline{1TT1})_{bal3}$。因此，任何有限長度的平衡三進位表示都只是近似值。

如何近似表示十進位 `0.8` 呢？利用 `0.8 = 1 - 0.2` 及平衡三進位的簡易負數轉換特性：
$$
\begin{flalign*}
1 - (0.1TT1...)_{bal3} &\approx (1.\overline{T11T}...)_{bal3} \\
(1.T11T)_{bal3} &= 1 - 1 \times 3^{-1} + 1 \times 3^{-2} + 1 \times 3^{-3} - 1 \times 3^{-4} \\
&= 1 - \frac{1}{3} + \frac{1}{9} + \frac{1}{27} - \frac{1}{81} \\
&= \frac{(81 - 27 + 9 + 3 - 1)}{81} = \frac{65}{81} \approx 0.8025 &&
\end{flalign*}
$$

上述近似計算，顯示平衡三進位中小數運算的對稱性。

接著，我們評估平衡三進位在計算負數時的效率。比較二進位 (採用二補數系統)與平衡三進位。以十進位數值 `114` 為例，求對應的負值：
* 二進位 (8 位元二補數)：給定 $(114)_{10} = (01110010)_2$
- 反轉所有位元: `10001101`
- 加 1: `10001110` (結果)
- 計算量: 8 次位元反轉 + 1 次 (可能帶進位的) 遞增操作
* 平衡三進位 (使用 5 trit 足以表示 $\pm 114$)：給定$(114)_{10} = (111T0)_{bal3}$
- 反轉所有非零 trit: $(TTT10)_{bal3}$ (結果)
- 計算量: 5 次 trit 反轉

若要量化這份差距，主要差異在於二進位的二補數所需的 `+1` 遞增操作。假設單個位元或 trit 的反轉時間相當，則平衡三進位的優勢在於省去加法步驟。n 位元二進位數加 1 時，發生進位傳播的平均位數期望值 $E(n)$ 為：
$$
E(n) = \sum_{k=0}^{n-1} \frac{k}{2^{k+1}} + \frac{n}{2^n} = 1 - \left(\frac{1}{2}\right)^n
$$

當 $n$ 增大時，期望值趨近於 $1$。這表示二進位取負運算平均會多出接近一次完整的加法器傳播延遲。因此，在轉換效率上，平衡三進位執行負數操作本質上更快。這個看似微小的延遲差異，在深度學習、大型語言模型訓練與推理等等需要巨量運算的情境中，可望累積顯著的效能提升。

### 現代應用：針對大型語言模型的壓縮、加速及節能
平衡三進位本身並非由歐拉數推導而來，但前述 radix economy 分析所揭示的 $e$ 最佳基底，為三進位在資訊編碼上的效率提供理論依據。近年 Microsoft 研究院發表的 [BitNet b1.58](https://arxiv.org/abs/2402.17764)，將大型語言模型 (LLM) 的權重限制為 $\{-1, 0, +1\}$ 三種狀態。名稱中的「1.58-bit」來自 $\log_2(3) \approx 1.58$，即三態所需的最小資訊量；但這是資訊論上的理論下限，實際部署時每個權重的儲存開銷仍取決於位元打包、索引結構與推論 kernel 的設計。

此方法的關鍵特質不在於「把數字硬塞進最少位元」，而在於運算邏輯的簡化。LLM 中佔主導地位的矩陣乘法 ($Y = WX$)，當權重 $W$ 僅含 $\{-1, 0, +1\}$ 時，浮點乘法被整數加減取代，$0$ 權重則直接跳過，從而降低計算延遲與記憶體頻寬壓力。相較於純粹的 1-bit 量化 ($\{-1, +1\}$)，保留 $0$ 能維持模型內在的稀疏性，減少資訊損失。根據該論文 (2024 年 2 月) 的實驗，在相同模型規模與訓練資料條件下，BitNet b1.58 可在延遲、記憶體用量、吞吐率與能耗上取得明顯改善。

這類運算 (僅需加/減/忽略) 也開啟專用硬體加速器 (ASIC 或 FPGA) 的設計空間，與前蘇聯「Сетунь」電子計算機的三進位硬體探索形成跨時代的呼應。平衡三進位兼具理論上接近 $e$ 的最佳整數基底效率，與實務上簡化運算、降低能耗的特性，使其在大型模型的壓縮、加速及節能上展現潛力。

延伸閱讀：
- [The Balanced Ternary Machines of Soviet Russia](https://dev.to/buntine/the-balanced-ternary-machines-of-soviet-russia)
- [Ternary Neural Networks for Resource-Efficient AI Applications](https://arxiv.org/abs/1609.00222)
- [Binarized Neural Networks: Training Deep Neural Networks with Weights and Activations Constrained to +1 or -1](https://arxiv.org/abs/1602.02830)
- [Bitnet.cpp: Efficient Edge Inference for Ternary LLMs](https://arxiv.org/abs/2502.11880)
- [Run DeepSeek R1 Dynamic 1.58-bit](https://unsloth.ai/blog/deepseekr1-dynamic)

## 歐拉數相關數學定理和方法
與歐拉數相關的數學定理和公式極為豐富，遍佈於數學與科學的多個領域，以下列舉其中幾個。

:::success
「永無止境地循環下去的數字，和讓人難以捉摸的虛數畫出簡潔的軌跡，在某一點落地。雖然沒有圓的出現, 但來自宇宙的 $\pi$ 飄然地來到 $e$ 的身旁，和害羞的 $i$ 握著手。他們的身體緊緊地靠在一起，屏住呼吸, 但有人加了 $1$ 以後，世界就毫無預警地發生了巨大的變化。一切都歸於 $0$。歐拉公式就像是暗夜中閃現的一道流星；也像是刻在漆黑的洞窟裡的一行詩句。」, 小川洋子《博士熱愛的算式》
:::

### [歐拉公式](https://en.wikipedia.org/wiki/Euler%27s_formula)與歐拉恆等式
$$
e^{i\pi} + 1 = 0
$$

上式稱為[歐拉恆等式](https://en.wikipedia.org/wiki/Euler%27s_identity) (Euler's identity)，將 5 個基本常數：$e$, $\pi$, $i$ (虛數單位), $1$ (乘法單位元)，和 $0$ (加法單位元)，以最簡潔的形式聯繫在一起，是歐拉公式 $e^{ix} = \cos x + i \sin x$ 在 $x = \pi$ 的特例。這不僅展現數學的對稱性與結構美，更是物理學、電路理論、波動方程、訊號處理等領域的根基。

為說明此公式的由來，先考慮一個粒子在實數線上運動，其位置為 $f(t)$，若假設粒子的速度等於當下位置，即：
$$
\frac{d}{dt} f(t) = f(t), \quad f(0) = 1
$$

這是一階線性常微分方程，其唯一解為 $f(t) = e^t$

此解代表一種「自我增長」的行為：當粒子的位置越大，其變化速率也越快，呈現典型的指數成長。

接著將此觀念推廣至複平面，設 $f(t) \in \mathbb{C}$ 表示粒子的位置，我們希望建構最簡單的圓周運動模型。假設：
- 粒子恆以單位速率運動
- 速度方向與位置向量常垂直 (即始終朝切線方向)
- 初始位置為 $f(0) = 1$

在複數表示中，向量旋轉 $90^\circ$ 相當於乘以虛數單位 $i$:

因此可建構以下運動方程：
$$
\frac{d}{dt} f(t) = i f(t), \quad f(0) = 1
$$

透過泰勒展開可驗證 $e^{it} = \cos t + i \sin t$。另一條推導路徑則從 [De Moivre 定理](https://en.wikipedia.org/wiki/De_Moivre%27s_formula) 與指數函數的極限定義出發。

複數乘法的幾何意義是「角度相加、模長相乘」。對單位圓上的複數 $\cos\theta + i\sin\theta$ (模長為 $1$)，反覆自乘 $n$ 次即得：
$$
(\cos\theta + i\sin\theta)^n = \cos n\theta + i\sin n\theta
$$

此為 De Moivre 定理。1714 年，Roger Cotes 已給出可解讀為 $\log(\cos x + i\sin x) = ix$ 的關係，但此式涉及複對數的多值性，需謹慎理解。

De Moivre 定理可自然地推出歐拉公式。令 $\theta = n\varphi$，則：
$$
\left(\cos\frac{\theta}{n} + i\sin\frac{\theta}{n}\right)^n = \cos\theta + i\sin\theta
$$

當 $n \to \infty$ 時，$\frac{\theta}{n} \to 0$。由 $\cos u = 1 + O(u^2)$ 與 $\sin u = u + O(u^3)$，底數滿足：
$$
\cos\frac{\theta}{n} + i\sin\frac{\theta}{n} = 1 + \frac{i\theta}{n} + O(1/n^2)
$$

因此其與 $1 + \frac{i\theta}{n}$ 的比值為 $1 + O(1/n^2)$，取 $n$ 次方後趨近 $1$，故：
$$
\lim_{n\to\infty} \left(1 + \frac{i\theta}{n}\right)^n = \cos\theta + i\sin\theta
$$

而前文已建立指數函數的極限定義 $\exp(z) = \lim_{n\to\infty}(1+z/n)^n$。以 $z = i\theta$ 代入，即得：
$$
e^{i\theta} = \cos\theta + i\sin\theta
$$

這正是歐拉公式的本體。此推導強調 De Moivre 定理所反映的複數乘法幾何意義，並以指數函數的極限定義銜接到 $e^{i\theta}$。

若令 $\theta = \pi$，由 $\cos\pi = -1$、$\sin\pi = 0$，則：
$$
e^{i\pi} = -1 \quad \Rightarrow \quad e^{i\pi} + 1 = 0
$$

這不僅是個形式優雅的恆等式，也將指數函數、三角函數與複數自然地統整為一體。

De Moivre 定理還揭示 $n$ 次單位根的結構。方程 $z^n = 1$ 的 $n$ 個複數解為：
$$
z_k = \cos\frac{2\pi k}{n} + i\sin\frac{2\pi k}{n} = e^{i \cdot 2\pi k/n}, \quad k = 0, 1, \ldots, n-1
$$

這 $n$ 個根均勻分佈在單位圓上，構成正 $n$ 邊形的頂點。單位根的概念貫穿代數、數論與訊號處理，例如離散傅立葉轉換 (DFT) 的核心運算 $W_N^k = e^{-i2\pi k/N}$ 正是 $N$ 次單位根。

在物理學中，$e^{it}$ 對應單位圓上的等速圓周運動，其在 $x$ 軸與 $y$ 軸上的投影分別為 $\cos t$ 與 $\sin t$。這類簡諧振動模型構成週期現象的基礎，如聲波、電磁波、光波與量子波動等皆可視為其延伸。

根據傅立葉分析，任何「行為良好」的週期函數 $f(t)$ (週期為 $T$)都可展開為一組離散頻率的複指數函數 $e^{i n \omega_0 t}$ (其中 $\omega_0 = 2\pi/T$ 是基本角頻率，而 $n$ 為整數) 的線性組合，稱為傅立葉級數 (Fourier series)：
$$
f(t) = \sum_{n=-\infty}^{\infty} c_n e^{i n \omega_0 t}
$$
其中係數 $c_n = \frac{1}{T} \int_{T} f(t) e^{-i n \omega_0 t} dt$ 代表函數在頻率 $n\omega_0$ 上的複數幅度。

即便是非週期訊號 $f(t)$ (或視為週期無限大的函數)，亦可藉由傅立葉轉換 (Fourier transform)，表示為連續頻率 $\omega$ 的複指數函數 $e^{i\omega t}$ 的積分形式 (一種連續的疊加)：
$$
f(t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} F(\omega) e^{i\omega t} \, d\omega \quad \text{(傅立葉逆轉換)}
$$
其中 $F(\omega)$ 稱為 $f(t)$ 的傅立葉頻譜，由以下傅立葉轉換計算得出：
$$
F(\omega) = \int_{-\infty}^{\infty} f(t) e^{-i\omega t} \, dt \quad \text{(傅立葉轉換)}
$$

在科學與工程中，面對複雜問題時，「分解」是種關鍵策略。傅立葉分析正是基於此思想，將複雜函數或訊號分解到一組更簡單、具有普適性的「基底函數」上。為了有效達成這種分解，我們借鑒向量空間的概念，特別是「正交性」和「內積」。


在三維歐幾里得空間中，一組標準正交基底為 $\mathbf{e}_1 = (1, 0, 0)$, $\mathbf{e}_2 = (0, 1, 0)$, $\mathbf{e}_3 = (0, 0, 1)$，它們滿足以下特性：
- 正交性 (orthogonality): 不同基底向量之間的內積 (點積)為零。例如：
$$
\langle \mathbf{e}_1, \mathbf{e}_2 \rangle = (1)(0) + (0)(1) + (0)(0) = 0
$$
- 正規性 (normality): 每個基底向量與自身的內積為 1 (即向量長度為 1)。例如：
$$
\langle \mathbf{e}_1, \mathbf{e}_1 \rangle = (1)(1) + (0)(0) + (0)(0) = 1
$$

這裡，二個向量 $\mathbf{u} = (u_1, u_2, u_3)$ 和 $\mathbf{v} = (v_1, v_2, v_3)$ 的內積定義為：
$$
\langle \mathbf{u}, \mathbf{v} \rangle = u_1 v_1 + u_2 v_2 + u_3 v_3
$$

任何一個三維向量 $\mathbf{v} = (v_1, v_2, v_3)$ 都可唯一地表示為這組正交基底的線性組合：
$$
\mathbf{v} = v_1 \mathbf{e}_1 + v_2 \mathbf{e}_2 + v_3 \mathbf{e}_3
$$
其中的係數 $v_k$ (即向量在 $\mathbf{e}_k$ 方向上的分量)可藉由內積計算得到：
$$
v_k = \langle \mathbf{v}, \mathbf{e}_k \rangle
$$
這是因為 $\langle \mathbf{v}, \mathbf{e}_k \rangle = \langle \sum_j v_j \mathbf{e}_j, \mathbf{e}_k \rangle = \sum_j v_j \langle \mathbf{e}_j, \mathbf{e}_k \rangle = v_k \langle \mathbf{e}_k, \mathbf{e}_k \rangle = v_k$ (利用基底的正交性 $\langle \mathbf{e}_j, \mathbf{e}_k \rangle = 0$ if $j \neq k$, 和正規性 $\langle \mathbf{e}_k, \mathbf{e}_k \rangle = 1$)。
因此，內積提供將向量「投影」到基底方向上，以獲得其分量 (係數)的方法。

傅立葉分析將這個概念從有限維的向量空間推廣到無限維的函數空間。我們可以將 (某些類型的)函數視為函數空間中的「向量」。為了定義函數空間中的「正交性」，我們需要推廣內積的概念。對於定義在 $(-\infty, \infty)$ 上的複數值函數 $f(t)$ 和 $g(t)$ (假設它們屬於適當的函數空間，如平方可積函數 $L^2(\mathbb{R})$)，它們的內積通常定義為：
$$
\langle f, g \rangle = \int_{-\infty}^{\infty} f(t) \overline{g(t)} \, dt
$$
其中 $\overline{g(t)}$ 表示 $g(t)$ 的複數共軛。

現在，考慮傅立葉分析中扮演基底角色的複指數函數。令 $f(t) = e^{i \omega_1 t}$ 和 $g(t) = e^{i \omega_2 t}$，計算它們的內積：
$$
\langle e^{i \omega_1 t}, e^{i \omega_2 t} \rangle = \int_{-\infty}^{\infty} e^{i \omega_1 t} \overline{e^{i \omega_2 t}} \, dt = \int_{-\infty}^{\infty} e^{i \omega_1 t} e^{-i \omega_2 t} \, dt = \int_{-\infty}^{\infty} e^{i (\omega_1 - \omega_2) t} \, dt
$$
這個積分在嚴格數學意義上不收斂，但在廣義函數 (或分佈)的框架下，其結果與狄拉克 δ 函數 (Dirac delta function) 相關：
$$
\int_{-\infty}^{\infty} e^{i (\omega_1 - \omega_2) t} \, dt = 2\pi \delta(\omega_1 - \omega_2)
$$
狄拉克 δ 函數的性質是：當 $\omega_1 \neq \omega_2$ 時 $\delta(\omega_1 - \omega_2) = 0$，而在 $\omega_1 = \omega_2$ 時為無窮大 (積分為 1)。這表明，不同頻率的複指數函數 $e^{i\omega_1 t}$ 與 $e^{i\omega_2 t}$ 在內積意義下是「正交」的 (它們的內積在 $\omega_1 \neq \omega_2$ 時為 0)。因此，$\{e^{i\omega t} \mid \omega \in \mathbb{R}\}$ 構成一組連續的「正交基底函數」。

藉由這組正交基底，我們就可以將任意 (行為良好的)函數 $f(t)$ 分解 (投影)到這些基底上。其在特定頻率 $\omega$ 上的「分量」或「投影係數」，正是由傅立葉轉換 $F(\omega)$ 給出：
$$
F(\omega) = \int_{-\infty}^{\infty} f(t) e^{-i\omega t} \, dt = \langle f(t), e^{i\omega t} \rangle
$$
$F(\omega)$ 是個複數，其大小 $|F(\omega)|$ 代表頻率 $\omega$ 在 $f(t)$ 中的強度 (幅度)，其相位 $\arg(F(\omega))$ 代表該頻率成分的相對相位。

而傅立葉逆轉換：
$$
f(t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} F(\omega) e^{i\omega t} \, d\omega
$$
則表示原始函數 $f(t)$ 可由其所有頻率分量 $e^{i\omega t}$ 按照頻譜 $F(\omega)$ 給出的權重 (複數幅度) 連續疊加 (積分) 而還原回來。

傅立葉分析的本質，就是將一個在時域 (time domain) 中描述的函數 $f(t)$，透過投影到正交的複指數基底 $e^{i\omega t}$ 上，轉換為在頻域 (frequency domain) 中描述的頻譜 $F(\omega)$。這種轉換讓我們能清晰地看到訊號由哪些頻率成分組成，及各個成分的強度與相位，這對於訊號處理、系統分析、物理建模等領域至關重要。這正是 $e$ (透過 $e^{i\omega t}$) 所構築的連續變化在頻率分析上的體現。

延伸閱讀: [圖解傅立葉分析](https://hackmd.io/@sysprog/fourier-transform)

### [質數定理](https://en.wikipedia.org/wiki/Prime_number_theorem)
> 「我認為質數的魅力在於無法說明它出現的秩序。每個質數都沒有因數，似乎隨意地嵌在數列之間。雖然數字越大，想要找出質數也就越困難，卻仍然無法依據任何已知規則來準確預測它們的出現。這種捉摸不定、難以駕馭的特性，對追求完美的博士而言，彷彿就是一種無法抗拒的誘惑。」

這段話出自小說《博士熱愛的算式》，凸顯質數既孤高又神秘，看似雜亂卻富含某種隱約的規律。人類對質數的探索已有二千多年，古希臘時期，歐幾里得以反證法證明質數的無窮多性，這是最早的論述，而既然知曉質數無窮存在，便產生新問題：前 $x$ 個自然數中究竟包含多少質數、它們又如何分佈？

十八、十九世紀，法國數學家 Adrien-Marie Legendre 和德國數學家高斯 (Gauss) 先後提出猜想：前 $x$ 個自然數裡的質數個數 $\pi(x)$ 大致與 $\frac{x}{\ln(x)}$ 或 $\frac{x}{\ln(x) - 1.08366}$ 相當，後者是 Legendre 公開的版本。高斯則提出對數積分
$$
\mathrm{Li}(x) = \int_{2}^{x} \frac{dt}{\ln t},
$$
用以更精細描述質數的實際增長。這些近似式不完美，但都強調質數分佈深受對數 $\ln(x)$ 影響。

歐拉曾為「質數數量」的研究奠下重要基礎。他考慮下列乘積：
$$
\prod_{p} \frac{p}{p - 1},
$$

其中 $p$ 逐一取自所有質數。他用不夠嚴謹但極具巧思的推導，將此「質數乘積」與「調和級數」 $\sum_{n=1}^\infty \tfrac{1}{n}$ 連結，顯示二者同樣發散，也再次印證質數不可能僅有限多個。進一步推廣後，歐拉定義
$$
\prod_{p} \frac{1}{1 - p^{-z}} = \sum_{n=1}^{\infty} \frac{1}{n^z},
$$
左邊針對全體質數的無窮乘積，右邊則針對全體自然數的無窮級數。當指數 $z$ 擴張至複數時，右邊級數成為黎曼 $\zeta$ 函數，左邊則保存質數的隱性編碼。這種「以質數乘積對應整個自然數系統」的理念在歐拉時代已現端倪，後經黎曼發展為複分析的有力工具，改變研究質數分佈的方式。

1859 年，黎曼在成功把 $\zeta$ 函數解析延拓至整個複平面 (除 $s = 1$ 處的簡單極點外)，將質數分佈與 $\zeta$ 函數零點緊密相連，首次以複分析方法探討實函數 $\pi(x)$。1896 年，法國數學家 Hadamard 與比利時數學家 de la Vallée Poussin 證明
$$
\lim_{x \to \infty}
\frac{\pi(x)}{x/\ln x} = 1,
$$
宣示質數定理成立。Chebyshev 曾顯示若上述比值的極限存在，必為 1，並由此推得[Bertrand-Chebyshev 定理](https://en.wikipedia.org/wiki/Bertrand%27s_postulate) (Bertrand Chebyshev theorem)：對任何自然數 $n$，區間 $[n,\,2n]$ 內必有質數。質數定理雖揭示整體分佈規律，許多微觀細節仍不乏謎團，例如，高斯的 $\mathrm{Li}(x)$ 在廣泛區域能很好貼近 $\pi(x)$，但二者的差 $\mathrm{Li}(x) - \pi(x)$ 會於極大範圍內多次翻轉正負，首次翻轉位置大到超越現今電腦極限，牽涉[斯奎斯數](https://en.wikipedia.org/wiki/Skewes%27s_number) (Skewes' number) 的研究領域。

經歷長期的純數論研究累積，質數的特性於 1970 年代末期發展成為現代公開金鑰密碼學 (public-key cryptography) 的關鍵。公元 1977 年，Ron Rivest、Adi Shamir 與 Leonard Adleman 在其發表的論文〈[A Method for Obtaining Digital Signatures and Public-Key Cryptosystems](https://dl.acm.org/doi/10.1145/359340.359342)〉中，提出以其姓氏首字母命名的 RSA 公鑰加密演算法。該演算法建基於選取二個極大的秘密質數 $p$ 與 $q$，計算其公開的乘積模數 $N=pq$。接著，選取一個與歐拉函數值 $\phi(N)=(p-1)(q-1)$ 互質的整數 $e$ 作為公共指數。數對 $(N, e)$ 構成公鑰並對外發布，而作為基礎的質因數 $p$ 與 $q$ 則必須保密 (用於產生私鑰)。RSA 演算法的安全性依賴於大整數[質因數分解問題](https://en.wikipedia.org/wiki/Integer_factorization) (integer factorization problem) 的計算複雜度，即從已知的公開模數 $N$ 反向分解出原始質因數 $p$ 與 $q$，在現有計算能力下被認為是不可行的 (computationally infeasible)。質數分解的困難性因此構成 RSA 演算法安全性的理論基石，並促使以質數理論為基礎的公鑰密碼學，在金融交易和雲端運算等多個關鍵領域得到廣泛應用。

質數分佈看似混沌，但質數定理揭示其宏觀規律：質數的平均密度可由以
$e$ 為底的自然對數 $1/\ln(x)$ 近似描述，將連續變化的標度 $e$ 引入離散的數論領域。此定理不僅成為連接離散數學與連續分析的橋樑，更將深刻的數論洞見轉化為資訊安全的關鍵方法，充分體現純粹理論與工程實踐的交織。

延伸閱讀：[Frans Oort 談質數](https://www.math.sinica.edu.tw/media/pdf/d421/42104.pdf)

### [常態分佈](https://en.wikipedia.org/wiki/Normal_distribution)
$$
f(x) = \frac{1}{\sqrt{2\pi} \cdot \sigma} \cdot e^{-\frac{(x - \mu)^2}{2\sigma^2}}
$$

其中 $\mu$ 為均值、$\sigma$ 為標準差、$\sigma^2$ 為變異數。此函數描述隨機變數在平均值附近出現的機率，並且依距離遠近呈對稱分佈。根據[中央極限定理](https://en.wikipedia.org/wiki/Central_limit_theorem)，只要有大量彼此獨立且同分佈的隨機變數，其總和在適當正規化後將趨近常態分佈。因此，常態分佈不僅是統計學的基礎模型，也廣泛應用於物理、生物、經濟、訊號處理等科學領域。

其中的指數項 $e^{-\frac{(x - \mu)^2}{2\sigma^2}}$ 起著關鍵作用。它不僅賦予整體分佈平滑的變化與對稱結構，也反映出「距離越遠，機率越小」的特性。

在物理學中，常態分佈的出現更具有深層意涵。例如，熱擴散方程的基本解即為高斯型 (即常態型) 分佈，這描述粒子在隨機熱運動下擴散的機率密度。該現象的隨機性來自於微觀層級上眾多獨立作用的累積，而其宏觀行為則收斂為常態分佈，這與中央極限定理完全吻合。事實上，在許多與隨機擾動相關的微分方程中，指數函數 $e^{-x^2}$ 經常自然浮現，作為穩定解或平衡態的形式。

接著思考常態分佈中為何會出現 $e$ 與 $x^2$，我們可從基本假設出發。設想一組由 $n$ 個獨立的隨機事件組成的系統，每個事件只有發生 (記作 $1$) 與不發生 (記作 $0$) 兩種結果。當 $n$ 很大時，整體行為會趨近於連續隨機變數，其總和的分佈逐漸接近常態分佈。此過程中，自然底數 $e$ 之所以出現，是因為其代表連續極限下的指數變化，與離散事件總數呈 $2^n$ 增長相呼應。

至於 $x^2$ 的來源，則可由空間對稱性與獨立性兩個條件導出。假設隨機變數在空間中呈旋轉對稱 (機率密度僅依據距離 $r$，與角度 $\theta$ 無關)，同時假設各軸方向上的誤差彼此獨立，那麼根據 Herschel 和 Maxwell 的分析，唯一能同時滿足這二個條件的聯合密度函數形式為：
$$
p(x, y) \propto e^{-c (x^2 + y^2)}
$$

其中 $c$ 為正常數。由於 $x^2 + y^2 = r^2$，密度僅取決於距原點的距離 $r$，體現旋轉對稱性。

這表示在二維或三維空間中，具有旋轉對稱性與軸向獨立性的隨機分佈，其機率密度會隨距離平方衰減，並自然形成以 $x^2$ 為變數的指數函數形式。

常態分佈中 $e^{-x^2/2}$ 的結構並非偶然，而是以下幾個自然假設的直接推論：
- 大量獨立隨機變數總和的分佈，根據中心極限定理，收斂為常態分佈
- 空間對稱性與軸向獨立性導致密度必然依據 $x^2$ (或更廣義的 $r^2$) 變化
- 指數衰減是維持總機率為 1 的最自然方式
- 底數 $e$ 為連續極限中描述增減變化最自然的常數

因此，$e^{-x^2/2}$ 不僅體現數學形式的優雅，也蘊含物理與統計行為背後的對稱和隨機特性，是自然與數學深刻結合的典範。

### [泊松分佈](https://en.wikipedia.org/wiki/Poisson_distribution)
描述在固定時間或空間範圍內，隨機事件發生的次數。其機率質量函數 (probability mass function) 為：
$$
P(k) = \frac{\lambda^k e^{-\lambda}}{k!}
$$

這裡的 $e^{-\lambda}$ (即 $P(0)$) 代表隨機事件「發生零次」的機率，類似於自然界中常見的衰減現象，例如前述的放射性衰變。除了提供描述隨時間 (或空間)隨機變化的數學形式外，$e^{-\lambda}$ 也確保機率總和為 1，是整體分佈「標準化」的關鍵。

泊松分佈可看作是[二項分佈](https://en.wikipedia.org/wiki/Binomial_distribution)在次數 $n \to \infty$、成功率 $p \to 0$ 且期望值 $\lambda = np$ 固定的極限形式。該極限推導自然導出 $e^{-\lambda}$ 的出現。

為什麼現實中常出現泊松分佈？回顧其定義：若某事件在單位時間內平均發生 $\lambda$ 次，且每次發生獨立、隨機、不可同時，於是在時間間隔內發生 $x$ 次的機率為：
$$
P(N = x) = \frac{\lambda^x e^{-\lambda}}{x!}
$$

這個形式與 $e^x$ 的泰勒展開驚人地相似：
$$
e^x = \sum_{n=0}^\infty \frac{x^n}{n!}
$$

二者的結合意義不只在於形式上的巧合。考慮將時間細分為 $n$ 段，每段事件發生的機率為 $p = \lambda/n$，則根據二項分佈：
$$
P(N = x) = C_n^x \left( \frac{\lambda}{n} \right)^x \left(1 - \frac{\lambda}{n} \right)^{n-x}
$$

當 $n \to \infty$ 時，上式逼近泊松分佈，也自然導出 $e^{-\lambda}$。這不僅是計算上的結果，更可視為平均每單位時間發生次數為 $\lambda$ 的一種極限模型。

進一步理解 $e^x$ 的意義，有助於揭示泊松分佈與 $e$ 的深層連結。若把時間看作連續且隨機事件可能不斷發生的場域，則 $e^x$ 表示某平均值為 $x$ 的複合增長。從極限定義出發：
$$
e^x = \lim_{n \to \infty} \left(1 + \frac{x}{n} \right)^n
$$

這可理解為：每一小段時間內有極小的發生機率，整體隨時間累積出平均次數為 $x$ 的事件，形成 $e^x$ 的指數形式。當所有乘項都是 $1$，表示從未發生；乘項中有個 $x/n$，則表示發生一次，依此類推。這與泊松分佈描述在眾多微小獨立機會下，最終觀測到特定事件次數的概念相符。

指紋識別晶片藉由數千個微型電容感測器生成指紋的拓撲圖，並針對每個畫素點進行比對，裡頭也能見到泊松分佈。

指紋可看作[手指的紋路](https://en.wikipedia.org/wiki/Fingerprint_scanner)，而常見的感測方式主要可分為光學式感應器與電容式感應器二種，前者透過光源與感光元件組合，儲存手指輪廓所反射的光線變化，後者則是藉由半導體晶片上的微型電容陣列，感測皮膚表面不同部位與感測板之間的電容差異。以電容式感應最為常見，iPhone 5s 所採用的正是這種技術。

電容式指紋辨識模組的基本原理是，對每個畫素點的電容感應顆粒預先充電到某一準位，當手指接觸感測器時，系統會偵測各畫素點的放電電流。由於指紋的脊 (凸起) 與峪 (凹下) 會產生不同的電容值，放電速率亦隨之改變，藉此可區分出凹凸紋理的位置，從而重建出完整的指紋圖像。這些圖像資料會儲存於快閃記憶體中作為後續比對使用，辨識流程則包含：從指紋圖中提取特徵點、建立特徵向量，並將其轉換為數位訊號進行比對。

電容式設計薄而輕巧，適合用於行動裝置上，其成本較高，且感應元件易受到汗水或油脂影響，為了解決這個問題，感應器上會額外加裝保護層，例如 iPhone 5s 採用藍寶石玻璃覆蓋其指紋感測器以增加耐久性。這些畫素點可視為獨立的伯努利試驗 (Bernoulli trial)，每個畫素點只有「是否誤配」的兩種結果。在試驗次數眾多、單次誤配機率極低的情況下，根據機率論的極限定理，這樣的重複伯努利試驗所導致的總誤配次數可近似為服從泊松分佈。

設 $\lambda$ 為雜訊造成的平均誤配次數，$X$ 為實際誤配的次數，則 $X$ 滿足：
$$
\mathbb{P}(X = k) = \frac{e^{-\lambda} \lambda^k}{k!}
$$

亦即，系統可藉此計算某次比對中出現恰好 $k$ 個錯誤對應的機率。當 $\lambda$ (平均噪點數) 超過某個閾值時，系統會判定這樣的比對極可能由雜訊造成，並觸發二次驗證機制。這也解釋為何在手指潮濕時，解鎖常常會失敗。

舉例來說，若系統估計正常狀況下雜訊造成的誤配次數為 $\lambda = 1.5$，並將比對閾值設為 $k = 5$，則系統會計算：
$$
\mathbb{P}(X \geq 5) = 1 - \sum_{k=0}^4 \frac{e^{-\lambda} \lambda^k}{k!}
$$

這樣的尾部機率可作為安全門檻的參考，決定是否強制啟用第二道驗證。$e$ 在這些生物識別技術中不只是數學符號，讓我們能以簡潔的數學形式，處理現實中充滿變數與雜訊的感測資料。

以下用金融案例，解說[泊松過程](https://en.wikipedia.org/wiki/Poisson_point_process)。

在日常生活中，我們經常會遇到排隊等候的情境，例如在銀行辦事時，客戶隨機到達櫃檯。那麼，若在某段時間內統計有多少人加入排隊？這類問題的數學建模可從最簡化的假設出發。

設 $N(t)$ 為在時間區間 $[0, t]$ 內到達的顧客數，且假設初始時無人排隊，即 $N(0) = 0$。我們對 $N(t)$ 所代表的隨機過程提出以下假設：
1. 任意短時間段中新增人數的機率只依賴時間長度，與目前時刻無關 (平穩性)
2. 不重疊的時間段內新增人數彼此獨立 (增量獨立性)
3. 時間段 $\Delta t$ 足夠小時，新增 1 人的機率近似為 $\lambda \Delta t$，其中 $\lambda > 0$ 為常數
4. 在 $\Delta t \to 0$ 的極限下，新增二人或以上的機率為高階無窮小，記為 $o(\Delta t)$

上述條件構成泊松過程的基礎假設。我們可數學化地寫為：
$$
\mathbb{P}\{N(t + \Delta t) - N(t) = 1\} = \lambda \Delta t + o(\Delta t),
$$

$$
\mathbb{P}\{N(t + \Delta t) - N(t) \geq 2\} = o(\Delta t).
$$

根據此，我們得到：
$$
\mathbb{P}\{N(t) = 0\} = 1 - \lambda t + o(t).
$$

設 $P(t) = \mathbb{P}\{N(t) = 0\}$，則：
$$
\begin{aligned}
P(t + \Delta t) &= \mathbb{P}\{N(t + \Delta t) = 0\} \\
&= \mathbb{P}\{N(t) = 0\} \cdot \mathbb{P}\{N(t + \Delta t) - N(t) = 0\} \\
&= P(t) (1 - \lambda \Delta t + o(\Delta t))
\end{aligned}
$$

兩邊減去 $P(t)$ 再同除以 $\Delta t$ 並令 $\Delta t \to 0$，得微分方程：
$$
P'(t) = -\lambda P(t)
$$

解得：
$$
P(t) = e^{-\lambda t}
$$

也就是 $\mathbb{P}\{N(t) = 0\} = e^{-\lambda t}$。$e$ 的出現揭示：即使泊松過程是個離散計數過程，其背後的變化機制仍可由連續變化的指數函數所支配。

若允許 $\lambda$ 隨時間變化，則稱為非齊次泊松過程；若 $\lambda(t)$ 本身也是隨機過程，則形成 [Cox 過程](https://en.wikipedia.org/wiki/Cox_process) (又稱 doubly stochastic Poisson process)。

泊松過程的應用之一是建構隨機跳變的模型。不同於布朗運動 (Brownian motion，連續變化)，泊松過程適用於處理突發事件，例如資產價格突升、信用違約事件等。

#### 資產價格跳變模型
在金融建模中，實際資產價格常會表現出「非連續、突變」的特徵，違反 [Black-Scholes 模型](https://en.wikipedia.org/wiki/Black%E2%80%93Scholes_model)中連續路徑的假設。此時，可將價格過程設為：
$$
S(t) = S(0) e^{\mu t + \sigma W(t)} \cdot \prod_{i=1}^{N(t)} J_i,
$$

其中 $N(t)$ 為泊松過程，$J_i$ 為跳變幅度的隨機變數。下圖中黑底圓圈 `1` 展現資產收益率不連續的變化:
> 1 minute USDJPY April 14, 2015。取自: Donnelly, B. (2019). The art of currency trading: a professional's guide to the foreign exchange market. John Wiley & Sons.

這種模型被稱為跳躍擴散模型 (jump-diffusion model)，由 Merton 最早提出。這裡的連續成分來自指數函數 $e^{\mu t}$，而跳躍次數則由泊松過程驅動，反映 $e$ 如何貫穿於連續與離散變化的混合結構中。

#### 信用風險建模與 hazard rate

在結構化信用風險模型中，違約事件常用泊松過程建模。若以 $\lambda(t)$ 表示[危險率](https://en.wikipedia.org/wiki/Failure_rate) (hazard rate)，即瞬時違約強度，則某公司在 $[0, t]$ 內發生違約的機率為：
$$
\mathbb{P}\{\text{default before } t\} = 1 - e^{-\int_0^t \lambda(s) ds}.
$$

這裡的 $e^{-\int \lambda}$ 結構亦再次彰顯 $e$ 在描述「累積風險」這類連續變化機制上的自然角色。此模型由 Jarrow-Turnbull (1995) 建立，是現今信用風險評價的其中一種主流架構。

#### 信用違約交換 (CDS) 與市場觀察

[信用違約交換](https://en.wikipedia.org/wiki/Credit_default_swap) (CDS, credit default swap) 報價可反推市場對違約率的預期。在最簡模型中：
$$
S = \lambda (1 - R)
$$

其中 $S$ 是 CDS 溢價，$R$ 為回收率。這可視為將市場報價轉化為違約機率估計的橋樑。

以下圖表顯示 2010 年歐洲主權債務危機時，各國 CDS 溢價快速上升，市場即透過泊松跳變預期模型來估算違約風險：
> [Sovereign Credit Default Swaps](https://en.wikipedia.org/wiki/Credit_default_swap) (2010 年歐債危機)

泊松過程因其簡潔的機率結構與高度可解性，成為排隊理論、風險管理與金融數學中不可或缺的工具，而與歐拉數 $e$ 的深刻連結，屢屢展現本文強調的觀念：離散事件的累積機率，往往可由連續變化的函數來精確描述。

### [秘書問題](https://en.wikipedia.org/wiki/Secretary_problem)與最佳停止理論
除了上述分佈，歐拉數 $e$ 在機率論的「最佳停止問題」(optimal stopping problem) 中也扮演決定性角色，最著名的案例便是「秘書問題」。

假設你要面試 $n$ 位應徵者。你的目標是從中選出最好的一位。規則如下：
1. 應徵者逐一前來面試，你必須在面試後立即決定是否錄取或拒絕
2. 錄取後面試終止；拒絕後不能反悔錄用已面試過的人
3. 你只能比較目前面試者與先前面試者的優劣，無法預知未來面試者的程度

最佳策略是：先面試前 $k$ 個人且一律不錄取，從第 $k+1$ 個人開始，只要遇到比前 $k$ 個人都優秀的應徵者，就立即錄取。$k$ 該如何選取？當 $n$ 很大時，最佳的 $k$ 值趨近於 $n/e$ (約佔總人數的 $37\%$)。使用此策略選到最強應徵者的機率，也同樣趨近於 $1/e$。

> 取自 [Wikipedia](https://en.wikipedia.org/wiki/Secretary_problem)

設成功錄取最佳者的事件依最佳者實際出現的位置 $i$ 拆解 ($i > k$，否則必被略過)。將事件分為位置事件與其下的條件事件：
- 最佳者恰好排在第 $i$ 個面試位置，機率為 $1/n$
- 在最佳者位於第 $i$ 位的條件下，從第 $k+1$ 個到第 $i-1$ 個皆未被錄取。這等價於前 $i-1$ 個人裡的最佳者落在前 $k$ 個觀察位置之中，否則策略會在到達 $i$ 之前就先選了某個次佳者。此條件機率為 $k/(i-1)$

對所有可能的 $i$ 求和，得到成功機率：
$$
P(k) = \sum_{i=k+1}^{n} \frac{1}{n} \cdot \frac{k}{i-1} = \frac{k}{n} \sum_{i=k+1}^{n} \frac{1}{i-1}
$$

右側為調和級數的部分和 (partial sum)。當 $n$ 足夠大時，以積分近似離散求和：
$$
P(k) \approx \frac{k}{n} \int_{k+1}^{n} \frac{1}{i-1} \, di = \frac{k}{n} \bigl[\ln(i-1)\bigr]_{k+1}^{n} = \frac{k}{n} \ln \frac{n-1}{k}
$$

由於 $n$ 很大時 $\ln \frac{n-1}{k} \approx \ln \frac{n}{k} = -\ln \frac{k}{n}$，令 $x = k/n$，問題化為求下列函數在 $(0, 1)$ 上的最大值：
$$
f(x) = -x \ln x
$$

對 $x$ 求導：
$$
f'(x) = -\ln x - 1
$$

令 $f'(x) = 0$ 解得 $\ln x = -1$，即 $x = 1/e$。又 $f''(x) = -1/x < 0$ 在 $(0, 1)$ 上恆成立，故此臨界點為唯一最大值。代回原式：
$$
f(1/e) = -\frac{1}{e} \ln \frac{1}{e} = \frac{1}{e} \approx 0.368
$$

於是 $k \approx n/e$ 時成功機率達到最大值，且該最大值本身也約為 $37\%$。換言之，最佳停止位置與最佳成功機率皆收斂於 $1/e$，是 $e$ 在最佳停止理論中的雙重出場。

此結論的直觀理解如下：若拒絕太少人 ($k$ 太小)，你可能還沒遇到真正優秀的人就先錄取平庸者；若拒絕太多人 ($k$ 太大)，優秀的人可能已經在前面被你拒絕，導致最後不得不錄取最後一人。$1/e$ 恰好是在「探索」(exploration) 與「利用」(exploitation) 之間取得平衡的數學分水嶺。這種策略廣泛應用於相親、買房尋價、乃至於電腦科學中的任務排程與資源分配。

### 隨機演算法與資料結構中的 $e$
$e$ 不僅出現在連續數學中，也頻繁浮現於離散的演算法分析。以下列舉幾個電腦科學中的經典情境。

#### 雜湊表的碰撞機率
將 $n$ 個鍵值隨機映射至 $n$ 個槽位的雜湊表中，某個特定槽位恰好為空的機率為 $(1 - 1/n)^n$，隨著 $n$ 增大，此值趨近 $1/e \approx 0.368$。因此，即使雜湊函數完全均勻，約有 $n/e$ 個槽位預期為空，$n(1 - 1/e)$ 個槽位被佔用，且平均每個被佔用的槽位含 $n / (n(1-1/e)) = 1/(1 - 1/e) \approx 1.58$ 個鍵值。此結果揭示理想雜湊的碰撞密度下界。

#### Bloom filter 的最佳雜湊函數數量
[Bloom filter](https://en.wikipedia.org/wiki/Bloom_filter) 使用 $k$ 個獨立雜湊函數與 $m$ 位元的位元陣列來判斷元素是否可能存在於集合中。插入 $n$ 個元素後，某個特定位元仍為 $0$ 的機率為：
$$
\left(1 - \frac{1}{m}\right)^{kn} \approx e^{-kn/m}
$$

誤判率 (false positive rate) 為所有 $k$ 個雜湊位元皆為 $1$ 的機率：
$$
\varepsilon \approx \left(1 - e^{-kn/m}\right)^k
$$

對 $k$ 求導並令其為零，可得最小化誤判率的最佳雜湊函數數量為：
$$
k^* = \frac{m}{n} \ln 2
$$

此時誤判率為 $\varepsilon^* = (1/2)^k = (1/2)^{(m/n)\ln 2}$。推導過程中，$\ln 2 = \log_e 2$ 的出現源自對 $e^{-kn/m}$ 求極值，而非任何人為的底數選擇。

#### 贈券蒐集者問題
[贈券蒐集者問題](https://en.wikipedia.org/wiki/Coupon_collector%27s_problem) (coupon collector's problem) 問：若有 $n$ 種不同的贈券，每次隨機獲得一種，要蒐集齊全所有種類，期望需要多少次？答案為：
$$
E[T] = n H_n = n \sum_{k=1}^{n} \frac{1}{k} \approx n \ln n + \gamma n
$$

其中 $H_n$ 為第 $n$ 個調和數，$\gamma \approx 0.5772$ 為 Euler-Mascheroni 常數 (注意：這是文章開頭區分的「歐拉常數」，而非歐拉數 $e$)。$\ln n$ 的出現源自調和級數 $H_n \approx \ln n$ 的漸近行為，而此漸近關係的根基正是 $\int_1^n \frac{1}{x} dx = \ln n$。此問題在分散式系統的負載平衡、隨機測試的覆蓋率分析中有直接應用。

#### Stirling 近似公式

當我們需要估算極大階乘 $n!$ 的增長時，[Stirling 公式](https://en.wikipedia.org/wiki/Stirling%27s_approximation)揭示 $e$ 在離散數學規律中的關鍵地位：
$$
n! \approx \sqrt{2\pi n} \left( \frac{n}{e} \right)^n
$$

或者取對數形式：
$$
\ln(n!) \approx n \ln n - n + \frac{1}{2} \ln(2\pi n)
$$

此公式在工程實務中極其重要。在本文「以 binary splitting 計算 $e$」的 C 實作中，我們正是利用此公式來精確預測：若要取得小數點後 $d$ 位精度的歐拉數，無窮級數 $\sum 1/k!$ 究竟需要展開到第幾項。具體而言，程式碼中的 `log_factorial` 函式便是 Stirling 公式的直接體現，它讓演算法得以在運算前就完成精確的資源分配，避免冗餘計算。

### 統計中的對數變換
前文已從多個角度闡述 $e$ 與 $\ln$ 在數學、物理和金融中的關鍵地位。然而在統計學與資料分析的實務中，「對變數取對數」幾乎是最常見的資料變換手法之一。這並非偶然，而是自然對數的數學性質與真實資料的生成機制深度契合的結果。

#### 穩定變異數
許多真實世界的量測資料 (如電力產量、股價、GDP) 具有一個共同特徵：隨著數值本身增大，其波動幅度也隨之增大。以美國每月電力生產資料為例，早期年份的產量低，波動幅度小；近年產量高，波動幅度明顯增大。這種「變異數隨平均值增大」的現象，在統計學中稱為[異質變異數](https://en.wikipedia.org/wiki/Heteroscedasticity) (heteroscedasticity)，違反許多統計模型 (如線性迴歸、時間序列分析) 要求誤差項具有同質變異數 (homoscedasticity) 的基本假設。

對數變換之所以能穩定變異數，可從以下推導理解。設時間序列 $Z_t$ 的標準差與其均值成正比：
$$
\sqrt{Var(Z_t)} = \mu_t \times \sigma, \quad \text{其中 } E(Z_t) = \mu_t
$$

由定義可推得：$Z_t = \mu_t \bigl(1 + \frac{Z_t - \mu_t}{\mu_t}\bigr)$。利用 $\ln(1+x) \approx x$ (當 $x$ 足夠小) 的近似性質：
$$
\ln(Z_t) \approx \ln(\mu_t) + \frac{Z_t - \mu_t}{\mu_t}
$$

於是 $E(\ln(Z_t)) \approx \ln(\mu_t)$，而 $Var(\ln(Z_t)) \approx \sigma^2$，為常數。原本隨均值膨脹的變異數，經對數變換後變為穩定常數，使資料適用於假設同質變異數的分析方法。

更一般地，當變異數 $\sigma^2$ 與均值 $E(y)$ 之間滿足不同的比例關係時，可選用不同的[變異數穩定化變換](https://en.wikipedia.org/wiki/Variance-stabilizing_transformation) (variance-stabilizing transformation)：

| $\sigma^2$ 與 $E(y)$ 的關係 | 適用變換 |
|---|---|
| $\sigma^2 \propto \text{常數}$ | $y' = y$ (無需變換) |
| $\sigma^2 \propto E(y)$ | $y' = \sqrt{y}$ (平方根，適用於 Poisson 計數資料) |
| $\sigma^2 \propto E(y)(1 - E(y))$ | $y' = \arcsin(\sqrt{y})$ (反正弦) |
| $\sigma^2 \propto (E(y))^2$ | $y' = \ln y$ (對數) |
| $\sigma^2 \propto (E(y))^3$ | $y' = y^{-1/2}$ (倒平方根) |
| $\sigma^2 \propto (E(y))^4$ | $y' = y^{-1}$ (倒數) |

對數變換對應的正是 $\sigma^2 \propto (E(y))^2$ 的情況，也就是變異數與均值平方成正比的資料。

#### 乘法機制與對數常態分佈
為何眾多自然與經濟資料呈右偏分佈，且取對數後近似常態？根本原因在於資料的生成機制是乘法性的，而非加法性的。

[中央極限定理](https://en.wikipedia.org/wiki/Central_limit_theorem)告訴我們：大量獨立隨機變數的加法和趨近常態分佈。然而，收入、企業營收、資產規模、生物族群、城市人口等變數的成長方式，本質上是乘法累積的，每期的變化率 (而非絕對變化量) 才是獨立隨機的。設 $X_T$ 為經 $T$ 期累積後的值：
$$
X_T = X_0 \cdot r_1 \cdot r_2 \cdots r_T
$$

其中 $r_t$ 為各期的獨立隨機成長因子。取對數：
$$
\ln X_T = \ln X_0 + \ln r_1 + \ln r_2 + \cdots + \ln r_T
$$

右側為獨立隨機變數的加法和，由中央極限定理可知 $\ln X_T$ 趨近常態分佈。因此 $X_T$ 本身服從[對數常態分佈](https://en.wikipedia.org/wiki/Log-normal_distribution) (log-normal distribution)，其機率密度函數為：
$$
f(x) = \frac{1}{x \sigma \sqrt{2\pi}} \exp\!\left(-\frac{(\ln x - \mu)^2}{2\sigma^2}\right), \quad x > 0
$$

對數常態分佈具有右偏、下界為零、長尾等特徵，與許多實證資料的形態高度吻合。取對數後，偏態消除，資料近似常態，使得依賴常態假設的統計推論 (如 $t$ 檢定、ANOVA、線性迴歸) 得以適用。

#### 增長率的近似與經濟學詮釋
自然對數的另一個實用性質是：對數的差分近似百分比變化。利用 $\ln(1+x) \approx x$ (當 $|x|$ 較小)：
$$
\ln(Z_t) - \ln(Z_{t-1}) = \ln\!\left(\frac{Z_t}{Z_{t-1}}\right) = \ln(1 + X_t) \approx X_t
$$

其中 $X_t = (Z_t - Z_{t-1})/Z_{t-1}$ 為增長率。這表示取對數後的差分直接近似於原始變數的百分比變動，使研究增長率分佈是否穩定的問題，轉化為研究對數差分是否服從常態分佈的問題。

在計量經濟學中，此性質賦予迴歸係數直觀的「彈性」(elasticity) 詮釋。考慮對數-對數模型：
$$
\ln Y = \beta_0 + \beta_1 \ln X + u
$$

係數 $\beta_1$ 的意義為：$X$ 每變動 $1\%$，$Y$ 平均變動 $\beta_1\%$。這正是經濟學中彈性的定義。類似地，在對數-水平模型 $\ln Y = \beta_0 + \beta_1 X + u$ 中，$\beta_1$ 表示 $X$ 每增加一個單位，$Y$ 近似變動 $100 \times \beta_1\%$。

人口成長提供一個具體的例子。若不考慮資源限制，生物族群的增長率 $r$ 為常數，人口 $P$ 滿足 $\frac{dP}{dt} = rP$，解為 $P(t) = P_0 e^{rt}$。取對數得 $\ln P(t) = \ln P_0 + rt$。若要比較兩個國家的增長率差異，可引入虛擬變數 $D$ (國家 A 為 $0$，國家 B 為 $1$) 並建構迴歸模型 $\ln P = \alpha + (r + \beta D) t$，其中 $\beta$ 表示國家 B 相對於國家 A 的增長率之差。若 $\beta > 0$，代表國家 B 的人口增長速率高於國家 A。在原始尺度上，兩國的人口比值 $P_B/P_A = e^{\beta t}$ 隨時間指數擴大，但在對數尺度上僅為斜率差異 $\beta t$，這正是對數變換將乘法結構線性化的體現。

#### Box-Cox 變換的特例
對數變換可視為 [Box-Cox 變換](https://en.wikipedia.org/wiki/Power_transform#Box%E2%80%93Cox_transformation)族的特例。Box-Cox 變換定義為：
$$
y_i^{(\lambda)} =
\begin{cases}
\dfrac{y_i^\lambda - 1}{\lambda}, & \lambda \neq 0 \\[6pt]
\ln y_i, & \lambda = 0
\end{cases}
$$

當 $\lambda = 0$ 時即為對數變換。$\lambda$ 由資料驅動估計：右偏 (長尾在右側) 資料通常對應較小的 $\lambda$ (含 $\lambda < 0$)，左偏資料通常對應較大的 $\lambda$。對數函數的凹性 (concavity) 使其能壓縮大值、展開小值的間距，從而緩解右偏。$|\lambda|$ 越大，對原始數值的調整越顯著。

然而，對離散計數資料 (如 Poisson 分佈的計數值) 直接取對數需格外謹慎。離散資料可能取 $0$，而 $\ln(0)$ 無定義，此時常用 $\ln(x+k)$ 的形式處理，但 $k$ 的選擇會引入額外偏誤。關於此議題，可參閱 O'Hara 與 Kotze 的研究 "Do not log-transform count data" (Methods in Ecology and Evolution, 2010)。

#### 最大概似估計與對數
在機率與統計推論中，[最大概似估計](https://en.wikipedia.org/wiki/Maximum_likelihood_estimation) (maximum likelihood estimation, MLE) 是最基本的參數估計方法。設觀測資料 $x_1, x_2, \ldots, x_n$ 來自參數為 $\theta$ 的機率模型，概似函數為：
$$
L(\theta) = \prod_{i=1}^n f(x_i; \theta)
$$

直接對連乘積求極值在計算上極為困難。取對數後，乘法變為加法：
$$
\ln L(\theta) = \sum_{i=1}^n \ln f(x_i; \theta)
$$

由於 $\ln$ 為嚴格遞增函數，最大化 $\ln L$ 與最大化 $L$ 等價。對 $\ln L$ 求導並令其為零：
$$
\frac{\partial \ln L}{\partial \theta_k} = 0, \quad k = 1, 2, \ldots, m
$$

即可求解參數估計值。這裡 $\ln$ 將乘法結構轉化為加法結構的性質，與前文 Napier 發明對數時的動機一脈相承，從 17 世紀簡化天文計算的工具，演變為現代統計推論的基礎運算。

#### 資訊理論中的 $\ln$ 與自然單位
Shannon 在公元 1948 年定義[資訊熵](https://en.wikipedia.org/wiki/Entropy_(information_theory))：
$$
H = -\sum_i p_i \log p_i
$$

對數的底數決定資訊量的單位：以 $\log_2$ 計量時，單位為 bit (binary digit)；以 $\ln$ (即 $\log_e$) 計量時，單位為 [nat](https://en.wikipedia.org/wiki/Nat_(unit)) (natural unit of information)；以 $\log_{10}$ 計量時，單位為 [ban](https://en.wikipedia.org/wiki/Ban_(unit)) (又稱 hartley)。各單位僅差一個常數因子，例如 $H_{\text{nat}} = H_{\text{bit}} \cdot \ln 2$。

此處 $\ln$ 的選擇不僅是為了計算便利，更與前文「基數經濟性」中 $e$ 為最佳基底的特性相契合：在自然單位 (nat) 下，資訊量的度量最能反映系統內在的自由度，且在處理連續機率分佈的變換時，雅可比行列式 (Jacobian) 的對數項直接以 $\ln$ 呈現，使理論框架保持一致。

為什麼非用對數不可？[Khinchin](https://en.wikipedia.org/wiki/Aleksandr_Khinchin) 證明，若要求不確定性度量 $H$ 滿足四條公理：
* 對機率 $p_i$ 連續
* 在均勻分佈時取最大值
* 增添機率為零的事件不改變 $H$ (可擴充性，即 $H(p_1, \ldots, p_n, 0) = H(p_1, \ldots, p_n)$)
* 對獨立系統具可加性 ($H(AB) = H(A) + H(B)$) — 則 $H$ 必為 $-C \sum p_i \log p_i$ 的形式，其中正常數 $C$ 僅決定單位

換言之，對數函數是唯一滿足上述公理的選擇，不存在以其他函數 (如冪函數或多項式) 替代的可能。此外，Shannon 的[信源編碼定理](https://en.wikipedia.org/wiki/Shannon%27s_source_coding_theorem)賦予 $H$ 一個具體的工程意義：對於任何無失真壓縮方案，平均每個符號至少需要 $H$ 個位元 (以 $\log_2$ 為底)。此結論的推導依賴 [Kraft-McMillan 不等式](https://en.wikipedia.org/wiki/Kraft%27s_inequality) $\sum_x 2^{-l_x} \le 1$ (其中 $l_x$ 為符號 $x$ 的碼字長度)：在此約束下，以 Lagrange 乘數法最小化平均碼長 $L = \sum_x p_x l_x$，可得實數最佳碼長 $l_x^* = -\log_2 p_x$，代回即得 $L^* = H$。實際的前綴碼要求碼長為整數，因此實務上的最佳平均碼長滿足 $H \le L < H + 1$；僅當所有 $p_x$ 恰為 2 的負冪時，才能以整數碼長精確達到 $L = H$。

上述公理化結論可用更基本的方式直接推導。設 $f(x)$ 表示發生機率為 $x$ 的事件所攜帶的資訊量 ($x \in (0, 1]$)，則兩個獨立事件同時發生的機率為 $x_1 x_2$ (機率相乘)，而其總資訊量等於個別資訊量之和 (資訊可加性)：
$$
f(x_1 x_2) = f(x_1) + f(x_2), \quad x_1, x_2 \in (0, 1]
$$

這是 Cauchy 函數方程在乘法形式下的變體。令 $x_1 = x_2 = 1$ 可得 $f(1) = 0$，即必然事件不攜帶資訊。為使後續微分論證成立，需先將 $f$ 的定義域從 $(0, 1]$ 擴充至全體正實數：對任意 $x > 1$，因 $1/x \in (0, 1)$，由 $f(x \cdot \frac{1}{x}) = f(x) + f(\frac{1}{x}) = f(1) = 0$ 可得 $f(x) = -f(1/x)$，故函數方程 $f(xy) = f(x) + f(y)$ 在 $(0, \infty)$ 上同樣成立。接著，對 $f(t)$ 求導：
$$
f'(t) = \lim_{dt \to 0} \frac{f(t + dt) - f(t)}{dt}
$$

利用 $t + dt = t \cdot \bigl(1 + \frac{dt}{t}\bigr)$ 及函數方程：
$$
f(t + dt) = f\!\left(t \cdot \frac{t + dt}{t}\right) = f(t) + f\!\left(1 + \frac{dt}{t}\right)
$$

代入可得：
$$
f'(t) = \lim_{dt \to 0} \frac{f\!\left(1 + \frac{dt}{t}\right)}{dt} = \frac{1}{t} \lim_{u \to 0} \frac{f(1 + u) - f(1)}{u} = \frac{f'(1)}{t}
$$

對兩端從 $1$ 至 $x$ 積分：
$$
f(x) = f'(1) \int_1^x \frac{1}{t}\,dt = f'(1) \ln x
$$

由於 $x \in (0, 1]$ 時 $\ln x \leq 0$，而資訊量應為非負，故 $f'(1) \leq 0$。令 $\gamma = -f'(1) > 0$，得 $f(x) = -\gamma \ln x$。透過換底公式 $\ln x = \log_a x \cdot \ln a$，可將結果改寫為 $f(x) = -\log_a x$ (其中 $a = e^{1/\gamma}$)。底數的選擇僅決定資訊量的單位 (bit、nat 或 ban)，而對數形式本身則是由獨立事件的機率乘法與資訊量加法這兩條前提所唯一決定。

上述推導也賦予對數一個直觀的操作意義：$\log_a x$ 可理解為「需要問多少個是非問題才能確定結果」。以猜 $1$ 至 $32$ 之間的整數為例，每次提問可排除一半的候選 (二分搜尋)，因此所需問題數為 $\log_2 32 = 5$。一般而言，若某事件的發生機率為 $p$，則確認該事件所需的最少是非問題數量級為 $-\log_2 p$，此即 Shannon 資訊量的主體概念。

在機器學習的理論推導中，偏好以 $e$ 為底有其結構性的原因。考慮 [KL 散度](https://en.wikipedia.org/wiki/Kullback%E2%80%93Leibler_divergence) (Kullback-Leibler divergence)：
$$
D_{\text{KL}}(p \| q) = \sum_i p_i \ln \frac{p_i}{q_i}
$$

若對數採用 $\ln$，則 KL 散度滿足以下性質：
* 非負性：$D_{\text{KL}}(p \| q) \geq 0$，等號成立若且唯若 $p = q$，此證明直接來自 $\ln x \leq x - 1$ (即 $e^x \geq 1 + x$ 的推論)
* 與 Fisher 資訊矩陣的關聯：在參數空間中，$D_{\text{KL}}$ 的二階展開為 $\frac{1}{2} \delta\theta^T I(\theta) \delta\theta$，其中 $I(\theta)$ 為 Fisher 資訊矩陣，此展開在 $\ln$ 下形式最簡
* 與最大概似估計的等價：最小化交叉熵 $H(p, q)$ 等價於最小化 $D_{\text{KL}}(p \| q)$ (因為 $H(p)$ 為常數)，而 MLE 恰好是此最佳化的解

以 $e$ 為底的選擇也使得資訊熵與統計力學中的 Boltzmann 熵 $S = k_B \ln W$ 形式一致，且在連續分佈的微分熵 $h(X) = -\int f(x) \ln f(x) \, dx$ 定義中，$\ln$ 與機率密度函數的積分自然相容，不會引入額外的換底常數。

對數在統計中的角色，歸結為三個數學性質的交匯：
* 將乘法轉換為加法 ($\ln(ab) = \ln a + \ln b$)，使乘法性的生成機制可套用加法性的定理 (如中央極限定理)
* 導數為 $1/x$，使得對數差分自然對應百分比變化，賦予迴歸係數直觀的彈性詮釋
* 作為 $e^x$ 的反函數，與指數族分佈 (常態、Poisson、指數等) 的充分統計量緊密相連，使最大概似估計的計算大幅簡化

這三者並非獨立的巧合，而是 $e$ 作為描述連續 變化之基石的不同面向。從 Napier 以對數簡化乘法運算，到現代統計學家以對數變換穩定變異數、線性化模型、簡化推論，$\ln$ 的角色始終如一：將複雜的乘法世界，映射至人類更擅長處理的加法世界。

## 橫跨連續現實與離散計算的橋樑
我們身處的世界，本質是連續的。時間猶如河流般流淌不息，汽車的速度、放射性元素的衰變速率、電容器兩端的電壓等現象，都在時時刻刻發生著平滑而連續的變化。在數學上，我們習慣使用連續變數的函數及微積分來描述這一切，其中，微分運算子 $D$ (即 $\frac{d}{dt}$ 或 $\frac{d}{dx}$) 用以描述連續的世界。

然而，數位電腦的世界運作方式截然不同。它們的內部時鐘以離散的節拍跳動，透過類比數位轉換器 (ADC) 在特定的時間點進行「取樣」(sampling)，並藉由數位類比轉換器 (DAC) 在特定時刻施加控制。電腦僅能在離散的時刻 $t_n$ (其間隔通常固定，記為 $h$ 或 $\Delta t$) 進行觀測與互動，對於二次取樣之間所發生的連續變化全然不知。

工程人員面臨的關鍵挑戰，便是如何運用這些離散的觀測結果，來理解、預測，甚至控制本質上連續運行的實體系統。Jack Crenshaw 博士在其於 [embedded.com](https://www.embedded.com/) 發表的經典系列文章中，提出極具洞察力的數學工具，巧妙地將其比作「羅塞達石碑」(Rosetta Stone)：
$$
z = e^{hD} \tag{8}
$$

此處 $D$ 代表連續領域中的微分運算子 ($\frac{d}{dt}$)，$z$ 代表離散時間域中的前移運算子 (shift operator，即 $z y_n = y_{n+1}$)，$h$ 是固定的取樣時間間隔。嚴格而言，此等式為形式冪級數關係，成立的前提是被操作的函數在取樣點鄰域內解析 (即泰勒級數收斂至函數本身)；在實務應用中，當 $h$ 足夠小且函數行為良好時，截斷展開即可得到精確度可控的數值近似。

「羅塞達石碑」源於其在破解古埃及象形文字上的關鍵作用。石碑上刻有同一段內容的三種不同文字 (古埃及象形文、埃及草書、古希臘文) ，使得學者能夠相互對照，最終解開失傳語言之謎，開啟現代埃及學的大門。Crenshaw 博士借此比喻，強調此方程如同石碑一般，扮演翻譯者的角色：它將描述連續系統動態的語言 (涉及微分運算子 $D$ 的微分方程) 與描述離散系統行為的語言 (涉及移位運算子 $z$ 的差分方程或數位濾波器) 聯繫起來，提供一座橫跨二個看似不同世界的數學橋樑，讓我們得以運用數位電腦的離散計算能力，來分析、模擬和駕馭真實世界中的連續動態系統。

「羅塞達石碑」的根基是泰勒級數，它根據函數 $f(x)$ 在某點 $x$ 的值及其所有導數，預測 $f(x + h)$：
$$
f(x + h) = f(x) + h \left[ \frac{df}{dx} \right]_x + \frac{h^2}{2!} \left[ \frac{d^2 f}{dx^2} \right]_x + \frac{h^3}{3!} \left[ \frac{d^3 f}{dx^3} \right]_x + \cdots \tag{9}
$$

其中 $[\cdot]_x$ 表示先取導數，再在 $x$ 求值。當 $h$ 很小時，高階項因階乘迅速衰減，級數變得實用。對於多項式，級數甚至精確終止。

$e$ 源於微分方程 $\frac{d}{dx} f(x) = f(x)$ (初值 $f(0) = 1$) ，其解為：
$$
e^x = 1 + x + \frac{x^2}{2!} + \frac{x^3}{3!} + \cdots $$

將 $x = hD$ 代入，泰勒級數可寫為：
$$
f(x + h) = e^{hD} f(x)$$

$e$ 將導數運算子編織成預測工具，奠定「羅塞達石碑」的連續基礎。

在離散世界，電腦處理序列 $y_n = f(t_n)$，其中 $t_n = t_0 + n h$。關鍵運算子包括：
- 移位運算子 $z$: $z y_n = y_{n+1}$ (前進一格)
- 延遲運算子 $z^{-1}$: $z^{-1} y_n = y_{n-1}$ (後退一格，適合即時系統)
- 前向差分 $\Delta$: $\Delta y_n = y_{n+1} - y_n$
- 後向差分 $\nabla$: $\nabla y_n = y_n - y_{n-1}$

這些運算子相互關聯：
$\Delta = z - 1$
$\nabla = 1 - z^{-1}$
$z = 1 + \Delta = \frac{1}{1 - \nabla}$
$\Delta = \nabla z$
這些關係為離散與連續的橋接提供操作基礎。

在連續世界，$f(t + h) = e^{hD} f(t)$ 描述函數隨時間推進。在離散世界，$y_{n+1} = z y_n$。若 $y_n = f(t_n)$，則 $y_{n+1} = f(t_n + h)$。

英國物理學家和電子工程師 [Oliver Heaviside](https://en.wikipedia.org/wiki/Oliver_Heaviside) 透過運算子微積分洞察，發現：
$$
D = \frac{1}{h} \ln(z) \tag{10}
$$

$e^{hD}$ 和 $\ln(z)$ 需透過冪級數展開為實用形式：
- $e^{hD} = 1 + h D + \frac{(h D)^2}{2!} + \frac{(h D)^3}{3!} + \cdots$
- 使用 $z = 1 + \Delta$：
$\ln(z) = \ln(1 + \Delta) = \Delta - \frac{\Delta^2}{2} + \frac{\Delta^3}{3} - \cdots$
- 使用 $z = \frac{1}{1 - \nabla}$：
$\ln(z) = -\ln(1 - \nabla) = \nabla + \frac{\nabla^2}{2} + \frac{\nabla^3}{3} + \cdots \tag{11}$
- 外插展開：
$z = \frac{1}{1 - \nabla} = 1 + \nabla + \nabla^2 + \nabla^3 + \cdots \tag{12}$

雙線性變換進一步實用化。設 $D = s$ (拉普拉斯域) ，則：
$$
s = \frac{1}{h} \ln(z) \approx \frac{2}{h} \left( \frac{z - 1}{z + 1} \right)
$$

這將連續傳遞函數轉換至離散域，廣泛應用於數位濾波器設計。

### 應用：汽車巡航控制
數值微分從離散資料估計導數。假設汽車巡航系統每 $0.1$ 秒取樣速度：
$$
v_0 = 60 \, \text{km/h}, v_1 = 61 \, \text{km/h}, v_2 = 61.8 \, \text{km/h}
$$

使用 $D v_n = \frac{1}{h} \ln(z) v_n$ 進行近似：
$$
D v_1 \approx \frac{v_1 - v_0}{h} = \frac{61 - 60}{0.1} = 10 \, \text{km/h/s}
$$

為提升精度，使用後向差分：
$$
D v_2 = \frac{1}{h} \left( \nabla + \frac{\nabla^2}{2} \right) v_2 \quad (13)
$$

其中
- $\nabla v_2 = 61.8 - 61 = 0.8$
- $\nabla^2 v_2 = 0.8 - 1 = -0.2$

$D v_2 \approx \frac{1}{0.1} (0.8 - 0.1) = 7 \, \text{km/h/s}$

微控制器透過環形緩衝區儲存歷史資料，即時計算加速度，調整油門以保持平順駕駛。

### 應用：音訊處理與 ODE 求解
數值積分求解常微分方程 (ODE) $\frac{dx}{dt} = f(t, x)$。

假設音頻取樣率 44100 Hz ($h \approx 22.68 \, \mu\text{s}$)，聲壓 $p(t)$ 需計算 $p_{n+1}$：
$$
p_{n+1} = p_n + h \left[ \frac{z - 1}{\ln(z)} \right] f_n \tag{14}
$$

[亞當斯預測法](https://en.wikipedia.org/wiki/Linear_multistep_method#Adams%E2%80%93Bashforth_methods)提供高精度：
$$
p_{n+1} = p_n + \frac{h}{24} [55 f_n - 59 f_{n-1} + 37 f_{n-2} - 9 f_{n-3}] \quad (15)
$$

假設 $p_0 = 0, f_0 = 1, f_1 = 1.1, f_2 = 1.2, f_3 = 1.3$，迭代計算提升精度，適用於 DSP 即時濾波。結合雙線性變換，可設計數位濾波器消除噪聲。

羅塞達石碑 $z = e^{hD}$ 是連續與離散世界的紐帶。從泰勒級數起步，透過運算子微積分，$e$ 孕育數值微分、外插與積分的實用演算法，成為數位訊號處理與控制系統的基石。該「石碑」讓數位世界得以精確地駕馭現實中的連續動態。

## 單一運算子產生所有初等函數
前述各節反覆強調 $\exp$ 與 $\ln$ 是連續數學的基本配對。乘法可化為加法 ($x \times y = e^{\ln x + \ln y}$)，三角函式可透過歐拉公式化為複指數，而羅塞達石碑 $z = e^{hD}$ 更將微分運算子與離散移位運算子橋接起來。下一個的問題是：這對 exp-log 配對究竟能走多遠？能否從中萃取出一個單一的原始運算子 (primitive operator)，使其在理論上足以產生所有[初等函數](https://en.wikipedia.org/wiki/Elementary_function)？

### 從 NAND 到 EML：離散與連續的功能完備
在數位邏輯中，NAND 閘是功能完備性的元件：僅靠 NAND 即可組合出任意布林操作。[Nand2Tetris](https://www.nand2tetris.org/) 課程具體示範這條 bootstrapping 鏈：從 NAND 出發，先自組合得到 NOT ($\mathrm{NAND}(a,a) = \overline{a}$)，再組合出 AND、OR、MUX，逐層堆疊至 ALU、CPU，最終構建完整的通用電腦。整個過程僅需一種閘，加上固定的接線規則。NAND 之所以能扮演此角色，關鍵在於二個性質：
- 透過輸入綁定實現一元反相運算: $\mathrm{NAND}(a,a) = \overline{a}$，從而取得 NOT
- 二輸入結構提供足夠的組合自由度，覆蓋所有布林操作

然而在連續數學中，一直缺乏對應的結果。計算 $\sin$、$\cos$、$\sqrt{\phantom{x}}$、$\log$ 等初等函數，過往總需要多種不同的運算。

2026 年 4 月，Jagiellonian University 理論物理研究所的 Andrzej Odrzywolek 在預印本 [All elementary functions from a single operator](https://arxiv.org/abs/2603.21852) 中發現連續數學的對應物。定義二元運算子 EML (Exp-Minus-Log)：
$$
\mathrm{eml}(x, y) = \exp(x) - \ln(y)
$$

搭配單一常數 $1$ (計算函式時另加輸入變數 $x, y, \ldots$ 作為終端符號)，即可透過反覆組合生成所有標準初等函數，涵蓋四則運算 ($+$、$-$、$\times$、$/$)、冪
運算、根號、三角與雙曲函式及其反函式，以及常數 $e$、$\pi$、$i$。

### 為何 exp 與 ln 的組合足以產生一切
EML 之所以可行，根源在於 $\exp$ 與 $\ln$ 這對互逆函式恰好橋接連續數學的三層結構：
- 加法 $\leftrightarrow$ 乘法：$\ln$ 將乘法降為加法 ($\ln(xy) = \ln x + \ln y$)，$\exp$ 反向提升 ($e^{a+b} = e^a \cdot e^b$)。該論文前言 Eq. (1) 即以此說
明 exp-log 可互相轉換加法與乘法；在純 EML 中，各運算則由窮舉搜尋得到的 EML tree 直接重建
- 實數 $\leftrightarrow$ 複數：歐拉公式 $e^{i\theta} = \cos\theta + i\sin\theta$ 使所有三角與雙曲函式歸結為複指數。EML 內部在複數域上操作，例如 $\cos x = (e^{ix} + e^{-ix})/2$，$\sin x = (e^{ix} - e^{-ix})/(2i)$，皆可由 $\exp$、算術與虛數單位 $i$ 組合而成
- 非交換性：$\mathrm{eml}(x,y) \ne \mathrm{eml}(y,x)$，因為 $\exp$ 作用於第一個引數而 $\ln$ 作用於第二個。該論文指出，非交換運算對於搜尋最小生成集至關重要，因為它同時提供表達式樹的成長能力與反演能力 (inversion capability)，使得從 $\exp$ 復原 $\ln$ 得以透過巢狀組合實現。

EML 的 bootstrapping 起點如下：
- $e^x = \mathrm{eml}(x, 1)$，因為 $\exp(x) - \ln(1) = e^x$
- $\ln x = \mathrm{eml}(1, \mathrm{eml}(\mathrm{eml}(1, x), 1))$，即 $e - \ln(e^e/x) = \ln x$

有了 $\exp$ 與 $\ln$ 之後，其餘運算皆可重建。論文 Table 4 顯示各運算的最短 EML 表示 (以直接搜尋所得的節點數計)：減法 $x - y$ 為 11，乘法 $x \times y$
為 17，除法 $x/y$ 為 17，加法 $x + y$ 為 19。值得注意的是，在 EML 體系中減法比加法更基本，乘除法比加法更短，這與傳統運算的層級不同。

也因此，一台只有二個按鍵 (EML 與 1) 的計算器，在功能上等價於完整的工程計算器。
> [source](https://x.com/greijukainen/status/2043754534615306601/photo/1)

### 統一的語法與發現過程
EML 表達式遵循極簡的 context-free 文法：
$$
S \to 1 \mid \mathrm{eml}(S, S)
$$

所有初等函數表達式皆映射為同構於 Catalan 結構的 full binary tree (每個節點恰有 0 或 2 個子節點)。原本結構各異的數學運算被嵌入同一個組合空間，正如 NAND 將所有布林電路統一為同質閘陣列。

Odrzywolek 的發現過程採用系統性的 ablation testing：從工程計算器的 36 個標準原始元件 (8 個常數、20 個函式、8 個二元運算) 出發，逐一移除元件並驗證剩餘集合是否仍能重建所有功能。這一過程產生一系列漸次精簡的計算器組態 (Calc 3、Calc 2、Calc 1、Calc 0)，最終收斂至 EML。關鍵洞察在於：所有最小組態都涉及互逆函式對搭配非交換運算，而 EML 恰好將這二個要素融合在單一運算子中。驗證結合 Mathematica 的符號化簡與多語言 (C、NumPy、PyTorch、mpmath) 的數值交叉驗證。

### 類比的邊界
NAND 與 EML 的類比在 functional completeness 層次上精確成立，但在工程可行性上有本質差異。NAND 的實用價值不僅來自 completeness，更來自數位邏輯的訊號再生能力 (signal regeneration) 與雜訊邊限 (noise margin)：每一級 NAND 閘都將輸入重新整形為乾淨的高/低電位，使誤差不會隨組合深度累積。Nand2Tetris 之所以可行，正是因為數位特性保證任意深度的組合仍然邏輯正確。EML 則內含 $\exp$ 與 $\ln$，計算成本高且高度非線性，數值誤差隨樹深迅速累積，以 EML 組合出簡單運算往往需要冗長的樹結構。論文 Table 4 顯示，即便是加法 $x + y$ 的最短 EML 表示也需要 19 個節點，而 $\pi$ 常數需要超過 53 個節點。

論文另一個較具實驗意味的部分是 symbolic regression 的應用。由於所有 EML 表達式被限制為固定深度的二元樹，問題轉化為在結構固定但連續參數可最佳化的計算圖上訓練。在論文的固定深度 master formula 中，深度 $n$ 的 full binary tree 有 $2^n - 1$ 個內部節點與 $2^n$ 個葉節點，每個葉節點的輸入以線性組合 $\alpha_i + \beta_i x + \gamma_i f$ 參數化 (其中 $f$ 為前一層 EML 的輸出)。如此一來，可直接使用 Adam 等梯度最佳化器對數值資料擬合，並在某些低深度情況 (depth $\le$ 4) 成功復原精確的閉式解。然而隨著樹深增加，節點數呈指數成長，搜尋空間迅速增長，加上 $\exp/\ln$ 帶來的梯度消失與爆炸，實際復原率從深度 2 的 100% 驟降至深度 3-4 約 25%、深度 5 不足 1%。

EML 的存在為 $e$ 提供新的詮釋：$\exp$ 與 $\ln$ 不僅是描述連續變化的工具，它們的組合在理論上足以產生連續數學的整個初等函數體系。正如 NAND 揭示所有布林邏輯其實是同一種閘的排列組合，EML 揭示所有初等函數其實是 $\exp$ 與 $\ln$ 的巢狀組合。

## 後記
回首探索之路，歐拉數 $e$ 從 17 世紀 Jacob Bernoulli 的複利利息問題中初露鋒芒，經由 Leonhard Euler 的系統化確立，最終成為橫跨微積分、機率統計、數值分析乃至人工智慧的關鍵常數。

筆者於二十年前 (2005 年前後) 拜讀 Jack Crenshaw 博士的一系列文章時深受啟發，驚嘆於 $e$ 如何在看似迥異的領域 (如放射性衰變與巡航控制系統) 中展現出結構上的必然。本文旨在彙整這些深刻的連結：從自然對數底數的誕生、數值穩定的 `expm1` 實作，到 BitNet 量化與大型語言模型的壓縮，得以見證一個常數如何定義連續變化的極致美感，並持續推動現代工程的邊界。期望這份整理能作為一座橋樑，助讀者領略隱藏在自然規律與數位計算背後的簡潔與優雅。