实数求根的一种快速实现方法

实数求根的一种快速实现方法

说明：本文主要根据如下链接的文章改写整理而成，框图及代码均为原文作者提供，版权归原文作者所有。有版权问题请联系博主。

http://ieeexplore.ieee.org/xpl/articleDetails.jsp?tp=&arnumber=5888664&contentType=Journals+%26+Magazines&queryText%3DMultiplier-Free+Divide%2C+Square+Root%2C+and+Log+Algorithms

嵌入式芯片运算处理单元一般都有加法器，有的还有专门的乘法器。但要实现求根的运算，特别是速度还要求比较高的时候，多是采用查表的方法。这里介绍一种实数求根的实现方法。这种方法在计算平方根的时候，甚至都不需要乘法运算，效率非常高。

给定一个正实数c，要计算其平方根，即sqrt(c)。基本的思路是定义两个序列，x(n)和c(n)，并且定义x(0)=0，c(0)=c。保证：

x(n+1)^2+c(n+1)=x(n)^2+c(n)=x(0)+c(0)=c (1)

这样，随着n的增加，c(n)逐步减小，直到0，x(n)则逐步增加，直到x(n)^2=c。

上述算法的实现框图如图1所示。参数c,h的含义已经介绍过了，Nbits是结果的精度要求，用来作为迭代的终止条件。比如要求除法结果的精度为小数点后8位，则Nbits选为8即可。若要求精确到64位，则Nbits = 64。这种算法不需要乘法运算的关键是h的选择。只要将h选为2的幂次方，则与h有关的乘法可用移位来实现，从而提高了算法的运算效率。




下面给出的是上述算法的matlab实现代码。

function [xn,k,Trajectory]=dcd_sqrt(c,h_init,Mb,Nu)

% [x,k,Trajectory]=dcd_sqrt(c,h_init,Mb,Nu)

% multiplication-free square root of c

% inspired from the algorithm proposed by Liu, Weaver and Zakharov

%

%

% h_init is the initial value of h (and should be a power of two)

% Mb is the number of bits used to code x

% Nu is the maximal number of successfull iterations.

%

% Trajectory variable is computed only for algorithm characterization.

%

% Francois Auger, august 2010

if (c<=0), error('c should be strictly positive'); end;



m=1; xn=0; cn=c; h=h_init;

hdiv2=h/2; htimes2=h*2; hsquare=h*h; % left or right shift



k=0; if (nargout>=3), Trajectory(:,1)=[xn;cn]; end;



while (m<=Mb)&(k<=Nu),



while (cn>=hsquare+xn*htimes2)||((cn+xn*htimes2>=hsquare)&&(xn<hdiv2))



if (cn>=hsquare+xn*htimes2) % left or right shift

cn=cn-hsquare-xn*htimes2; xn=xn+h; % left or right shift

else

cn=cn-hsquare+xn*htimes2; xn=xn-h; % left or right shift

end; % end if (cn>hsquare+xn*htimes2)



k=k+1; if (nargout>=3), Trajectory(:,k+1)=[xn;cn]; end;



end; % end while (cn>=hsquare+xn*htimes2)|((cn+xn*htimes2>=hsquare)&(xn<hdiv2))



h=h/2; hdiv2=hdiv2/2; htimes2=htimes2/2; hsquare=hsquare/4; m=m+1; % right shift

end; % end while (m<=Mb)&(k<=Nu)