Verilog -- 求两数最大公因数(gcd)和最小公倍数(lcm)

Verilog – 求两数最大公因数和最小公倍数

@(verilog)

1. 原理简介

1.1 辗转相除法求公因数

求最大公因数的常用算法为辗转相除法，又被称为欧几里德（Euclidean）算法， 是求最大公约数的算法。辗转相除法首次出现于欧几里得的《几何原本》（第VII卷，命题i和ii）中，而在中国则可以追溯至东汉出现的《九章算术》。

两个数的最大公约数是指能同时整除它们的最大正整数。辗转相除法的基本原理是：两个数的最大公约数等于它们中较小的数和两数之差的最大公约数。
例如，252和105的最大公约数是21（252 = 21 × 12；105 = 21 × 5）；因为252 − 105 = 147，所以147和105的最大公约数也是21。在这个过程中，较大的数缩小了，所以继续进行同样的计算可以不断缩小这两个数直至其中一个变成零。这时，所剩下的还没有变成零的数就是两数的最大公约数。

可以简单证明一下：
设$A$和$B$的最大公因数为$gcd=X$,且$A>B$,可以将他们表示为：
$$\begin{gathered} A=mX,m\in Z \\ B=nX,n\in Z \end{gathered}$$
其中，$m,n$互质。计$A,B$的差为$C$：
$$C=(m-n)X$$
首先，可以证明，$m-n$与$m,n$都互质，因为$m-n$如果跟$n$有公因数，则$m,n$必也有公因式，违反了条件。因此互质数的和差与原来两数都互质。
其次，由于$m-n$与$n$互质，因此$B,C$的最大公因式就是$X$,因为没办法再从$m-n$与$n$中提出一个因式放到$X$里。
所以，$A,B$的最大公约数就是$B,C$的最大公约数。
由此，就可以将$A,B$“辗转”地相互做差，这样不断地进行，最后一定会出现要做差的两个数都为$X$的情形，因为做差做到最后，$(m^{\prime }-n^{\prime })$与$n^{\prime }$一定会变成最小的互质数对$(1,1)$,也就是经过若干次辗转后较小的那个数与他们的差相等，此时就得到了最大公约数。

事实上，上面的是辗转相减法，为了加速相减的过程，实际上可以使用除法加速辗转的过程。（除法就是若干次减法）

1.2 最小公倍数求法

最小公倍数实际上等于两数之积除以他们的最大公约数。为什么呢？还是可以从上面的假设进行推导：
首先我们知道，互质数的最小公倍数一定是他们的积。而要求$A,B$的最小公倍数，因为它们已经有了最大的公约数$X$,所以其实只要求互质的两个数$m,n$的最小公倍数，所以$A,B$的最小公倍数$lcm=mnX=AB/X$。

2. 代码实现

module lcm_gcd
#(
parameter DATAWIDTH=8
)
(
input   [DATAWIDTH-1:0]       A,
input   [DATAWIDTH-1:0]       B,
input                         en,
input                         clk,
input                         rstn,

output  wire [DATAWIDTH*2-1:0]lcm_out,
output  wire [DATAWIDTH-1:0]  gcd_out,
output  wire                  vld_out,
output  wire                  ready
);
reg [DATAWIDTH-1:0]     gcd,a_buf,b_buf;
reg [DATAWIDTH*2-1:0]   lcm,mul_buf;
reg [1:0]               current_state,next_state;
reg done;
parameter IDLE=2'b00,
          SUB  =2'b01,
		  DONE  =2'b10;
        
assign ready=(current_state==IDLE)?1'b1:1'b0;

always@(posedge clk or negedge rstn)
  if(!rstn) current_state <= IDLE;
  else current_state <= next_state;

always@(*) begin
  next_state = 2'bx;
  case(current_state)
    IDLE:if(en)  next_state <= SUB;
         else    next_state <= IDLE;	
    SUB:  if(a_buf!=b_buf) next_state <= SUB;
         else             next_state <= DONE;
    DONE:  next_state <= IDLE;
  endcase
end


  
always@(posedge clk or negedge rstn)
     if(!rstn)begin
		gcd    <= 0;
        lcm    <= 0;
		a_buf  <= 0;
		b_buf  <= 0;	
        mul_buf<= 0;	
        done   <= 0;    
      end
     else begin
	    case(current_state)
        IDLE:begin
            if(en)begin
		       a_buf   <= A;
			   b_buf   <= B;
               mul_buf <= A*B;              
             end
             done <= 0;
           end
	    SUB:if(a_buf != b_buf)begin
		     if(a_buf > b_buf)begin
			    a_buf <= a_buf-b_buf;
				b_buf <= b_buf;
			 end
			 else begin
			    b_buf <= b_buf-a_buf;
				a_buf <= a_buf;			 
              end
           end		
		DONE:begin
            gcd <= b_buf;
            lcm <= mul_buf / b_buf; 
            done <= 1;
          end

	    endcase
      end	 

assign gcd_out = gcd;
assign lcm_out = lcm;
assign vld_out = done;

endmodule

testbench:

`timescale 1ns/1ps



module lcm_gcd_tb();

parameter DATAWIDTH = 8;
reg   [DATAWIDTH-1:0]       A;
reg   [DATAWIDTH-1:0]       B;
reg                         en;
reg                         clk;
reg                         rstn;

wire [DATAWIDTH*2-1:0]      lcm_out;
wire [DATAWIDTH-1:0]        gcd_out;
wire                        vld_out;
wire                        ready;


reg [DATAWIDTH*2-1:0]      lcm_ref;
reg [DATAWIDTH-1:0]        gcd_ref;

always #1 clk = ~clk;

function automatic integer gcd(input integer A,input integer B);
    if(A>B) 
      return gcd(B,A-B);
    else if (A<B) 
      return gcd(A,B-A);
    else 
      return A;
endfunction


integer i;
initial begin
  clk = 1;
  rstn = 1;
  en = 0;
  #2 rstn = 0; #2 rstn = 1;
  repeat(2) @(posedge clk);

  for(i=0;i<20;i=i+1) begin
    en <= 1;
    A <= $urandom()%20+1; // no zero input
    B <= $urandom()%20+1;
    gcd_ref = gcd(A,B);
    lcm_ref = A*B/gcd_ref;
    @(posedge vld_out);
   end

end


initial begin
    $fsdbDumpvars();
    $fsdbDumpMDA();
    $dumpvars();
    #1000 $finish;
  end


lcm_gcd #(
    .DATAWIDTH                      ( DATAWIDTH                     ))
U_LCM_GCD_0(
    .A                              ( A                             ),
    .B                              ( B                             ),
    .en                             ( en                            ),
    .clk                            ( clk                           ),
    .rstn                           ( rstn                          ),
    .lcm_out                        ( lcm_out                       ),
    .gcd_out                        ( gcd_out                       ),
    .vld_out                        ( vld_out                       ),
    .ready                          ( ready                         )
);



endmodule

上面代码里用到的除法如果无法综合或者延时过大可以考虑使用时序逻辑实现，此外，该代码无法处理输入为0的情况，而且其实输入为0时求公因数也无意义，如有需要可自行修改。（在最开始判断一下是否为0应该就行，懒得改了）
下面是仿真波形：