牛顿迭代法求近似倒数及其应用
这是继上篇软浮点（传送门：软件模拟实现IEEE-754单精度浮点数运算）后进一步的研究，主要是为了找到除法运算一种更快的方法。
全部的代码依然在代码仓库 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$ 的倒数了，不过对于初值的估计是一个问题。因为如果估计偏差过大，可能会导致误差不断变大而无法在较短的步数内收敛。同时，这个式子对于浮点数来说实现起来完全没什么问题，不过如果要定点化就稍微需要改一下了。
浮点数实现
对着上述公式，很容易就可以写出浮点数的版本：
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位定点小数规范如下：
回到牛顿迭代上来，我们可以把迭代的每一步拆开来：
t = a * x;
t = 2 - t;
x = x * t;
-
对于 t = a * x，可以确定变量 a 是Q31的，而 x 必然是Q32的，那么其结果是Q1.63，所以实际上可以取其高32位，而舍弃低32位，从而变成Q1.31。
-
对于 t = 2 - t，由于此时变量 t 是Q1.31的，而 2 可以理解为是Q2.31，那么可以确定的是 2 - t 之后的结果依然是Q1.31的。
-
对于 x = x * t，显然 x 是Q32，t 是Q1.31，相乘是Q1.63，但是由于下一轮迭代的 x 需要Q32，所以可以使结果右移31位，获取Q1.63的Q63部分。
所以写成实际的代码应该是这样的：
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位从表中获取初值，具体方法如下：
x = initial_table_8x8[(a >> 28) & 7] << 24;
如果是估算4位的数据，那么就需要获取第30-27位，方法如下：
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}$，这样对于各种情况，倒数的估计会更加准确。
具体生成表的代码如下：
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字节。至于表是不是越大越好，后面会进行分析。
完整代码
这里给出实际可以使用的代码：
const uint8_t initreciptable8x8[8] = {
0xf1, 0xd8, 0xc3, 0xb2, 0xa4, 0x98, 0x8d, 0x84
};
const uint8_t initreciptable16x8[16] = {
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[n >> 27 & 15] << 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}$ 个，完全可以使用穷举测试。
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("[eq] %d, [maxpos] %d, [maxneg] %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:
[eq] 874319370, [maxpos] 0, [maxneg] -3
[-1] 1126936446, [-2] 145419076, [-3] 808756
16x8:
[eq] 971865634, [maxpos] 0, [maxneg] -3
[-1] 1050310821, [-2] 125307190, [-3] 3
可以看到8x8的表已经非常不错了，误差为-3的概率只有0.038%。而16x8的表主要进一步减少了-3的个数，不过居然有3个漏网之鱼，这就让我好奇这3个数究竟是什么。然后就用暴力法再搜一下：
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:
[eq] 970775900, [maxpos] 0, [maxneg] -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。
用代码来表示就是这样：
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位可以得到所需的近似结果。
实现的代码如下：
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默认使用的随机数算法，因为简单所以拿来用了。
uint32_t lcg_seed;
uint32_t lcg_rand()
{
lcg_seed = lcg_seed * 134775813 + 1;
return lcg_seed;
}
至于随机数的初值，在main函数的开头用 lcg_seed = time(NULL); 就可以了。
在之前的 div_approx() 函数基础上，再使用如下代码进行一次除法的误差分析：
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;
}
随机测试代码如下：
void random_test_error()
{
const int LOOP_MAX = 10000000;
int n[8] = {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[diff]++;
}
for (int i=0; i<8; i++) {
printf("[-%d]%7.4f%%, ", i, 100.0 * n[i] / 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[i];
printf("average fix times = %.3f\n", 1.0 * total_fix / LOOP_MAX);
}
以及main函数：
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$ 可以得出如下代码：
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 一书中提到的两次短除法的操作也是比较费事的，而本文所述的牛顿法就可以用在这种情况，获得更快的速度。
部分代码如下：
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）：
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$，所以下面的代码就很好理解了：
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的误差，所以实际结果需要进行修正。于是结果精确的代码如下，适当修改一下参数还可以顺便获得余数：
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节所述的方法求乘法逆元，然后直接相乘就行了。用定点小数再实现整除确实没有必要，只是为了说明这种可行性。
参考资料
-
https://blog.segger.com/algorithms-for-division-part-3-using-multiplication/
-
https://blog.segger.com/algorithms-for-division-part-4-using-newtons-method/
-
Hacker's Delight，论坛搜 算法心得：高效算法的奥秘 可以找到
[复制链接]
发表于 2023-10-2 00:14
|
发表于 2023-10-4 23:47
发表于 2023-10-2 05:16
