owned this note
owned this note
Published
Linked with GitHub
# Fixed Point Arithmetic
為了避免使用浮點運算單元 (FPU)，我們可以在 Monte-Carlo Tree Search 中以定點數來替代浮點數運算。以 32 位元的整數為例，其表示方式可以寫成:
$$X = \sum_{i = 0}^{31} q_i \cdot 2^{i} , q_i \in \{0, 1\}$$
其中，$q_i$ 為整數 $X$ 在第 $i$ 位的二進制位元。例如，十進制數 $9_{10}$ 的二進制表示為:
$$9_{10} = 1001_{2} = 1 \cdot 2^3 + 0 \cdot 2^2 + 0 \cdot 2^1 + 1 \cdot 2^0$$
同樣地，定點數也可以用二進制方式表示，若是我們將整數的二進制表示向右位移16位元我們可以將值域在 $[0, 2^{32} - 1]$ 的無號整數表示位移至無號定點數二進制表示 $[2^{-16}, 2^{16} - 2^{-16}]$ 的值域，整數表示裡 $i$ 的值域為 $[0,31]$ ，在定點數表示中，前面例子指數 $i$ 的範圍為 $[-16, 15]$。 這相當於我們將整數表示的每個位元都乘上一個位移量 $2^{-16}$，因此我們可以將前面的整數表式轉換成定點數表示:
\begin{aligned}
2^{-16} \times X &= 2^{-16} \times \sum_{i = 0}^{31} q_i \cdot 2^{i} , q_i \in \{0, 1\} \\
&= \sum_{i = 0}^{31} q_i \cdot 2^{i - 16} , q_i \in \{0, 1\}
\\ \\
2^{-16} \times 9_{10} &= 0.0001373291015625\\
& = 1 \cdot 2^{-13} + 0 \cdot 2^{-14} + 0 \cdot 2^{-15} + 1 \cdot 2^{-16} \\
&= 1001_{2} \times 2^{-16}
\end{aligned}
我們可以使用 $UQm.n$ ([Q notation](https://en.wikipedia.org/wiki/Q_(number_format))) 表示無號定點數的格式，其中 $m$ 為整數部分的位元數，$n$ 為小數部分的位元數，總共使用的位元數為 $m + n$（基於 Texas Instruments 的定義）。若為有號數，需額外增加一位元作為符號位。
例如，$UQ16.16$ 格式的定點數可以表示為：
$$X = \sum_{i = -16}^{15} q_i \cdot 2^{i} , q_i \in \{0, 1\}$$
$1.625_{10}$ 我們可以使用定點數 $UQ1.3$ 表示為:
$$1.625_{10} = 1101 = 1 \cdot 2^0 + 1 \cdot 2^{-1} + 1 \cdot 2^{-2} + 1 \cdot 2^{-3}$$
## 四則運算
### 加減法
在二進制的定點數的加法及減法計算會與十進制整數相同，
* Integer Addtion and Substraction:
\begin{aligned}
9_{10} + 1_{10} &= 1001_{2} + 0001_{2} \\
&= \{1 \cdot 2^3 + 0 \cdot 2^2 + 0 \cdot 2^1 + 1 \cdot 2^0 \}+ \{1 \cdot 2^0 \} \\
&= {1 \cdot 2^3 + 0 \cdot 2^2 + 1 \cdot 2^1 + 0 \cdot 2^0} \\
&= 1010_2 = 10_{10}
\end{aligned}
\begin{aligned}
9_{10} - 1_{10} &= 1001_{2} - 0001_{2} \\
&= \{1 \cdot 2^3 + 0 \cdot 2^2 + 0 \cdot 2^1 + 1 \cdot 2^0 \} - \{1 \cdot 2^0 \} \\
&= {1 \cdot 2^3 + 0 \cdot 2^2 + 0 \cdot 2^1 + 0 \cdot 2^0} \\
&= 1000_2 = 8_{10}
\end{aligned}
* Unsigned Fixed Point Addtion and Substraction:
\begin{aligned}
0.625_{10} + 0.375_{10} &= 0101_{2} + 0011_{2} \\
&= \{0 \cdot 2^0 + 1 \cdot 2^{-1} + 0 \cdot 2^{-2} + 1 \cdot 2^{-3} \} \\
&+ \{0 \cdot 2^0 + 0 \cdot 2^{-1} + 1 \cdot 2^{-2} + 1 \cdot 2^{-3} \} \\
&= {1 \cdot 2^0} \\
&= 1000_2 = 1_{10}
\end{aligned}
\begin{aligned}
0.625_{10} - 0.375_{10} &= 0101_{2} - 0011_{2} \\
&= \{0 \cdot 2^0 + 1 \cdot 2^{-1} + 0 \cdot 2^{-2} + 1 \cdot 2^{-3} \} \\
&-\{0 \cdot 2^0 + 0 \cdot 2^{-1} + 1 \cdot 2^{-2} + 1 \cdot 2^{-3} \} \\
&= {1 \cdot 2^{-2}} \\
&= 0010_2 = 0.25_{10}
\end{aligned}
### 乘法
但是在定點數的乘法會需要額外的位移來幫助我們得到正確的運算，以 $1.375_{10} \times 1.125_{10}$ 為例子，在 $UQ1.7$ 下的二進制運算為:
\begin{aligned}
1.375_{10} \times 1.125_{10}
&= (176_{10} \times 2^{-7}) \times (144_{10} \times 2^{-7})\\
&= 25344_{10} \times 2^{-14}\\
&= \{1 \cdot 2^{7} + 1 \cdot 2^{5} + 1 \cdot 2^{4}\} \times \{1 \cdot 2^{7} + 1 \cdot 2^{4}\} \times 2^{-14}\\
&= \{1 \cdot 2^{14} + 1 \cdot 2^{13} + 1 \cdot 2^{9} + 1 \cdot 2^{8}\} \times 2^{-14}\\
&= 0110\_0011\_0000\_0000_2 \times 2^{-14} \\
\end{aligned}
正確的運算結果應為 $1.546875_{10} = 1100\_0110_{2}$ ，由於電腦運算時仍然是使用整數系統計算，所以上面的結果我們必須向右位移7個位元才能得到正確的數值。在實作定點數的乘法實必需考量到整數溢位(overflow)的可能性，如同上述例子，若是我們在一個8位元的定點數 $UQ1.7$ 下進行 ${1011\_0000}_{2} \times 1001\_0000_{2}$ 運算，因為整數溢位的關係導致最終結果可能為 0 的答案。
### 除法
除法的運算也會需要額外的位移幫助得到正確的結果，以 $3_{10} \div 2_{10}$ 為例子，在 $UQ2.6$ 下的二進制運算為:
\begin{aligned}
3_{10} \div 2_{10}
&= (192_{10} \times 2^{-6}) \div (128_{10} \times 2^{-6})\\
&= 1.5_{10}\\
&= \{[1 \cdot 2^{7} + 1 \cdot 2^{6}] \times 2^{-6}\} \div \{[1 \cdot 2^{7}] \times 2^{-6}\} \\
&= \{{1100\_0000}_{2} \times 2^{-6}\} \div \{1000\_0000_{2} \times 2^{-6}\} \\
&= 0000\_0001_2 (因為整數運算的關係，小數部分被截取)\\
\end{aligned}
正確的運算結果應為 $1.5_{10} = 0110\_0000_{2}$，與上面得到得結果大相逕庭，這是因為運算結果的值域為 $[0, 2^8 - 1]$ ，所以除法計算中指數小於 0 的部分都會被截取，如同上述的運算應為 $2^0 + 2^{-1}$，但是小數部分被截取得到只有整數的 $2^0$。為了避免小數部分因為整數運算關係被截取，可以先將被除數透過位移的方式將其提升至更高的值域，最後得到結果可以落在我們所期望的範圍:
\begin{aligned}
3_{10} \div 2_{10}
&= (192_{10} \times 2^{-6}) \div (128_{10} \times 2^{-6})\\
&= 1.5_{10}\\
&= \{[1 \cdot 2^{7} + 1 \cdot 2^{6}] \times 2^{-6}\} \div \{[1 \cdot 2^{7}] \times 2^{-6}\} \\
&= \{[(1 \cdot 2^{7} + 1 \cdot 2^{6}) \times {\color{red}{2^{6}}}] \div [1 \cdot 2^{7}]\} \times ({\color{red}{2^{-6}}})\\
&= \{2^{6} + 2^{5}\} \times 2^{-6} \\
&= 0110\_000_2 \times 2^{-6}\\
&= 1.5_{10}
\end{aligned}
這種方式會需要額外的位元儲存被除數，以上述為例在 $UQ2.6$ 下我們會需要 2(整數)+6(小數) = 8 個位元表示一個定點數，但是在除法的計算過程中如果把被除數向左位移6個位元的話我們會需要總共 14 位元來幫助我們運算，雖然我們可以透過對除數使用向右位移來減少運算過程中額外位元的使用(如上述運算，我們可以將除數向右位移 7 位元，此時的運算不需要額外的位元)，但是只靠單純的整數除法仍然無法完全避免額外的位元使用，想像一下在 $UQ2.6$ 下若除數為 $2.015625_{10}$ = $1000\_0001_{2}$ ，可以發現沒有空間可以向右位移，必須利用額外的位元幫助運算。
## Pseduo Division
[Pseudo Division](https://archived.hpcalc.org/laporte/Meggitt_62.pdf) 是一種模擬除法逐位元計算的演算法，IBM 研究員 J.E. Meggitt 改良了一個中古世紀時人們計算對數的演算法 [Radix Method](https://archived.hpcalc.org/laporte/The%20Radix%20Method.htm)，模擬除法的運算期間將除數進行調整，可以幫助我們獲得 $\log{x}$, $\sqrt{x}$, $\tan^{-1}{x}$ 的近似值，由於其中的運算只需要加法和位移，不需要使用到浮點數運算，因此適合實作在計算機當中，世界上第一台手持式計算機 [HP-35](https://en.wikipedia.org/wiki/HP-35) 內就是實作這種演算法計算 $\ln{x}$。
>The method used is essentially Brigg’s method modified
這是 Meggitt 在論文中闡述演算法核心概念源自於16、17世紀的數學家 [Henry Briggs](https://en.wikipedia.org/wiki/Henry_Briggs_(mathematician)) 的方法，他改良了原本對數，使用以10為底的方式建立後來人們常用的對數表。由於論文中的這句話導致我一開始先入為主地認為 Meggitt 使用的方法就是 Briggs 建立對數表的使用方法，因此花了兩三天研究 [《Arithmetica Logarithmica》](https://www.17centurymaths.com/contents/albriggs.html): [Chapter6](https://www.17centurymaths.com/contents/briggsagain/ALfrevCh6.pdf), [Chapter7](https://www.17centurymaths.com/contents/briggsagain/ALfrevCh7.pdf), [Chapter8](https://www.17centurymaths.com/contents/briggsagain/ALfrevCh8.pdf)，後來發現其實 Meggitt 所描述的方法是在 [Chapter14](https://www.17centurymaths.com/contents/briggsagain/ALfrevCh14.pdf) 介紹的方法，不過因為這個意外我也學習了很多以前如何計算根號和對數的方式，其實相當有趣，關於這部分我放到 [ToDo]。 Briggs 在他的書中並沒有為他的方法取名字，但是根據推測第一個提出這個概念的作者為 [William Oughtred](https://archived.hpcalc.org/laporte/The%20Radix%20Method.htm)，其想法其實很簡單:
$$\log{(a\times b)} = \log{(a)} + \log{(b)}$$
這個我們從高中就可以學到的對數性質，可以將 N 分解成若干個 $N_0(1 + q_1)(1 + q_2)...$ 的乘積，而 $\log{N} = \log{N_0} + \log{(1 + q_1)} + \log{(1 + q_2)} + ...。$
\begin{aligned}
1716 &= 1000 \times 1.1 \times 1.2 \times 1.3\\
\log{1716} &= 3 + \log{1.1} + \log{1.2} + \log{1.3}
\end{aligned}
這裡我們會使用現代數學符號來表示 Briggs 在書中使用的例子來說明如何找到任一數的對數，接下來的說明引用自[《A reconstruction of the tables of Briggs’ Arithmetica logarithmica》: Section 3.7](https://locomat.loria.fr/briggs1624/briggs1624doc.pdf)，以 $x = 3041.851528656$ 為例子，另外定義 $\nu{(x)}$ 為截取 $x$ 至第一個有效位元(one significant digit)，舉例來說，$\nu{(0.0002800160 . . .)} = 0.0002$。首先，我們先將整數部分提取出來 $x_{0} = 3041$，接著依序建立一個數列:
\begin{aligned}
r_1 &= x - x_0 = 0.851528656,\ q_1 = \nu(\frac{r_1 }{x - r_1}) = \nu(\frac{0.851528656}{3041}) = \nu(0.000280016) = 0.0002 \\
r_2 &= x - x_0(1 + q_1) = 0.243328656 ,\ q_2 = \nu(\frac{r_2 }{x - r_2}) = \nu(\frac{0.243328656}{3,041.6082})= 0.00008 \\
r_3 &= x - x_0(1 + q_1)(1 + q_2) = 0
\end{aligned}
由上面可以得到 $x = (3041)\times(1.0002)\times(1.00008) = 3041.851528656$，因此可以得到其對數為 $\log{x} = \log{3041.851528656} = \log{3041} + \log{1.0002} + \log{1.00008}$，至此我們可以透過事先計算的對數表組合出任意一個數的對數。
Meggitt 在 [Pseudo Division and Pseudo Multiplication Processes](https://archived.hpcalc.org/laporte/Meggitt_62.pdf) 這篇論文中先介紹一般計算除法的演算法，如下所示:
> 論文中所討論的演算法皆為10進制表示，後續會接著討論關於二進制的處理方式
<figure>
<img
src="https://hackmd.io/_uploads/H1GOMBm4ye.png"
>
<figcaption>Figure 1.</figcaption>
</figure>
$A$ 和 $B$ 皆為一個儲存 $n$ 位元 $10$ 進制的暫存器，分別儲存被除數 $y$ 和除數 $x$。首先，將 $y$ 不斷的減去 $x$ 直到 $A$ 為負數，使用一個計數器(counter)儲存減去的次數，使用另一個暫存器 $Q$ 儲存目前計數器的內容，此時所記錄的是由左向右第 $j$ 位的商數。接著將 A 和 Q 中的內容向左位移一位(乘 10)，對 $j$ 加 1，重複前面的過程計算第 $j + 1$ 位的商數。這個過程會重複執行 n 次，此時 $Q$ 的內容為 $\frac{y}{x}$ 的商數。這個方法暫存器最多會只需要 $n + 1$ 位元即可完成運算，與前一小節所提的方法相比使用更少的位元，但是要注意的是可能會需要使用迴圈的方式實作，迴圈的成本也需要加入實作考量。
<figure>
<img
src="https://hackmd.io/_uploads/r1hSMrmNkx.png"
>
<figcaption>Figure 2.</figcaption>
</figure>
### Logarithm
為了計算對數，透過修改前面所提的 radix method:
$$ y + x = x \prod_{j=0}(1 + 10^{-j})^{q_j}$$
對兩邊取 $log$ 後:
$$\log{(1 + \frac{y}{x})} = \sum_{j}{q_j\log{(1+10^{-j})}}$$
這裡的 $q_j$ 表示在第 $j$ 位元的表示，因此 $0\le q_j < 10$，計算 $q_j$ 的過程會類似前面所介紹的除法，所以才會稱作 pseudo division。若是已知 $q_0 ... q_{j - 1}$，那我們會希望使下面的式子越小越好:
$$ y - x(\prod_{k=0}^{j - 1}(1 + 10^{-k})^{q_k} - 1) $$
接著要找下一位數 $q_j$，可以得到等式:
\begin{aligned}
y_a^{\ (j)} &= y - x\{\ [ \ \prod_{k=0}^{j - 1}(1 + 10^{-k})^{q_k} \ ](1 + 10^{-j})^a - 1 \ \} \\
\text{for}\ a &= 0, 1, 2 ... q_{j} \ ,\ \text{then} \ y_{q_j}^{(j)} \ge 0 >y_{q_{j + 1}}^{(j)}
\end{aligned}
利用上述的不等式可以得知，我們要尋找的 $q_j$ 會使得 $y_a^{\ (j)}$ 盡量靠近 $0$ 但是仍然為正數。接著使用 $x_a^{{j}}$ 簡化等式，可以得到遞迴式:
\begin{aligned}
x_a^{{j}} &=x[\ \prod_{k = 0}^{j - 1}(1 + 10^{-k})^{q_k} \ ](1 + 10^{-j})^a \\
y_{a + 1}^{{(j + 1)}} &= y_{a}^{{(j)}} - 10^{-j}x_a^{(j)}\\
x_{a + 1}^{{(j + 1)}} &= x_{a}^{(j)} + 10^{-j}x_a^{-j}
\end{aligned}
利用這個遞迴式可以不斷的逼近目標 $y$，可以觀察這個遞迴式，如果將 $x_a^{(j)}$ 的係數 $10^{-j}$ 移除，那這個遞迴式就是前面所提的普通除法，再次呼應為何這是 pseudo division。為了保持精準度，接著更進一步地簡化遞迴式:
\begin{aligned}
z_a^{(j)} &= 10^jy_a^{(j)} \\
z_{a + 1}^{(j)} &= z_a^{(j)} - x_a^{(j)} \\
x_{a + 1}^{{(j + 1)}} &= x_{a}^{(j)} + 10^{-j}x_a^{(j)} \\
\text{then,} \ \ z_{q_j}^{(j)} &\ge 0 > z_{q_{j+1}}^{(j)}
\end{aligned}
這個遞迴式與普通除法幾乎相同，只差在每次更新時都會加上 $10^{-j}x_a^{-j}$。
假設已知 $q_j$ ，則計算 $q_{j+1}$ 的初始條件為:
\begin{aligned}
z_0^{(j+1)} &= 10^{j+1}y_0^{(j+1)} = 10^{j+1}y_{q_j}^{(j+1)} = 10z_{q_j}\\
x_0^{(j+1)} &= x_{q_j}^{(j)}
\end{aligned}
接著透過不斷地進行減法運算，以提取 $q_{j+1}$，這個過程就跟除法運算相同。整個 pseudo division 的運算不斷地提取 $q_j$ 直到第 $n - 1$ 個位元。一旦我們得到 $q_0, q_1, q_2 ...$ 後，可以使用事先計算的對數表得知所求的對數值，可以得到等式:
$$\log{(1 + \frac{y}{x})} = \sum_{j = 0}{q_j\log{(1+10^{-j})}}$$
* Register lengths:
$Q$ 會記錄 $n$ 位元的運算結果，而暫存器 $B$ 儲存 $x$ 可以從前面的遞迴式得知在進行加法的過程可能會導致進位，所以會需要 $n+1$ 的位元數，$A$ 則是儲存 $y$ 在進行偽除法的過程中，有可能因為 $q_j = 0$ 而必須多借一位，所以會需要 $n + 2$ 位進行偽除法運算。
#### Binary format:
這裡會介紹一篇由浮點數之父 [William Kahan](https://en.wikipedia.org/wiki/William_Kahan) 撰寫的一篇文章 [Pseudo-Division Algorithms for Floating-Point Logarithms and Exponentials](https://ieeemilestones.ethw.org/w/images/3/30/Wk_pseudo_division_log_exp_may01.pdf) ，主要的內容只有6頁，但是卻可以言簡意賅的闡述從推導到實際演算法的過程。
> What makes the pseudo-division algorithms so attractive is their close resemblance to the simplest binary multiplication and division algorithms
這段話可以得知為何人們一開始會選擇使用 pseudo-division 的原因，因為整個演算法只需要單純的加法和位移操作即可實作，因此文中也有提到像是 8087 數值協同處理器也有使用此演算法。
>Algorithms like these are used by the Intel 8087 family of numeric coprocessors.
雖然這篇文章的內容是應用在浮點數， 但是對於 pseudo division 的核心想法同樣適用在定點數上。
首先定義等式:
$$ y = \frac{1}{(1 + \frac{q_1}{2})(1 + \frac{q_2}{4})(1 + \frac{q_3}{8})(...)(1 + \frac{q_k}{2^k})(1 + ...)} \ , \ q_j \in \{0, 1\}$$
$y$ 的最小值為:
\begin{aligned}
\frac{1}{\eta(1)} &= 0.41942244\\
\frac{1}{\eta(t)} &=\frac{1}{(1 + \frac{t}{2})(1 + \frac{t}{4})(1 + \frac{t}{8})(...)(1 + \frac{t}{2^k})(1 + ...)}
\end{aligned}
這個無窮乘積需要使用 [Euler function](https://en.wikipedia.org/wiki/Euler_function) (不是 [Euler's totient function](https://en.wikipedia.org/wiki/Euler%27s_totient_function)) 才能計算[[註](https://math.stackexchange.com/questions/1924882/evaluating-prod-n-1-infty-left1-frac12n-right)]，可以得到 [link](https://www.wolframalpha.com/input?i=product+1%2F%281%2B2%5E%28-j%29%29%2C+j%3D1+to+infinity) 與論文中相同的答案。可以得到 $y$ 會介於 $\frac{1}{\eta{(1)}}$ 和 $1$ 之間。對於同一 $y$ 可以有不唯一的表式方式:
\begin{aligned}
\frac{2}{3} &= \frac{1}{1 + \frac{1}{2}}\\
&= \frac{1}{(1 + 2^{-2})(1 + 2^{-3})(1 + 2^{-4})(1 + 2^{-8})(1 + 2^{-16})(1 + ...)} \\
&= (0.666666666821...) \times (1 + ...) \ , \ \text{calculate with W11 calculator} \\
&= \frac{1}{(1 + 2^{-2})(1 + 2^{-3})(1 + 2^{-4})(1 + 2^{-8})(1 + 2^{-16}) \prod_{k = -33}^{-51}(1 + 2^k)}\\
&= (0.666666666666666962750...) \ , \ \text{calculate with W11 calculator}
\end{aligned}
對於單精度的浮點數至小數點後第7位即足夠($log_{10}{2^{24}}\approx 7.224719$)，雙精度則會是小數點後15位($log_{10}{2^{53}}\approx 15.954590$)，因為這個不唯一的性質使得我們可以利用 pseudo-division 演算法產生 y 的 pseudo-quotient bits $q_k$，對 y 取 $log$ 後可以得到:
$$\log{y} = -q_1\log{\frac{3}{2}} - q_2\log{\frac{5}{4}} - ... - q_k\log{(1 + 2^{-k})}-...$$
對於一個浮點數 Y，想計算其 $\log{(Y)}$，可以利用提取 $2^n$ 將其值域轉換至 $y$ 的區間:
$$\frac{1}{\eta{(1)}} < \frac{1}{2} \le y = \frac{Y}{2^n} \le 1$$
可以知道"正" [單精度的浮點數](https://en.wikipedia.org/wiki/Single-precision_floating-point_format) 表示為:
\begin{aligned}
\text{Normalized:} \ \ Y &= (2^0 + \sum_{i=-1}^{-23}q_i2^i) \times 2^{e-127} \ ,\ 2^{-126}\le Y\le (2 - 2^{-23})\times 2^{127}\\
\text{Denormalized:} \ \ Y &= (\sum_{i=-1}^{-23}q_i2^i) \times 2^{-126} \ , \ 2^{-149} \le Y \le (1 - 2^{-23})\times 2^{-126}\\
\end{aligned}
只考慮正數部分，對於 Normorlized 的部分可以將 $1$ 化為 $2^0$ ，可以得到:
$$\frac{1}{2} \le y = \frac{Y}{2^{n}} = \sum_{i=1}^{24}q_i2^{-i} \le 1\ ,\ \text{for} \ \ n\in[-126, 128]$$
對於 Denormorlized 的部分同樣可以得到:
$$\frac{1}{2} \le y = \frac{Y}{2^{n}} = \sum_{i=1}^{23}q_i2^{-i} \le 1\ ,\ \text{for} \ \ n\in[-149, -126]$$
接著可以利用 pseudo-division 得到 $\log{(Y)} = n + \log{(y)}$ 。 $\log{(y)}$ 可以利用對數表事先計算 $\log{(1+2^{-k})}$ 的對應值，透過查表的方式進行計算。然而對數表能夠表示是有限的，小於對數表的位元必須要做出補償。假設我們要對小於 $\log{(1+2^{-L})}$ 的部分進行補償，則等式可以寫成如下:
\begin{aligned}
y_L&=y(1+2^{-1})(1+2^{-2})(1+2^{-3})(...)(1+2^{-L}) \\
&=\frac{1}{(1+\frac{q_{L+1}}{2^{L+1}})(1+\frac{q_{L+2}}{2^{L+2}})(1+\frac{q_{L+3}}{2^{L+3}})(...)} \\
\end{aligned}
可以知道 $\frac{1}{\eta{(2^{-L})}} \le y_L \le 1$，而當 $L$ 夠大時 $y_L$ 會非常接近 1， $\log{(y_L)}$ 可以同樣地使用多項式逼近:
$$ \log{Y} = n - \sum_{1}^{L}q_k\log{(1+2^{-k})} + \log{(y_L)}$$
同時可以得到不等式:
$$0\le y_L = y - \sum_{1}^{L}q_k\log{(1+2^{-k})} = \sum_{L+1}^{\infty}q_k\log{(1+2^{-k})} < \log{(\eta{(2^{-L})})}$$
當 $L$ 夠大時，此時 $y_L$ 會非常的小。
#### Cancellation
主要有兩種捨入誤差的會使得精準度因此下降，分別是災難性抵銷(cancellation)和累積捨入誤差(accumulation)。
考慮 $Y$ 在 $\log{(Y)}$ 時僅超出 $1$ 一點點，也就是 $Y = 1 + \xi$ 當 $\xi$ 非常小時， $\log{(Y)} = \frac{\xi}{\ln{2}}$ 會非常接近。當 $n=1$ 時，對於 $y = \frac{Y}{2^n}$ 會使 $\log{(y)} =-(1 - \frac{\xi}{\ln{(2)}})$ 同樣的接近。
>But if $\xi$ is tiny enough, a few orders of magnitude bigger than rounding errors in $\log{(y)}$, then $log{(y)}$ will round to $-(1-\frac{\eta}{\ln{(2)}})$ instead, where $\eta$ matches in at most its first few sig. bits even though those sig. bits lie far to the right of the binary point, and then the value computed for $\log{Y} = 1 + \log{(y)}$ will be $\frac{\eta}{\ln{(2)}}$ instead of $\frac{\xi}{\ln{(2)}}$ differing in all but the first few sig. bits.
這段話主要是在闡述當 $\xi$ 夠小時 $\log{(y)}$ 會非常接近 $-(1-\frac{\eta}{\ln(2)})$ ，所以可以使用 $\frac{\eta}{\ln{(2)}}$ 取代 $\frac{\xi}{\ln{(2)}}$ 計算 $\log{Y} = 1 -(1-\frac{\eta}{\ln(2)})$。以10進制捨入至第五位元為例子，假設 $\log{(Y)} = 0.00014426$ ;而 $\log{(y)} = -0.99985574$ 會捨入至 $-0.99986$，因此會得到計算 $\log{(Y)} = 1 - 0.99986 = 0.00014000$，可以發現雖然絕對誤差(absolute error)非常的微小，但是相對誤差(relative error)卻省略了一半的有效位元(significant bits)，若是計算像是 $\frac{\log{(Y)}}{(Y - 1)}$ 等運算可能會出現意想不到的數值。
作者因此提出限制 $y = \frac{Y}{2^n}$ 範圍，使其接近 1 ，避免前面所述會可能在 $\log{(Y)}$ 時僅超出 $1$ 外一點點會造成的災難性抵銷。這裡作者提出兩種解決方案。
1. 將 $y = \frac{Y}{2^n}$ 一開始就限制在 $\frac{1}{\sqrt{2}} < y < \sqrt{2}$ 。如果此時 $y < 1$ 那就沒有任何問題，可以繼續進行後續操作;否則用 $y' = \frac{1}{y}$ 取代 $y$ ，接著計算 $y' = 1 - \frac{(y-1)}{y} \rightarrow 1 - y' =\frac{(y-1)}{y}$，所以可以使用 $-\log{(y')}$ 取代 $\log{(y)}$，可以得到當 $y > 1$ 時 $1 < y' = \frac{1}{y} < \frac{1}{\sqrt{2}}$，所以當 $Y$ 非常靠近 1 、而 $n=0$ 時不會再有嚴重的抵銷發生。
> The Intel 8087 family of numeric coprocessors does this to keep the error in their log functions below two units in the last (64th) bit.
Intel 8087 的數值協同處理器內就是採取這種方式使 $log$ 函式內的誤差保持在最後一位元(64th)。但也代表需要額外的成本處理 $1-y'$。
2. 第二種方法透過事先建立的對數表:
$$-\log{(1-2^{-k})}$$
這種方式因為不需要除法運算所以會比前一種還要的快。首先利用 $y=\frac{Y}{2^n}$ 得到 $\frac{2}{3} < y < \frac{4}{3}$ ，接著計算 $y-1$，可以得到 $M$:
$$-2^{-M} < y-1 \le -2^{-M-1} \text{ or } 2^{-M-1} < y-1 \le 2^{-M}$$
假設 $M = 2$。如果 $y < 1$ 維持原本的偽除法;否則當 $y>1$ 時，首先用 $y'$ 取代 $y-1$ :
\begin{aligned}
&\frac{1}{8} < y - 1 < \frac{1}{4}\\
&\frac{9}{8} < y < \frac{5}{4} \\
&y'-1 = (1 - 2^{-M})y-1 =(y-1)-2^{-M}y \\
&y'-1 = (\frac{3}{4})y-1
\end{aligned}
可以發現 $y'$ 會因此介於 $(\frac{27}{32},1)$ 之間，正好同時落在 $\frac{1}{\eta{(\frac{1}{4})}} = 0.786417$ 和 1 這個區間，這也是當偽除法從 $q_3$ 開始計算的區間，表示前面的 pseudo qoutient bit 為 0 所以可以略過不需計算，最後計算:
$$\log{(y)} = -\log{(1-2^{-M})} + \log{(y')}$$
這部分的敘述強烈建議閱讀論文中的第三小節。
#### Fixed point
可以發現不管是浮點數或是定點數，都可以透過提取出 $2^{n}$ 將目標 $y$ 轉換至 $(\frac{1}{\eta{(1)}}, 1]$ 的值域當中。除了前面提到的浮點數，若是定點數 $Z$ 使用 $M$ 個位元表示小數如下:
\begin{aligned}
Z &= \sum_{i = 0}^{31} q_i 2^{i-M} \ , \ q_i \in \{0, 1\} \\
\frac{1}{\eta{(1)}} &< \frac{1}{2} \le y = \frac{Z}{2^n} \le 1
\end{aligned}
仍然可以提取出 $2^n$ 使得 $y$ 落在 $[\frac{1}{\eta{(1)}}, \ 1]$。因此對於浮點數和定點數可得:
\begin{aligned}
y &= \frac{1}{(1 + \frac{q_1}{2})(1 + \frac{q_2}{4})(1 + \frac{q_3}{8})(...)(1 + \frac{q_k}{2^k})(1 + ...)} \ \\ \\
y_k&=y(1+\frac{q_1}{2})(1+\frac{q_2}{4})(1+\frac{q_3}{8})(1+...)(1+\frac{q_k}{2^k}) \\
&= \frac{1}{(1+\frac{q_{k+1}}{2^{k+1}})(1+\frac{q_{k+2}}{2^{k+2}})(1+\frac{q_{k+3}}{2^{k+3}})(1+...)} \ge \frac{1}{\eta(2^{-k})} \ \\
&, \ \text{for} \ k = 1,2,3,4...,L \ \ \text{in turn}
\end{aligned}
對於每個 pseudo-quotient bit $q_k \in \{0,1\}$ ，可以透過下面的不等式來決定:
\begin{aligned}
\frac{1}{(1+2^{-k})} < y_{k-1} \le 1 \ &\text{then} \ q_k = 0 \ \text{and} \ y_k = y_{k-1} \\
&\text{so} \ \frac{1}{\eta{(2^{-k})}} < \frac{1}{1+2^{-k}} < y_k \le 1; \\
\frac{1}{\eta{(2^{1-k})}} < y_{k-1} < \frac{1}{\eta{(2^{-k})}} \ &\text{then} \ q_k = 1 \ \text{and} \ y_k = y_{k-1}(1+2^{-k}) \\
&\text{so} \ \frac{1}{\eta{(2^{-k})}} < y_{k} < 1.
\end{aligned}
這個不等式可以理解為，為了找出符合所求 $y$ 的 $\frac{q_k}{1+2^{(-k)}}$ 乘積組合，對於 $k=1,2,3...,L$ 利用 $y_{k-1}(1+2^{(-k)}) \le 1$ 找出所求的乘積組合，所以當 $y_{k-1}(1+2^{(-k)}) > 1$ 時則超出我們所求的範圍，且 $y_{k-1}$ 同時也必須落在一開始所限制 $[\frac{1}{\eta{(t)}}, 1]$ 區間中，所以才有了第一個不等式 $\frac{1}{(1+2^{-k})} < y_{k-1} \le 1$。
從上面的不等式可以得到下面的遞迴式:
\begin{aligned}
y_{k - 1}(1+2^{-k}) > 1 \ &\text{then} \ q_k = 0 \ and \ y_k = y_{k-1} \\
&\text{else} \ q_k = 1 \ and \ y_k = y_{k-1}(1+2^{-k}).
\end{aligned}
這個遞迴式為整個演算法的核心。但是他仍然有兩個問題，第一個問題是初始化條件沒有還未定義，第二個問題是不斷的累加的誤差
。接著引入下面縮放變數 $z_k$:
\begin{aligned}
z_k &= 2^k(1 - y_k) \\
y_k &= 1 - 2^{-k}z_k \\
1 - 2^{-k}z_k &= 1 - 2^{-(k-1)}z_{k - 1}(1 + 2^{-k}) \\
z_k &= 2z_{k - 1}+ 2^{1-k}z_{k - 1}
\end{aligned}
引入 $q_k$ 可得:
\begin{aligned}
\ y_k &= y_{k-1}(1+2^{-k}q_k) \ge q_k , \ q_k \in \{0,1\}\\
z_k &= 2z_{k - 1}+ 2^{1-k}q_{k}z_{k - 1} \ge q_k \\
z_k &= 2z_{k - 1}+ 2^{1-k}q_{k}z_{k - 1} - q_k\ge 0 \\
z_k &= 2z_{k - 1} (1 + 2^{-k}q_{k}) - q_k\ge 0
\end{aligned}
可以發現捨入誤差會發生在運算 $2^{(-k)}z_k$ 被移出暫存器時，而後期移出的位元比早期移出的位數更右的位置，因此它們不會累積。
z-遞迴應該如何開始？回想一下，偽除法可行的條件是 $\frac{1}{\eta(2^{-k})} < y_k \le 1$ 對每個 $k$ 成立，如果第一個 $k$ 成立則對每一個 $k$ 都成立。則對於 $z_k$:
\begin{aligned}
&1-\frac{1}{\eta(2^{-k})} >1 - y_k \ge 1-1 \\
&1-\frac{1}{\eta(2^{-k})} >1 - y_k \ge 0 \\
&0 \le 2^k(1 - y_k) < 2^k(1-\frac{1}{\eta(2^{-k})}) \\
&0 \le z_k < 2^k(1-\frac{1}{\eta(2^{-k})}) < 1 \\
&0 \le z_k < \zeta{(2^{-k})} < 1 \ ,\ \zeta{(t)} = \frac{(1-\frac{1}{\eta{(t)}})}{t} \\
\end{aligned}
從前面的章節可以知道如果要避免 cancellation 則 $y$ 會落入下面三個區間中:
$$\frac{1}{\sqrt{2}} < y \le 1 \text{ or } \frac{2}{3} < y \le 1 \text{ or } \frac{27}{32} \le y < 1$$
接著引入 $z = 1 -y$ :
$$0 \le z < \frac{1}{2 + \sqrt{2}} \text{ or } 0 \le z < \frac{1}{3} \text{ or } 0 < y \le \frac{5}{32}$$
可以很輕易地知道從浮點數中可以輕易地提取出其指數部分(對於定點數也可以透過 `clz` 函數得到)為整數 $K$ 會使得 $\frac{1}{4} \le z_K = 2^Kz < \frac{1}{2}$，這就是 $z_k$-recurrence 對於 $z_k$ 的初始值。但是必須先考慮兩個條件;首先，如果 $q_k = 0$ 則不斷地運算 $K$ 至 $K+1$ 直到 $q_k=1$，會得到第一個非零的 pseudo-quotient bit $q_K$，並且保存下來。若是 $z$ 落在前兩個區間中，則 $K$ 會從 $K \ge 1$ 開始，在最後一個區間則是 $K\ge2$。$k=K$ 的第二個考量點是，當 $K\ge L$ 或 $z=0$ 時偽除法的運算都可以省略，當 $z=0$ 時 $\log{(y)} = 0$;否則將 $K$ 加 $1$ 得到 $\frac{1}{2} \le z_K < 1$ 然後使用前面的方式逼近 $\log{(1-2^{-K}z_K)}$。
至此已經介紹完關於定點數計算 $\log$，文章內同時有介紹如何計算 $2^x$，有興趣的讀者可以詳細閱讀原文。
#### Implementation
由前一小節得到的遞迴式:
\begin{aligned}
y_{k - 1}(1+2^{-k}) > 1 \ &\text{then} \ q_k = 0 \ and \ y_k = y_{k-1} \\
&\text{else} \ q_k = 1 \ and \ y_k = y_{k-1}(1+2^{-k}). \ \ \ \ \ \ \ \ \ \text{(1)}
\end{aligned}
可以實作出32位元 `UQ16.16` 的定點數 $log_2$ 演算法，為了將定點數 $X = \sum_{i = 0}^{31} q_i 2^{i-16} \ , \ q_i \in \{0, 1\}$ 轉換至 $y$ 值域 $(\frac{1}{\eta{(1)}}, 1]$ 中，使用 `__builtin_clz` 得到最高項的係數 `expo = 32 - lz`，可以將前面的定點數表示為 $X =2^{expo} \times y$，第 14 行的計算是就是將得到的 $y$ 位移成 `UQ1.16`。第 17 ~ 35 行為上述遞迴式的實作，當等式成立時 pseudo-qoutient bit $q_k = 1$ 為:
\begin{aligned}
\log{(X)} &= \log{(2^{expo} \times y)} = lz + \log{(y)} \\
&=lz -q_1\log{1 + 2^{-1}} - q_2\log{1 + 2^{-2}} - ... - q_k\log{(1 + 2^{-k})}-...
\end{aligned}
下面的實作透過變數 `r` 儲存當 $q_k = 1$ 時對應的 $\log{(1 + 2^{-k})}$ 。以 $\log{(1 + 2^{-1})}$ 為例子:
\begin{aligned}
\log_2{1+2^{-1}} &= 0.584962500721 \\
0.584962500721 \times 2^{16} &= 38336.102447 \\
38336_{10} &= 95C0_{16}
\end{aligned}
必須把事先計算的對數值轉換成定點數的格式，所以此時必須乘 $2^{16}$ 可以得到 `0x95C0`。
由前面可以知道如果要對小於 $\log{(1+2^{-L})}$ 的部分進行補償，可以在結果的最後加上 $\log{(y_L)} = \log{(\eta{(2^{-L})})}$ ， 若是以 $UQ16.16$ 為例子此時的 $L = 16$ ，使用 [link](https://www.wolframalpha.com/input?i=product+1%2F%281%2B2%5E%28-j-16%29%29%2C+j%3D1+to+infinity&lang=zh) 計算得到 $\eta{(2^{-L})} \sim 0.999984741366...$ ，由下面式子計算:
\begin{aligned}
\eta{(2^{-L})} &= 0.999984741366... \\
\log_2{(\eta{(2^{-L})})} &\simeq \log{(0.999984741366)}= -0.00002201372 \\
-0.00002201372 \times 2^{16} &= -1.44269115392
\end{aligned}
因此所得到的補償值會是 $-2 < \log{(y_L)} = \log{(\eta{(2^{-L})})} \simeq -1.44269115392 < -1$ ，對應下面實作最後將所得的 r 減去 2。
最終，將前面得到的 `expo` 轉換成 $UQ16.16$ 計算 `(expo << FSHIFT32) - r` 得到輸出的答案。下面是對應的實作:
```c=
#define FSHIFT32 16
#define FIXED32_1 (1 << FSHIFT32)
uint32_t log2(uint32_t x)
{
if (x == 1)
return 0;
if (x == 1) return 0;
uint32_t z, t, one = FIXED32_1; // UQ1.16
uint32_t r = 0;
uint32_t lz, expo;
/* split: x = a * 2**expo, 0 < a < 1 */
lz = __builtin_clz(x);
expo = 32 - lz;
z = (uint32_t) ((x << lz) >> (32 - FSHIFT32)); // UQ1.16
/* try factors (1+2**(-i)) with i = 1, ..., 16 */
if ((t = z + (z >> 1)) < one) { z = t; r += 0x095c0; }
if ((t = z + (z >> 2)) < one) { z = t; r += 0x0526a; }
if ((t = z + (z >> 3)) < one) { z = t; r += 0x02b80; }
if ((t = z + (z >> 4)) < one) { z = t; r += 0x01664; }
if ((t = z + (z >> 5)) < one) { z = t; r += 0x00b5d; }
if ((t = z + (z >> 6)) < one) { z = t; r += 0x005ba; }
if ((t = z + (z >> 7)) < one) { z = t; r += 0x002e0; }
if ((t = z + (z >> 8)) < one) { z = t; r += 0x00171; }
if ((t = z + (z >> 9)) < one) { z = t; r += 0x000b8; }
if ((t = z + (z >> 10)) < one) { z = t; r += 0x00052; }
if ((t = z + (z >> 11)) < one) { z = t; r += 0x0002e; }
if ((t = z + (z >> 12)) < one) { z = t; r += 0x00017; }
if ((t = z + (z >> 13)) < one) { z = t; r += 0x0000c; }
if ((t = z + (z >> 14)) < one) { z = t; r += 0x00006; }
if ((t = z + (z >> 15)) < one) { z = t; r += 0x00003; }
if ((z + (z >> 16)) < one) { /* z = t;*/ r += 0x00001; }
r-=2;
r = (expo << FSHIFT32) - r;
return r;
}
```
原本的遞迴式可以透過迭代加上查找表的方式實作，但是直接展開可以利用編譯器的最佳化選項改善，另外最後一個判斷式內的 `z = t` 其實是冗餘，展開後可以直接省略。這裡要注意的是這個函式輸入為32位元的無號整數，而輸出則是 $UQ16.16$ 的定點數。
#### Correct cumulative error
接著引入前面所提及的縮放變數 $z_k = 2^k(1 - y_k)$ 和 $q_k$ 可以得到不等式
$$z_k = 2z_{k - 1} (1 + 2^{-k}q_{k}) - q_k\ge 0$$
由這個不等式可以得到遞迴式:
\begin{aligned}
z_k = 2z_{k - 1} (1 + 2^{-k}q_{k}) \ge q_k &\text{ then } q_k = 0 \text{ and } z_k = 2z_{k - 1} = 0\\
& \text{ else } q_k = 1 \text{ and } z_k = 2z_{k - 1} (1 + 2^{-k}) - 1 \\
\end{aligned}
下面第 13~31 行判斷式內為 $z_{k - 1}(1 + 2^{-k})$ 的計算，若是 $q_k = 1$ ，則將 $z_{k - 1}(1 + 2^{-k})$ 減去 1 並且加上對應的 $\log{(1 + 2^{-k})}$。
若是要實作 $UQ16.16$ 定點數的 $\log_2$ 演算法，則 $z$ 的初始值由前面的說明可以知道會落在下面其中一個區間:
$$0 \le z < \frac{1}{2 + \sqrt{2}} \text{ or } 0 \le z < \frac{1}{3} \text{ or } 0 < z \le \frac{5}{32}$$
若是 $z$ 落在前兩個區間則 $k\ge2$ ， 反之則從 $k \ge 1$ 開始，因為具體實作時需要額外增加 `if-else` 判斷式，因此下面實作並未加入此判斷。
* 實作
```c=
#define FSHIFT32 16
#define FIXED32_1 (1 << FSHIFT32)
uint32_t log2(uint32_t x)
{
if (x == 1)
return 0;
const uint32_t one = FIXED32_1 >> 1; // UQ16.16
uint32_t z;
uint32_t r = 0; // UQ16.16
uint32_t lz;
/* split: x = a * 2**expo, 0 < a < 1 */
lz = __builtin_clz(x);
z = ((x << lz) >> (32 - (FSHIFT32 - 1)));
z = one - z;
z = z << 1; if (z + (z >> 1) >= one) { z = z + (z >> 1) - one; r += 0x095c0; }
z = z << 1; if (z + (z >> 2) >= one) { z = z + (z >> 2) - one; r += 0x0526a; }
z = z << 1; if (z + (z >> 3) >= one) { z = z + (z >> 3) - one; r += 0x02b80; }
z = z << 1; if (z + (z >> 4) >= one) { z = z + (z >> 4) - one; r += 0x01664; }
z = z << 1; if (z + (z >> 5) >= one) { z = z + (z >> 5) - one; r += 0x00b5d; }
z = z << 1; if (z + (z >> 6) >= one) { z = z + (z >> 6) - one; r += 0x005ba; }
z = z << 1; if (z + (z >> 7) >= one) { z = z + (z >> 7) - one; r += 0x002e0; }
z = z << 1; if (z + (z >> 8) >= one) { z = z + (z >> 8) - one; r += 0x00171; }
z = z << 1; if (z + (z >> 9) >= one) { z = z + (z >> 9) - one; r += 0x000b8; }
z = z << 1; if (z + (z >> 10) >= one) { z = z + (z >> 10) - one; r += 0x00052; }
z = z << 1; if (z + (z >> 11) >= one) { z = z + (z >> 11) - one; r += 0x0002e; }
z = z << 1; if (z + (z >> 12) >= one) { z = z + (z >> 12) - one; r += 0x00017; }
z = z << 1; if (z + (z >> 13) >= one) { z = z + (z >> 13) - one; r += 0x0000c; }
z = z << 1; if (z + (z >> 14) >= one) { z = z + (z >> 14) - one; r += 0x00006; }
z = z << 1; if (z + (z >> 15) >= one) { z = z + (z >> 15) - one; r += 0x00003; }
z = z << 1; if (z + (z >> 16) >= one) { r += 0x00001; }
r = ((32 - lz) << FSHIFT32) - r;
return r;
}
```
#### Implementation comparision
這裡分別會比較幾個實作，
1. y-recurrence log2
2. compensated log2
3. cumulative error correction log2
1.是一開始 $y$-遞迴式的實作; 2.為經過補償過後的版本; 3.則是利用縮放變數 $z_k$ 減少累積誤差實作。測試的輸入為 $[1,10^7]$，以及所比較的對象(reference)會是單精度版本的 `log2f`。
首先這三個實作的平均及最大誤差為:
| | y-recurrence | compensated | cumulative error correction |
|:----------:|:------------:|:-----------:|:---------------------------:|
| avg. error | 0.000075 | 0.000074 | 0.000078 |
| max. error | 0.000242 | 0.00021172 | 0.000206 |
由前面的表格可以觀察到，額外增加的補償值在平均與最大誤差上與原先的 `y-recurrence` 相比分別改善了 $1.33\%$ 和 $12.51\%$。可是對於減少累績誤差的版本在平均誤差卻沒有改善，而最大誤差則比原先好 $14.88\%$。
下圖為在輸入為 $[1,10^7]$ 下與 `log2f` 的誤差值。
下面三張圖是使用 `__rdtscp` 後每個方法跟 `log2f` 的 `cycles` 數比較。
很可惜的是三種方法都沒有比原本的 `log2f` 還要來的好，由原本的實作可以觀察到， `digit by digit` 的方法會有一個根本上的問題，那就是其精確度會跟 `branch` 的數量成正比關係，因此若是有速度上的需可以透過減少 `branch` 數降低精確度來得到更好的計算速度。
### Squre Root
由 [從 √2 的存在談開平方根的快速運算](https://hackmd.io/@sysprog/sqrt#Digit-by-digit-Method) 可以知道對於開根號有非常多種的方式可以實作，這裡會接著介紹如合從 Meggitt 的 10 進制 pseudo division 計算開根號並推廣至 2 進制的實作，核心想法與 [Digit-by-digit calculation](https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Digit-by-digit_calculation) 類似。
首先定義等式:
\begin{aligned}
y &= x[\sum_{j=0}q_j10^{-j}]^2 \\
\sqrt{\frac{y}{x}} &= \sum_{j=0}q_j10^{-j}
\end{aligned}
假設對於 pseudo qoutient bit $q_0, ... q_{j-1}$ 都已經確定，如果要確定 $q_j$ 可以有下面的表達式:
$$y_{a}^{(j)} = y - x[\sum_{k=0}^{j-1}q_{k}10^{-k} + a10^{-j}]^2\ , \ \text{for } a=0,1,2,...$$
此等式的目標是要盡可能的使 $y_a^{(j)}$ 靠近 0 並保持為正數。將平方展開後會有下面的規律:
\begin{aligned}
(\sum_{i=1}^{n}d_i)^2 &=(d_1 + d_2 + d_3 + ... + d_n)^2 \\
&= d_1^2 + 2d_1d_2 + d_2^2 + 2d_1d_3 + d_3^2 + ... + 2\sum_{i=1}^{n-1}(d_i)d_n + d_n^2 \\
&= d_1^2 + (2d_1 + d_2)d_2 + [2(d_1 + d_2)+d_3]d_3 + ... [2(\sum_{i=1}^{n-1}d_i) + d_n]d_n
\end{aligned}
利用上面的性質可以得到下面:
\begin{aligned}
y &= x[\sum_{j=0}q_j10^{-j}]^2 \\
&= x\{q_0^210^{2(0)} + 2q_0q_110^{0+(-1)} + q_1^210^{-2} + [2(q_010^0 + q_110^{-1})+q_210^{-2}]q_210^{-2} + \\
& ... + [2\sum_{i=0}^{j-1}q_i10^{-i}+q_j10^{-j}]q_j10^{-j} \} \\
&= 2x\{\frac{1}{2}q_0^210^{2(0)} + q_0q_110^{0+(-1)} + \frac{1}{2}q_1^210^{-2} + [(q_010^0 + q_110^{-1})+\frac{1}{2}q_210^{-2}]q_210^{-2} + \\
& ... + [\sum_{i=0}^{j-1}q_i10^{-i}+\frac{1}{2}q_j10^{-j}]q_j10^{-j} \} \\
\end{aligned}
由上面的展開式可以知道對於 $q_j10^{-j}$ 的展開部份為 $2x\{[\sum_{i=0}^{j-1}q_i10^{-i}+\frac{1}{2}q_j10^{-j}]q_j10^{-j}\}$，這部份可以改寫:
\begin{aligned}
2x10^{-j}[\sum_{i=0}^{j-1}q_iq_j10^{-i}+\frac{1}{2}q_j^210^{-j}] = 2x10^{-j}[\sum_{i=0}^{j-1}q_i'10^{-i}+\frac{1}{2}q_j'10^{-j}]
\end{aligned}
$$x_a^{(j)} = 2x[\sum_{i=0}^{j-1}q_{i}10^{-i} + a10^{-j}]+x10^{-j}$$
由此定義遞迴式:
\begin{aligned}
x_{a+1}^{(j)} &= x_{a}^{(j)} + 2x10^{-j} \\
y_{a+1}^{(j)} &= y_{a}^{(j)} - 10^{-j}x_a^{(j)}
\end{aligned}
跟之前一樣簡化遞迴式:
\begin{aligned}
z_a^{(j)} &= 10^jy_a^{(j)}\\
z_{a+1}^{(j)} &= z_{a}^{(j)} - x_a^{(j)} \\
x_{a+1}^{(j)} &= x_{a}^{(j)} + 2x10^{-j} \\
\end{aligned}
遞迴式會一直計算 $z_a^{(j)}$ 直到負數，得到係數 $q_j$ 有 $z_{q_j}^{(j)} \ge 0 > z_{q_j + 1}^{(j)}$ 。
對於初始條件可以知道，當從第 $j$ 個位元要計算第 $j + 1$ 位元時，必須將被除數進位也就是乘上 10:
$$z_0^{(j+1)} = 10z_{q_j}^{(j)}$$
另外可以知道，
$$x_0^{(j+1)} = x_{q_j}^{(j)} - x10^{-j} + x10^{-j-1}$$
由上面性質可以得到，當需要計算 $q_{j+1}$ 時，必須減去 $0.9 \times 10^{-j}x$ 。最後利用下面的初始值，
$$y_0^{(0)} = y \ , \ x_0^{(0)} = x$$
完成整個演算法。
#### Binary
接著將演算法推廣至2進制，假設定點數:
\begin{aligned}
y &= x[\sum_{j=0}q_j2^{-j}]^2 \ ,\ \text{for } q_j \in \{0,1\}\\
\sqrt{\frac{y}{x}} &= \sum_{j=0}q_j2^{-j}
\end{aligned}
$$x_a^{(j)} = 2x[\sum_{i=0}^{j-1}q_{i}2^{-i} + a2^{-j}]+x2^{-j}$$
定義遞迴式:
\begin{aligned}
x_{a+1}^{(j)} &= x_{a}^{(j)} + 2x2^{-j} \\
&= x_{a}^{(j)} + x2^{1-j} \\
y_{a+1}^{(j)} &= y_{a}^{(j)} - 2^{-j}x_a^{(j)}
\end{aligned}
簡化遞迴式:
\begin{aligned}
z_a^{(j)} &= 2^jy_a^{(j)}\\
z_{a+1}^{(j)} &= z_{a}^{(j)} - x_a^{(j)} \\
x_{a+1}^{(j)} &= x_{a}^{(j)} + x2^{1-j} \\
\end{aligned}
初始條件為，
$$y_0^{(0)} = y \ , \ x_0^{(0)} = x$$
對於初始條件可以知道，當從第 $j$ 個位元要計算第 $j + 1$ 位元時，必須將被除數進位也就是乘上 2:
$$z_0^{(j+1)} = 2z_{q_j}^{(j)}$$
由前面 10 進制的介紹對於 $x_{q_j}^{(j)}$ 可以得到，
\begin{aligned}
x_0^{(j+1)} &= x_{q_j}^{(j)} - x2^{-j} + x2^{-j-1} \\
&= x_{q_j}^{(j)} + x2^{-j}(-1 + 2^{-1}) \\
&= x_{q_j}^{(j)} - x2^{-j-1} \\
\end{aligned}
#### Implementation
為了將輸入改寫成下面的形式:
$$y = x[\sum_{j=0}q_j2^{-j}]^2 \ ,\ \text{for } q_j \in \{0,1\}$$
需要將 $2^K$ 提取出來使定點數 $Y$ 可以有:
\begin{aligned}
2^0 \le y = \frac{Y}{2^K} &= \frac{Y}{(2^k)^2} = \sum_{i=0}(q_i2^{-i})^2 < 2^2 \\
1 \le \sqrt{y} = \sqrt{\frac{Y}{2^K}} &= \frac{\sqrt{Y}}{2^k} = \sum_{i=0}(q_i2^{-i}) < 2
\end{aligned}
可以發現下面第 7 行就是為了將 $Y$ 轉換至 $[1,4)$ 的區間當中，第 8 行則是相當於提取出 $2^K$。
從前面的說明可以知道如果 $q_k = 1$，則關於更新 $x_{q_k}^{(j)}$ 的遞迴式如 (1)，當 $q_k = 0$ 仍然滿足 $x_{0}^{(j)} = x_{0}^{(j)}$; 而從第 $j$ 位計算 $j + 1$ 時
\begin{aligned}
x_{q_k}^{(j)} &= x_{0}^{(j)} + x2^{1-j}q_k \ \ \ \ \ \ \text{(1)}\\
x_0^{(j+1)} &= x_{q_j}^{(j)} - x2^{-j-1} \ \ \ \ \ \ \ \text{(2)}\\
y_0^{(j+1)} &= 2y_{q_k}^{(j)} = 2(y_{q_k}^{(j)} - 2^{-j}x_{q_k}^{(j)})\ \ \ \ \ \ \text{(3)}\\
\end{aligned}
\begin{aligned}
x_{q_k}^{(j+1)} &= x_{q_{j}}^{(j)} + x2^{-j-1}(-1 + 2q_k) \ \ \ \ \ \ \text{(1)}\\
\end{aligned}
由上面的說明可以知道初始條件為，
$$y_0^{(0)} = y = \frac{Y}{2^K} \ , \ x_1^{(0)} = 1$$
利用 Figure 2. 的流程圖模擬一次計算過程，以定點數 $Y=\text{0x1e400}$ 為例子，一開始利用 $K = 16$ 將 $Y$ 的值域轉換至 $[1,4)$，如下面的式子:
\begin{aligned}
Y &=\text{0x1e400} = 2^{16} + 2^{15} + 2^{14} + 2^{13} + 2^{10} \\
y &=\frac{Y}{2^{16}}=\frac{Y}{(2^{8})^2} = 2^{0} + 2^{-1} + 2^{-2} + 2^{-3} + 2^{-6}\\
\end{aligned}
此時的 $x$ 按照我們的定義應該是 $2^{16}$ ，但是為了實作上方便可以將 $2^{16}$ 從 $\sqrt{\frac{y}{x}}$ 提取出來，最後再將結果乘上 $2^8$，而此時的 $x=1$。 因此可以得到紀錄被除數的暫存器 $A=y = 2^{0} + 2^{-1} + 2^{-2} + 2^{-3} + 2^{-6}$ ，儲存除數暫存器的 $B =x= 1$，$j$ 為紀錄目前正在計算的位元， CNT 為表示 $q_k$ 的係數，Q 為儲存結果的暫存器。計算過程如下面的表格 (利用不同顏色幫助區分計算不同位元係數，可以參考 Meggitt 論文中的 Table 3 理解):
| j | B | A | CNT | Q |
|:---:| --------------------------------------------:| --------------------------------------------:|:---:|:-------:|
| 0 | $2^0$ | $2^0+2^{-1}+2^{-2}+2^{-3}+2^{-6}$ | | |
| | $+) \ \ 2^{1}$ | $-) \ \ 2^0$ | | |
| | = $2^1+2^0$ | $= \ \ 2^{-1}+2^{-2}+2^{-3}+2^{-6}$ | 1 | |
| | $-) \ \ 2^{-1}$ | $\text{lsh 1)}$ | | |
| 1 | $\color\red{= \ \ 2^1+2^{-1}}$ | $\color\red{= \ \ 2^0+2^{-1}+2^{-2}+2^{-5}}$ | | 1 |
| | $\color\red{-) \ \ 2^{-2}}$ | $\color\red{\text{lsh 1)}}$ | | |
| 2 | $= \ \ 2^1+2^{-2}$ | $= \ \ 2^1+2^{0}+2^{-1}+2^{-4}$ | | 10 |
| | $+) \ \ 2^{-1}$ | $-) \ \ 2^1+2^{-2}$ | | |
| | $= \ \ 2^1+2^{-1}+2^{-2}$ | $= \ \ 2^0+2^{-2}+2^{-4}$ | 1 | |
| | $-) \ \ 2^{-3}$ | $\text{lsh 1)}$ | | |
| 3 | $\color\red{= \ \ 2^1+2^{-1}+2^{-3}}$ | $\color\red{= \ \ 2^1+2^{-1}+2^{-3}}$ | | 101 |
| | $\color\red{+) \ \ 2^{-2}}$ | $\color\red{-) \ \ 2^1+2^{-1}+2^{-3}}$ | | |
| | $\color\red{= \ \ 2^1+2^{-1}+2^{-2}+2^{-3}}$ | $\color\red{= \ \ 0}$ | 1 | |
| | ... | $\color\red{\text{lsh 1)}}$ | | |
| 4 | ... | ... | | 1011 |
| | ... | ... | | 0x16000 |
由於計算至第 3 位元時就已經計算完成，此時 A 內的結果為 0 ，之後不管暫存器 B 內的值為多少，結果都為 $A > B$ ，所以 Q 會一直左移直到計算結束。
由前面的等式 $1 \le y = \frac{Y}{2^K} = \sum_{i=0}^{n}(q_i2^{-i}) < 4$ 可以知道要表示 $y$ 會需要至少 $n+1$ 個位元，而每次計算一個位元後 y 都必須向左進一位，因此暫存器 A 會需要至少 $n+2$ 儲存 $y_{q_k}^{(j)}$。
下面程式碼是對應的實作:
```c=
uint32_t UQ16P16_pseudo_division_sqrt(uint32_t n)
{
const int lz = __builtin_clz(n);
const uint32_t one = 1UL << 29;
uint32_t y = n;
uint32_t x = one, d = one << 2;
uint32_t r = 0, count = 0;
y = (y << lz >> 1) >> (lz & 1);
int expo = (FSHIFT - lz + (lz & 1)) / 2;
while (count < (FSHIFT) ) {
d >>= 1;
if (y >= x) {
y -= (x);
x = x + (d);
r += (one >> count >> 14);
}
y <<=1;
x-= (d >> 2);
count++;
}
return (expo >= 0) ? r << (expo) : r >> (-expo);
}
```
這裡附上 [Digit-by-digit calculation](https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Digit-by-digit_calculation) 方式的實作，詳細推導可以參考 wiki 內的說明。
```c=
uint32_t UQ16P16_DBD_sqrt(uint32_t n)
{
uint32_t Y = n;
const int lz = __builtin_clz(n);
uint32_t x = 0;
uint32_t d = 1UL << 30;
Y = (Y << lz) >> (lz & 1);
const int expo = (int) (FSHIFT - 1 - lz) / 2;
if (Y >= x + d) { Y -= (x + d); x = x + (d << 1); } x >>= 1; d >>= 2;
if (Y >= x + d) { Y -= (x + d); x = x + (d << 1); } x >>= 1; d >>= 2;
if (Y >= x + d) { Y -= (x + d); x = x + (d << 1); } x >>= 1; d >>= 2;
if (Y >= x + d) { Y -= (x + d); x = x + (d << 1); } x >>= 1; d >>= 2;
if (Y >= x + d) { Y -= (x + d); x = x + (d << 1); } x >>= 1; d >>= 2;
if (Y >= x + d) { Y -= (x + d); x = x + (d << 1); } x >>= 1; d >>= 2;
if (Y >= x + d) { Y -= (x + d); x = x + (d << 1); } x >>= 1; d >>= 2;
if (Y >= x + d) { Y -= (x + d); x = x + (d << 1); } x >>= 1; d >>= 2;
if (Y >= x + d) { Y -= (x + d); x = x + (d << 1); } x >>= 1; d >>= 2;
if (Y >= x + d) { Y -= (x + d); x = x + (d << 1); } x >>= 1; d >>= 2;
if (Y >= x + d) { Y -= (x + d); x = x + (d << 1); } x >>= 1; d >>= 2;
if (Y >= x + d) { Y -= (x + d); x = x + (d << 1); } x >>= 1; d >>= 2;
if (Y >= x + d) { Y -= (x + d); x = x + (d << 1); } x >>= 1; d >>= 2;
if (Y >= x + d) { Y -= (x + d); x = x + (d << 1); } x >>= 1; d >>= 2;
if (Y >= x + d) { Y -= (x + d); x = x + (d << 1); } x >>= 1; d >>= 2;
if (Y >= x + d) { Y -= (x + d); x = x + (d << 1); }
return (expo >= 0) ? x << (expo) : x >> (-expo);
}
```
#### Error comparision
## Conclusion
在定點數的計算 log 和 開根號的實作當中，可以發現都必須事先使用 `clz` 取得輸入 $Y$ 的 $2^k$ 有效最大值，並且位移至合適的位置開始計算，而從浮點數的格式可以發現，浮點數可以從 expo 的部份直接得到 $Y$ 的 $2^k$ ， mantissa 的部份也剛好對齊不必在演算法開始及最後時位移，從這裡發現浮點數在格式上的設計剛好可以符合演算法的需求，以減少額外的計算開銷，實際研究過一次後對浮點數的設計有更多的體悟。
## Reference
1. [Pseudo Division](https://archived.hpcalc.org/laporte/Meggitt_62.pdf)
2. [stackoverflow: How do I compute a fixed point binary logarithm?](https://stackoverflow.com/questions/32159300/compute-logarithmic-expression-without-floating-point-arithmetics-or-log)
3. [stackoverflow: Compute logarithmic expression without floating point arithmetics or log](https://stackoverflow.com/questions/76908525/how-do-i-compute-a-fixed-point-binary-logarithm)
4. [DIGIT BY DIGIT Methods](https://archived.hpcalc.org/laporte/digit_by_digit.htm)
5. [Pseudo-Division Algorithms for Floating-Point Logarithms and Exponentials](https://ieeemilestones.ethw.org/w/images/3/30/Wk_pseudo_division_log_exp_may01.pdf)
6. [Methods of computing square roots](https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Binary_numeral_system_(base_2))
7. [Computing floating-point logarithms with fixed-point operations
](https://hal.science/hal-01227877/document)
8. [The Difference Method of Henry Briggs](https://archived.hpcalc.org/laporte/The%20method%20of%20Henry%20briggs.htm)
9. [Henry Briggs and the HP35](https://archived.hpcalc.org/laporte/Briggs%20and%20the%20HP35.htm)
10. [A reconstruction of the tables of Briggs’ Arithmetica logarithmica](https://locomat.loria.fr/briggs1624/briggs1624doc.pdf)
11. [《HPM 通訊》第 13 卷第 12 期 : 數學史融入教學─以對數表為例](https://hpmsociety.tw/wp-content/uploads/hpmletters/hpm13(12).pdf) )
12. [The Radix Method](https://archived.hpcalc.org/laporte/The%20Radix%20Method.htm)
13. https://www.wolframalpha.com/input?i=product+%281%2B2%5E%28-j%29%29%2C+j%3D1+to+infinity
14. [evaluating infinity production](https://math.stackexchange.com/questions/1924882/evaluating-prod-n-1-infty-left1-frac12n-right)
Google Drive
Gist
Clipboard
Sign in with Wallet
