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\),可以將他們表示為：
其中，\(m,n\)互質。計\(A,B\)的差為\(C\)：
首先，可以證明，\(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'-n')\)與\(n'\)一定會變成最小的互質數對\((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應該就行，懶得改了）
下面是仿真波形：
