一、对整型数据快速开根号快速算法
应用场合：嵌入式某些需要开根号情况，同时被开根数为整型，比如AD值，由于嵌入式对实时性较高，如果对精确度不是很高，比如精确到对100A电流只需要精确到0.01A时，可以使用以下方法。
适用被开根数：32位 整型
提高精度方法：先放大，计算出整型基础考虑定标位置 sqrt(NUM)=sqrt(NUM*10000)*0.01，当 NUM*10000<0xffff ffff
下面给出的方面本人在DSP2808上测试了，时间效果不错，只需要20%的原使用库sqrt函数执行时间。
===================================================================
这一部分转载自http://blog.sina.com.cn/s/blog_6e350d880100okvn.html
作者：未知
因为工作的需要，要在单片机上实现开根号的操作。目前开平方的方法大部分是用牛顿
迭代法。我在查了一些资料以后找到了一个比牛顿迭代法更加快速的方法。不敢独享，介
绍给大家，希望会有些帮助。
1.原理
因为排版的原因，用pow(X,Y)表示X的Y次幂，用B[0]，B[1]，...，B[m-1]表示一个序列，
其中[x]为下标。
假设：
B[x],b[x]都是二进制序列,取值0或1。
M = B[m-1]*pow(2,m-1) + B[m-2]*pow(2,m-2) + ... + B[1]*pow(2,1) +B[0]*pow(2,0)
N = b[n-1]*pow(2,n-1) + b[n-2]*pow(2,n-2) + ... + b[1]*pow(2,1) +n[0]*pow(2,0)
pow(N,2) = M
(1) N的最高位b[n-1]可以根据M的最高位B[m-1]直接求得。
设 m 已知,因为 pow(2, m-1) <= M <= pow(2, m)，所以 pow(2, (m-1)/2) <= N <=pow(2, m/2)
如果 m 是奇数，设m=2*k+1,
那么 pow(2,k) <= N < pow(2, 1/2+k) < pow(2, k+1),n-1=k, n=k+1=(m+1)/2
如果 m 是偶数，设m=2k,
那么 pow(2,k) > N >= pow(2, k-1/2) > pow(2, k-1),n-1=k-1,n=k=m/2
所以b[n-1]完全由B[m-1]决定。
余数 M[1] = M - b[n-1]*pow(2, 2*n-2)
(2) N的次高位b[n-2]可以采用试探法来确定。
因为b[n-1]=1，假设b[n-2]=1，则 pow(b[n-1]*pow(2,n-1) + b[n-1]*pow(2,n-2) ,2)
= b[n-1]*pow(2,2*n-2) + (b[n-1]*pow(2,2*n-2) + b[n-2]*pow(2,2*n-4)),
然后比较余数M[1]是否大于等于 (pow(2,2)*b[n-1] + b[n-2]) * pow(2,2*n-4)。这种比较只须根据B[m-1]、B[m-2]、...、B[2*n-4]便可做出判断，其余低位不做比较。
若 M[1] >= (pow(2,2)*b[n-1] + b[n-2]) * pow(2,2*n-4), 则假设有效，b[n-2] =1；
余数 M[2] = M[1] - pow(pow(2,n-1)*b[n-1] + pow(2,n-2)*b[n-2], 2) = M[1] -(pow(2,2)+1)*pow(2,2*n-4)；
若 M[1] < (pow(2,2)*b[n-1] + b[n-2]) * pow(2,2*n-4), 则假设无效，b[n-2] =0；余数 M[2] = M[1]。
(3) 同理，可以从高位到低位逐位求出M的平方根N的各位。
使用这种算法计算32位数的平方根时最多只须比较16次，而且每次比较时不必把M的各位逐
一比较，尤其是开始时比较的位数很少，所以消耗的时间远低于牛顿迭代法。
2. 流程图
(制作中，稍候再上)
3. 实现代码
这里给出实现32位无符号整数开方得到16位无符号整数的C语言代码。
----------------------------------------------------------------------------------------------------------------------
-
unsigned int sqrt_16(unsigned long M)
-
{
-
unsigned int N, i;
-
unsigned long tmp, ttp; // 结果、循环计数
-
if (M == 0) // 被开方数，开方结果也为0
-
return 0;
-
N = 0;
-
tmp = (M >> 30); // 获取最高位：B[m-1]
-
M <<= 2;
-
if (tmp > 1) // 最高位为1
-
{
-
N ++; // 结果当前位为1，否则为默认的0
-
tmp -= N;
-
}
-
for (i=15; i>0; i--) // 求剩余的15位
-
{
-
N <<= 1; // 左移一位
-
tmp <<= 2;
-
tmp += (M >> 30); // 假设
-
ttp = N;
-
ttp = (ttp<<1)+1;
-
M <<= 2;
-
if (tmp >= ttp) // 假设成立
-
{
-
tmp -= ttp;
-
N ++;
-
}
-
}
-
return N;
-
}
===================================================================
4.测试
芯片 DSP2808
主频 100MHZ
测试数据范围：0-2000,000
测试程序附上
-
#include "IQmathLib.h"
-
#include "CMMN_SciA.h"
-
#include <math.h>
-
#include "CommonUseFun.h"
-
-
extern UINT32 Get_Timer0_Timecnt();
-
#define CalTimeUse(TimeCntUse,TimerCnt0,TimerCnt1) do \
-
{\
-
if (TimerCnt0>=TimerCnt1)\
-
{\
-
TimeCntUse=TimerCnt0-TimerCnt1;\
-
}\
-
else\
-
{\
-
TimeCntUse=TimerCnt0+100000uL-TimerCnt1;\
-
}\
-
} while (0)
-
UINT32 NumIn=0;
-
UINT16 StepInTimeCnt=0;
-
-
void Sqrt_Test()
-
{
-
UINT32 TimerCnt[3]={0};
-
UINT32 sqrt_num[3]={0};
-
UINT16 TimeCntUse[2]={0};
-
++StepInTimeCnt;
-
if (StepInTimeCnt<10)//考虑串口发送的累计，10ms计算一次
-
{
-
return;
-
}
-
-
StepInTimeCnt=0;
-
if (NumIn>2000000uL) //4096*4096=16777216
-
{
-
return;
-
}
-
-
NumIn+=10;
-
-
TimerCnt[0]=Get_Timer0_Timecnt();
-
sqrt_num[0]=(UINT32)(sqrt(NumIn)*100);
-
TimerCnt[1]=Get_Timer0_Timecnt();
-
sqrt_num[1]=(UINT32)(sqrt_16(NumIn*100)*10);//有bug,10->10uL
-
TimerCnt[2]=Get_Timer0_Timecnt();
-
-
CalTimeUse(TimeCntUse[0],TimerCnt[0],TimerCnt[1]);
-
CalTimeUse(TimeCntUse[1],TimerCnt[1],TimerCnt[2]);
-
send_uint32data_by_char_in_QueueBuf(NumIn);
-
send_1char_in_QueueBuf('\t');
-
send_uint32data_by_char_in_QueueBuf(sqrt_num[0]);
-
send_1char_in_QueueBuf('\t');
-
send_uint32data_by_char_in_QueueBuf(sqrt_num[1]);
-
send_1char_in_QueueBuf('\t');
-
send_uidata_by_char_in_QueueBuf(TimeCntUse[0]);
-
-
send_uidata_by_char_in_QueueBuf(TimeCntUse[1]);
-
if (ABS(sqrt_num[0],sqrt_num[1])>=100)
-
{
-
send_string_in_QueueBuf("\tNE\t");
-
}
-
else
-
{
-
send_string_in_QueueBuf("\t\t");
-
}
-
-
if (TimeCntUse[1]<(TimeCntUse[0]-500))//差值大于5us
-
{
-
send_string_in_QueueBuf("S\n");
-
}
-
else if(TimeCntUse[1]>TimeCntUse[0])
-
{
-
send_string_in_QueueBuf("L\n");
-
}
-
else
-
{
-
send_string_in_QueueBuf("E\n");
-
}
-
-
}
测试结果
注:
第1列 Num
第2列 库sqrt()函数计算后放大100结果
第3列 sqrt_16()函数计算后放大100倍计算结果
第4列 库sqrt()函数计算TimeCnt 精确到0.01
第5列 sqrt_16函数计算TimeCnt 精确到0.1
第6列 库sqrt()函数和sqrt_16函数计算结果差值 NE-相差大 无-相近
第7列 库sqrt()函数和sqrt_16函数计算所花时间，S计算时间远小
-
680 2607 2600 4241 826 S
-
690 2626 2620 4239 823 S
-
700 2645 2640 4241 826 S
-
710 2664 2660 4243 823 S
-
720 2683 2680 4241 823 S
-
730 2701 2700 4241 820 S
-
-
10680 10334 10330 4239 823 S
-
10690 10339 10330 4239 823 S
-
10700 10344 10340 4237 823 S
-
10710 10348 10340 4239 823 S
-
10720 10353 10350 4243 820 S
-
10730 10358 10350 4239 820 S
-
10740 10363 10360 4241 823 S
-
10750 10368 10360 4239 823 S
-
10760 10373 10370 4239 820 S
-
-
10810 10397 10390 4243 817 S
-
10820 10401 10400 4243 826 S
-
10830 10406 10400 4243 826 S
-
10840 10411 10410 4239 823 S
但是也出现了如下的问题
-
429500 65536 65530 4225 811 S
-
429510 65537 65530 4227 811 S
-
429520 65537 65530 4225 811 S
-
429530 65538 65530 4227 811 S
-
429540 65539 65530 4227 811 S
-
429550 65540 4 4227 811 NE S
-
429560 65540 4 4227 811 NE S
-
429570 65541 4 4227 811 NE S
-
429580 65542 4 4227 811 NE S
-
429590 65543 4 4227 811 NE S
-
429600 65543 4 4227 811 NE S
-
429610 65544 4 4227 811 NE S
-
429620 65545 4 4229 811 NE S
经分析，原因是测试程序中
sqrt_num
[1
]=(UINT32
)(sqrt_16
(NumIn
*100
)*10
)存在问题，
应修改为sqrt_num
[1
]=(UINT32
)(sqrt_16
(NumIn
*100
)*10
uL)
结论：
使用上述的
unsigned int sqrt_16(unsigned long M)函数，计算时间较库函数时间短很多
库sqrt()函数计算时间：42.2~42.3us
sqrt_16()函数计算时间：8.1~8.2us
可见sqrt_16()对长整型数据计算快速性很好。
但是需要对其优化才能在精度上达到要求。
阅读(4982) | 评论(0) | 转发(0) |
