postergre除法 raptor除法取余数

文章目录

• 1 基于恢复余数（Restoring）算法的除法器

• 1.1 原理
• 1.1 Verilog 实现

• 2 基于不恢复余数（Non-Restoring）算法的除法器

• 2.1 原理
• 2.1 Verilog 实现
• 2.2 用于有符号数的改进

• 3 基于级数展开算法的除法器

• 4.2 Verilog 实现

• 4 基于 Newton-Raphson 算法的除法器

• 4.1 原理

参考《基于FPGA的数字信号处理（第二版）》，对其中的除法器算法进行仿真。

1 基于恢复余数（Restoring）算法的除法器

Restoring 算法无法直接用于与有符号数，对于有符号数需要先转换为无符号数。然后根据除数与被除数的符号判断商与余数的符号。

1.1 原理

被除数（Dividend）5621 / 除数（Divisor）3 = 商（Quotient）1873 … 余数（Remainder）2

1873
     _______
  3 |  5621
    /  3000
       ____
       2621
       2400  
       ____
        221
        210
        ___
         11
          9
          _
          2

$r_i = r_{i+1}-q_i·D·10^i$

$q_i$通过部分余数的正负确定，如果余数为负表示$q_i$较大，需要将其恢复到部分余数为正以获取正确的$q_i$，这便是“恢复”的含义。

用恢复余数算法表示上式计算过程：

	余数	是否恢复
	5621-3×10^3×1=2621
	5621-3×10^3×2=-379	恢复为1，余2621
	2621-3×10^2×1=2321
	2621-3×10^2×2=2021
……	……
	2621-3×10^2×8=221
	2621-3×10^2×9=-79	恢复为8，余221
……	……
	221-3×10^1×7=11
	221-3×10^1×8=-19	恢复为7，余11
……	……
	11-3×10^0×3=2
${q_3q_2q_1q_0}$	11-3×10^0×4=-1	${q_3q_2q_1q_0}$恢复为3，余2

得 ${q_3q_2q_1q_0}$

以上是十进制算法，对于二进制：

$r_i = r_{i+1}-q_i·D·2^i$

$q_i$只能取0或1，所以每位只需一次运算即可判断是否恢复。

此外还需确定各参数的字长：假定除数D为 nbit，商Q为 nbit，则 $Q\le2^n-1$。由于R<D，所以R的位宽可以设置为 nbit。则被除数Y:
$Y=D·Q+R < (2^n-1)D+D=2^nD$
意味着Y的高 nbit 小于D。

基于 restoring 算法的无符号数除法运算具有如下特征：

1. 各参数字长：被除数为 2nbit，除数、商和余数为 nbit。
2. 被除数的高 nbit 需小于 除数。
3. 算法需迭代 n次。

算法流程：

基本运算单元（移位和减法）硬件架构：

延时分析：
上述硬件架构运算 latency 为2，对于 nbit 商，需要 2n latency。去掉寄存输出可降低延迟。

1.1 Verilog 实现

流水线实现方式比较简单，这里实现运算单元分时复用。

`timescale 1ns/1ps

module div_restoring_pip (
    input             I_sys_clk  ,
    input             I_reset_n  ,
    input             I_valid    ,
    input      [15:0] I_dividend ,
    input      [7:0]  I_divisor  ,
    output reg        O_valid    ,
    output reg [7:0]  O_quotient ,
    output reg [7:0]  O_remainder
);
//--- internal signal definitions ---
//=== parameter definitions ===			
                                   
//=== reg definitions ===  
reg        R_valid           ;			
reg [7:0]  R_divisor         ;
reg [3:0]  R_index_cnt       ;
reg [15:0] R_remainder_t     ;
reg [15:0] R_di              ;
reg [7:0]  R_quotient        ;
reg        R_restoring_valid ;


//=== wire definitions === 	
wire        W_restoring_valid;
wire        W_quotient       ;			
wire [15:0] W_remainder      ;                                   
                                    
//--- Main body of code ---   
// 非流水线处理，需要 blocking
always @(posedge I_sys_clk or negedge I_reset_n)
begin
    if(~I_reset_n) 
    begin
        R_valid    <= 1'b0;
        R_divisor  <= 1'b0;
    end 
    else
    begin
        R_valid    <= I_valid;
        if (I_valid)
        begin
            R_divisor  <= I_divisor ;
        end
    end
end

// 共计算8位的商
always @(posedge I_sys_clk or negedge I_reset_n)
begin
    if(~I_reset_n) 
    begin
        R_index_cnt <= 4'd0;
    end 
    else
    begin
        if (I_valid)
        begin
            R_index_cnt <= 4'd8;
        end 
        else if (W_restoring_valid)
        begin
            R_index_cnt <= R_index_cnt - 4'd1;
        end
    end
end

// 得到下一次计算的输入
always @(posedge I_sys_clk or negedge I_reset_n)
begin
    if(~I_reset_n) 
    begin
        R_remainder_t <= 16'd0;
        R_di <= 16'd0;
    end 
    else
    begin
        if (I_valid)
        begin
            R_remainder_t <= I_dividend ;
            R_di <= I_divisor << 3'd7;
        end
        else if (W_restoring_valid)
        begin
            R_remainder_t <= W_remainder ;
            R_di <= R_divisor << (R_index_cnt-2);
        end
    end
end

// 保存每位的商
always @(posedge I_sys_clk or negedge I_reset_n)
begin
    if(~I_reset_n) 
    begin
        R_quotient <= 8'd0;
    end 
    else
    begin
        if (I_valid)
        begin
            R_quotient <= 8'd0;
        end
        else if (W_restoring_valid)
        begin
            R_quotient <= {R_quotient[6:0], W_quotient};
        end
    end
end

// 输出
always @(posedge I_sys_clk or negedge I_reset_n)
begin
    if(~I_reset_n) 
    begin
        O_valid     <= 1'd0;
        O_quotient  <= 8'd0;
        O_remainder <= 8'd0;
    end 
    else
    begin
        if (W_restoring_valid && (R_index_cnt == 4'd1))
        begin
            O_valid     <= 1'd1;
            O_quotient  <= {R_quotient[6:0], W_quotient};
            O_remainder <= W_remainder;
        end
        else 
        begin
            O_valid     <= 1'd0;
            O_quotient  <= 8'd0;
            O_remainder <= 8'd0;
        end
    end
end

always @(posedge I_sys_clk or negedge I_reset_n)
begin
    if(~I_reset_n) 
    begin
        R_restoring_valid <= 1'b0;
    end 
    else
    begin
        R_restoring_valid <= W_restoring_valid;
    end
end

restoring restoring_u (
    .I_sys_clk(I_sys_clk),
    .I_reset_n(I_reset_n),
    .I_valid  (R_valid | (R_restoring_valid && (R_index_cnt != 4'd0))),
    .I_R      (R_remainder_t),
    .I_Di     (R_di),
    .O_valid  (W_restoring_valid),
    .O_Q      (W_quotient),
    .O_R      (W_remainder)  
);

endmodule


module restoring (
    input             I_sys_clk,
    input             I_reset_n,
    input             I_valid,
    input      [15:0] I_R,
    input      [15:0] I_Di,
    output reg        O_valid,
    output reg        O_Q,
    output reg [15:0] O_R
);
//--- internal signal definitions ---
//=== parameter definitions ===			
                                
//=== reg definitions ===  		
reg        R_valid;		
reg [15:0] R_d1;
reg [15:0] R_d2; 

//=== wire definitions === 				
                                    
                                    
//--- Main body of code ---   
always @(posedge I_sys_clk or negedge I_reset_n)
begin
    if(~I_reset_n) 
    begin
        R_valid <= 1'b0;
        R_d1 <= 16'd0;
        R_d2 <= 16'd0;
    end 
    else
    begin
        R_valid <= I_valid;
        R_d1 <= I_R;
        R_d2 <= I_R - I_Di;
    end
end

always @(posedge I_sys_clk or negedge I_reset_n)
begin
    if(~I_reset_n) 
    begin
        O_valid <= 1'b0;
        O_Q <= 16'd0;
        O_R <= 16'd0;
    end 
    else
    begin
        O_valid <= R_valid;
        if (R_valid)
        begin
            O_Q <= ~R_d2[15];
            O_R <= R_d2[15] ? R_d1 : R_d2;
        end
        else 
        begin
            O_Q <= 16'd0;
            O_R <= 16'd0;
        end
    end
end

endmodule

仿真结果：

2 基于不恢复余数（Non-Restoring）算法的除法器

2.1 原理

恢复余数第k位迭代：
$r_k = r_{k+1}-q_k·D·2^k$
若 $r_k<0$ 表明 $q_k=0$，使得 $ r_k = r_{k+1} > 0 $ ，下一位（k-1位）迭代变为：
$r_{k-1} = r_{k}-D·2^{k-1} = r_{K+1}-D·2^{k-1}$

如果不恢复余数，$r_k$ 保留为负值$r_{k+1}-D·2^k $，下一位（k-1位）迭代为：
$r_{k-1} = r_{k}-D·2^{k-1} = r_{K+1}-D·2^{k}-D·2^{k-1}= r_{K+1}-3·D·2^{k-1}$

上式的结果是错的，如果将k-1位迭代改为：
$r_{k-1} = r_{k}+D·2^{k-1} = r_{K+1}-D·2^{k}+D·2^{k-1}= r_{K+1}-D·2^{k-1}$

这个结果是对的了。所以 Non-Restoring 算法流程为：

Non-Restoring 算法的相比 Restoring 算法的特征：

1. $q_i$ 取决于 $r_i$ 的正负。而不是预先假定为1，再根据 $r_i$
2. $r_i$
3. $r_0$ 的正负决定是否需要余数校正。若为负，最终余数为 $r_0+D$。

基本迭代单元硬件架构：

以 33/7 为例，Non-Restoring 算法迭代过程：

k	余数
3	33
2	33-1×7×2^2=5
1	5-1×7×2^1=-9
0	-9+1×7×2^1=-2
最终计算结果为 Q=3’b100=4，余数为-2。需要对余数进行校正，最终余数为 。

2.1 Verilog 实现

注意：实现有待优化，复杂的组合逻辑需要插入pipeline才能用于高速处理。

`timescale 1ns/1ps


module divider (
    input             I_sys_clk  ,
    input             I_reset_n  ,
    input             I_valid    ,
    input      [15:0] I_dividend ,
    input      [7:0]  I_divisor  ,
    output reg        O_valid    ,
    output reg [7:0]  O_quotient ,
    output reg [7:0]  O_remainder
);

//--- internal signal definitions ---
//=== parameter definitions ===			
                                  
//=== reg definitions ===  				
reg [3:0]  R_index_cnt;
reg        R_valid    ;
reg [15:0] R_remainder;
reg [7:0]  R_divisor  ;
reg [15:0] R_remainder_t;
reg [7:0]  R_quotient_t;

//=== wire definitions === 				                                  

                                    
//--- Main body of code ---  
always @(posedge I_sys_clk or negedge I_reset_n)
begin
    if(~I_reset_n) 
    begin
        R_valid     <= 1'b0;
        R_remainder <= 16'd0;
        R_divisor   <= 8'd0;
    end 
    else
    begin
        R_valid     <= I_valid;
        if (I_valid)
        begin
            R_remainder <= I_dividend;
            R_divisor   <= I_divisor ;
        end
    end
end

always @(posedge I_sys_clk or negedge I_reset_n)
begin
    if(~I_reset_n) 
    begin
        R_index_cnt <= 4'd0;
    end 
    else
    begin
        if (R_valid)
        begin
            R_index_cnt <= 4'd8;
        end
        else if (|R_index_cnt)
        begin
            R_index_cnt <= R_index_cnt - 4'd1;
        end
    end
end

always @(posedge I_sys_clk or negedge I_reset_n)
begin
    if(~I_reset_n) 
    begin
        R_remainder_t <= 16'd0;
    end 
    else
    begin
        if (I_valid)
        begin
            R_remainder_t <= I_dividend;
        end
        else if (|R_index_cnt)
        begin
            R_remainder_t <= R_remainder_t[15] ? (R_remainder_t + (R_divisor << (R_index_cnt-2))) : (R_remainder_t - (R_divisor << (R_index_cnt-2))) ;
        end
    end
end

always @(posedge I_sys_clk or negedge I_reset_n)
begin
    if(~I_reset_n) 
    begin
        R_quotient_t <= 8'd0;
    end 
    else
    begin
        if (I_valid)
        begin
            R_quotient_t <= 8'd0;
        end
        else if (|R_index_cnt && (R_index_cnt <= 4'd7))
        begin
            R_quotient_t <= {R_quotient_t[6:0], ~R_remainder_t[15]};
        end
    end
end


always @(posedge I_sys_clk or negedge I_reset_n)
begin
    if(~I_reset_n) 
    begin
        O_valid     <= 1'd0;
        O_quotient  <= 8'd0;
        O_remainder <= 16'd0;
    end 
    else
    begin
        if (R_index_cnt == 4'd1)
        begin
            O_valid     <= 1'd1;
            O_quotient  <= {R_quotient_t[6:0], ~R_remainder_t[15]};
            O_remainder <= R_remainder_t[15] ? (R_remainder_t + R_divisor) : R_remainder_t;
        end
        else 
        begin
            O_valid     <= 1'd0;
            O_quotient  <= 8'd0;
            O_remainder <= 16'd0;
        end
    end
end

endmodule

仿真结果：

2.2 用于有符号数的改进

3 基于级数展开算法的除法器

将D归一化为[0.5,1)，定义：
$f(x)=\frac{1}{D}=\frac{1}{1+x}$
根据Taylor级数展开：
$f(x)=\frac{1}{1+x}=1-x+x^2-x^3+x^4-...$
进一步改写为：
$\frac{1}{D}=(1-x)(1+x^2)(1+x^4)(1+x^8)(1+x^{16})(1+x^{32})...$

又有：
$1+x^n=2-(1-x^n)$

则 1/D 的求解过程为：

1. 制作查找表，存储$(1-x)(1+x^2)(1+x^4)$
2. $1-x^8=[(1-x)(1+x^2)(1+x^4)](1+x)$，进行1次乘法运算
3. $1+x^8=2-(1-x^8)$，进行一次加法运算
4. $1-x^{16}=(1-x^8)(1+x^8)$，进行1次乘法运算
5. $1+x^{16}=2-(1-x^{16})$，进行一次加法运算
6. $1-x^{32}=(1-x^{16})(1+x^{16})$，进行1次乘法运算
7. $1+x^{32}=2-(1-x^{32})$，进行一次加法运算
   … 直到达到需要的精度。

步骤1查找表的指标：深度和宽度。宽度取决于数据位宽，即系统精度。深度取决于地址位宽，即除数D的位宽。为节省内存可取D的高m-1位。

D为[0.5,1)，则符号位S为0。

若要求最大误差为 $\epsilon$

$m \ge log_2(\frac{1}{\epsilon})+2$

因式求解需 3次乘法和 3次加法。1/D求解中因式相乘需3次乘法。最终 Y/D 1次乘法。共7次乘法，3次加法。

4.2 Verilog 实现

摆烂了 😃

4 基于 Newton-Raphson 算法的除法器

4.1 原理

基本思想：将非线性方程线性化，以线性方程的解逼近非线性方程的解。
构造：
$f(x)=\frac{1}{x}-D$

求解其解1/D，迭代过程：

构造函数，然后基于 Newton-Raphson 算法求解，与级数展开算法本质一致，都是求1/D，只是出发点不同。