牛顿迭代法求近似倒数及其应用
本帖最后由 DEATHTOUCH 于 2023-10-27 00:41 编辑# 牛顿迭代法求近似倒数及其应用
这是继上篇软浮点（传送门：[软件模拟实现IEEE-754单精度浮点数运算](https://www.52pojie.cn/thread-1830228-1-1.html)）后进一步的研究，主要是为了找到除法运算一种更快的方法。
全部的代码依然在代码仓库 <https://gitee.com/peazomboss/softfloat32> 可以找到，放在了文件夹 _approxrecip_ 里面，同时软浮点除法部分的代码也已经更新。
有关牛顿法的介绍，可以看文末参考资料 *Hacker's Delight* 的附录B，或者上网搜索。
## 公式简介
首先看一下牛顿迭代公式：
$x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}$
对于求一个数 $a$ 的倒数，我们可以用公式 $f(x) = \frac{1}{x} - a$ 计算，令 $f(x) = 0$，其根就是 $a$ 的倒数。将该式子带入上述迭代公式，化简得到：
$x_{n+1} = x_n(2-ax_n)$
用这个式子不断迭代就可以得到数 $a$ 的倒数了，不过对于初值的估计是一个问题。因为如果估计偏差过大，可能会导致误差不断变大而无法在较短的步数内收敛。同时，这个式子对于浮点数来说实现起来完全没什么问题，不过如果要定点化就稍微需要改一下了。
## 浮点数实现
对着上述公式，很容易就可以写出浮点数的版本：
```cpp
float newton_recip(float a)
{
float x = get_initial(a); // 获取初值，可以用查表法
for (int i=0; i<N; i++) // 可以展开循环，N不用很大
x = x * (2 - a*x);
return x;
}
```
对于浮点数的初值，可以获取其指数，然后根据指数大小来查表估计。表生成方法就是遍历不同指数的浮点数，然后提前计算其倒数即可。
不过由于现代支持浮点数操作的CPU通常有较快的运算，所以这个方法并不怎么实用，实际上直接用 `1.0 / a` 就可以得到倒数了。因此牛顿法求倒数主要是给定点小数来使用的。
## 定点小数实现
假设有两个16位的定点小数相除，一种可行的方案是将被除数扩展到32位。对于一个32位数来说，如果除数是16位的，那么商必然是16位或17位的数。只有被除数小于除数，那么其扩展到32位后做除法的结果才不至于是17位。也就是说为了确保除法的结果是不会溢出的，我们需要确保被除数小于除数。
对于32位的数来说，如果按照这个方法，需要扩展到64位去除。但是64位数的除法不一定是所有32位的CPU都支持的，那么一个可行的方法就是先求除数的倒数，然后再用被除数去乘这个倒数，之后修正一下商即可。
基于此，可以定义32位定点小数规范如下：
- Q1.31：表示1位整数和31位尾数，范围是1.0-2.0，其中2.0取不到。
- Q32：表示32位尾数，所以范围是0-1.0，其中1.0取不到。
回到牛顿迭代上来，我们可以把迭代的每一步拆开来：
```cpp
t = a * x;
t = 2 - t;
x = x * t;
```
1. 对于 `t = a * x`，可以确定变量 a 是Q31的，而 x 必然是Q32的，那么其结果是Q1.63，所以实际上可以取其高32位，而舍弃低32位，从而变成Q1.31。
2. 对于 `t = 2 - t`，由于此时变量 t 是Q1.31的，而 2 可以理解为是Q2.31，那么可以确定的是 2 - t 之后的结果依然是Q1.31的。
3. 对于 `x = x * t`，显然 x 是Q32，t 是Q1.31，相乘是Q1.63，但是由于下一轮迭代的 x 需要Q32，所以可以使结果右移31位，获取Q1.63的Q63部分。
所以写成实际的代码应该是这样的：
```cpp
uint32_t fixedpoint_newton_recip(uint32_t a)
{
uint32_t x = get_initial(a); // 获取初值，可以用查表法
for (int i=0; i<N; i++)
{
uint32_t t = ((uint64_t)a * x) >> 32;
t = -t;
x = ((uint64_t)x * t) >> 31;
}
return x;
}
```
这里值得一提的是迭代过程中的一句 `t = -t;`，看似是取负数，实际上可以理解为 `t = 2 - t`，这是因为 2 是用Q2.31表示的，表示为 10.000...，这样做减法实际上就是取负数，或者说是补码。
不过事情没有那么简单，因为如果简单取负数，会发现一个奇怪的问题，即迭代后的估计值偏大。这很不正常，因为根据牛顿法的原理，对于 $f(x) = \frac{1}{x^2}-a$ 这样的函数来说，其在坐标上画出来是一个凹函数，所以最终的估计值理论上是偏小的。
但是仔细分析一下就会发现，这是由于截断造成的误差。在代码中我们求出了 a 和 x 的乘积后，只是简单右移，忽略了低32位的结果。尽管这个值很小，但是仍然会造成1位的误差，也就是偏小了1，此时取补码显然是偏大了1，再去相乘，最终的结果自然也会偏大1或2。
考虑到求补码实际上是按位取反再加一的操作，所以完全可以直接按位取反，也就是求反码即可。这样确保了最终结果存在的误差只能是偏小的。所以迭代步骤中的 `t = -t;` 一句就需要修改为 `t = ~t;`，这样可以让结果更加符合预期，避免造成奇怪的问题。
### 初值估计
和浮点数一样，估算初值是一个重要的操作。由于我们已经定义了Q1.31的格式，所以这实际上不难实现。我们可以根据参数 a 的取值范围，选择一个合适的初值。因为 a 的最高位必定是1，所以可以对后面的3-4位进行估计。经过我的实测，3位就可以满足要求的，如果要求更高的精度，那么4位也是可以的。
然后就是表的宽度了，一般情况下对于32位来说，如果初值能达到4位的精度，那么只要3次迭代就可以达到接近32位的精度了。为了能够利用到牛顿法的平方收敛性，最好将初值的精度提高一点，使用8位就足够了，用16位不如把表先扩大一点。
使用查表法，获取 a 的第30-28位，然后根据这3位从表中获取初值，具体方法如下：
```cpp
x = initial_table_8x8[(a >> 28) & 7] << 24;
```
如果是估算4位的数据，那么就需要获取第30-27位，方法如下：
```cpp
x = initial_table_16x8[(a >> 27) & 15] << 24;
```
这里假设表是Q8格式的，需要左移将其扩展到Q32格式。
生成表的方法也不难，因为 a 是Q1.31的，倒数是Q32的，所以根据 a 值提前计算倒数即可。a 的高4位一定是8~15，为了一定精度，可以取 a 的高5位进行计算，多出一位的作用是取平均值。比如说获取到高4位是8，其取值是 $[1,\frac{9}{8})$，则倒数取值是 $(\frac{8}{9},1)$。但是如果按照取值是1，那么其倒数就只能是近似于0.9999...了；而取值接近 $\frac{9}{8}$ 的时候，这个倒数就会有较大的误差。所以可以取一个折中的方案，也就是 $\frac{17}{16}$，这样对于各种情况，倒数的估计会更加准确。
具体生成表的代码如下：
```cpp
void gentable8x8()
{
printf("TABLE 8x8:\n");
for (int i=8; i<=15; i++) {
uint32_t u = (i << 28) + 0x8000000;
uint32_t x = (uint32_t)0x7fffffff / (u >> 24);
uint16_t w = x;
x >>= 16;
if (w >= 0x8000)
x++;
printf("0x%x, ", x);
}
putchar(10);
}
void gentable16x8()
{
printf("TABLE 16x8:\n");
for (int i=0; i<=15; i++) {
uint32_t u = 0x80000000 + (i << 27) + 0x4000000;
uint32_t x = (uint32_t)0x7fffffff / (u >> 24);
uint16_t w = x;
x >>= 16;
if (w >= 0x8000)
x++;
printf("0x%x, ", x);
if ((i+1) % 8 == 0)
putchar(10);
}
putchar(10);
}
int main()
{
gentable8x8();
gentable16x8();
}
```
由于实际上就用了5-6位，所以完全可以把这个数当成Q1.7来看待，那么使用Q1.31表示的1来除以它，就可以得到Q24的结果了。对于这个结果，保留其高8位作为Q16的表，同时使用其低16位用来判断是否舍入。比如当 i = 8 的时候，其结果应该是0xf0f0f0，那么按照直接截断法，应该保留0xf0作为Q8，而将其舍入到0xf1则可以提高一定的精度（虽然不多）。
上述代码运行输出的内容如下：
```
TABLE 8x8:
0xf1, 0xd8, 0xc3, 0xb2, 0xa4, 0x98, 0x8d, 0x84,
TABLE 16x8:
0xf8, 0xea, 0xdd, 0xd2, 0xc8, 0xbf, 0xb6, 0xae,
0xa7, 0xa1, 0x9b, 0x95, 0x90, 0x8b, 0x86, 0x82,
```
这样初值的查找表就构造完成了，第一种表共8个元素，每个元素1字节，一个只有8字节；第二种表则多了一倍，16个元素，占16字节。至于表是不是越大越好，后面会进行分析。
### 完整代码
这里给出实际可以使用的代码：
```cpp
const uint8_t initreciptable8x8 = {
0xf1, 0xd8, 0xc3, 0xb2, 0xa4, 0x98, 0x8d, 0x84
};
const uint8_t initreciptable16x8 = {
0xf8, 0xea, 0xdd, 0xd2, 0xc8, 0xbf, 0xb6, 0xae,
0xa7, 0xa1, 0x9b, 0x95, 0x90, 0x8b, 0x86, 0x82
};
uint32_t getapproxrecip8x8(uint32_t n)
{
uint32_t x, t;
x = initreciptable8x8[(n >> 28) & 7] << 24; // 查表估计初值
// 第一轮迭代
t = ((uint64_t)x * n) >> 32;
t = ~t;
x = ((uint64_t)x * t) >> 31;
// 第二轮迭代
t = ((uint64_t)x * n) >> 32;
t = ~t;
x = ((uint64_t)x * t) >> 31;
// 第三轮迭代
t = ((uint64_t)x * n) >> 32;
t = ~t;
x = ((uint64_t)x * t) >> 31;
return x;
}
uint32_t getapproxrecip16x8(uint32_t n)
{
uint32_t x, t;
x = initreciptable16x8 << 24;
t = ((uint64_t)x * n) >> 32;
t = ~t;
x = ((uint64_t)x * t) >> 31;
t = ((uint64_t)x * n) >> 32;
t = ~t;
x = ((uint64_t)x * t) >> 31;
t = ((uint64_t)x * n) >> 32;
t = ~t;
x = ((uint64_t)x * t) >> 31;
return x;
}
```
这里提供了两种代码，主要差别是初值表的不同。二者都是迭代了3次，因为对于这样的查找表，其估算精度实际上有大约4位，3次迭代则可以达到约32位。函数 `getapproxrecip8x8` 使用了表 `initreciptable8x8`，精度稍微差一些，但是表节省了8字节的空间；而函数 `getapproxrecip16x8` 则用了更大的表，有更高的精度，最终修正误差也会更快一点。至于更大的表其实没有必要了。
## 误差分析
对于上述提供的代码，需要进行测试分析实际产生的误差。对于32位数来说，可以使用CPU自带的64位除法来计算其倒数，然后和上述方法产生的结果进行比对，即可得到实际的误差。然后再可以使用随机数进行测试，得到误差的分布，从而计算出大致的误差范围。
对于一个Q1.31的倒数来说，可以使用随机数覆盖测试，不过由于数据量不大，只有 $2^{31}$ 个，完全可以使用穷举测试。
```cpp
uint32_t (*getapproxrecip)(uint32_t); // 求近似倒数的函数指针
void bruteerrortest()
{
int maxpos, maxneg, eq, neg1, neg2, neg3;
maxpos= 0; mindiff = 0; eq = 0; neg1 = 0; neg2 = 0; neg3 = 0;
for (uint32_t a = 0x80000000; (int32_t)a < 0; a++) { // 遍历
uint32_t recip1 = getapproxrecip(a); // 求近似倒数
uint32_t recip2 = 0x7fffffffffffffff / a; // 求实际倒数
int diff = recip1 - recip2; // 求误差
if (diff == 0) // 误差为0
eq++;
else if (diff == -1) // 误差为-1
neg1++;
else if (diff == -2) // 误差为-2
neg2++;
else if (diff == -3) // 误差为-3
neg3++;
if (diff > maxpos) // 最大正误差
maxpos = diff;
if (diff < maxneg) // 最大负误差
maxneg = diff;
}
printf(" %d, %d, %d\n", eq, maxpos, maxneg);
printf("[-1] %d, [-2] %d, [-3] %d\n", neg1, neg2, neg3);
}
void brutetest()
{
getapproxrecip = getapproxrecip8x8; // 测试8x8表的误差
printf("8x8:\n");
bruteerrortest();
getapproxrecip = getapproxrecip16x8; // 测试16x8表的误差
printf("16x8:\n");
bruteerrortest();
}
```
以下是暴力测试的结果：
```
8x8:
874319370, 0, -3
[-1] 1126936446, [-2] 145419076, [-3] 808756
16x8:
971865634, 0, -3
[-1] 1050310821, [-2] 125307190, [-3] 3
```
可以看到8x8的表已经非常不错了，误差为-3的概率只有0.038%。而16x8的表主要进一步减少了-3的个数，不过居然有3个漏网之鱼，这就让我好奇这3个数究竟是什么。然后就用暴力法再搜一下：
```cpp
void findneg3()
{
for (uint32_t a = 0x80000000; (int32_t)a < 0; a++) {
uint32_t recip1 = getapproxrecip16x8(a);
uint32_t recip2 = 0x7fffffffffffffff / a;
int diff = recip1 - recip2;
if (diff == -3)
printf("%x\n", a);
}
}
```
发现是0x80083b6a、0x80083f4c和0x8011120c这三个数，这就好办了，因为显然这三个数都是通过查表的第一个元素0xf8来计算的，那是不是就说明0xf8不是很准确呢？
把0xf8改成0xf7，发现误差数据变多了，于是改成0xf9，发现找不到误差是-3的值了。
以下是修改之后的暴力测试结果：
```
16x8:
970775900, 0, -2
[-1] 1050401445, [-2] 126306303, [-3] 0
```
不过这样改了意义不大，因为仔细看会发现虽然没有了误差是-3的，但是误差是-1和-2的情况有所增加，总体上是更差的。而且更重要的是，对于实际的除法应用，这点差别更是可以忽略不计，所以表的内容有个大概就行了，没必要为了微弱的提升去找到一张所谓最好的表。
总体看来，使用8x8和16x8的表没有明显的差距，所以实际上根据情况选择即可。对于代码体积要求高的情况下可以选择8x8的表，而16x8的表可以减少约8.5%的误差，不过对于实际情况并没有明显的改善。
## 定点小数除法
设变量 a 是被除数，变量 b 是除数，现在要计算商 q = a / b 以及余数 r = a % b，最简单的办法就是用语言提供的除法指令，但是不要忘了我们现在用的是定点小数，直接使用除法会造成大量的误差（想一下 `0x1234 / 0x2345` 的结果）。所以下面介绍一般的方法，以及利用上述牛顿法求得的近似倒数来计算定点小数除法的方法。
### 一般的方法
对于定点小数的除法，为了保证精度，一种简单有效的方案就是长除法。长除法可以理解为被除数的位宽是除数的两倍，比如64位数除以32位数，结果保留到32位数。Intel 的CPU一直都是支持长除法的，比如8086的除法指令 `div`，如果参数是8位寄存器，那么被除数是 `ax` 寄存器；如果参数是16位寄存器，那么被除数是由 `ax` 和 `dx` 组成的32位数，其中 `ax` 是低16位。同样的，对于386以后的32位CPU，则是 `eax` 和 `edx`；对于x86-64则是 `rax` 和 `rdx`。不过即使 Intel 的CPU支持长除法，但是对于绝大多数语言来说，并不能直接生成这样的代码，所以一般情况下使用长除法就需要手动写汇编来实现了。
注意到这样一个细节：对于除法 `0x8000 / 0x80`，其结果是0x100，无法用一个8位寄存器保存。扩展一下，可以确定对于位宽2W的被除数去除以一个位宽W的除数，如果被除数的高W位值大于等于除数值，那么其商必须用W+1位来保存。
所以为了避免除完之后结果溢出，需要确保被除数的高位是小于除数的。这实际上不难实现。根据之前的定义，Q1.31的数 a 去除以同样是Q1.31的数 b，使用长除法的条件是如果 a 小于 b，则扩展 a 到Q1.63，否则扩展到Q2.62。
用代码来表示就是这样：
```cpp
uint32_t ulong_div(uint32_t a, uint32_t b)
{
uint64_t a64;
if (a < b)
a64 = (uint64_t)a << 32;
else
a64 = (uint64_t)a << 31;
return a64 / b;
}
```
这样就实现了32位定点小数的除法了。如果硬件不支持一个64位数去除以一个32位数，也可以使用软件的方法实现，这里一个常用的方法是恢复余数法，还有一个就是用两次短除法实现一次长除法。这在 **Hacker's Delight** 有所提到，有兴趣可以去阅读。
### 近似倒数法
该方法的基本原理是对于Q1.31的被除数和除数，按照之前的方法对除数求倒数，其结果就是Q32的。因此计算 $a / b$，就可以计算 $a * \frac{1}{b}$，由于 $a$ 和 $b$ 是Q1.31格式，那么 $\frac{1}{b}$ 自然就是Q32了。
对于Q1.31和Q32的数相乘，显然结果是Q1.63，由于假定了 $a$ 都是小于 $b$ 的，所以这个结果是小于1的，即结果的最高位是0，于是右移31位可以得到所需的近似结果。
实现的代码如下：
```cpp
uint32_t div_approx(uint32_t a, uint32_t b)
{
uint32_t recip_b = getapproxrecip(b);
uint64_t quo = (uint64_t)a * recip_b;
return quo >> 31;
}
```
当然这个除法是存在误差的，需要去修正这个误差。
### 误差修正
由于之前已经测试出了0x80000000到0xffffffff的倒数误差分布情况，所以可以确定倒数的误差在 $[-3,0]$ 的范围内，因为最大可能出现的误差是-3，在乘2的情况下可以达到-6，又因为右移31位可能产生的截断误差导致结果偏小的情况，所以最大误差可以达到-7，这个值刚好是3位的范围内。当然如果两个数的低位的0数量越多，精度就越高，但是通常误差不可消除，所以需要进行修正。
在修正误差之前，我们先看看实际的误差分布情况。
由于两个数相除的组合情况非常多，如果还像之前暴力穷举的方法那肯定不行了（一般的电脑是吃不消的），所以需要使用随机数来进行操作。由于C语言的 `rand()` 函数的最大值是32767，不能满足我们所需要的Q1.31的随机数的生成，所以需要先写一个随机数生成器（简称RNG）。
这里写一个简单的用线性同余法（Linear Congruential Method）实现的随机数生成器（线性同余生成器，LCG，Linear Congruential Generator，不是52的LCG哦），这个随机数生成算法是Delphi默认使用的随机数算法，因为简单所以拿来用了。
```cpp
uint32_t lcg_seed;
uint32_t lcg_rand()
{
lcg_seed = lcg_seed * 134775813 + 1;
return lcg_seed;
}
```
至于随机数的初值，在main函数的开头用 `lcg_seed = time(NULL);` 就可以了。
在之前的 `div_approx()` 函数基础上，再使用如下代码进行一次除法的误差分析：
```cpp
static int check_div_diff(uint32_t a, uint32_t b)
{
uint32_t q = ((uint64_t)a << 32) / b;
uint32_t q2 = div_approx(a, b);
return q2 - q;
}
```
随机测试代码如下：
```cpp
void random_test_error()
{
const int LOOP_MAX = 10000000;
int n = {0};
for (int i=0; i<LOOP_MAX; i++) {
int a = (lcg_rand() % 0x80000000 + 0x80000000);
int b = (lcg_rand() % 0x80000000 + 0x80000000);
if (a >= b) // 保证a小于b
a >>= 1;
int diff = check_div_diff(a, b);
diff = -diff;
n++;
}
for (int i=0; i<8; i++) {
printf("[-%d]%7.4f%%, ", i, 100.0 * n / LOOP_MAX);
if ((i+1) % 4 == 0)
putchar(10);
}
uint64_t total_fix = 0; // 总修正次数
for (int i=1; i<8; i++)
total_fix += i * n;
printf("average fix times = %.3f\n", 1.0 * total_fix / LOOP_MAX);
}
```
以及main函数：
```cpp
int main()
{
lcg_seed = time(NULL);
getapproxrecip = getapproxrecip8x8;
printf("8x8:\n");
random_test_error();
getapproxrecip = getapproxrecip16x8;
printf("16x8:\n");
random_test_error();
}
```
其中一次运行的结果如下：
```
8x8:
[-0] 6.7992%, [-1]33.7324%, [-2]35.0911%, [-3]17.2700%,
[-4] 6.0442%, [-5] 1.0029%, [-6] 0.0588%, [-7] 0.0014%,
average fix times = 1.853
16x8:
[-0] 8.1363%, [-1]36.1467%, [-2]33.7419%, [-3]15.8187%,
[-4] 5.3604%, [-5] 0.7756%, [-6] 0.0204%, [-7] 0.0000%,
average fix times = 1.765
```
可以看到误差在-1和-2的概率是最大的，其它的概率就非常低。这里给出了平均修正次数的估算，指的是对于随机的情况下，误差需要几次修正才能得到正确的结果。这样来看似乎使用16x8的表可以获得更少的平均修正次数，但是二者差的并不多。
这也说明即使是8x8和16x8的表的差距都不过大约5%，更不用说对于16x8修改一下表的内容能有多少差距，实际上连0.5%都达不到。
因为我们保证了商是偏小的，所以按照公式 $q * b \le a$ 可以得出如下代码：
```cpp
uint32_t div_exact(uint32_t a, uint32_t b)
{
uint32_t q = div_approx(a, b);
uint64_t rem = ((uint64_t)a << 32) - ((uint64_t)q * b);
while (rem >= b) {
q++;
rem -= b;
}
return q;
}
```
这个代码在近似商的基础上修正了结果，可以得到精确的定点小数Q1.31除法的商。按照之前所述，while内语句执行的平均次数小于2，while的判断次数平均小于3，总体开销不大。
## 用于软浮点除法
熟悉 _IEEE 754_ 单精度浮点数的应该知道，其尾数有23位，算上规格化省去的一位，实际上是24位，按照定点小数Q1.23的格式存储。而在计算除法的过程中，需要对Q1.23的定点小数进行除法运算，一种简单的实现方式就是长除法，但是长除法通常比较耗时，即使使用 _Hacker's Delight_ 一书中提到的两次短除法的操作也是比较费事的，而本文所述的牛顿法就可以用在这种情况，获得更快的速度。
部分代码如下：
```cpp
if (sig1 < sig2) { // sig1和sig2是Q1.23格式
exp--; // 这里为了防止溢出做的判断
sig1 <<= 8; // 因为除法需要小数除以大数
} // 因此被除数必须小于除数
else {
sig1 <<= 7;
}
sig2 <<= 8; // 除数扩展到Q1.31
sig = (uint64_t)sig1 * approx_recip(sig2) >> 31;
uint64_t rem = (uint64_t)sig1 << 32;
rem -= (uint64_t)sig * sig2; // r = a - q * b
while (rem >= sig2) { // 修正商，对于Q1.23最多4次
sig++;
rem -= sig2;
}
```
完整代码可以在代码仓库找到，对于软浮点有兴趣可以看我之前有关软浮点的介绍，内容也是在不断更新的。
## 用于整数除法
对于定点数除法，通常使用CPU自带的32位硬件除法指令就可以完成，但是可能存在有些CPU没有硬件除法指令的情况（虽然不多见了，可能一些比较老的没有）。上述定点小数除法也可以用在定点整数除法。
例如对于 a / b，如果 a < b，则返回0，否则执行定点小数除法，求出 1 / b，然后再用 a 去乘这个倒数，不过由于定点数和定点小数的区别，所以需要将 b 扩展到Q1.31的数。
于是需要求数 b 的前导0数量，具体方法如下（来自*Hacker's Delight*）：
```cpp
int count_leading_zero(uint32_t x)
{
if (x == 0)
return 32;
int n = 1;
if ((x >> 16) == 0) { n += 16; x <<= 16; }
if ((x >> 24) == 0) { n += 8; x <<= 8; }
if ((x >> 28) == 0) { n += 4; x <<= 4; }
if ((x >> 30) == 0) { n += 2; x <<= 2; }
n -= (x >> 31);
return n;
}
```
由于 b 扩展到了Q1.31，所以实际上 b 的位宽从 $w_b$ 扩大到了 32，而扩大的位数 $k$ 就是前导0的数量。因此结果需要右移一定的量才能得到实际的宽度。由于 b 的倒数是Q32的格式，所以和 a 相乘后得到的结果的宽度是 $w_a+32$ 位。按照除法的原理，实际上 a / b 的宽度应该是 $w_a - w_b + 1$，而这里得到的比实际的偏大了 $w_b + 31$ 位，所以需要右移这么多，又因为 $w_b = 32-k$，所以移位值就是 $63-k$，所以下面的代码就很好理解了：
```cpp
uint32_t udiv_approx(uint32_t a, uint32_t b)
{
if (a < b)
return 0;
int lz = count_leading_zero(b);
uint32_t recip_b = getapproxrecip8x8(b << lz);
return ((uint64_t)a * recip_b) >> (63 - lz);
}
```
当然由于近似倒数存在最大-3的误差，所以实际结果需要进行修正。于是结果精确的代码如下，适当修改一下参数还可以顺便获得余数：
```cpp
uint32_t udiv(uint32_t a, uint32_t b)
{
if (a < b)
return 0;
int lz = count_leading_zero(b);
uint32_t recip_b = getapproxrecip8x8(b << lz);
uint32_t quo = ((uint64_t)a * recip_b) >> (63 - lz);
uint32_t rem = a - quo * b;
while (rem >= b) {
quo++;
rem -= b;
}
return quo;
}
```
不过对于整数除法还有一个方法，就是按照 Hacker's Delight 第10.16.2节所述的方法求乘法逆元，然后直接相乘就行了。用定点小数再实现整除确实没有必要，只是为了说明这种可行性。
## 参考资料
1. <https://blog.segger.com/algorithms-for-division-part-3-using-multiplication/>
2. <https://blog.segger.com/algorithms-for-division-part-4-using-newtons-method/>
3. Hacker's Delight，论坛搜 **算法心得：高效算法的奥秘** 可以找到 本帖最后由 DEATHTOUCH 于 2023-10-6 00:41 编辑
# 更新内容
## 迭代过程针对部分32位CPU的优化
这一部分主要是针对如下代码在部分32位CPU下的优化：
```cpp
uint32_t getapproxrecip8x8(uint32_t n)
{
uint32_t x, t;
x = initreciptable8x8[(n >> 28) & 7] << 24;
t = ((uint64_t)x * n) >> 32;
t = ~t;
x = ((uint64_t)x * t) >> 31;
t = ((uint64_t)x * n) >> 32;
t = ~t;
x = ((uint64_t)x * t) >> 31;
t = ((uint64_t)x * n) >> 32;
t = ~t;
x = ((uint64_t)x * t) >> 31;
return x;
}
```
仔细观察代码中的 `x = ((uint64_t)x * t) >> 31;` 这一部分，对于大部分32位CPU来说这句代码的执行需要较大的开销。对于乘法，有些32位CPU不一定支持32位相乘得到64位的操作，但是可以使用一定的算法计算高32位结果。除了乘法指令，右移31位一般需要的指令超过两条（除了Intel，刚好是两条），因为一个64位结果是保存在两个不同的寄存器的。而后面要讲的优化代码对于Intel的CPU来说可能效果不大。
考虑到牛顿迭代法的二次收敛性，其实一开始的两次迭代从4位到8位再到16位精度，实际上不需要右移31位来确保精度。所以完全可以在前两次迭代使用高32位（第63位到32位）的结果，然后最后一次结果才使用完整的第62到31位结果（即右移31位）。
改进后的代码如下（误差和原来一样）：
```cpp
uint32_t getapproxrecip8x8_32(uint32_t n)
{
uint32_t x, t;
x = initreciptable8x8[(n >> 28) & 7] << 24;
t = ((uint64_t)x * n) >> 32;
t = -t;
x = ((uint64_t)x * t) >> 32;
x = x << 1;
t = ((uint64_t)x * n) >> 32;
t = -t;
x = ((uint64_t)x * t) >> 32;
x = x << 1;
t = ((uint64_t)x * n) >> 32;
t = ~t;
x = ((uint64_t)x * t) >> 31;
return x;
}
```
这段代码对于编译器来说会生成更短的代码，可以提高执行速度。尽管按照之前的分析使用 `t = -t;` 会造成误差，但是这里用了高32位结果，忽略了1位，所以影响不大，而用 `t = ~t;` 甚至会过于偏小了，导致平均误差偏大。但是最后一次迭代必须确保结果不偏大，所以就需要使用 `t = ~t;` 了。
如果CPU不支持结果是64位的乘法，那么可以用如下代码（不过会导致最大误差达到-4）：
```cpp
uint32_t getapproxrecip8x8_32(uint32_t n)
{
uint32_t x, t;
x = initreciptable8x8[(n >> 28) & 7] << 24;
t = -mulhu(x, n);
x = mulhu(x, t) << 1;
t = -mulhu(x, n);
x = mulhu(x, t) << 1;
t = ~mulhu(x, n);
x = mulhu(x, t) << 1;
return x;
}
```
这段用到的 `mulhu` 函数的作用是计算两个32位无符号整数的乘积，返回乘积的高32位结果。具体代码的实现，可以在 Hacker's Delight 第 8.2 节找到。
上述代码都经过穷举测试，应该是效果比较好的了。
## 计算64位定点小数的近似倒数
要计算Q1.63的近似倒数，只需在Q1.31近似倒数的基础上再进行一次或两次迭代就行了，一种可以用于64位CPU和大多数32位CPU的方法如下：
```cpp
uint64_t getapproxrecip64_1(uint64_t n)
{
uint64_t x, t;
x = (uint64_t)getapproxrecip16x8(n >> 32) << 32;
t = ~mulhu64(x, n);
x = mulhu64(x, t) << 1;
return x;
}
uint64_t getapproxrecip64_2(uint64_t n)
{
uint64_t x, t;
x = (uint64_t)getapproxrecip8x8(n >> 32) << 32;
t = ~mulhu64(x, n);
x = mulhu64(x, t) << 1;
t = ~mulhu64(x, n);
x = mulhu64(x, t) << 1;
return x;
}
```
注意一次迭代为了精度，使用的是 `getapproxrecip16x8` 来获取32位近似倒数，而两次迭代则无所谓了。
其中的 `mulhu64` 的实现如下：
```cpp
uint64_t mulhu64(uint64_t u, uint64_t v)
{
uint64_t w0,w1,w2,u0,u1,v0,v1,t;
u0 = u & 0xffffffff;
u1 = u >> 32;
v0 = v & 0xffffffff;
v1 = v >> 32;
w0 = u0 * v0;
t = u1 * v0 + (w0 >> 32);
w1 = t & 0xffffffff;
w2 = t >> 32;
w1 = u0 * v1 + w1;
return u1 * v1 + w2 + (w1 >> 32);
}
```
这样就能得到精度尚可的Q64小数结果。
为了分析误差，需要计算精确的倒数，这对于x86-64架构来说丝毫没有难度，然而对于大多数CPU来说确存在问题。为了简化操作，就只选择使用x86-64来进行分析了，否则不光实现复杂，耗时也会很长。
以下代码可以获得精确的Q1.63小数的Q64倒数：
```cpp
__attribute__((noinline)) uint64_t exactrecip64(uint64_t n)
{
asm("movq $0xffffffffffffffff, %rax\n");
asm("movq $0x7fffffffffffffff, %rdx\n");
#ifdef _WIN64
asm("divq %rcx\n");
#else
asm("divq %rdi\n");
#endif
}
```
之后用随机数进行误差分析，可以得到大致的误差范围。
一次迭代的倒数误差在 $[-12,0]$ 之间，两次迭代则缩小到了 $[-3,0]$ 之间。
## 计算64位定点小数除法
在得到倒数以后，显然可以去计算定点小数的除法了，为了和精确结果比较，所以需要计算精确的结果，和之前一样，用x86-64汇编来实现：
```cpp
__attribute__ ((noinline,sysv_abi))
uint64_t exactdiv64(uint64_t u, uint64_t v)
{
asm("movq %rdi, %rdx\n");
asm("xorq %rax, %rax\n");
asm("divq %rsi\n");
}
```
为了不针对不同平台写两份汇编，就用了 `sysv_abi` 这个调用约定，这样前两个参数就是 `rdi` 和 `rsi` 了。
为了节约篇幅，下面就直接贴出计算除法以及修正商的代码：
```cpp
uint64_t recipdiv64(uint64_t u, uint64_t v)
{
uint64_t r1, r0, q;
q = mulhu64(u, getapproxrecip64_1(v)) << 1;
r1 = u - mulhu64(q, v) - 1;
r0 = -(q * v);
while ((r1 != 0) || (r0 >= v)) {
q++;
if (r0 < v)
r1--;
r0 -= v;
}
return q;
}
```
由于没有128位结果的乘法以及128位的加减法，所以需要拆开来处理，这个实际上不难实现。上述代码中变量 `r1` 在32位CPU实现可以用 `uint32_t` 类型，其它就没什么了。
有了64位定点小数的实现，实际上实现一个 `double` 类型浮点数的模拟也就不是什么难事了。 说到牛顿迭代法，我就想到那个神奇的「快速求解平方根倒数算法」：
```c
float q_rsqrt(float number)
{
long i;
float x2, y;
const float threehalfs = 1.5F;
x2 = number * 0.5F;
y= number;
i= * ( long * ) &y; // evil floating point bit level hacking
i= 0x5f3759df - ( i >> 1 ); // what the fuck?
y= * ( float * ) &i;
y= y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
// y= y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed
return y;
}
```
爱飞的猫 发表于 2023-10-2 05:16
说到牛顿迭代法，我就想到那个神奇的「快速求解平方根倒数算法」：
```c
没看错的话这是john carmak的神作吧{:301_997:} 爱飞的猫 发表于 2023-10-2 05:16
说到牛顿迭代法，我就想到那个神奇的「快速求解平方根倒数算法」：
```c
这个实在是太经典了，网上资料铺天盖地，算是大部分人提到牛顿法一下就能想到的，相比之下我这个定点数的牛顿法就相当少见了。 hrh123 发表于 2023-10-2 09:24
没看错的话这是john carmak的神作吧
这个代码最早就是Quake3的源码里发现的，不过据研究这个方法在之前就有人内部使用这个代码了。这在网上可以找到许多，Hacker's Delight也有所涉及。 DEATHTOUCH 发表于 2023-10-2 12:13
这个代码最早就是Quake3的源码里发现的，不过据研究这个方法在之前就有人内部使用这个代码了。这在网上可 ...
不过这个magic number应该不是最优的吧 hrh123 发表于 2023-10-5 15:11
不过这个magic number应该不是最优的吧
是的，后来Lomont博士研究出来了更好的0x5f375a86，还有人研究出了double类型的0x5fe6ec85e7de30da。 DEATHTOUCH 发表于 2023-10-5 15:37
是的，后来Lomont博士研究出来了更好的0x5f375a86，还有人研究出了double类型的0x5fe6ec85e7de30da。
好像说最佳值是0x5f37642f,只不过在迭代次数较少时表现不佳
双精度的话0x5fe6ec85e7de30da也是lomont博士搞出来的,但后来的研究表明，代入0x5fe6eb50c7aa19f9的结果精确度更高,也有说0x5FE6EB50C7B537AA和0x5FE6EB50C7B537A9的 本帖最后由 DEATHTOUCH 于 2023-10-5 15:58 编辑
hrh123 发表于 2023-10-5 15:43
好像说最佳值是0x5f37642f,只不过在迭代次数较少时表现不佳
双精度的话0x5fe6ec85e7de30da也是lomont博 ...
按照书中所说，其实只要穷举测试每一个可能的结果，然后判断误差就行了。
这样的话应该就可以针对不同的迭代次数需求找出各自的最优解。