Kuekua's blog

YesterDay you said tomorrow!

平方根运算的软件与硬件的加速计算

• 分类: Algorithm Optimization • 标签: Algorithm, Algorithm Optimization

1. 平方根运算软件算法

1.1 二分法

利用二分进行开平方的思想很简单：假定中值为最终解。假定下限为0，上限为x，然后求中值；然后比较中值的平方和x的大小，并根据大小修改下限或者上限；重新计算中值，开始新的循环，直到前后两次中值的距离小于给定的精度为止。**需要注意的一点是，如果x小于1，我们需要将上限置为1。

代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
float SqrtByBisection(float n)
{
	float low,up,mid,last; 
	low=0,up=(n<1?1:n); 
	mid=(low+up)/2; 
	do
	{
		if(mid*mid>n)
			up=mid; 
		else 
			low=mid;
		last=mid;
		mid=(up+low)/2; 
	}while(fabsf(mid-last) > eps);  //eps为精度，一般为0.000001
	return mid; 
}

这里有一点需要特别注意：在精度判别时不能利用上下限而要利用前后两次mid值，否则可能会陷入死循环！这是因为由于精度问题，在循环过程中可能会产生mid值和up或low中的一个相同。这种情况下，后面的计算都不会再改变mid值，因而在达不到精度内时就陷入死循环。

1.2 牛顿迭代法

将中值替换为切线方程的零根作为最终解，原理可以利用下图解释：

这里写图片描述

具体代码：

1
2
3
4
5
6
7
8
9
10
11
12
float SqrtByNewton(float x)
{
	int temp = (((*(int *)&x)&0xff7fffff)>>1)+(64<<23);
	float val=*(float*)&temp;       //初始值
	float last;
	do
	{
		last = val;
		val =(val + x/val) / 2;     //利用一阶梯度更新结果
	}while(fabsf(val-last) > eps);  //eps为精度，一般为0.000001
	return val;
}

注意到代码中初始值的选择，解释如下：
IEEE浮点标准用的形式来表示一个数，将该数存入float类型之后变为：

现在需要对这个浮点数进行开方，我们看看各部分都会大致发生什么变化。指数E肯定会除以2,127保持不变，m需要进行开方。由于指数部分是浮点数的大头，所以对指数的修改最容易使初始值接近精确值。幸运的是，对指数的开平方我们只需要除以2即可，也即右移一位。但是由于E+127可能是奇数，右移一位会修改指数，我们将先将指数的最低位清零，这就是& 0xff7fffff的目的。然后将该转换后的整数右移一位，也即将指数除以2，同时尾数也除以2（其实只是尾数的小数部分除以2）。由于右移也会将127除以2，所以我们还需要补偿一个64，这就是最后还需要加一个(64<<23)的原因。

这里大家可能会有疑问，最后为什么加(64<<23)而不是(63<<23)，还有能不能不将指数最后一位清零？答案是都可以，但是速度都没有我上面写的快。这说明我上面的估计更接近精确值。下面简单分析一下原因。首先假设e为偶数，不妨设e=2n，开方之后e则应该变为n,127保持不变，我们看看上述代码会变为啥。e+127是奇数，会清零，这等价于e+126，右移一位变为n+63，加上补偿的64，指数为n+127，正是所需！再假设e为奇数，不妨设e=2n+1，开方之后e应该变为n+1（不精确），127保持不变，我们看看上述代码会变为啥。e+127是偶数等于2n+128，右移一位变为n+64，加上补偿的64，指数为n+1+127，也是所需！这确实说明上述的估计比其他方法更精确一些，因而速度也更快一些。

1.3 卡马克算法—快速平方根倒数算法

此法实质上是进行了1~2次的牛顿迭代法求开方倒数，求平方根的倒数，实际就是求方程$$ 1/(x^2)-a=0 $$的解。将该方程按牛顿迭代法的公式展开为：$$ x_{(n+1)}=1/2x_n(3-ax_nx_n) $$。这个方法的一个关键地方还在于起始值的选取，算法选取一个“magic number”，使算法的初始解非常接近精确解！

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
float SqrtByCarmack( float number )  
{  
    int i;  
    float x2, y;  
    const float threehalfs = 1.5F;  
  
    x2 = number * 0.5F;  
    y  = number;  
    i  = * ( long * ) &y;       
    i  = 0x5f375a86 - ( i >> 1 );   //得到初始化值
    y  = * ( float * ) &i;  
    //第一次迭代，可以根据精度要求重复此过程，一般一次迭代精度就够用了
    y  = y * ( threehalfs - ( x2 * y * y ) );   
    return number*y;  
}

2. 整数平方根———硬件加速算法

2.1 手工开平方算法法

先以10进制为例，解释手动开平方算法：

首先将数字每两位分成一段。如：745836942。就分成：
7|45|83|69|42，即从个位开始分。共分成五段，第一段是7。对第一段的数字7开方取整，可得a=2。此时，要在2后面接一个数字b，并在7后面加上下一段的数45，使产生的两位数ab的平方不大于745。
我们知道，数ab实际值表示为10a+b，其平方是$$100a^2+20ab+b^2$$。我们可以暂时忽略$$b^2$$，而产生一个“试商”b。即$$b = (745 - 100a^2) / 20a = (745 - 10022) / (20 * 2) = 8.625 ≈ 8$$(向下取整)。但是，我们发现$$28^2 = 784>745$$，于是这个试商需要减少为7。然而，当a=0时，上述求试商的方法不在适用，但我们可以直接取下一段的两位数开方。如√45≈6。求出试商后，用$$745-ab^2$$得到新的“第一段”的数。

具体过程：

取出第一段的数mn，开方得到a，然后接上第二段的数pq，用$$mnpq-100a^2$$得到“余数”x， $$x/20a$$得到试商b，然后调整b（当$$ 20ab+bb>x $$时b需要减少），调整后，将 $$x-20ab-b*b$$作为新的余数x’，ab就是第二轮开平方的结果。由于前面已经将$$100a^2$$减掉，所以后面每次都不用再减去$$100a^2$$。重复步骤，直到开方完毕，或达到要求的精度为止。最后得到的a就是平方根。

如求$$745836942$$的平方根：
$$7|45|83|69|42$$
$$①a=sqrt(7)≈2，b=(745-400)/40≈8，28^2=784>745 ==> b=7 ==> 27^2=729<745，745-729 = 16$$

$$②a=27，b=1683/540≈3，27203+3*3=1629<1683，1683-1629 = 54$$

$$③a=273，b=5469/5460≈1，273201+1*1=5461<5469 ==""> b=1， 5469-5461=8$$

当运用到二进制时，数ab实际值表示为a<<1+b，其平方是$$a^2<<2+ab<<2+b^2$$。我们可以暂时忽略$$b^2$$，而产生一个“试商”b，依此类推下去。以下代码为求一个32位数的平方根，算法将原始值两两分组进行计算，注意：根应至多为16位且每次计算的b值只能是0或1：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
unsigned int sqrt_16_1(unsigned int M) 
{
	unsigned int N,N2, i,j; 
	unsigned int tmp; // 每轮计算的残差 
	if (M == 0)      // 被开方数为0，开方结果也为0 
		return 0;
	N = 0;
	N2 = 0;
	tmp = 0;
	for (i=16; i>0; i-=1) //结果为16位 
	{ 
		N <<= 1;        // 左移一位 
		N2 = N << 1;   //N2代表a<<2
		tmp <<= 2;   // tmp储存的是每次计算完剩下的残差
		tmp += ((M&0xc0000000) >> 30);   //再加上此次计算的2bit的值
		M <<= 2; 
		if( tmp >= N2 + 1 )              //试值ab<<2+bi^2与残差的比较，此时设bi=1
		{
			tmp -= N2 + 1;               //计算此轮的残差
			N++;                         //确定bi=1
		}
	} 
	return N; 
}

将原始值分类大小由2变成4，那么每次要计算2位二进制数值，可能情况有00,01,10,11，00情况时由于不更新结果，直接移动两位即可，所以可不进行处理，内部循环因此为1~3！
也可把分类大小变成8,16…只要是2的倍数均可，算法流程类似！

具体代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
unsigned int sqrt_16_2(unsigned int M) 
{
	unsigned int N,N2, i,j; 
	if (M == 0)                 // 被开方数，开方结果也为0 
		return 0;
	N = 0;
	N2 = 0;
	tmp = 0;
	for (i=16; i>0; i-=2)        // 求剩余的15位，每次2位 
	{ 
		N <<= 2;                 // 左移两位 
		N2 = N << 1;             // N2代表2a
		tmp <<= 4;               // tmp储存的是每次查表完剩下的残差
		tmp += (M >> 28);        // 再加上下一个4bit的值
		for (j=1;j<4;j++)
		 { 
			 ttp[j] = N2*j + j*j; //求2ab+bi^2，bi可为01,10,11
		}
		M <<= 4;
		for (j=3;j>0;j--)
		{
			if (tmp >= ttp[j])     //试值ab<<2+bi^2与残差的比较，bi可为01,10,11
			{
				tmp -= ttp[j]; 
				N+=j;
				break;            //bi只可取一个值
			}
		}
	} 
	return N; 
 }//32->16
© 2017 Kuekua Wu
博客全站共85.5k字