基于Verilog的Mitchell近似算法

---

参考文献：
[1] 苏剑林. (Dec. 14, 2020). 《Mitchell近似：乘法变为加法，误差不超过1/9 》
[2] 《Computer Multiplication and Division Using Binary Logarithms》

前言

  本文章基于1962年的论文《Computer Multiplication and Division Using Binary Logarithms》，作者是John N. Mitchell，该算法在二进制下，可以通过加法来近似两个数的乘法，最大误差不超过1/9。本文前文引用来自参考文献[1]。

---

一、快速对数与指数

  说到乘法变加法，大家应该很自然能联想到对数和指数，即
$$pq=a^s,s=lo{g}_{a}p+lo{g}_{a}q\text{\hspace{1em}(1)}$$
  本文是基于二进制的，所以$a=2$。问题在于上式虽然确实将乘法转为加法了，但是对数$lo{g}_{2}p$, $lo{g}_{2}q$和转换后的指数$2s$算起来都不是一件容易的事情。所以要利用这个等式做乘法，关键在于要实现快速对数和指数运算。
  对于十进制的非负数p，我们假设它的二进制表示为
$$z_n{z}_{n-1}\cdots z_1{z}_0.z_{-1}{z}_{-2}\cdots z_{-(m-1)}{z}_{-m}\text{\hspace{1em}(2)}$$

  其中$z_n=1$且各个$z_i\in 0,1$，那么我们就有
$$p=2^n+\sum _{i=-m}^{n-1}{z}_i{2}_i=2^n(1+\sum _{i=-m}^{n-1}{z}_i{2}^{i-n})\text{\hspace{1em}(3)}$$

  记$x={\sum }_{i=-m}^{n-1}{z}_i{2}^{i-n}$，我们得到
$$lo{g}_{2}p=n+lo{g}_2(1+x)\text{\hspace{1em}(4)}$$

  在这里，Mitchell做了一个相当果断而美妙的近似，那就是log2(1+x)≈x（此处会带来误差，后续会进行分析），于是得到
$$lo{g}_{2}p\approx n+x\text{\hspace{1em}(5)}$$

  这个结果妙在何处呢？首先n是一个整数，等于p的二进制的整数部分的位数减去1，它转换为二进制自然自然也是整数；那x呢？根据x的定义我们不难看出，实际上x的二进制表示就是
$$0.z_{n-1}\cdots z_1{z}_0{z}_{-1}{z}_{-2}\cdots z_{-(m-1)}{z}_{-m}\text{\hspace{1em}(6)}$$

  也就是说x其实就是上述近似的小数部分，它的二进制表示只不过是p的二进制表示的简单变形（小数点平移）。
  综上所述，我们得到对数的Mitchell近似算法：

1、输入十进制的p；
2、将p转换为二进制数znzn−1⋯z1z0.z−1z−2⋯z−(m−1)z−m，这里zn=1；
3、将n转换为二进制数ykyk−1⋯y1y0；
4、那么log2p的二进制近似为ykyk−1⋯y1y0.zn−1⋯z1z0z−1z−2⋯z−(m−1)z−m。

  将上述过程逆过来，就得到了指数的Mitchell近似算法：

1、输入二进制的znzn−1⋯z1z0.z−1z−2⋯z−(m−1)z−m；
2、将znzn−1⋯z1z0转换为十进制数n；
3、那么它的指数的二进制近似为1z−1z−2⋯z−(n−1)z−n.z−(n+1)z−(n+2)⋯z−(m−1)z−m；
4、将上述结果转换为十进制。

  所以，在二进制下对数和指数的近似计算就只是数数和拼接操作而已！

---

二、近似乘法过程以及例子

  假设x和y为输入的两个乘法，z=x*y。对x和y进行近似可得
$$x=2^{k_1}(1+m_1)->lo{g}_{2}x=k_1+lo{g}_2(1+m_1)->lo{g}_{2}x\approx k_1+m_1$$
$$y=2^{k_2}(1+m_2)->lo{g}_{2}y=k_2+lo{g}_2(1+m_2)->lo{g}_{2}y\approx k_2+m_2$$
$$lo{g}_{2}z=lo{g}_{2}xy=lo{g}_{2}x+lo{g}_{2}y\approx k_1+m_1+k_2+m_2$$
  对$lo{g}_{2}z$进行逆变换可得
$$z=2^{k_1+k_2}(1+m_1+m_2)$$
  以上便是近似算法的全过程，无非就是先将乘数进行对数相加，相加再进行逆变换恢复结果。

  有了快速的（近似的）对数和指数算法，我们就可以乘法了，现在我们就用这个思路来算一个具体例子，由此加深我们对上述过程理解。

  我们要算的是$12.3\times 4.56$，计算流程如下表（近似指数是对求和的结果算的）：

p,q / 十进制	12.3	4.56
p,q / 二进制	1100.010011	100.1000111
n / 十进制	3	2
n / 二进制	11	10
x / 二进制	0.100010011	0.001000111
n+x / 二进制	11.10001001	10.001000111
求和 / 二进制	101.10101101
n / 二进制	101
n / 十进制	5
x / 二进制	0.10101101
近似指数 / 二进制	110101.101
近似指数 / 十进制	53.625
精确乘积 / 十进制	56.088
相对误差 / 十进制	4.39%

  其中p,q二进制表示为无限循环小数，这里只截断了有限位，如果保留精确的二进制表示，那么最终结果是53.68，跟上表的53.625相差不大。可以看到，整个过程的主要计算量有两部分：求和、十进制与二进制之间的转换。然而，尽管十进制与二进制之间的转换对于我们人来说是计算量比较大的操作，但对于计算机来说，它本身就是二进制存储的，我们可以认为两者之间的转换计算量可以忽略不计；又或者说，只要我们还是以十进制为输出，那么不管什么算法这部分计算量都是必然存在的，因此我们不把它算进去算法的复杂度中。所以，整个算法的计算量就只有求和这一步，实现了将乘法（近似地）转化为加法运算.

---

四、最大误差分析

  标题写到，这个算法的误差不超1/9，现在我们就来证明这一点。证明需要全部转换到十进制来理解，Mitchell近似实际上是用了如下近似式：
$$lo{g}_2{2}^n(1+x)\approx n+x,2^{n+x}\approx 2^n(1+x)\text{\hspace{1em}(7)}$$
  其中n是整数，而x∈[0,1)，所以分析误差也是从这两个近似入手。

  假设两个数为$p=2^{n1}(1+x_1),q=2^{n2}(1+x_2)$，那么根据近似就有$lo{g}_{2}p+lo{g}_{2}q\approx n_1+n_2+x_1+x_2$，可见要分两种情况讨论：第一种情况$x_1+x_2<1$，那么近似指数的结果就是$2_{n1+n2}(1+x_1+x_2)$，因此近似程度就是
$$\frac{2^{n1+n2}(1+x_1+x_2)}{2^{n1}(1+x1)\times 2^{n2}(1+x2)}=\frac{1+x1+x2}{1+x1+x2+x1x2}\text{\hspace{1em}(8)}$$

  第二种情况$x_1+x_2\ge 1$，这时候近似指数的结果就是$2^{n1+n2+1}(x1+x2)$，因此近似程度就是
$$\frac{2^{n1+n2+1}(x_1+x_2)}{2^{n1}(1+x_1)\times 2^{n2}(1+x_2)}=\frac{2(x_1+x_2)}{1+x_1+x_2+x_1{x}_2}\text{\hspace{1em}(9)}$$

  可以按部就班地证明，上面两种情况的最小值，都是在$x_1=x_2=0.5$时取到，其结果为8/9，所以最大的相对误差为1/9（如果是除法变成减法，那么它的最大误差是12.5%）。因为按部就班的证明有点繁琐，我们这里就不重复了，而是直接用软件画出它的等高线图，看出误差最大的地方应该是中心处：

作图代码：

import numpy as np
import matplotlib.pyplot as plt

x = np.arange(0, 1, 0.001)
y = np.arange(0, 1, 0.001)

X, Y = np.meshgrid(x, y)
Z1 = (1 + X + Y) / (1 + X + Y + X * Y)
Z2 = 2 * (X + Y) / (1 + X + Y + X * Y)
Z = (X + Y < 1) * Z1 + (X + Y >= 1) * Z2
plt.figure(figsize=(7, 6))
contourf = plt.contourf(X, Y, Z)
plt.contour(X, Y, Z)
plt.colorbar(contourf)
plt.show()

  其实这个误差本质上取决于$lo{g}_2(1+x)\approx x$的近似程度，我们知道x是自然对数ln(1+x)的一阶泰勒展开式，而e=2.71828…更加接近于3，所以如果计算机使用3进制的话，本算法会有更高的平均精度。事实上，确实有一些理论分析表明，对于计算机来说其实最理想是e进制，而3比2更接近于e，所以三进制计算机在不少方面会比二进制更优，我国和前苏联都曾研究过三进制计算机，不过由于二进制实现起来更加容易，所以当前还是二进制计算机的天下了。

三、 如何处理负数

  上文讲述的是对于非负数十进制p和q，其实对于负数而言，该过程是相似的。总所周知，数据在计算机中是以补码的方式进行存储的。对于负数来说，可以转换成原码相乘，符号位再进行亦或即可。

---

五、Verilog实现

以下代码实现的是int8近似乘法。
代码如下（仅作参考）：

//This code implements Mitchell approximation calculation for signed numbers p and q.
module Mitchell(
    input [7:0] a,      //有符号数字
    input [7:0] b,  
    // output sign,    
    output [15:0] out
    );
    //reg define
    reg  [ 2:0] n_0,n_1;
    //wire define
    wire [ 7:0] x_0 , x_1;
    wire [15:0] n_x_0,n_x_1;  //n+x
    wire [15:0] sum;          //n_x_0 + n_x_1
    wire [ 7:0] sum_n;        //sum[15:8]
    wire [15:0] sum_x;        //sum[ 7:0]
    wire [ 7:0] p_c,q_c;
    wire [15:0] c;

    // complement signed numbers
    complement complement_m0(
        .in           (a        ),
        .out          (p_c      )
    );
    complement complement_m1(
        .in           (b        ),
        .out          (q_c      )
    );
    complement #(
        .WIDTH(16)
    )
    complement_m2(
        .in           (c        ),
        .out          (out      )
    );
   always @(*)begin
       if(p_c[6]) n_0 = 'd6;
       else if(p_c[5]) n_0 = 'd5;
       else if(p_c[4]) n_0 = 'd4;
       else if(p_c[3]) n_0 = 'd3;
       else if(p_c[2]) n_0 = 'd2;
       else if(p_c[1]) n_0 = 'd1;
       else if(p_c[0]) n_0 = 'd0;
       else
           n_0 = 'd7; 
   end
    always @(*)begin
        if(q_c[6]) n_1 = 'd6;
        else if(q_c[5]) n_1 = 'd5;
        else if(q_c[4]) n_1 = 'd4;
        else if(q_c[3]) n_1 = 'd3;
        else if(q_c[2]) n_1 = 'd2;
        else if(q_c[1]) n_1 = 'd1;
        else if(q_c[0]) n_1 = 'd0;
        else
            n_1 = 'd7; 
    end

    assign x_0   = {1'b0,p_c[6:0]}<<(8-n_0);
    assign x_1   = {1'b0,q_c[6:0]}<<(8-n_1);
    assign n_x_0 = {n_0,x_0};
    assign n_x_1 = {n_1,x_1};
    assign sum   = n_x_0 + n_x_1;
    assign sum_n = sum[15:8];
    assign sum_x = (sum_n>8)? {8'b0,sum[7:0]}<<(sum_n-8):{8'b0,sum[7:0]}>>(8-sum_n) ;
    assign c[14:0] = ((&n_0) | (&n_1))? 'd0:(sum_x | (16'd1 << (sum_n)));
    assign c[15]   =  ((&n_0) | (&n_1))? 1'b0:p_c[7]^q_c[7];
endmodule

module complement#(
    parameter WIDTH = 8
)
(
    input  [WIDTH-1:0] in,
    output [WIDTH-1:0] out
);
    assign out[WIDTH-2:0] = (in[WIDTH-1])?((~in[WIDTH-2:0]) + 1'b1):in[WIDTH-2:0];
    assign out[WIDTH-1] = in[WIDTH-1];
endmodule //complement

测试文件

module tb(
    );
    reg signed [7:0] p,q;
    wire signed [15:0] out;
    wire signed [15:0] ref;
    integer i;
    wire signed [15:0] diff;

    initial begin
        p=0;q=110;
        for(i=0;i<10;i=i+1)begin
            #30 p<=$random%128;
            q<=$random%128;
        end
        #10  p<=0; q<=24;
        #10  p<=1; q<=24;
    end
    assign ref = p*q;
    assign diff = (ref > out)?(ref - out): (out - ref);
    Mitchell Mitchell(
    .a         (p    ),
    .b         (q    ),
    .out       (out    )
    );
endmodule

---

六、仿真

  可见仿真结果和C语言得到的结果一致。

---

总结

  本文基于verilog实现了Mitchell近似算法，由于后续会对近似算法进行相关研究，而Mitchell近似是很多近似算法的基础，本文主要以实现算法逻辑为目的，体会其基本原理。