C语言除法算法和取模运算的实现(多种算法,多种思路)

对计算机来说,除法与求模是整数算术运算中最复杂的运算。相对其他运算(如加法与减法)来说,这两种算法的执行速度非常慢。例如,ARM 硬件上不支持除法指令,编译器调用 C 库函数来实现除法运算。直接利用 C 库函数中的标准整数除法程序要花费 20~100 个周期,消耗较多资源。

在非嵌入式领域,因为 CPU 运算速度快、存储器容量大,所以执行除法运算和求模运算消耗的这些资源对计算机来说不算什么。但是在嵌入式领域,消耗大量资源带来的影响不言而喻。因此,从理论上讲,我们应该在程序表达式中尽量减少对除法运算与求模运算的使用,尽量使用其他方法来代替除法与求模运算。例如,对于下面的示例代码:
if (x/y>z)
{
    // ...
}
我们可以将其修改成如下形式:
if (((y>0)&&(x>y*z))||((y<0)&&(x<y*z)))
{
    // ...
}
这样就简单地避免了一些除法运算。同时,也可以在表达式中通过合并除法的方式来减少除法运算,下面通过示例来讲解。对于如下代码:
double x=a/b/c;
double y=a/b+c/b;
根据数学结合原则,上面的代码可以通过合并的方式减少代码中的除法运算,修改后的代码如下:
double x=a/(b*c);
double y=(a+c)/b;
同样,对于求模运算,也可以采用相应的方法来代替,如下面的示例代码:

a=a%8;

可以修改为:

a=a&7;

对于下面的表达式:

x=(x+y)%z;

可以通过如下方式来避免使用模操作符:
x+=y;
while(x>=z)
{
    x-=z;
}
通过上面的阐述,相信大家对如何减少使用除法与模运算有了初步了解。下面将详细讨论如何优化除法运算与求模运算。

用倒数相乘来实现除法运算

何为倒数相乘?其实很简单,它的核心思想就是利用乘法来代替实现除法运算。例如,在 IA-32 处理器中,乘法指令的运算速度比除法指令要快 4~6 倍。因此,在某些情况下尽量使用乘法指令来代替除法指令。

那么,我们该如何利用乘法来代替实现除法运算呢?原理就是被除数乘以除数的倒数,用公式表现为:

x/y=x*(1/y)

例如,计算 10/5,可以根据公式 x/y=x*(1/y) 这样来计算:

10/5=10*(1/5)=10*0.2=2

在实际应用中,一些编译器也正是基于这个原理才得以将除法运算转换为乘法运算的。现在我们来看一个除法运算示例:
#include <stdio.h>
int main(void)
{
    int x = 3/2;
    float y = 3.0/2.0;
    printf("3/2 = %d\r\n3.0/2.0 = %1.1f\n",x,y);
    return 0;
}
运算结果为:
3/2 = 1
3.0/2.0 = 1.5

通过该除法运算示例可以看出,很明显没能充分考虑到浮点类型。另外,在 C 语言中,一般情况下 1 除以任何数其结果皆为 0。那么怎样才能解决这个问题呢?编译器采用了一种称为“定点运算”(fixed-point arithmetic)的方法。

那么何为定点运算,定点运算有什么特点呢?

前面已经阐述过,由于计算机表示实数时为了在固定位数内能表示尽量精确的实数值,分配给表示小数部分的位数并不是固定的,也就是说“小数点是浮动的”,因此计算机表示的实数数据类型也称为浮点数。

相对于“小数点是浮动的”来讲,定点运算根据字面意思来理解就是“小数点是固定的”。有了定点运算,表示小数时不再用阶码(exponent component,即小数点在浮点数据类型中的位置),而是要保持小数点的位置固定不变。这和硬件浮点数机制截然不同,硬件浮点数机制是由硬件负责向整数部分和小数部分分配可用的位数。有了这种机制,浮点数就可以表示很大范围的数——从极小的数(在 0~1 的实数)到极大的数(在小数点前有数十个 0)。这种小数的定点表示法有很多优点,尤其能极大地提高效率。当然,作为代价,同样也必须承受随之而来的精度上的损失。

对于定点数表示法(fixed-point),相信大家并不陌生。所谓定点格式,即约定机器中所有数据的小数点位置是固定不变的。在计算机中通常采用两种简单的约定:将小数点的位置固定在数据的最高位之前(即定点小数),或者固定在最低位之后(即定点整数)。

其中,定点小数是纯小数,约定的小数点位置在符号位之后、有效数值部分的最高位之前。若数据 x 的形式为 x=x0x1x2…xn(其中 x0 为符号位,x1,…,xn 是数值的有效部分,也称为尾数,x1 为最高有效位),则在计算机中的表示形式为:


一般说来,如果最末位 xn=1,前面各位都为 0,则数的绝对值最小,即 |x|min=2-n;如果各位均为 1,则数的绝对值最大,即 |x|max=1-2-n。因此定点小数的表示范围是:


定点整数是纯整数,约定的小数点位置在有效数值部分最低位之后。若数据 x 的形式为 x=x0x1x2…xn(其中 x0 为符号位,x1,…,xn 是尾数,xn 为最低有效位),则在计算机中的表示形式为:



由此可知,定点整数的表示范围是:


当数据小于定点数能表示的最小值时,计算机将它作 0 处理,称为下溢;当数据大于定点数能表示的最大值时,计算机将无法表示,称为上溢,上溢和下溢统称为溢出

当计算机采用定点数表示时,对于既有整数又有小数的原始数据,需要设定一个比例因子,数据按该比例缩小成定点小数或扩大成定点整数再参加运算。在运算结果中,根据比例因子,将数据还原成实际数值。若比例因子选择不当,往往会使运算结果产生溢出或降低数据的有效精度。

使用牛顿迭代法求除数的倒数

在上一小节,我们阐述了如何使用倒数相乘(x/y=x*(1/y))的方法来实现除法运算。然而,对于如何能够快速有效地取倒数,牛顿迭代法(Newton’s method)是最佳方案。

对于牛顿迭代法,相信学过高等数学的读者并不陌生,它又称为牛顿-拉夫逊方法(Newton-Raphson method),是牛顿在 17 世纪提出的一种在实数域和复数域上近似求解方程的方法,它将非线性方程线性化,从而得到迭代序列的一种方法。

对于方程 f(x)=0,设 x0 为它的一个近似根,则函数 f(x) 在 x0 附近截断高次项可用一阶泰勒多项式展开为如下形式:


这样,由式(1)我们可以将 f(x)=0 转化为如下形式:


 
在这里,我们设 f′(x)≠0,则有:


 
取 x 作为原方程新的近似根 x1,再代入方程,如此反复,于是就产生了迭代公式:


 
有了迭代公式(4)之后,现在我们继续来看如何用牛顿迭代公式来求倒数,即求除数 a 的倒数 1/a。

这里我们设:


 
式中 x 为 a 的倒数,方程 f(x)=0 为一非线性方程。现在把 f(x)=0 代入牛顿迭代序列式(4)中,就可以得出求倒数的公式,如下所示:


 
在式(5)中,xn 为第 n 次迭代的近似根。

如式(5)所示,用牛顿迭代法求倒数,每次迭代需要一次减法与两次乘法,所用的迭代次数决定最终的计算速度和精度。迭代次数越多,则精度越高。但迭代次数越多,速度也越慢,因此实际运用时应综合考虑速度和精度两方面的因素,选择合适的迭代次数。

其实,牛顿迭代法在程序中应用得非常广泛,如最常用的开方、开方求倒数等。在 QuakeⅢ 源码中,在 game/code/q_math.c 文件中就有一个函数 Q_rsqrt,它的作用是将一个数开平方后取倒,其运行效率也非常高。其函数实现为:
float Q_rsqrt(float number)
{
    long i;
    float x2, y;
    const float threehalfs = 1.5F;
    x2 = number * 0.5F;
    y  = number;
    i  = * (long * ) &y;
    i  = 0x5f3759df - (i >> 1);
    y  = * ( float * ) &i;
    // 第一次迭代
    y  = y * ( threehalfs - ( x2 * y * y ));
    // 第二迭代
    // y  = y * ( threehalfs - ( x2 * y * y ) );
    return y;
}
从代码中可以看出,程序首先猜测出一个接近 1.0/sqrt(number) 的近似值,然后两次使用牛顿迭代法进行迭代(实际只需要使用一次)。这里需要特别注意的是 0x5f3759df 这个值,因为通过执行语句“0x5f3759df-(i>>1)”,得出的值出人意料地接近 1/sqrt(number) 的值,因此,我们只需要一次迭代就可以求得近似解,或许这就是数学的神奇。

用减法运算来实现整数除法运算

我们知道,减法运算比除法运算要快得多。因此,对整数除法运算来说,如果知道被除数是除数很小的倍数,那么可以使用减法运算来代替除法运算。例如,对于下面的示例代码:
unsigned int x=300;
unsigned int y=100;
unsigned int z=x/y;
我们可以将“z=x/y”表达式修改成如下形式:
unsigned int x=300;
unsigned int y=100;
unsigned int z=0;
while (x>=y)
{
    x-=y;
    ++z;
}
这里使用减法来代替除法运算,虽然代码看起来不是很直观,但是在运行效率上确实要快许多。当然,具体效率也要取决于被除数与除数的倍数。如果倍数比较大,那么相应的循环次数就会增多,采取这种方法就得不偿失了。

用移位运算实现乘除法运算

用移位运算来实现乘除法运算的方法,相信大家并不陌生,实际上有很多 C 编译器都能够自动地做好这个优化。通常,如果需要乘以或除以 2n,都可以用移位的方法代替。例如:
a=a*2;
b=b/2;
可以修改为如下形式:
a=a<<1;
b=b>>1;
其中,除以 2 等价于右移 1 位,乘以 2 等价于左移 1 位。同理,除以 4 等价于右移 2 位,乘以 4 等价于左移 2 位;除以 8 等价于右移 3 位,乘以 8 等价于左移 3 位,以此类推。

其实,利用上面的原理,只要是乘以或除以一个整数,均可以用移位运算的方法来得到结果,例如:
a=a*5;
可以将其分解为 a*(4+1),即 a*4+a*1。由此,我们就可以很简单地得到下面的程序表达式:

a=(a<<2)+a

尽量将浮点除法转化为相应的整数除法运算

有时候,如果不能够在代码中避免除法运算,那么尽量使除数和被除数是无符号类型的整数。实际上,有符号的除法运算执行起来比无符号的除法运算更加慢,因为有符号的除法运算要先取得除数和被除数的绝对值,再调用无符号除法运算,最后再确定结果的符号。

同时,对于浮点除法运算,可以先将浮点除法运算转化为相应的整数除法运算,最后对结果进行相应处理。例如,可以将浮点除法运算的分子和分母同时放大相同的倍数,就可以将浮点除法运算转换成相同功能的整数除法运算。