作者论坛账号：DEATHTOUCH

这是继上篇软浮点（传送门：软件模拟实现IEEE-754单精度浮点数运算 https://www.52pojie.cn/thread-1830228-1-1.html）后进一步的研究，主要是为了找到除法运算一种更快的方法。

全部的代码依然在代码仓库 https://gitee.com/peazomboss/softfloat32 可以找到，放在了文件夹 approxrecip 里面，同时软浮点除法部分的代码也已经更新。

公式简介

首先看一下牛顿迭代公式：

��+1=��−�(��)�′(��)

对于求一个数 � 的倒数，我们可以用公式 �(�)=1�2−� 计算，令 �(�)=0，其根就是 � 的倒数。将该式子带入上述迭代公式，化简得到：

��+1=��(2−���)

用这个式子不断迭代就可以得到数 � 的倒数了，不过对于初值的估计是一个问题。因为如果估计偏差过大，可能会导致误差不断变大而无法在较短的步数内收敛。同时，这个式子对于浮点数来说实现起来完全没什么问题，不过如果要定点化就稍微需要改一下了。

浮点数实现

对着上述公式，很容易就可以写出浮点数的版本：

 复制代码 隐藏代码
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取不到。

回到牛顿迭代上来，我们可以把迭代的每一步拆开来：

 复制代码 隐藏代码
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部分。

所以写成实际的代码应该是这样的：

 复制代码 隐藏代码
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...，这样做减法实际上就是取负数，或者说是补码。

不过事情没有那么简单，因为如果简单取负数，会发现一个奇怪的问题，即迭代后的估计值偏大。这很不正常，因为根据牛顿法的原理，对于 �(�)=1�2−� 这样的函数来说，其在坐标上画出来是一个凹函数，所以最终的估计值理论上是偏小的。

但是仔细分析一下就会发现，这是由于截断造成的误差。在代码中我们求出了 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,98)，则倒数取值是 (89,1)。但是如果按照取值是1，那么其倒数就只能是近似于0.9999...了；而取值接近 98 的时候，这个倒数就会有较大的误差。所以可以取一个折中的方案，也就是 1716，这样对于各种情况，倒数的估计会更加准确。

具体生成表的代码如下：

 复制代码 隐藏代码
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的倒数来说，可以使用随机数覆盖测试，不过由于数据量不大，只有 231 个，完全可以使用穷举测试。

 复制代码 隐藏代码
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的。因此计算 �/�，就可以计算 �∗1�，由于 � 和 � 是Q1.31格式，那么 1� 自然就是Q32了。

对于Q1.31和Q32的数相乘，显然结果是Q1.63，由于假定了 � 都是小于 � 的，所以这个结果是小于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%都达不到。

因为我们保证了商是偏小的，所以按照公式 �∗�≤� 可以得出如下代码：

 复制代码 隐藏代码
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 的位宽从 �� 扩大到了 32，而扩大的位数 � 就是前导0的数量。因此结果需要右移一定的量才能得到实际的宽度。由于 b 的倒数是Q32的格式，所以和 a 相乘后得到的结果的宽度是 ��+32 位。按照除法的原理，实际上 a / b 的宽度应该是 ��−��+1，而这里得到的比实际的偏大了 ��+31 位，所以需要右移这么多，又因为 ��=32−�，所以移位值就是 63−�，所以下面的代码就很好理解了：

 复制代码 隐藏代码
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;
}

迭代过程针对部分32位CPU的优化

这一部分主要是针对如下代码在部分32位CPU下的优化：

 复制代码 隐藏代码
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位）。

改进后的代码如下（误差和原来一样）：

 复制代码 隐藏代码
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）：

 复制代码 隐藏代码
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的方法如下：

 复制代码 隐藏代码
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 的实现如下：

 复制代码 隐藏代码
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倒数：

 复制代码 隐藏代码
__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汇编来实现：

 复制代码 隐藏代码
__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 了。

为了节约篇幅，下面就直接贴出计算除法以及修正商的代码：

 复制代码 隐藏代码
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 类型浮点数的模拟也就不是什么难事了。

参考资料

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，论坛搜 算法心得：高效算法的奥秘 可以找到

-官方论坛

www.52pojie.cn