[Leka 0]

Не ту ссылку дал, перепутал с делением.

Уоррен "Алгоритмические трюки для программистов", там есть "аппаратный" алгоритм на Си:

unsigned x, y;
unsigned sqrt(){
  y = 0;
  unsigned m = 1 << 30;    
  while( m ){
    unsigned b =  y | m;
    y >>= 1;
    if( x >= b ){ 
        x -= b; 
        y |= m; 
    }             
    m >>= 2;
  }    
}

Изменено 2 мая, 2013 пользователем Leka

[Timmy 1]

Там код оптимизирован под ECP2, чтобы его собрать в общем виде, надо сменить архитектуры add_sub, add_mux на common.

А вообще этот ужас надо переписать по нормальному non-restoring алгоритму, но руки пока не дошли, и так работает и в 200МГц укладывается.

sqrt.7z

[Timmy 1]

Вот, всё-таки дошли руки, всё очень просто.

library ieee;
use ieee.std_logic_1164.all;
use ieee.std_logic_arith.all;
use ieee.std_logic_unsigned.all;
use std_lib.all;
entity sqrt_e is
port(
	clk, rst, start:std_logic;
	di:std_logic_vector;
	do:out std_logic_vector;
	so:out std_logic
);
end entity;

architecture arch of sqrt_e is
   signal sqrt:std_logic_vector(do'length-1 downto 0);
   signal sqrt_ext:std_logic_vector(sqrt'length downto 0);
   signal rm:std_logic_vector(sqrt'length+2 downto 0);
   signal sum, a1, a2:unsigned(sqrt'length+2 downto 0);
   signal d:unsigned(di'length-1 downto 0);
   signal sum_h:std_logic_vector(0 downto 0);
begin
sqrt_ext <= sqrt&(not sum_h);
sum_h(0) <= sum(sum'high);
a1 <= trunc_high(sum, 2) & part_high(d, 2);
a2 <= unsigned(sqrt_ext&sum_h&"1");
process(clk) is begin
   if rising_edge(clk) then
       sqrt <= trunc_high(sqrt_ext, 1);
       sum <= iif(sum_h(0)='1', a1+a2, a1-a2);
       d <= trunc_high(d, 2)&"00";
       if start = '1' then
           sum <= (others=>'1');
           sqrt <= (others=>'0');
           d <= unsigned(di);
       end if;
       do <= std_logic_vector(round(unsigned(sqrt_ext)));
   end if;
end process;
end architecture;

И подправленный тестбенч, теперь на 2 такта дольше.

library ieee;
library std;
use ieee.std_logic_1164.all;
use ieee.std_logic_arith.all;
use work.std_lib.all;
entity sqrt_tb is
end entity;

architecture tb of sqrt_tb is
constant dd:integer:=10;
signal d:unsigned(dd-1 downto 0);
signal do:unsigned(dd-1 downto 0);
signal dop:std_logic_vector(dd-1 downto 0);
signal clk:std_logic:='0';
signal start:std_logic;
signal test,testm,testp,testo:unsigned(dd*2-1 downto 0);
signal error:boolean:=false;
begin
main:process is 
function sqr(v:unsigned)return unsigned is begin
	return v*v;
end function;
variable test_h,test_l,orig:unsigned(dd*2+1 downto 0);
begin
for i in 0 to 2**dd-1 loop
	d <= conv_unsigned(i, d'length);
	start <= '1';
	wait until rising_edge(clk);
	start <= '0';
	for i in 0 to dd+2 loop
		wait until rising_edge(clk);
	end loop;
	test_l := sqr((do-1)&"1");
	test_h := sqr(do&"1");
	orig := align_high(d, orig'length);
	if i/=0 then
		if test_l>orig or test_h < orig then	
			error <= true;
		end if;
	else
		if do/=0 then error <= true; end if;
	end if;
end loop;
std.env.stop(2);
end process;
test <= do*do;
testm <= (do-1)*(do-1);
testp <= (do+1)*(do+1);
--mov_high_v(testo, d);
testo(testo'high downto testo'high - d'length+1) <= d;
testo(testo'high - d'length downto 0) <= (others=>'0');

clk <= not clk after 5 ns;
UUT:entity work.sqrt_e port map(clk=>clk, start=>start, di=>std_logic_vector(d),do=>dop,rst=>'0');
do <= unsigned(dop);
end architecture;

 

[Maverick_ 15]

Не ту ссылку дал, перепутал с делением.

Уоррен "Алгоритмические трюки для программистов", там есть "аппаратный" алгоритм на Си:

unsigned x, y;
unsigned sqrt(){
  y = 0;
  unsigned m = 1 << 30;    
  while( m ){
    unsigned b =  y | m;
    y >>= 1;
    if( x >= b ){ 
        x -= b; 
        y |= m; 
    }             
    m >>= 2;
  }    
}

Вы мне не расскажите цифра 30 для m это ограничение или как?

Алгоритм похож на мой найденный :) и приведено для него описание (последнее)

[Leka 0]

цифра 30 для m это ограничение или как?

Это упрощение для 32-разрядного процессора, для FPGA такого ограничения нет.

 

[Maverick_ 15]

Если кто-то перепишит состояние calc (описание в сообщении), чтобы увеличить быстродействие (понимаю, что при этом увеличится число тактов на вычисление квадратного корня и скорее всего понадобиться большее кол-во состояний автомата)... Буду благодарен...

 

PS Вижу, что нужно

op <= op - (res+one);

разнести вычисление по времени...

[Leka 0]

вместо op <= op - (res + one);

можно op <= op - (res | one); //или "OR" (не знаю, как надо в VHDL)

в Спартанах это один уровень логики займет, не считая переноса.

[Maverick_ 15]

вместо op <= op - (res + one);

можно op <= op - (res | one); //или "OR" (не знаю, как надо в VHDL)

в Спартанах это один уровень логики займет, не считая переноса.

заменил в 2-х местах (в тех строчках - в коментариях "старое" описание), увеличил быстродействие на ~10 МГц (только здесь Generic ( b = 32), в ранее приведенном 16 - пишу чтобы учитывали в тетсбенче). Сравнение проводил для 32 битных входных данных. ПЛИС - Spartan 6.

 

-- algorithm in C code
------------
--int sqrt(int num) {
--    int op = num;
--    int res = 0;
--    int one = 1 << 30; // The second-to-top bit is set: 1L<<30 for long
--
--    // "one" starts at the highest power of four <= the argument.
--    while (one > op)
--        one >>= 2;
--  
--    while (one != 0) {
--        if (op >= res + one) {
--            op -= res + one;
--            res += 2 * one;
--        }
--        res >>= 1;
--        one >>= 2;
--    }
--    return res;
--}
-----------



library IEEE;
use IEEE.STD_LOGIC_1164.ALL;
use IEEE.NUMERIC_STD.ALL;

entity SQRT is
    Generic ( b  : natural range 4 to 32 := 32 );
    Port ( clk    : in   STD_LOGIC;
           start  : in   STD_LOGIC;
           value  : in   STD_LOGIC_VECTOR (31 downto 0);
           result : out  STD_LOGIC_VECTOR (15 downto 0);
           busy   : out  STD_LOGIC);
end SQRT;

architecture Behave of SQRT is
signal op  : unsigned(b-1 downto 0); --
signal res : unsigned(b-1 downto 0); --
signal one : unsigned(b-1 downto 0); --

signal bits : integer range b downto 0;
type aaa is (idle, shift, calc, done);
signal z : aaa;

begin
   process begin
      wait until rising_edge(clk);
      case z is
         when idle =>
            if (start='1') then
               z <= shift;
               busy <= '1';
            end if;
            one <= to_unsigned(2**(b-2),b);
            op  <= unsigned(value);
            res <= (others=>'0');

         when shift =>
            if (one > op) then
               one <= one/4;
            else
               z   <= calc;
            end if;

         when calc =>
            if (one /= 0) then
             if (op >= res+one) then
                 op   <= op - (res or one); --op   <= op - (res+one);
                  res  <= res/2 or one; --res  <= res/2 + one;
               else
                  res  <= res/2;
               end if;
               one <= one/4;
            else
               z    <= done;
            end if;
            
         when done =>
            busy <= '0';
            if (start='0') then
               z <= idle;
            end if;
      end case;
   end process;
  
   result <= std_logic_vector(res(result'range));
end;

[Leka 0]

Быстродействие должно ограничиваться только этим:

if (op >= (res or one)) then

op <= op - (res or one);

 

Можно попробовать оптимизировать описание этого кусочка по результатам синтеза, например на Верилоге (VHDL не знаю):

wire[32:0] tmp = op - (res | one);

...

if(!tmp[32]) op <= tmp;

 

Д/б >200МГц

[Maverick_ 15]

signal bits : integer range b downto 0;

можно удалить...

в последнем приведенном описании...

[eugen_pcad_ru 0]

Тут уже много рекомендаций:)

Если ни одна не устроит, то была где-то древняя реализация на АHDL. Если нужно, пишите в личку, поищу.

[Maverick_ 15]

Всем спасибо.

[Maverick_ 15]

добавлю еще алгоритм

[ukpyr 0]

еще вариация (16 бит, но расширяется до 32 или 64 бита):

uint8_t sqrt16(uint16_t n) {
    uint8_t  guess = 0x00;
    uint8_t  bit   = 0x80;
    
    do {
        uint8_t  guess1 = guess | bit;
        if (n > (guess1 * guess1)) guess = guess1;
        bit >>= 1;
    } while (bit);
    
    return guess;
}

Изменено 21 июля, 2013 пользователем ukpyr

[Kuzmi4 0]

А как кто решает проблемы с дробными частями, ну что там после запятой ?

Потому как, например, вышеприведённый пример вычисления корня даёт нам для входного значения "35" - результат "5", то есть 0.916... остались неучтёнными..