前言
在定点数开方算法文章的第1部分中，我们介绍了定点数开根号的几种算法原理，包括牛顿-拉弗森法和逐比特确认法。
在定点数开方算法文章的第2部分中，我们具体推导了牛顿-拉弗森法的定点数推广和实现。
本文继续写逐比特确认法的实现细节。
逐比特确认法
确认根位宽
在第一部分里我们证明了根的位宽 $B_r$：
$$ B_r = \begin{cases} \frac{B_N}{2} &\text{if } B_N=2k, k\in \mathbb{Z} \\ \frac{B_N+1}{2} &\text{if } B_N=2k+1, k \in \mathbb{Z} \end{cases} $$
这里不作赘述。
迭代公式
不使用余数的迭代公式是最直接的，同样在第1部分已经有过论述和证明。
确认比特位的公式是：
$$ r_{n-1} = \begin{cases} 1 &\text{if } N\ge{(R_n+2^{n-1})^2} \\ 0 &\text{if } N\lt{(R_n+2^{n-1})^2} \end{cases} $$
已经求得的根的迭代公式则是：
$$ R_{n-1} = R_{n} + r_{n-1}2^{n-1} $$
定点数推广
逐比特确认法求的是正整数的根，而我们需要的是求实数 $V$ 的根 $B$。
这一步我们回顾第2部分牛顿法的定点数推广，有下面的结论：
$n$ 比特精度定点数 $V_{int}$ 的实数根 $R$ 的定义：
$$ R = \frac{\sqrt{V_{int}}}{2^{\frac{n}{2}}} $$
实数根 $R$ 的 $n$ 比特精度定点数形式 $R_{int}$ 定义：
$$ R_{int} = \sqrt{V_{int}} \cdot 2^{\frac{n}{2}} $$
我们使用逐比特确认法迭代得到 $\sqrt{V_{int}}$ 的值后，还要经过一个缩放因子 $2^{\frac{n}{2}}$ 才能得到定点数形式的根 $R_{int}$ 。
考虑到精度问题（同样在第2部分有过论述，我们使用 Q31.32 迭代，最终得到
$$ \sqrt{V_{q32}} = \sqrt{V} \cdot 2^{16} $$
具体证明过程见第2部分。
代码实现
下述代码的
Fix32类型是 Q15.16 格式。
// BitByBitSquareRoot performs fixed-point arithmetic square root of a Fix32 number in bit-by-bit method
func (f Fix32) BitByBitSquareRoot() (Fix32, error) {
if f < 0 {
return 0, ErrNaN
}
if f == 0 {
return 0, nil
}
var (
x int64
v = int64(f) << 16 // scale up to Q31.32
bitwidth = bits.Len64(uint64(v))
rootBitwidth int
)
if bitwidth&1 > 0 {
rootBitwidth = (bitwidth + 1) / 2
} else {
rootBitwidth = bitwidth / 2
}
for i := rootBitwidth-1; i >= 0; i-- {
intermediate := x | (1 << i)
if v >= intermediate*intermediate {
x |= 1 << i
}
}
return Fix32(x), nil
}
关于余数法
第 $X_{n-1}$ 比特的确认公式是：
$$ X_{n-1}= \begin{cases} 1 &\text{if } X_n \ge{R_n2^{n} + 2^{2n-2}} \\ 0 &\text{if } X_n \lt{R_n2^{n} + 2^{2n-2}} \end{cases} $$
余数的迭代公式为：
$$ X_{n-1} = X_n - 2R_n\cdot r_{n-1} \cdot 2^{n-1} - r_{n-1}^2 \cdot 2^{2n-2} $$
根的迭代公式
$$ R_{n-1} = R_{n} + r_{n-1}2^{n-1} $$
余数法的实现会比直接比较复杂很多，还要考虑到余数法公式中 $2^{2n-2}$ 这样的式子转为代码时的复杂性，性能也没有明显优势，这里不再展开。