码迷,mamicode.com
首页 > 编程语言 > 详细

大整数算法[08] Comba乘法(原理)

时间:2015-04-12 19:05:44      阅读:292      评论:0      收藏:0      [点我收藏+]

标签:

        引子

         原本打算一篇文章讲完,后来发现篇幅会很大,所以拆成两部分,先讲原理,再讲实现。实现的话相对复杂,要用到内联汇编,要考虑不同平台等等。

         在大整数计算中,乘法是非常重要的,因为在公钥密码学中模幂运算要频繁使用乘法,所以乘法的性能会直接影响到模幂运算的效率。下面将会介绍两种乘法:基线乘法和 Comba 乘法,尽管他们的原理和计算看起来十分类似,而且算法的时间复杂度都是 O(n^2),但是他们的效率差别是很大的。

 

        基线乘法 (Baseline Multiplication)

         基线乘法是 Bignum Math 一书中的一个术语,这里直接拿来用。实际上基线乘法就是通常所说的笔算乘法。先来看看计算 123 * 456 用笔算乘法应怎么做:

                                                               1          2          3

                                                      x      4           5          6

                            -----------------------------------------------

                                                             7           3           8

                                                6           1           5

                                  4           9            2

                            ------------------------------------------------

                                  5            6            0           8           8

          在笔算算法中,第二个数的每一位分别与第一个数的每一个位相乘,并且将进位传递到下一位,比如 3 * 6 = 18 保留本位 8,进位 1 传递到下一位。按照笔算的原理,很容易得到基线乘法的算法思路:

          设大整数 x 和 y 分别有 x->used 和 y->used 个数位,进位 c,相乘结果放在 z 中,算法的步骤是:

 

          1. 对于下标 i,从 0 到 x->used + y->used - 1 循环,将 z 中每一个数位清 0。其实就是把 z 设为 0。

          2. 对于下标 i,从 0 到 y->used - 1 循环,执行如下几个操作:

                2.1  c = 0

                2.2  对于下标 j,从 0 到 x->used - 1循环计算:(注:i,j 都是从 0 开始的,r 是一个双精度变量)

                         r  = z(i + j) + x(i) * y(j) + c。其中 z(i + j),x(i),y(j) 分别表示 z 的第 i + j + 1个数位,x 的第 i + 1 个数位,y 的第 j +1 个数位。

                         z(i + j) = r & BN_MASK          //取 r 的低半部分作为 z(i + j) 的本位。

                         c = r >> biL                                  //取 r 的高半部分作为进位 c。

                 2.3  z(i + x->used) = c

             3. 返回结果 z。

 

             根据上面的算法思路,便可以写出基线乘法,C 代码如下:

int bn_baseline_mul(bignum * z, const bignum *x, const bignum *y)
{
    int ret;
    bn_udbl r;
    size_t i, j, c;
    size_t *px, *py, *pz;
    bignum ta[1], tb[1];

    bn_init(ta);
    bn_init(tb);

    if(x == z)
    {
        BN_CHECK(bn_copy(ta, x));
        x = ta;
    }
    if(y == z)
    {
        BN_CHECK(bn_copy(tb, y));
        y = tb;
    }

    BN_CHECK(bn_grow(z, x->used + y->used));
    BN_CHECK(bn_set_word(z, 0));
    z->used = x->used + y->used;

    for(i = 0; i < y->used; i++)
    {
        c = 0;

        px = x->dp;
        py = y->dp + i;
        pz = z->dp + i;

        for(j = 0; j < x->used; j++)
        {
            r = (*pz) + (bn_udbl)(*px++) * (*py) + c;
            *pz++ = (bn_digit)r;
            c = (bn_digit)(r >> biL);
        }

        *pz = c;
    }

    z->sign = x->sign * y->sign;
    bn_clamp(z);

clean:
    bn_free(ta);
    bn_free(tb);

    return ret;
}

          上面的代码中 ta,tb 临时变量是为了处理类似 x = x * y,y = x * y,x = x * x 之类的输入和输出是同一个变量的情况,因为在计算中,z 会先被置 0,如果不用临时变量存储 x 或 y,则 x 或 y 的值就会被改变。要注意的是这个乘法并没有考虑到没有双精度变量的情况,如果在 64 位的环境下,没有 128 bit 的双精度类型,则处理起来要复杂很多,这里暂时不讲,下一篇会慢慢讨论。

 

         Comba 乘法 (Comba  Multiplication)原理

        Comba 乘法以(在密码学方面)不太出名的 Paul G. Comba 得名。上面的笔算乘法,虽然比较简单, 但是有个很大的问题:在 O(n^2) 的复杂度上进行计算和向上传递进位,看看前面的那个竖式,每计算一次单精度乘法都要计算和传递进位,这样的话就使得嵌套循环的顺序性很强,难以并行展开和实现。Comba 乘法则无需进行嵌套进位来计算乘法,所以虽然其时间复杂度和基线乘法一样,但是速度会快很多。还是以计算 123 * 456 为例:

                                                                  1            2            3

                                                       x        4             5            6

                                -----------------------------------------------

                                                                 6           12          18

                                       5          10         15

                             4           8          12

                               ------------------------------------------------

向上传递进位:    4          13          28         27         18

                                       4          13          28         28           8

                                       4          13          30          8     

                             4       16         0

                                       5          6

                            0         5

                       ------------------------------------------------------

                                        5           6             0            8             8

            和普通的笔算乘法很类似,只是每一次单精度乘法只是单纯计算乘法,不计算进位,进位留到每一列累加后进行。所以原来需要 n * n 次进位,现在最多只需要 2n 次即可。

            以上就是 Comba 乘法的原理,不过这里有个比较严重的问题:如何保证累加后结果不溢出。上面的例子,假设单精度数 1  位数,双精度是两位数,那万一累加后的结果超过两位数则么办?那没办法,只能用三精度变量了。在大整数算法中,单精度能表示的最大整数是 2^n - 1(n 是单精度变量的比特数),用三个单精度变量 c2,c1,c0 连在一起作为一个三精度变量(高位在左,低位在右),则 c2 || c1 || c0 能表示的最大整数是 2^(3n) - 1,最多能存放 (2^(3n) - 1) / ((2^n - 1)^2) 个单精度乘积结果。当 n = 32 时,能够存放 4294967298 个单精度乘积结果;当 n = 64 时,能够存放约 1.845 * 10^19 个单精度乘积结果,而我一开始规定 bignum 不能超过 25600 个数位,这样使用三精度变量就可以保证累加结果不会溢出了。

              有了上面的铺垫,下面就把 Comba 乘法的思路列出来:

              计算 c = a * b,c0,c1,c2 为单精度变量。

 

              1. 增加 c 到所需要的精度,并且令 c = 0,c->used = a->used + b->used。

              2. 令 c2 = c1 = c0 = 0。

              3. 对于 i 从 0 到 c->used - 1 循环,执行如下操作:

                   3.1  ty = BN_MIN(i, b->used - 1)

                   3.2  tx = i - ty

                   3.3  k = BN_MIN(a->used - tx, ty + 1)

                   3.4  三精度变量右移一个数位:(c2 || c1 || c0)  =  (0 || c2 || c1)

                   3.5  对于 j 从 0 到 k - 1 之间执行循环,计算:

                            (c2 || c1 || c0) = (c2 || c1 || c0) + a(tx + j) * b(ty - j)

                   3.6  c(i) = c0

                4. 压缩多余位,返回 c。      

                

                BN_MIN 是宏定义:#define BN_MIN(x, y)                  (((x) < (y)) ? (x) : (y))          

 

                以上就是 Comba 乘法的思路,先计算每一列,然后求和累加,将累加结果求余数得到本位,进位传递到下一列。

                第三步中,分别计算 tx, ty 和 k 的值,用于控制列的计算,光这么说比较抽象,所以还是举个例子比较直观。假设要计算 12345 * 678。

 

         下标 i 的值:                                                 6            5               4                 3                2                 1               0

                                                        --------------------------------------------------------------------------------------------------

                                                                                                                       1                 2                3                 4               5

                                                                                                 x                                                          6                 7               8

                                                        --------------------------------------------------------------------------------------------------

                                                                                                                        8               16              24              32            40

                                                                                                       7             14              21              28              35

                                                                                        6           12            18             24              30

             来自上一列累加产生的进位:          2             4              6                8                7                 4               0

                                                       --------------------------------------------------------------------------------------------------

           结果:                                                              8              3              6               9                9                 1               0  

            

           过程:a->used = 5, b->used = 3

           i = 0,ty = 0,tx = i - ty = 0,a->used - tx = 5,ty + 1 = 1,k = 1,循环计算第一列:5 * 8 = 40

           i = 1,ty = 1,tx = i - ty = 0,a->used - tx = 5,ty + 1 = 2,k = 2,循环计算第二列:5 * 7 = 35,4 * 8 = 32

           i = 2,ty = 2,tx = i - ty = 0,a->used - tx = 5,ty + 1 = 3,k = 3,循环计算第三列:5 * 6 = 30,4 * 7 = 28,3 * 8 = 24

           i = 3,ty = 2,tx = i - ty = 1,a->used - tx = 4,ty + 1 = 3,k = 3,循环计算第四列:4 * 6 = 24,3 * 7 = 21,2 * 8 = 16

           i = 4,ty = 2,tx = i - ty = 2,a->used - tx = 3,ty + 1 = 3,k = 3,循环计算第五列:3 * 6 = 18,2 * 7 = 14,1 * 8 = 8

           i = 5,ty = 2,tx = i - ty = 3,a->used - tx = 2,ty + 1 = 3,k = 2,循环计算第六列:2 * 6 = 12,1 * 7 = 7

           i = 6,ty = 2,tx = i - ty = 4,a->used - tx = 1,ty + 1 = 3,k = 1,循环计算第七列:1 * 6

 

           可以看出,Comba 乘法每一次都是计算一列,所以如果想进行并行计算的话,那会方便很多:同时计算所有列并累加,最后进行一次进位传递。

           前面提到基线乘法和 Comba 乘法的时间复杂度是一样的,这是对于单精度乘法而言,Comba 乘法之所以要比基线乘法快,主要是因为减少了进位传递的次数,所以减少了加法的计算量。

 

        总结

         篇幅好像有点长了,所以暂时写到这里。说了这么多,其实关键就是:都是使用笔算乘法,基线乘法计算每一行的结果,在累加的同时计算进位,Comba 乘法则按列计算,先累加,然后传递进位,减少了计算量,所以速度快。这次先把 Comba 乘法的原理讲清楚,下一篇文章讲讲讲如何实现 Comba 乘法,包括在没有双精度的环境下和在 x86 环境下使用内联汇编优化速度。

 

 

   【回到本系列目录】 

 

 

版权声明
原创博文,转载必须包含本声明,保持本文完整,并以超链接形式注明作者Starrybird和本文原始地址:http://www.cnblogs.com/starrybird/p/4419444.html

 

大整数算法[08] Comba乘法(原理)

标签:

原文地址:http://www.cnblogs.com/starrybird/p/4419444.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!