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

Allevi · 发表于 2019-5-24 10:22

EDA365欢迎您登录！

您需要登录才可以下载或查看，没有帐号？注册

x

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

ut 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位整数.

helendcany · 发表于 2019-5-24 18:23

谢谢分享

帐号		自动登录	找回密码
密码			注册

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

EDA365欢迎您登录！

浏览过的版块

推荐内容 /1