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
default:next_state<= IDLE;
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应该就行，懒得改了）
下面是仿真波形：
