牛顿迭代法求近似倒数及其应用
2023-10-7 15:41:40 Author: mp.weixin.qq.com(查看原文) 阅读量:0 收藏

作者论坛账号:DEATHTOUCH

由于本文使用了大量的数学公式,公众号无法正常显示,建议点击左下角原文进入论坛阅读。

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

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

公式简介

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

+1=()()

对于求一个数  的倒数,我们可以用公式 ()=12 计算,令 ()=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...,这样做减法实际上就是取负数,或者说是补码。

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

但是仔细分析一下就会发现,这是由于截断造成的误差。在代码中我们求出了 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


文章来源: https://mp.weixin.qq.com/s?__biz=MjM5Mjc3MDM2Mw==&mid=2651139766&idx=1&sn=0c43c9cb4d80e57a6d010489fdab30e6&chksm=bd50bee28a2737f40fec5cc6cfb7111d8a6ce1a91063fb9bbaffd9e3dbcebdd76e36170a58c0&scene=58&subscene=0#rd
如有侵权请联系:admin#unsafe.sh