zoukankan      html  css  js  c++  java
  • 大整数算法[09] Comba乘法(原理)

            引子

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

             在大整数计算中,乘法是非常重要的,因为在公钥密码学中模幂运算要频繁使用乘法,所以乘法的性能会直接影响到模幂运算的效率。下面将会介绍两种乘法:基线乘法和 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(j) * y(i) + c。其中 z(i + j),x(j),y(i) 分别表示 z 的第 i + j + 1个数位,x 的第 j + 1 个数位,y 的第 i +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。

             Index 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 = 6

     

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

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

            总结

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

       【回到本系列目录】 

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

     

  • 相关阅读:
    WPF 关于拖拽打开文件的注意事项
    asp.net core 3.1中对Mongodb BsonDocument的序列化和反序列化支持
    用百度webuploader分片上传大文件
    多线程学习笔记
    web.config数据库连接字符串加密
    Visual Studio 2010 常用快捷方式
    Team Foundation Server 2013 日常使用使用手册(四)分支与合并
    Team Foundation Server 2013 日常使用使用手册(三)上传新工程、创建任务、创建bug、设置预警
    Team Foundation Server 2013 日常使用使用手册(二)修改、签入、撤销、回滚、对比代码变更
    Team Foundation Server 2013 日常使用使用手册(一)-本地连接TFS、查看任务
  • 原文地址:https://www.cnblogs.com/starrybird/p/4419444.html
Copyright © 2011-2022 走看看