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


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

本文回顧 $e$ 的起源，一探數學家如何定義它、發現它的奇妙性質 (例如 $e^x$ 的導數仍是自身)，及它何以稱為「自然對數」的底。讀者也將理解 $e$ 在碳 14 測定法、微積分、金融、機率統計、物理領域所奠定的根基，並探討其在嵌入式系統、人工智慧等資訊系統中的應用。讀者最終會體會到，$e$ 如何優雅地幫助我們理解這個世界的運作。

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

:::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$，歐拉常數的定義是調和級數與自然對數的差值。論及物理與工程系統時，[Euler number](https://en.wikipedia.org/wiki/Euler_number_(physics)) (縮寫為 $Eu$) 與數學中 Euler's number 完全不同，前者出現在流體力學與應用物理中，後者則源自數學分析。

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

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

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

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

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

再進一步，若銀行每 4 個月 (即 $\tfrac{1}{3}$ 年) 給付一次利息 (下圖紅色、紫色圓圈):

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

你或許會貪婪地推論：倘若銀行每分每秒都給付利息，是否會使存款金額趨近無窮大呢？試算如下：
| $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.7182817813$ 元，並非原本貪婪設想的無窮大。換言之，上述計算得到的數值相當接近歐拉數 $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 首先研究複利頻率對本金增長的影響，發現上述極限的存在。然而，真正將此數深入研究並證明其無理性質的則是瑞士數學家 Leonhard Euler，他於公元 1748 年的著作《[Introductio in Analysin Infinitorum](https://en.wikipedia.org/wiki/Introductio_in_analysin_infinitorum)》詳述該數的性質，後人遂以 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!}
$$

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


在物理中的[波茲曼常數](https://en.wikipedia.org/wiki/Boltzmann_constant) 描述氣體分子在重力場中或一般熱平衡狀態下的能量分佈，表現為 $e^{-E/kT}$，其中 $E$ 為粒子能量，$k$ 為玻爾茲曼常數，$T$ 為絕對溫度。統計力學中的[馬克士威-波茲曼分布](https://en.wikipedia.org/wiki/Maxwell%E2%80%93Boltzmann_distribution)描述氣體分子速度的機率分佈，其中指數函數的形式反映高能分子出現的可能性隨速度平方指數遞減的特性。隨機運動方面，[布朗運動](https://en.wikipedia.org/wiki/Brownian_motion)是微小顆粒在流體中因為分子碰撞而產生的隨機運動，其數學模型依賴於指數衰減函數來描述粒子擴散的機率分佈。

在核物理學中，[放射性衰變](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 (Carbon-14) ，也具備同樣的衰變行為。

> 圖片出處: [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}$ 定年法，是種利用碳的放射性同位素 $^{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$) 而逐漸減少。$^{14}\text{C}$ 的半衰期 (Half-life) 約為 5730 年，意即每經過這段時間，樣品中 $^{14}\text{C}$ 的數量會衰減為原來的一半，此方法可有效測定約 5 萬年內的樣本年代。

$^{14}\text{C}$ 定年法自公元 1949 年由 Arnold 與 Libby 證實可行性後，其在定年方面的實用性及優越性已被廣泛肯定，對考古學產生革命性影響，為遺址和文化演進提供精確的時間框架。在地質學等領域，它也廣泛應用於測定晚第四紀 (如末次冰期結束) 以來的地質事件年代。

半衰期是描述衰變快慢的常用概念，符號為 $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$ 這個微小的時間間隔內，發生衰變的原子核數量 $dN$ (注意 $dN$ 在這裡是負值，因為數量減少) 會大約等於「總原子核數量 $N$」乘以「單一原子核在此時間內衰變的機率 $p$」，亦即:
$$
dN \approx - 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}$ 或更複雜的運算子) 的特徵函數和特徵值緊密相關。

特徵函數通常代表系統的一種「固有」行為模式或狀態。例如，在振動系統中，它們可能代表特定的自然振動頻率；在量子力學中，定態薛丁格方程式 $H\psi = E\psi$ 其實就是個典型的特徵值問題：能量運算子 $H$ 作用在系統的穩定狀態波函數 $\psi$ (特徵函數) 上，得到能量值 $E$ (特徵值) 乘以該波函數。理解這種關係，是理解量子世界本質的關鍵。

回到上述碳-14 衰變問題。$\frac{dN}{dt} = -\lambda N$ 這個看似單純的微分方程描述的是「變化率與目前數量成正比」的過程。這不僅適用於放射性衰變，也同樣適用於銀行存款的連續複利計算、物體的冷卻過程 (即牛頓冷卻定律) 等眾多自然與經濟現象。指數函數 $e^{kt}$ 正是描述這類過程最自然、最貼切的數學語言。

微積分的建立，正是透過對這種連續變化的嚴謹研究與描述。$e$ 並非僅是個數學常數，而是貫穿於自然規律之中的根本法則。愛因斯坦曾提出一個深刻的疑問：「數學，它明明是人類心智的創造，與實際經驗並無直接關聯，那它究竟是如何能夠如此精準地描述我們所處的現實世界？」(How can it be that mathematics, being after all a product of human thought independent of experience, is so admirably adapted to the objects of reality?)

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

## 自然底數的由來
$e$，又稱為「自然底數」，是自然對數 (以 $e$ 為底的對數) 所對應的基礎常數。

雖然在人口成長與利息複利等應用問題中，可用指數函數進行推演，但 $e$ 並非直接來自這些應用，而是微積分領域中的一項關鍵發現。先前我們藉由極限的方式給出 $e$ 的定義，不過，這並非科學史上研究自然對數的起點。

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

公元 1614 年，John Napier 發表對數與對數表，而直到公元 1637 年，法國數學家笛卡兒才引入現代指數記號系統，也就是說，對數的出現，竟然早於指數的形式化定義。真正釐清二者之間的關係，則已是百餘年後的公元 1770 年，由歐拉首次明確指出：「對數乃指數的反運算」。

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

對數的誕生，汲取積化和差的主體想法，卻將其精煉為更通用、更易操作的數學工具。只要藉由 $x \cdot y = 10^{\log x + \log y}$，即可將乘法化為單純的加法與一次反查，操作上極為高效。流程如下：
- 查對數表二次
- 加總對數值
- 查一次反對數表，即得乘積

一張對數表足以大幅加速乘除與開根號等運算，其後演化出的對數尺，更是工程師必備工具，猶如現今理工科學生獲贈的科學計算器。
然而，當時的對數計算主要用[常用對數](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 是等比例遞減。接下來想要找一個數學函數去描述這兩個質點位置的關係。先寫成簡單的表格。
| 時間 | $\beta$ 位置 | $\beta$ 速度 | log($\beta$ 位置) = B位置 | log($\beta$ 位置) 變化=B速度 |
| -------- | ------| - | - | - |
| 0 | $\alpha =x = sin(\theta)$ | 速度等比例遞減 | $log(x)=y$ | 等速 |
| 1 | $\gamma = x - \frac{x}{2} = sin(\theta')$ | $\frac{x}{2}$ | $log(\frac{x}{2})$ | $log(\frac{x}{2}) - log(x)$ |
| 2 | $\delta = x - \frac{x}{2} - \frac{x}{4}$ | $\frac{x}{4}$ | $log(x/4)$ |$log(\frac{x}{2}) -log(\frac{x}{4})$ |
| 3 | $\epsilon = x - \frac{x}{2} - \frac{x}{4} - \frac{x}{8}$ | $\frac{x}{8}$ | $log(x/4)$ | $log(\frac{x}{4}) - log(\frac{x}{8})$ |

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

接下來看這段在[〈指數與對數發展史簡介〉](https://www.sec.ntnu.edu.tw/uploads/asset/data/66cf6cde0e0b30672f4bad39/1983-056-10_58-73_.pdf)中的描述。
> 〈指數與對數發展史簡介〉乙、對數的發明 二、Napier對數
> Napier 引進對數的方法並不像前面所說的那麼代數化，他的方法是幾何式的。下面我們來介紹他的定義方法 : 取二線 $\overrightarrow{AB}$ 與 $\overrightarrow{CD}$，令 $\overrightarrow{AB}$ 之長為 $10^7$ 。動點 $P$ 沿 $\overrightarrow{AB}$ 由 $A$ 點移動至 $B$ 點，動點 $Q$ 沿 $\overrightarrow{CD}$ 由 $C$ 點移動至 $D$ 點。$P$ 點的初速度為 $10^7$ ，而在位置 $P$ 時的速度為 $\overrightarrow{PB}$ 。 $Q$ 點 $\overrightarrow{CD}$ 上保持等速，其初速度也是 $10^7$，假設兩個動點同時出發，某段時間後兩動點分別到達 $P$ 與 $Q$ 的位置，他定義 $\overrightarrow{CD}$ 為 $\overrightarrow{PB}$ 的對數。

示意圖 : 

分析上面的話得到下面的資訊 :
- $\overrightarrow{AB}$ 長度 $10^7$

如果 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)$，表示等速。但當時就是為了想要用數學描述這兩者間的關係而有這個對數的數學式子。

圖中展示 Napier 定義對數的運動學模型：想像二個質點沿線段運動，上方 Moving $b$ Particle，等速運動，其軌跡對應函數是 $y = \log(x)$，而下方 Moving $\beta$ Particle，其軌跡對應函數則是 $x = \sin(\theta)$，皆由左至右移動，對應比例常數 $R = \text{Sinus Tontus} = 10,\!000,\!000$。每個對應位置以垂直線段連結，展現在不同變化速率下的對數關係。上方標記點 $A–I$ 與下方的 $\alpha–\omega$ 構成一對一映射，呈現二條運動軌跡在時間參數下的幾何對應。該模型巧妙地將等差運動與等比運動聯繫起來，反映他將對數視為描述運動變化的比例關係，並與當時天文計算所需之三角函數計算緊密結合。

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

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

數學家在編撰對數表的過程中，面對繁瑣的計算步驟，意識到若採用以下形式作為底數：
$$
\Bigl(1 + \frac{1}{n}\Bigr)^n, \quad n \gg 0
$$
可顯著簡化計算。Napier 採用的底數近似於 $\Bigl(1 - \frac{1}{10^7}\Bigr)^{10^7}$ (此值近似 $1/e$)，而 Bürgi 則選擇 $\Bigl(1 + \frac{1}{10^4}\Bigr)^{10^4}$ (此值近似 $e$)。儘管當時尚未明確建立「底數」的概念，他們已運用接近 $e$ 或其倒數的值，儘管 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 |

我們從中發現規律：若底數形如 $b^n$，且 $n$ 為一大整數，則底 $a$ 越接近 1，真數間距就越均勻，計算上也更合理。若我們選擇 $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 正是採用該方法。要計算出 $\log_a x = 0.000N$ 對應的真數 $(1.0001)^N$，需要進行 $N$ 次乘法。因此，要填滿表格至 $\log_a x = 1.0000$ (即 $N=10000$)，總共需要的乘法次數約為 $\sum_{N=1}^{10000} N = \frac{10000 \cdot (10000 + 1)}{2} = 50,\!005,\!000$，已是數千萬次的規模，全靠手工運算。這樣的設計讓早期的對數表建構者得以確保高精度，進而利用以下換底公式推導出常用對數：
$$
\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$ 計算到所需精度。

換言之，構造高精度對數表的訣竅，在於選擇形如 $(1 + \frac{1}{n})^n$ 的底數，並讓 $n$ 足夠大，如 $10000$ 或以上，如此方能兼顧計算便利與精度要求。

接著我們用現代的數學語言來選擇 $a = 1 + c$，其中 $c$ 為極小的正數 (例如 $c = 10^{-7}$)，並觀察這樣的底數下，指數變化對應的函數值是否能逼近等差關係。同理，也可考慮 $a = 1 - c$ 的情況，以對稱地探討正負微小變動對函數行為的影響。

為了觀察當底數趨近於 $1$ 時，指數函數是否能展現近似等差的變化，我們考慮底數為 $a = 1 \pm c$ 的二種情況，其中 $c$ 是極小正數，例如 $c = 10^{-7}$。分別計算下列：
- 若 $a = 1 + c$：
$$
\begin{aligned}
(1 + c)^0 &\approx 1 + 0 \cdot c = 1 \\
(1 + c)^1 &\approx 1 + 1 \cdot c \\
(1 + c)^2 &\approx 1 + 2 \cdot c \\
(1 + c)^3 &\approx 1 + 3 \cdot c
\end{aligned}
$$
- 若 $a = 1 - c$：
$$
\begin{aligned}
(1 - c)^0 &\approx 1 - 0 \cdot c = 1 \\
(1 - c)^1 &\approx 1 - 1 \cdot c \\
(1 - c)^2 &\approx 1 - 2 \cdot c \\
(1 - c)^3 &\approx 1 - 3 \cdot c
\end{aligned}
$$

這些近似來自以下二項式展開：
$$
(1 \pm c)^k = 1 \pm kc + \frac{k(k - 1)}{2} c^2 + \cdots \approx 1 \pm kc
$$

其中忽略高階項是因為 $c \ll 1$。此處我們明確觀察到：儘管左側為等比級數，其近似值卻對應到右側的等差級數。再次，我們得知，將等比變化近似對應為等差變化，正是對數設計的主體想法。

設 $x = a^b$。為了將等比關係 $x$ 映射到近似等差的值 $y$，我們嘗試探索一種關係，使得 $y$ 的微小變化對應於 $b$ 的微小變化，且與 $c$ 相關。考慮令 $b = y / (\pm c)$，如此 $y$ 的單位變化就對應 $b$ 的 $1/(\pm c)$ 變化。代入 $x = (1 \pm c)^b$：
$$
x = \left( (1 \pm c)^{\frac{1}{\pm c}} \right)^y
$$

在極限情況下，當 $c \to 0$ 時，基底 $(1 \pm c)^{1/(\pm c)}$ 趨近於 $e$ (利用 $\lim_{n\to\infty} (1 + 1/n)^n = e$ 以及 $\lim_{n\to\infty} (1 - 1/n)^{-n} = e$)。因此，極限關係為 $x = e^{\pm y}$，其反函數即為 $y = \pm \ln x$。

為了展現上述極限行為，我們可用先前案例中的極小值 $c=10^{-7}$ (即 $\frac{1}{c} = 10^7$) 來近似，分別對應於以下函數：
- 對於 $a = 1 + 10^{-7}$，關係近似 $x = (1.0000001)^{10^7 y} = [(1+10^{-7})^{10^7}]^y \approx [e]^y = e^y$。其反函數近似為 $y \approx \ln x$
- 對於 $a = 1 - 10^{-7}$，關係近似 $x = (0.9999999)^{10^7 y} = [(1-10^{-7})^{10^7}]^y \approx [e^{-1}]^y = e^{-y}$。其反函數近似為 $y \approx -\ln x$

這些近似函數幾乎與自然對數函數 $y = \ln x$ 以及 $y = -\ln x$ 重合，如圖一：
> 圖一: $y = \ln{x}$

$x = 0.9999999^{10^7 y}$ 幾乎與 $y = -\ln{x}$ 重合，如圖二：
> 圖二: $y = -\ln{x}$

為更嚴謹地刻劃這樣的關係，以下探討函數 $f(x)$ 的變化率。設 $x = a^b$，則對應的值為 $y = f(x)$。以差商形式近似導數：
$$
f'(x) \approx \frac{f(x) - f(x / a)}{x - x / a} = \frac{\Delta y}{x (1 - 1 / a)}
$$

對於 $a = 1 \pm c$，有：
$$
f'(x) \approx \frac{\pm c}{x (1 - 1 / (1 \pm c))} = \frac{\pm c}{x \left(1 - \frac{1}{1 \pm c} \right)}
$$

簡化得：
$$
f'(x) \approx (1 \pm c) \cdot \frac{1}{x}.
$$

當 $c$ 趨近於 $0$，有 $1 \pm c \approx 1$，因此導數趨近於：
$$
f'(x) \approx \frac{1}{x}
$$

接下來對 $x$ 積分可得函數形式：
$$
f(x) = \int \frac{1}{x} \, dx = \ln x + C
$$

為滿足初始條件 (例如 $f(1) = 0$) ，我們取 $C = 0$，因此：
$$
f(x) = \ln x
$$

這條推導路徑顯示，當底數趨近 $1$ 時，指數函數之變化率呈現出對應等差關係的特徵，進而引出自然對數 $\ln x$ 作為其反函數的自然形式。

換言之，我們已成功從「將幾何級數近似為算術級數」的出發點，建構出自然對數的定義基礎。這也說明 John Napier 當初選擇底數極接近 $1$ 的實際考量。他在尚未發展出微積分與極限概念的時代，透過大量近似運算與觀察，已隱然掌握 $e$ 所代表的極限形式，並以此建構對數表，簡化複雜的連乘與比例運算。

Napier 在微積分尚未成形的時代，獨力以二十年反覆運算與觀察，建立早期的對數表，並在過程中觸及 $e$ 的極限形式。雖他未能系統化地提出其理論，卻已隱然揭示 $e$ 的數學意義，這也是為何 $e$ 有時稱為「Napier 常數」，儘管不精確。他過世後出版的《De arte logistica》(直譯為「計算之術」) 展現其關注如何透過數學工具實際提升計算效率，與現代數學強調結構與理論的取向相當不同。

在十七世紀微積分萌芽的年代，介於 Napier 與牛頓、萊布尼茲之間的數學家，逐步將 $e$ 納入連續變化與極限計算的架構中，使其從早期的計算工具演變為現代分析、機率與數值方法的關鍵常數。自然對數的誕生，反映科學的計算需求，也是人類對「連乘極限」與「變化率」本質的探索。
> 延伸閱讀: [對數與約翰．納皮爾 (John Napier)](https://episte.math.ntu.edu.tw/articles/mm/mm_03_4_07/)

費馬 (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) 的結果是含虛數的對數，進而在三角函數與所謂的「橢圓對數」之間建立關聯。需要注意的是，John 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 的研究中，我們不難發現，無窮級數的展開對許多函數而言都至關重要。儘管主要是使用冪級數展開，但更為人所知的泰勒展開和牛頓展開，在本質上只差若干系數或運算步驟，後者則在隨後的研究中獲得更廣泛的應用。

延伸閱讀:
* [自然對數漫談](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$，則 $f'(x) = f(x)$。

為何要如此定義呢？

在探究指數函數 (例如以 2 為底的 $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$，也可透過另一個等價的極限來定義：
$$
a = \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$，使得 $\frac{d}{dx} e^x = e^x$，這也是所謂的「自然」。

因此，$e$ 的出現本質上是為了處理「指數函數的導數究竟該怎麼算」這個問題。它之所以與許多自然及社會現象相關，並非 $e$ 本身有何特殊魔力，而是因這些現象普遍遵循「成長 (或衰減) 率與目前數值呈正比」的數學規律。

至於「為何用 $e$ 這個底數？」則是因為 $e^x$ 在計算導數或積分的時候最為方便，對任何其他底數，都能轉寫成與 $e$ 相關的形式。$e$ 因其「導數等於函數自身」的獨特性，成為所有指數函數的基礎，於是 $e$ 不僅是個數值，而是讓指數函數在微分與積分運算中最自然、最簡潔的底數。

延伸閱讀:
* [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)

## 證明極限存在
牛頓在發展微積分的過程中，研究對數函數的性質，牛頓利用他發展的廣義二項式定理，推導出 $\ln(1+x)$ 的冪級數展開式 (大約在公元 1665 年，但較晚發表)：
$$
\ln(1+x) = x - \frac{x^2}{2} + \frac{x^3}{3} - \frac{x^4}{4} + \cdots = \sum_{n=1}^\infty (-1)^{n-1} \frac{x^n}{n}
$$

這個結果是藉由對 $\frac{1}{1+x} = 1 - x + x^2 - x^3 + \cdots$ 逐項積分而得到。

雖然牛頓主要處理的是對數函數，但他意識到對數函數存在一個反函數：若 $y = \ln x$，於是 $x$ 就是 $y$ 的某種「反對數」函數。藉由對 $\ln(1+x)$ 級數的探索，牛頓間接得到該反函數的級數形式，也就是我們今日所知的指數函數 $\operatorname{Exp}(x)$ 的級數展開：
$$
\operatorname{Exp}(x) = 1 + x + \frac{x^2}{2!} + \frac{x^3}{3!} + \cdots
$$

然而，在牛頓的時代，數學家著重於函數本身的級數表示及其作為對數反函數的角色，但不會特別關注或定義這個級數在 $x=1$ 處的值 (即 $e$) 作為一個獨立的、基礎性的數學常數。他為指數函數奠定級數基礎，但 $e$ 這個數本身尚未系統性討論。

真正將自然常數 $e$ 和指數函數 $e^x$ 進行系統研究並賦予其現代形式者，是瑞士數學家歐拉，其貢獻在於：
1. 命名與符號：歐拉首次使用字母 $e$ 來表示這個常數，約在公元 1727 或 1728 年
2. 確立級數定義：歐拉運用牛頓等人得到的指數函數級數展開。他將指數函數 $\operatorname{Exp}(x)$ (他寫作 $e^x$) 的定義牢固地建構於冪級數之上：
$$
e^x = \operatorname{Exp}(x) := \sum_{n=0}^\infty \frac{x^n}{n!} = 1 + \frac{x}{1!} + \frac{x^2}{2!} + \frac{x^3}{3!} + \cdots
$$
該級數定義是從指數函數最根本的性質 (即它是微分方程 $\frac{dy}{dx}=y$ 且滿足 $y(0)=1$ 的唯一解) 出發，假設解可表示為冪級數而導出) 。
3. 統一不同定義：歐拉證明基於級數定義的 $e^x$ 函數，其在 $x=1$ 處的值：
$$
e = e^1 = \sum_{n=0}^\infty \frac{1}{n!} = 1 + 1 + \frac{1}{2!} + \frac{1}{3!} + \cdots
$$
確實等於先前由連續複利極限所定義的常數：
$$
\sum_{n=0}^\infty \frac{1}{n!} = \lim_{n \to \infty} \left(1 + \frac{1}{n}\right)^n
$$

歐拉的級數收斂很快，便於精確計算 $e$ 的值，且該級數形式在數學分析中極易處理 (如求導數函數、積分、證明性質等) 。因此，在現代數學分析中，通常將 $e^x = \sum \frac{x^n}{n!}$ 作為指數函數的基本定義。

以下證明極限 $\lim_{n \to \infty} \left(1 + \frac{1}{n}\right)^n$ 確實存在，正因極限存在，我們才能將 $e$ 定義為常數。

**證明思路**
觀察 $\left(1 + \frac{1}{n}\right)^n$，隨著 $n$ 越大，其當前項相比前項越大，換言之，有序數列 $\lbrace \left(1 + \frac{1}{i}\right)^i \rbrace$ 是遞增數列。若再說明該有序數列有一上界存在，則可以使用單調遞增收斂定理說明該有序數列收斂至一最小上界。

首先，證明一個不等式。對於任意滿足 $b > a > 0$ 的實數 $a$ 和 $b$，我們希望證明：
$$
b^n - a^n < (b - a)\, n\, b^{n-1}
$$

該不等式成立的原因如下：
$$
\begin{align}
b^n - a^n
&= b^n\left[1-\left(\frac{a}{b}\right)^n\right] \\
&= b^n\left[1-\left(\frac{a}{b}\right)\right]\left(1+\left(\frac{a}{b}\right)+\left(\frac{a}{b}\right)^2+\dots+\left(\frac{a}{b}\right)^{n-1}\right) \\
&= b\left[1-\left(\frac{a}{b}\right)\right]b^{n-1}\left(1+\left(\frac{a}{b}\right)+\left(\frac{a}{b}\right)^2+\dots+\left(\frac{a}{b}\right)^{n-1}\right) \\
&= (b - a)\left(b^{n-1} + b^{n-2}a + b^{n-3}a^2 + \cdots + b\,a^{n-2} + a^{n-1}\right)
\end{align}
$$
其中，第二個等式成立使用到因式分解 $(1-x^n)=(1-x)(1+x+x^2+\dots+x^{n-1})$。隨後我們分段分析該等式左側的乘數 $\left(b^{n-1} + b^{n-2}a + b^{n-3}a^2 + \cdots + b\,a^{n-2} + a^{n-1}\right)$ 由於 $a < b$，則，觀察乘數中每一項的不等式可以得到以下不等式關係
$$
\begin{align}
&ab^{n-2} > bb^{n-2} = b^{n-1} \\
&a^2b^{n-3} > b^2b^{n-3} = b^{n-1} \\
&\vdots \\
&a^{n-1} > b^{n-1}
\end{align}
$$
由於 $b > a > 0$ 皆為正數，故不等式加總不會變號，則乘數 $\left(b^{n-1} + b^{n-2}a + b^{n-3}a^2 + \cdots + b\,a^{n-2} + a^{n-1}\right)$ 有以下不等式關係
$$
\left(b^{n-1} + b^{n-2}a + b^{n-3}a^2 + \cdots + b\,a^{n-2} + a^{n-1}\right) <nb^{n-1}
$$

兩側不等式乘上被乘數 $(b-a)$ 則得到
$$
b^n-a^n=(b - a)\left(b^{n-1} + b^{n-2}a + b^{n-3}a^2 + \cdots + b\,a^{n-2} + a^{n-1}\right) < (b-a)nb^{n-1}
$$
移項並整理可得：
$$
a^n > b^{n-1} \left[ b - (b - a)n \right] \tag{5}
$$

>:warning: 注意
>$b - (b - a)n$ 是特意構造的，目的是帶入特定的代數 $b>a>0$ ，以使得該表達式 $b - (b - a)n=1$ 這能讓我們將不等式再化簡。

接下來，構造一個遞增數列。令整數 $n > 1$，並指定
$$
a = 1 + \frac{1}{n}, \quad b = 1 + \frac{1}{n - 1}
$$

由於 $b > a > 0$，可將其代入不等式 (5)，先處理 $(b - a)n$ 則：
$$
(b - a)n = \left[\left(1 + \frac{1}{n - 1}\right) - \left(1 + \frac{1}{n}\right)\right]n= \frac{1}{n - 1}
$$

所以 $b-(b - a)n = 1$，最終得到：
$$
\left(1 + \frac{1}{n}\right)^n > \left(1 + \frac{1}{n - 1}\right)^{n - 1}
$$

這表示 $\left(1 + \frac{1}{n}\right)^n$ 隨著 $n$ 單調遞增。

最後，證明數列有[上界](https://en.wikipedia.org/wiki/Upper_and_lower_bounds) (upper bound)。指定
$$
a = 1, \quad b = 1 + \frac{1}{2(n - 1)}
$$

將其代入不等式 (5)，先處理 $(b - a)n$ 則：
$$
(b - a)n = \frac{n}{2(n - 1)}
$$

所以 $b -(b - a)n = \frac{1}{2}$

最終可得不等式：
$$
2 > \left(1 + \frac{1}{2(n - 1)}\right)^{n - 1}
$$

兩邊平方，得到：
$$
4 > \left(1 + \frac{1}{2(n - 1)}\right)^{2(n - 1)}
$$

由於這對任意 $n > 1$ 成立，可知數列 $x_k = \left(1 + \frac{1}{k}\right)^k$ 中，指數為偶數 $k=2(n-1)$ 的子序列 $x_{2(n-1)}$ 均小於 $4$。又因已知整個數列 $x_m = \left(1 + \frac{1}{m}\right)^m$ 是單調遞增，對於任何 $m$，總能找到一個偶數 $k > m$，使得 $x_m < x_k < 4$。因此，數列 $\left(1 + \frac{1}{m}\right)^m$ 存在上界 (例如 $4$)。

綜合上述結果，我們已證明：
- 有序數列 $\lbrace\left(1 + \frac{1}{n}\right)^n\rbrace_n$ 單調遞增；
- 且存在上界。

上述關係為[單調有界數列收斂原理](https://en.wikipedia.org/wiki/Monotone_convergence_theorem) (Monotone convergence theorem)的[充分條件](https://en.wikipedia.org/wiki/Necessity_and_sufficiency) (sufficiency)，則根據單調有界數列收斂原理，該數列必然收斂，因此極限存在且等於一最小上界，則令該最小上界為 $e = \lim_{n \to \infty} \left(1 + \frac{1}{n}\right)^n$。

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

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

歐拉數的泰勒展開如下：
$$
\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}
$$

根據指數函數的圖形：

我們可推想，若有某個 $a$ 讓在 $x = 0$ 處的切線斜率剛好是 $1$，於是該函數的導數就會處處等於函數本身。這個 $a$ 就定義為 $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.00236299$。

可發現，當 $x$ 越接近 $1$ 時，這種方法計算 $\ln{x}$ 的誤差越小。若把 $y = \ln{x}$ 與 $y = (\sqrt[2^k]{x} - 1) \times 2^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.00766%。不僅展示 $e$ 與錯位排列機率之間的關聯，也說明 Minecraft 作為計算模擬與實驗平台的潛力，將抽象的數學概念轉化為遊戲中高度視覺化的趣味過程。
展示影片: [Approximation of e in Minecraft](https://youtu.be/YH-wx5YQftY) 及[播放清單](https://www.youtube.com/playlist?list=PLscpLh9rN1Rf-dqLO4r3GAwvm1O_xL7D1)。

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

該級數具有快速收斂性，若期望小數點後達 $d$ 位正確，我們需確保截尾誤差 $\sum_{k>n} 1/k! < 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$ 位精度為例，計算所需項數為約 250 萬。我們採取 [GMP](https://gmplib.org/) 多精度整數 (`mpz_t`) 與浮點數 (`mpf_t`) 為基礎進行實作，主要流程如下：
* 使用遞迴函式 `calc(mpz_t a, mpz_t b, mpz_t p, mpz_t q)` 計算 $\sum_{k=a}^{b-1} 1/k!$ 的精確分數表示 $p/q$ (其中 `p`, `q` 為輸出參數)
* 遞迴分解：將計算區間 $[a, b)$ 對半分割為 $[a, m)$ 與 $[m, b)$，其中 $m = \lfloor(a+b)/2\rfloor$
* 遞迴計算：
* 呼叫 `calc(a, m, p_L, q_L)` 計算左半區間和 $p_L/q_L$ (此處 $p_L$, $q_L$ 對應程式碼中的輸出參數 `p`, `q1`)
* 呼叫 `calc(m, b, p_R, q_R)` 計算右半區間和 $p_R/q_R$ (此處 $p_R$, $q_R$ 對應程式碼中的輸出參數 `p2`, `q`)
* 合併結果：將左右兩區間的和相加 $\frac{p_L}{q_L} + \frac{p_R}{q_R} = \frac{p_L q_R + p_R q_L}{q_L q_R}$。程式碼中透過以下操作將結果存回輸出參數 `p` 與 `q` (此處 `p` 重用 $p_L$ 的儲存空間，`q` 重用 $q_R$ 的儲存空間，而 `q_L` 和 `p_R` 則用區域變數 `q1` 和 `p2`)：
1. 計算新的分子：$p_{\text{new}} = p_L \cdot q_R + p_R \cdot q_L$ (對應程式碼 `mpz_mul(p, p, q); mpz_add(p, p, p2)`)
2. 計算新的分母：$q_{\text{new}} = q_L \cdot q_R$ (對應程式碼 `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 記憶體。對一億位數則需約 33 MB。參考程式碼如下:

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <limits.h>
#include <gmp.h>
#include <math.h>
#include <time.h>
#include <stdint.h>

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

/* 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 log factorial using Stirling's approximation */
static long double log_factorial(uint64_t n)
{
    if (n <= 1) return 0;
    return (logl(2.0 * PI) / 2.0 + logl(n) / 2.0 + n * logl(n / E)) / logl(10.0);
}

/* 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 %llu digits\n\n", digits);
    fprintf(fp, "e = ");

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

    /* Write first line with integer part and decimal point */
    fprintf(fp, "%c.", str[0]);
    for (uint64_t i = 1; i <= digits; i++) {
        fprintf(fp, "%c", str[i]);
        if (i % 10 == 0 && i != digits) {
            fprintf(fp, " ");
        }
        if (i % 50 == 0) {
            fprintf(fp, "   [%6llu]\n", i);
            if (i < digits) {
                fprintf(fp, "      ");
            }
        }
    }
    if (digits % 50 != 0) {
        fprintf(fp, "\n");
    }

    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");
}

/* Recursive calculation with memory optimization */
static void calc(mpz_t a, mpz_t b, mpz_t p, mpz_t q)
{
    mpz_t m, p2, q1;
    mpz_inits(m, p2, q1, NULL);

    mpz_sub(m, b, a);
    if (mpz_cmp_ui(m, 0) <= 0) {
        mpz_set_ui(p, mpz_cmp_ui(a, 0) == 0 ? 1 : 0);
        mpz_set_ui(q, 1);
    } else if (mpz_cmp_ui(m, 1) == 0) {
        if (mpz_cmp_ui(a, 0) == 0) {
            mpz_set_ui(p, 2);
            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 %llu digits\n", digits);
    printf("Terms: %llu, Precision: %llu bits\n", num_terms, precision);

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

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

    mpz_set_ui(z, 0);
    mpz_set_ui(n, num_terms);
    calc(z, 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(z, 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 %llu\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/)

## 快速計算指數函數
標準函式庫中的 [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) 方法，透過梯度下降 (gradient descent) 來擬合多項式係數。
> 利用三次多項式進行曲線擬合，參見 [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;

經由上述程式碼可獲得多項式的係數。其產生的多項式在區間 $[0.5, 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) 的浮點數，能表示的最大指數為 $710$。因此，我們可建立一個範圍為 $[-710, 710]$ 的表格，以提供快速精確的 $\exp(i)$ 計算結果：

#define TABLE_MIN -710
#define TABLE_MAX 710
#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;

    /* Minimax 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));

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

    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)

## 深受 $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) 的資訊量
- 支援三進位指令集與可變運算長度 (例如 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$，而是將可能的狀態擴充為 `0`, `1`, `2`，在 [balanced ternary](https://en.wikipedia.org/wiki/Balanced_ternary) (平衡三進位) 中，就是 `-1`, `0`, `+1` 等三個可能狀態，又可簡寫為 `-`, `0`, `+`。
這裡的 "ternary"指「三元的」或「以三為基底的」(base-3)。標準三進位使用 $({0, 1, 2})$，而平衡三進位則使用對稱的 ${-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 it's own value and each successive symbol has it's value multiplied by the base, B, raised to the power of it's 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...$。這意味著從理論上講，以 $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$。這表示二進位取負運算平均會多出接近一次完整的加法器傳播延遲。因此，在轉換效率上，平衡三進位執行負數操作本質上更快。這個看似微小的延遲差異，在深度學習、大型語言模型訓練與推理等等需要巨量運算的場景中，可望累積顯著的性能提升。

### 現代應用：針對大型語言模型的壓縮和節能設計
平衡三進位的思想在現代人工智慧領域，特別是大型語言模型 (LLM) 的技術演化上，重新獲得關注。2023 年 10 月，Microsoft 研究院發表〈[BitNet: Scaling 1-bit Transformers for Large Language Models](https://arxiv.org/abs/2310.11453)〉，初步探索極低位元量化的可行性。隨後在 2024 年 2 月，該團隊發表影響更廣的〈[The Era of 1-bit LLMs: All Large Language Models are in 1.58 Bits](https://arxiv.org/abs/2402.17764)〉，正式提出 BitNet b1.58 方法，並釋出對應的開放原始碼專案 [BitNet](https://github.com/microsoft/BitNet)。此方法正是採用平衡三進位 $({ -1, 0, +1 })$ 來表示大語言模型的權重 (weight)。這種方法之所以被稱為「1.58-bit 量化」，是因為表示三種狀態所需的最少資訊量為 $\log_2(3) \approx 1.58$ 位元，巧妙地介於純粹的 1-bit $({-1, +1})$ 與 2-bit 之間。

BitNet b1.58 的特點及其體現的平衡三進位優勢，首先在於其三元權重 (Ternary Weights) 的量化策略。此方法將傳統以浮點數 (如 FP16 或 [BFloat16](https://en.wikipedia.org/wiki/Bfloat16_floating-point_format)) 表示的模型權重，透過縮放與閾值判斷，映射到 $({-1, 0, +1})$ 這三個離散值之一，並透過 Straight-Through Estimator (STE) 手法來訓練非可微的量化函數，訓練過程中會保留高精度「主權重」用以累積梯度，而前向與反向傳播則使用低位元權重。相較於純粹的 1-bit 量化僅使用 $({-1, +1})$，BitNet b1.58 最關鍵的優勢在於保留零值 (preserving zeros and sparsity)。由於大型語言模型訓練後，常含有大量接近零的權重，強制將其映射至 `-1` 或 `+1` 會引入顯著誤差，並破壞模型內在的稀疏性；而平衡三進位允許這些權重被精確量化為 $0$，從而有效減少資訊損失，維持模型性能。

這種三元量化直接帶來顯著的計算效率提升和能耗降低。在 LLM 中佔據主導地位的矩陣乘法 ($Y = WX$)，當權重 $W$ 僅包含 $({-1, 0, +1})$ 時，其運算過程得以極大簡化：原本的浮點乘法運算被完全取代，權重為 $0$ 的部分直接忽略，權重為 $+1$ 和 $-1$ 的部分則分別轉化為高效的整數加法和減法。這種運算邏輯的簡化顯著降低計算延遲，同時，由於權重儲存從 16 位元大幅壓縮至約 1.58 位元，讀取權重所需的記憶體頻寬壓力也大為減輕，進而降低整體能耗，這對於記憶體受限的推理場景及大規模部署尤其重要。

該手法與 2016 年的 [BinaryConnect](https://arxiv.org/abs/1511.00363) 和 [IBM TrueNorth](https://research.ibm.com/publications/truenorth-design-and-tool-flow-of-a-65-mw-1-million-neuron-programmable-neurosynaptic-chip) 等低精度類神經網路概念相呼應。這些研究早已證明：即使將權重量化為 1 位元，神經網路仍能正常訓練與推論。不過 BitNet b1.58 採用的平衡三進位在訊號能量與資訊表示上找到理想平衡點。此外，雖然是高度量化，$({-1, 0, +1})$ 這三個值具備潛在的表達能力 (potential expressiveness)，可直觀對應神經連結的抑制 ($-1$)、無作用 ($0$) 與激勵 ($+1$) 三種狀態，其熵值為 $\log_2(3) \approx 1.58$，比 2 進位的 1-bit 多，但遠低於 FP16 或 8-bit。

從神經科學的視角來看，平衡三進位的 $({-1, 0, +1})$ 權重也有其生物直觀性：例如 `-1` 可視為抑制性神經元 (inhibitory neuron) 的輸出訊號，`+1` 則代表興奮性輸出，而 `0` 則表示未觸發。實際神經元的行為當然更為複雜，它們不依賴時鐘同步運作，並可整合訊號強度與時間累積效應 (如多個突觸電位的整合) 。脈衝的頻率與時間分佈也能影響神經訊號強度，這與我們從視覺感知中所見到的亮度與色彩變化不無關聯。儘管目前的類神經網路模型尚無法完全模擬這些生物特性，然而像 Intel 的 [Loihi 2](https://www.intel.com/content/www/us/en/research/neuromorphic-computing.html) 等神經形態晶片 (neuromorphic chips) 正逐步接近這種非同步、類比式的訊號處理方式。

根據 BitNet b1.58 的研究論文，其效能表現相當亮眼：在 30 億參數規模以上時，其模型性能已能接近 FP16 的基準模型；而相較於同等規模的 Llama 模型，700 億參數的模型能達成約 4 倍的推理速度提升和 7 倍的記憶體節省，並支援更大的批次處理量以提高吞吐率。

如此簡化的運算 (僅需加/減/忽略) 也開啟硬體實作的潛力，非常適合設計專用的硬體加速器 (ASIC 或 FPGA) ，猶如前蘇聯科學家在「Сетунь」電子計算機所投入的三進位硬體探索。

這些研究成果表明，源於對數值系統效率 ($e$) 思考且具備優雅對稱性的平衡三進位，在現代人工智慧的脈絡下，不僅是理論上的最佳整數基底探討，包含對 $0$ 的處理和計算簡化在內的特性，已將平衡三進位轉化為極具潛力的實用技術，為大型模型的壓縮、加速及節能提供全新的途徑。

延伸閱讀：
- [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
$$

該恆等式將 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
$$

此為複數微分方程，其解為 $f(t) = e^{it}$

透過泰勒展開或解析法可知：
$$
e^{it} = \cos t + i \sin t
$$

這正是歐拉公式的本體。若令 $t = \pi$，則：
$$
e^{i\pi} = -1 \quad \Rightarrow \quad e^{i\pi} + 1 = 0
$$

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

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

根據[傅立葉分析](https://en.wikipedia.org/wiki/Fourier_analysis)，任何「行為良好」的週期函數 $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{(傅立葉轉換)}
$$

由此可見，複指數函數 $e^{i\omega t}$ (結合 $e$, $i$, 以及與週期 $\pi$ 相關的 $\omega$) 不僅是數學分析中的強大工具，更是描述波動、振動、訊號傳播等眾多涉及週期性或頻率特性的自然現象之基本數學語言。

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


在三維歐幾里得空間中，一組標準正交基底為 $\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$ 函數解析延拓至整個複平面，將質數分布與 $\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 - a)^2}{2\sigma^2}}
$$

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

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

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

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

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

其中 $r = \sqrt{x^2 + y^2}$ 為點到原點的距離，$c$ 為常數。

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

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

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

### [泊松分布](https://en.wikipedia.org/wiki/Poisson_distribution)
描述在固定時間或空間範圍內，隨機事件發生的次數。其機率密度函數為：
$$
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}
$$

兩邊同除以 $\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)。

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

- [ ] 資產價格跳變模型

在金融建模中，實際資產價格常會表現出「非連續、突變」的特徵，違反 [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/Hazard_ratio) (hazard ratio)，則某公司在 $[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$ 的深刻連結，屢屢展現本文強調的觀念：離散事件的累積機率，往往可由連續變化的函數來精確描述。

## 橫跨連續現實與離散計算的橋樑
我們身處的世界，本質是連續的。時間猶如河流般流淌不息，汽車的速度、放射性元素的衰變速率、電容器兩端的電壓等現象，都在時時刻刻發生著平滑而連續的變化。在數學上，我們習慣使用連續變數的函數及微積分來描述這一切，其中，微分運算子 $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$ 是固定的採樣時間間隔。

「羅塞達石碑」源於其在破解古埃及象形文字上的關鍵作用。石碑上刻有同一段內容的三種不同文字 (古埃及象形文、埃及草書、古希臘文) ，使得學者能夠相互對照，最終解開失傳語言之謎，開啟了現代埃及學的大門。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
$$

[亞當斯預測法](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$ 孕育出數值微分、外插與積分的實用演算法，作為數位訊號處理與控制系統的基石。該「石碑」讓數位世界與現實和諧共舞。

延伸閱讀:
* [Two Different Worlds](https://www.embedded.com/two-different-worlds/)
* [More on the Rosetta Stone](https://www.embedded.com/more-on-the-rosetta-stone/)
* [Taking the last lap around the Rosetta Stone](https://www.embedded.com/taking-the-last-lap-around-the-rosetta-stone/)

## 後記
整整二十年前，筆者因拜讀 Jack Crenshaw 博士於 [embedded.com](https://www.embedded.com/) 發表的一系列文章而深受啟發，複習歐拉數 $e$ 的極限定義 $\lim_{n\to\infty}(1+1/n)^n$、其獨特的導數性質 $(d/dx\,e^x = e^x)$、優雅的級數表達 $\sum 1/n!$ 乃至於連結數學基石的歐拉公式 $(e^{i\pi}+1=0)$ 等抽象概念，如何在工程領域 (特別是數值計算、控制與訊號處理) 中展現出如此具體的威力與結構上的必然性。這促使筆者在後續的開發與研究中，持續關注 $e$ 作為描述連續增長、衰減等基礎現象的核心數學語言。本文彙整對 $e$ 的多方面探索，涵蓋從 Napier 與 Bürgi 的早期對數工作、歐拉的系統性確立，到其在微積分中的基礎地位、在機率統計 (常態與泊松分佈) 中的關鍵作用，以及在現代計算機上高效近似的方法。期望本文能與讀者一同領略 $e$ 這個常數的深刻內涵，並理解它如何作為自然規律、數學理論與工程應用之間一座重要的橋樑。