前言
把第一份工作后写的定点数运动-碰撞检测库重新拿出来看了下，打算手搓个 go 的定点数。
搓完四则运算就来到了第一个坑点，开方。
手动开方法常用的有两种：
- 牛顿-拉弗森法
- 逐比特确认法
完成这篇博客花费的时间长过我的预期，所以分成两部分完成。
牛顿-拉弗森法
我们求一个数 $a$ 的平方根，可以看作求解方程 $x^2-a=0$ 的正根。我们可以画出函数 图像 $f(x)=x^2-a$ 。
牛顿-拉弗森法的基本思路是可以任取一点 $x_1$ 作为初始值，当 $f(x_1)=0$ 时 $x_1$ 就是 $a$ 的根，如果不是，那么从 $f(x)=x^2-a$ 上这一点做切线，切线与 x 轴的交点会更接近 $a$ 的根。
我们先列出这个函数的导函数：$f'(x)=2x$
还有切线方程：$m(x-x_n)=y-y_n$ ，其中 n 是迭代轮数，$x_n$ 表示第 n 轮迭代的 x 。 $m$ 是原函数在 $x_n$ 的导数。
变形可以得到求切线零点的方程：$m(x-x_n)+y_n=0$
同时，我们知道这个方程中 $x$ 、 $m$ 、 $y_n$ 的的含义：
- $x$ 是我们下一轮迭代的取值，即 $x_{n+1}$
- $m$ 是 $f'(x_n)$
- $y_n$ 是 $f(x_n)$
替换进方程后可得到：
$$ f'(x_n)(x-x_n)+f(x_n)=0 $$
接下来整理下这个方程：
$$ \begin{aligned} f'(x_n)(x_{n+1}-x_n)+f(x_n)&=0 \\ x_{n+1}-x_n &= -\frac{f(x_n)}{f'(x_n)} \\ x_{n+1} &= x_n - \frac{f(x_n)}{f'(x_n)} \end{aligned} $$
就得到了牛顿-拉弗森迭代法公式。
求给定数字 $a$ 的平方根根的场景，我们有：
$$ \begin{aligned} f(x) &= x_n^2-a \\ f'(x) &= 2x_n \\ \end{aligned} $$
可以进一步代入并化简这条公式
$$ \begin{aligned} x_{n+1} &= x_n - \frac{f(x_n)}{f'(x_n)} \\ &= x_n - \frac{x_n^2-a}{2x_n} \\ &= \frac{2x_n^2 - (x_n^2-a)}{2x_n} \\ &= \frac{x_n^2+a}{2x_n} \\ &= \frac{1}{2}\frac{x_n^2+a}{x_n} \\ &= \frac{1}{2}(x_n+\frac{a}{x_n}) \\ \end{aligned} $$
在迭代结束后，我们就有了$x_n\approx\sqrt{a}$ 。
NOTE：对于开根号的这个特定场景，$x^2-a=0$ 函数图像的顶点一定在 $x=0$ 只要避开 $x=0$ 的 猜测值，牛顿法就不会发散。不过选择良好的初始值还是有益的，可以加快收敛速度。
牛顿-拉弗森法求商
牛顿-拉弗森法可以用于实现快速除法运算。
对于除法我们可以列出这样的方程：$q=\frac{n}{d}$
其中：
- q 是商
- n 是被除数
- d 是除数
我们知道除法可以表示为被除数和除数的倒数的乘积，由此得到方程：$q=n\frac{1}{d}$
只要能求得 $\frac{1}{d}$ 的近似值，就可以用乘法运算求得 $q$ 的近似值。
我们设 $x=\frac{1}{d}$ 可以得到 $d=\frac{1}{x}$ ，写成关于 $x$ 的函数就是：$f(x)=d-\frac{1}{x}$ ，我们求函数的零点，即 $d-\frac{1}{x}=0$
这个函数的导函数是 $f'(x)=\frac{1}{x^2}$
如果直接使用牛顿-拉弗森迭代法，会发现存在一个递归定义：$\frac{1}{x}$ 依然是除法。
这里需要将牛顿-拉弗森迭代公式展开化简：
$$ \begin{aligned} x_{n+1} &= x_n-\frac{f(x_n)}{f'(x_n)} \\ x_{n+1} &= x_n-\frac{d-\frac{1}{x_n}}{\frac{1}{x_n^2}} \\ x_{n+1} &= x_n-(d-\frac{1}{x_n})x_n^2 \\ x_{n+1} &= x_n-x_n^2d+x_n \\ x_{n+1} &= 2x_n-x_n^2d \\ x_{n+1} &= x_n(2-dx_n) \end{aligned} $$
化简后的牛顿-拉弗森求商迭代公式就没有除法了。
我们得到迭代后的 $x_n$ 是原除数的倒数 $\frac{1}{d}$ 的近似，此时就可以用乘法计算 出原除法的商了：
$$ \begin{aligned} q&=n\frac{1}{d} \\ &=nx_n \end{aligned} $$
逐比特确认法
逐比特确认法是一种类似长除法(竖式除法开根号)的思路，从高位开始逐个比特确认。
例如，设被开方数为 $N$ ，$N$ 的比特数 $B_N=8$ ，也就是 $N$ 的最高位是 8 比特。
从根的最高位第4个比特开始判断是否为1，当第4个比特为是1时，根是 $2^3=8$ 。 如果 $N>2^3$ ，那么根的第4个比特位就是 $1$ ，反之则是 $0$ 。
确认根的比特数
那么根的最高位怎么确认呢？最直觉告诉我们，正整数的根肯定小于原数，我们从 $N$ 的最高位开始肯定没问题，但这样会引入无效计算。
我们定义任意正整数最小需要 $B_X$ 位表示，可得
$$ B_X=\lfloor{\log_{2}{X}}\rfloor+1 $$
可设 $B_r$ 为 $N$ 的根的比特数：
$$ \begin{aligned} B_r&=\lfloor\log_{2}{\sqrt{N}}\rfloor+1 \\ &=\lfloor\log_{2}{N^{\frac{1}{2}}}\rfloor+1 \\ &=\lfloor\frac{1}{2}\log_{2}{N}\rfloor+1 \\ \end{aligned} $$
$N$ 的最高位是 $B_N$ 位时，$N$ 满足：
$$ 2^{B_N-1}\le{N}\lt{2^{B_N}} $$
将不等式的每一项都求2为底的对数，得到：
$$ B_N-1\le{\log_{2}{N}}\lt{B_N} $$
所有项全部除以2，得到：
$$ \frac{B_N-1}{2}\le{\frac{1}{2}\log_{2}{N}}\lt{\frac{B_N}{2}} $$
由于 $B_r$ 的定义不能直接代入，我们先设 $x=\frac{1}{2}\log_{2}{N}$ 代入不等式考察。
$$ \frac{B_N-1}{2}\le{x}\lt{\frac{B_N}{2}} $$
现在我们开始考虑 $B_N$：
偶数：
如果 $B_N$ 是偶数，我们可以定义 $B_N=2k$
$$ \frac{2k-1}{2}\le{x}\lt{\frac{2k}{2}} $$
化简得到
$$ k-\frac{1}{2}\le{x}\lt{k} $$
基于 $x$ 的范围我们向下取整，得到 $\lfloor{x}\rfloor=k-1$ ，代入 $B_r$ 定义可得 $B_r=k-1+1=k$
奇数：
如果 $B_N$ 是奇数，我们定义 $B_N=2k+1$
$$ \frac{2k+1-1}{2}\le{x}\lt{\frac{2k+1}{2}} $$
化简得到
$$ k\le{x}\lt{k+\frac{1}{2}} $$
我们向下取整，得到 $\lfloor{x}\rfloor=k$ 。 代入 $B_r$ 定义可得，当 $B_N$ 是奇数时， $B_r=k+1$ ，代回推导 $B_r$ 和 $B_N$ 的关系：
$$ \begin{aligned} 2k+1&=B_N \\ 2k &=B_N-1 \\ k &=\frac{B_N-1}{2} \\ k+1 &=\frac{B_N-1}{2}+1 \\ B_r &=\frac{B_N-1+2}{2} \\ B_r &=\frac{B_N+1}{2} \end{aligned} $$
得出结论：
$$ \begin{cases} B_r = \frac{B_N}{2} &\text{if } B_N=2k, k\in \mathbb{Z} \\ B_r = \frac{B_N+1}{2} &\text{if } B_N=2k+1, k \in \mathbb{Z} \end{cases} $$
逐比特迭代
定义根的第 $n$ 个比特位为 $r_n$，$n\ge{0}$，每一轮迭代确认未被确认的最高位。 $R_n$ 为已经确认的根。
$$ \begin{cases} r_{n-1} = 1 &\text{if } N\ge{(R_n+2^{n-1})^2} \\ r_{n-1} = 0 &\text{if } N\lt{(R_n+2^{n-1})^2} \end{cases} $$
确认 $r_n$ 后加入已确认的根 $R_n$ ： $R_{n-1} = R_{n} + r_{n-1}2^{n-1}$
然后开始下一轮迭代。
逐比特迭代（余数）
我们也可以使用余数迭代，避免重复计算 $(R_n+2^{n-1})^2$，我们定义余数 $X_n$ 替换原式子中的 $N$。
$$ \begin{cases} r_{n-1} = 1 &\text{if } R_n^2+X_n\ge{(R_n+2^{n-1})^2} \\ r_{n-1} = 0 &\text{if } R_n^2+X_n\lt{(R_n+2^{n-1})^2} \end{cases} $$
将 $(R_n+2^{n-1})^2$ 展开：
$$ \begin{aligned} (R_n+2^{n-1})^2 &= R_n^2 + 2R_n2^{n-1} + 2^{2n-2} \end{aligned} $$
代回原不等式，将不等式化简：
$$ \begin{aligned} R_n^2+X_n &\ge{(R_n+2^{n-1})^2} \\ R_n^2+X_n &\ge{R_n^2 + 2R_n2^{n-1} + 2^{2n-2}} \\ X_n &\ge{2R_n2^{n-1} + 2^{2n-2}} \\ X_n &\ge{R_n2^{n} + 2^{2n-2}} \end{aligned} $$
余数更新公式的定义同理。
确认到第 $n$ 比特位时，余数 $X_n$ 的定义是：$X_n = N-R_n^2$ 。
我们把 $N$ 替换为 $R_n^2+X_n$ 再考虑新迭代出的比特位，得到余数的迭代公式。
$$ \begin{aligned} X_{n-1} &= N-(R_n+r_{n-1}2^{n-1})^2 \\ X_{n-1} &= R_n^2+X_n-(R_n+r_{n-1}2^{n-1})^2 \\ X_{n-1} &= R_n^2+X_n-(R_n^2 + 2 \cdot R_n \cdot r_{n-1} \cdot 2^{n-1} + r_{n-1}^2\cdot2^{2n-2}) \\ X_{n-1} &= R_n^2+X_n-R_n^2 - 2R_n\cdot r_{n-1} \cdot 2^{n-1} - r_{n-1}^2 \cdot 2^{2n-2} \\ X_{n-1} &= X_n - 2R_n\cdot r_{n-1} \cdot 2^{n-1} - r_{n-1}^2 \cdot 2^{2n-2} \\ \end{aligned} $$
根的迭代公式则没有变化。
$$ \begin{aligned} R_{n-1} &= R_{n} + r_{n-1}2^{n-1} \end{aligned} $$