FPGA开方实现一种简单的算法

eeleader

FPGA开方实现一种简单的算法 [复制链接]

理论基础:

   我们先来看看10进制下是如何手工计算开方的。
先看下面两个算式，
x = 10*p + q  (1)
公式(1)左右平方之后得：
x^2 = 100*p^2 + 20pq + q^2 (2)
现在假设我们知道x^2和p，希望求出q来，求出了q也就求出了x^2的开方x了。
我们把公式(2)改写为如下格式：
q = (x^2 - 100*p^2)/(20*p+q) (3)

这个算式左右都有q，因此无法直接计算出q来，因此手工的开方算法和手工除法算法一样有一步需要猜值。

我们来一个手工计算的例子：计算1234567890的开方

首先我们把这个数两位两位一组分开，计算出最高位为3。也就是(3)中的p，最下面一行的334为余数，也就是公式(3)中的(x^2 - 100*p^2)近似值
    3
  ---------------
 / 12 34 56 78 90
    9
  ---------------
 /  3 34

下面我们要找到一个0-9的数q使它最接近满足公式(3)。我们先把p乘以20写在334左边：
                           3  q
                         ---------------
                        / 12 34 56 78 90
                           9
                         ---------------
(20*3+q)*q      /  3 34

我们看到q为5时(60+q)*q的值最接近334，而且不超过334。于是我们得到：
      3  5
    ---------------
   / 12 34 56 78 90
      9
    ---------------
65 /  3 34
      3 25
    ---------------
         9 56

接下来就是重复上面的步骤了，这里就不再啰嗦了。

这个手工算法其实和10进制关系不大，因此我们可以很容易的把它改为二进制，改为二进制之后，公式(3)就变成了：
q = (x^2 - 4*p^2)/(4*p+q) (4)

我们来看一个例子，计算100(二进制1100100)的开方：
       1  0  1  0
      -----------
     / 1 10 01 00
       1
      -----------
 100 / 0 10
       0 00
      -----------
1001 /   10 01
         10 01
      -----------
          0 00

这里每一步不再是把p乘以20了，而是把p乘以4，也就是把p右移两位，而由于q的值只能为0或者1，所以我们只需要判断余数(x^2 - 4*p^2)和(4*p+1)的大小关系，如果余数大于等于(4*p+q)那么该上一个1，否则该上一个0。

下面给出完成的C语言程序，其中root表示p，rem表示每步计算之后的余数，divisor表示(4*p+1)，通过a>>30取a的最高 2位，通过a<<=2将计算后的最高2位剔除。其中root的两次<<1相当于4*p。程序完全是按照手工计算改写的，应该不难理解。
unsigned short sqrt(unsigned long a){
  unsigned long rem = 0;
  unsigned long root = 0;
  unsigned long divisor = 0;
  for(int i=0; i<16; ++i){
    root <<= 1;
    rem = ((rem << 2) + (a >> 30));
    a <<= 2;
    divisor = (root<<1) + 1;
    if(divisor <= rem){
      rem -= divisor;
      root++;
    }
  }
  return (unsigned short)(root);
}

根据上面算法:设计一个36二进制数的开方实现程序如下:

library ieee;
use ieee.std_logic_1164.all;
use ieee.std_logic_arith.all;
use ieee.std_logic_unsigned.all;
entity Svg_Sqrt is
 port
   (
----系统时钟32MHZ, 16MHZ经过PLL 2倍频，复位信号低电平有效
     In_rst       : in std_logic;
     In_clkmain   : in std_logic; 
     In_En:in std_logic;
     In_x2:in std_logic_vector(35 downto 0);
 
     Out_x:out std_logic_vector(17 downto 0)
  --   Test1: out std_logic_vector(35 downto 0);
  --   Test2: out std_logic_vector(35 downto 0);
  --   Temp : buffer std_logic_vector(7 downto 0);
  --   Cnt: buffer integer range 0 to 17
    );
end Svg_Sqrt;

architecture Arch_Svg_Sqrt of Svg_Sqrt is
signal Temp : std_logic_vector(7 downto 0);
signal m2: std_logic_vector(35 downto 0);
signal quot:std_logic_vector(35 downto 0);
signal root: std_logic_vector(35 downto 0);
signal divisor: std_logic_vector(35 downto 0);
signal cnt  : integer range 0 to 17 ;

begin
process(In_rst,In_clkmain)
begin
  if (In_rst='0') then
           root<=(others=>'0');
           quot<=(others=>'0');
           divisor<=(others=>'0');
           m2<=(others=>'0');
           Temp<=x"00";
           Cnt<=0;
  elsif In_clkmain'event and In_clkmain='1'  then
     if (Temp=x"00") then
        if (In_En='1') then
           m2<=In_x2;
           root<=(others=>'0');
           quot<=(others=>'0');
           divisor<=(others=>'0');
           Temp<=x"01";
        end if;
     elsif (Temp=x"01") then
           root<=root(34 downto 0) & "0";
           quot<=((quot(33 downto 0) & "00")+((x"00000000" & "00") &  m2(35 downto 34)));
       --    Test2<=((quot(33 downto 0) & "00")+((x"00000000" & "00") &  m2(35 downto 34)));   
           Temp<=x"02";
     elsif (Temp=x"02") then
           divisor<=(root(34 downto 0) & "0")+'1';
           m2<=m2(33 downto 0) & "00";  
           Temp<=x"03";
        --   Test1<=(root(34 downto 0) & "0")+'1';
     elsif (Temp=x"03") then
           if (divisor<=quot) then
               quot<=quot-divisor;
               root<=root+'1';
           end if;
           Temp<=x"04";
     elsif (Temp=x"04") then
           if (Cnt=17) then
               Cnt<=0;
               Out_x<=root(17 downto 0);
               Temp<=x"00";
           else
               Cnt<=Cnt+1;
               Temp<=x"01";
           end if;
     else
           Temp<=x"00";
     end if;
 end  if;
end process;
end Arch_Svg_Sqrt;

36位整数开方结果为18位整数.

zhuoyue

不错~~~很详细~

常见泽1