大整数算法[13] 单数位乘法

时间  2019-11-11
标签 整数 算法 单数 乘法 繁體版
原文   原文链接

        ★ 引子html

        最近在折腾 wxWidgets，同时拖延症又犯了，因此中断了很久。此次来说讲单数位乘法，前面讲到 Comba 和 Karatsuba 乘法，这两个算法适合用来处理比较大的整数，可是对于一个大整数和一个单精度数相乘，其效果反而会很差，由于计算量过多。实际上单数位乘法只是基线乘法的一个特例，不存在嵌套循环进位，所以能够经过优化减小计算量。另外与完整的乘法不一样的是，单数位乘法不须要什么临时变量存储和内存分配（目标精度增长除外）。git

 

        ★ 算法思路算法

        单数位乘法相似于计算 1234567890 * 8 这种计算，第二个数只有一位，在大整数中就是一个单精度变量，只须要执行 O(n) 次单精度乘法就能够完成主要的计算。每一次单精度乘法计算完后，执行进位传递。具体的实现思路以下：
优化

      

         计算 z = x * y，其中 z 和 x 是 bignum，y 是无符号的单精度数。spa

         1. olduse = z->used       记录当前 z 使用了多少数位，用于辅助后面的高位清零。htm

         2. 目标精度增长 1 ： z->used = x->used +1blog

         3. z->sign = x->sign （由于 y 为无符号数，因此结果符号只和 x 有关）ip

         4. 初始化进位： u = 0内存

         5. 对于 i 从 0 到 x->used - 1 之间进行循环：get

              5.1  r = u + x(i) * y                //r 是双精度变量

              5.2  z(i) = r & BN_MASK0     //取低半部分做为本位  （32 bit 下，BN_MASK0 = 0xFFFFFFFF，其余状况同理，为 2^n - 1）

              5.3  u = r >> biL                  //取高半部分做为进位

          6. z(x->used) = u                      //传递最后一个进位

          7. 剩余的高位清零

          8. 压缩多余位

 

          这里 y 是无符号数，若是想计算和 -y 的乘积，在计算完后将 z 的符号取反便可。

 

        ★ 实现

         和 Comba 乘法同样，先给出整体实现，具体细节后面讲。考虑到计算效率和可移植性的问题，第五步的循环关键代码仍是写在宏定义里面，而后按照 1，4，8，16 和 32 的步进展开乘法器，减小循环控制的开销。

int bn_mul_word(bignum *z, const bignum *x, const bn_digit y)
{
    int ret;
    size_t i, olduse;
    bn_digit u, *px, *pz;

    olduse = z->used;
    z->sign = x->sign;
    BN_CHECK(bn_grow(z, x->used + 1));

    u = 0;
    px = x->dp;
    pz = z->dp;

    for(i = x->used; i >= 32; i -= 32)
    {
        MULADDC_WORD_INIT
        MULADDC_WORD_CORE    MULADDC_WORD_CORE
        MULADDC_WORD_CORE    MULADDC_WORD_CORE
        MULADDC_WORD_CORE    MULADDC_WORD_CORE
        MULADDC_WORD_CORE    MULADDC_WORD_CORE

        MULADDC_WORD_CORE    MULADDC_WORD_CORE
        MULADDC_WORD_CORE    MULADDC_WORD_CORE
        MULADDC_WORD_CORE    MULADDC_WORD_CORE
        MULADDC_WORD_CORE    MULADDC_WORD_CORE

        MULADDC_WORD_CORE    MULADDC_WORD_CORE
        MULADDC_WORD_CORE    MULADDC_WORD_CORE
        MULADDC_WORD_CORE    MULADDC_WORD_CORE
        MULADDC_WORD_CORE    MULADDC_WORD_CORE

        MULADDC_WORD_CORE    MULADDC_WORD_CORE
        MULADDC_WORD_CORE    MULADDC_WORD_CORE
        MULADDC_WORD_CORE    MULADDC_WORD_CORE
        MULADDC_WORD_CORE    MULADDC_WORD_CORE
        MULADDC_WORD_STOP
    }
    for(; i >= 16; i -= 16)
    {
        MULADDC_WORD_INIT
        MULADDC_WORD_CORE    MULADDC_WORD_CORE
        MULADDC_WORD_CORE    MULADDC_WORD_CORE
        MULADDC_WORD_CORE    MULADDC_WORD_CORE
        MULADDC_WORD_CORE    MULADDC_WORD_CORE

        MULADDC_WORD_CORE    MULADDC_WORD_CORE
        MULADDC_WORD_CORE    MULADDC_WORD_CORE
        MULADDC_WORD_CORE    MULADDC_WORD_CORE
        MULADDC_WORD_CORE    MULADDC_WORD_CORE
        MULADDC_WORD_STOP
    }
    for(; i >= 8; i -= 8)
    {
        MULADDC_WORD_INIT
        MULADDC_WORD_CORE    MULADDC_WORD_CORE
        MULADDC_WORD_CORE    MULADDC_WORD_CORE
        MULADDC_WORD_CORE    MULADDC_WORD_CORE
        MULADDC_WORD_CORE    MULADDC_WORD_CORE
        MULADDC_WORD_STOP
    }
    for(; i >= 4; i -= 4)
    {
        MULADDC_WORD_INIT
        MULADDC_WORD_CORE    MULADDC_WORD_CORE
        MULADDC_WORD_CORE    MULADDC_WORD_CORE
        MULADDC_WORD_STOP
    }
    for(; i > 0; i--)
    {
        MULADDC_WORD_INIT
        MULADDC_WORD_CORE
        MULADDC_WORD_STOP
    }

    *pz++ = u;

    for(i = x->used + 1; i < olduse; i++)
        *pz++ = 0;

    z->used = x->used + 1;
    bn_clamp(z);

clean:

    return ret;
}

         

          以上是单数位乘法的整体实现，关键的地方都在宏定义中，下面将讲讲不一样环境下的实现方式。

 

        ★ 单双精度变量都有的状况

          此状况下， bn_digit 和 bn_udbl 同时有定义，最容易实现的。三个宏的定义以下：

#define MULADDC_WORD_INIT                          \
{                                                  \
    bn_udbl r;                                     \

#define MULADDC_WORD_CORE                          \
	                                               \
    r = u + (bn_udbl)(*px++) * y;                  \
    *pz++ = (bn_digit)r;                           \
    u = (bn_digit)(r >> biL);                      \

#define MULADDC_WORD_STOP                          \
}
 

           这种状况彻底是按照思路实现的，具体原理就很少说了。

 

        ★ 只有单精度变量的状况

        若是遇到这种状况，则 bn_udbl 无定义，单精度乘法须要转换成 4 个 半精度的乘法来计算，相对比较复杂。具体的实现原理和 Comba 的乘法器相似，参考此处：http://www.cnblogs.com/starrybird/p/4441022.html 。 具体的宏实现以下：

#define MULADDC_WORD_INIT                           \
{                                                   \
    bn_digit a0, a1, b0, b1;                        \
    bn_digit t0, t1, r0, r1;                        \

#define MULADDC_WORD_CORE                           \
	                                                \
    a0 = (*px << biLH) >> biLH;                     \
    b0 = (  y << biLH) >> biLH;                     \
    a1 = *px++ >> biLH;                             \
    b1 =  y    >> biLH;                             \
    r0 = a0 * b0;                                   \
    r1 = a1 * b1;                                   \
    t0 = a1 * b0;                                   \
    t1 = a0 * b1;                                   \
    r1 += (t0 >> biLH);                             \
    r1 += (t1 >> biLH);                             \
    t0 <<= biLH;                                    \
    t1 <<= biLH;                                    \
    r0 += t0;                                       \
    r1 += (r0 < t0);                                \
    r0 += t1;                                       \
    r1 += (r0 < t1);                                \
    r0 += u;                                        \
    r1 += (r0 < u);                                 \
    *pz++ = r0;                                     \
    u = r1;                                         \

#define MULADDC_WORD_STOP                           \
}

 

        ★ 使用内联汇编的状况

        C 的内联汇编细节就很少说了，若是你不会能够跳过。

 

        VC x86:

#define MULADDC_WORD_INIT                           \
{                                                   \
    __asm   mov   esi, px                           \
    __asm   mov   edi, pz                           \
    __asm   mov   ecx, u                            \

#define MULADDC_WORD_CORE                           \
	                                                \
    __asm   lodsd                                   \
    __asm   mul   y                                 \
    __asm   add   eax, ecx                          \
    __asm   adc   edx, 0                            \
    __asm   mov   ecx, edx                          \
    __asm   stosd                                   \

#define MULADDC_WORD_STOP                           \
	                                                \
    __asm   mov   px,  esi                          \
    __asm   mov   pz,  edi                          \
    __asm   mov   u,   ecx                          \
}

#endif

 

        GCC x86:

#define MULADDC_WORD_INIT                           \
{                                                   \
    asm                                             \
    (                                               \
       "movl %3, %%esi       \n\t"                  \
       "movl %4, %%edi       \n\t"                  \
       "movl %5, %%ecx       \n\t"                  \

#define MULADDC_WORD_CORE                           \
	                                                \
       "lodsl                \n\t"                  \
       "mull %6              \n\t"                  \
       "addl %%ecx, %%eax    \n\t"                  \
       "adcl $0, %%edx       \n\t"                  \
       "movl %%edx, %%ecx    \n\t"                  \
       "stosl                \n\t"                  \

#define MULADDC_WORD_STOP                           \
                                                    \
       "movl %%esi, %0       \n\t"                  \
       "movl %%edi, %1       \n\t"                  \
       "movl %%ecx, %2       \n\t"                  \
       :"=m"(px),"=m"(pz),"=m"(u)                   \
       :"m"(px),"m"(pz),"m"(u),"m"(y)               \
       :"%eax","%ecx","%edx","%esi","%edi"          \
    );                                              \
}

 

        GCC x64:

#define MULADDC_WORD_INIT                           \
{                                                   \
    asm                                             \
    (                                               \
       "movq %3, %%rsi       \n\t"                  \
       "movq %4, %%rdi       \n\t"                  \
       "movq %5, %%rcx       \n\t"                  \

#define MULADDC_WORD_CORE                           \
	                                                \
       "lodsq                \n\t"                  \
       "mulq %6              \n\t"                  \
       "addq %%rcx, %%rax    \n\t"                  \
       "adcq $0, %%rdx       \n\t"                  \
       "movq %%rdx, %%rcx    \n\t"                  \
       "stosq                \n\t"                  \

#define MULADDC_WORD_STOP                           \
                                                    \
       "movq %%rsi, %0       \n\t"                  \
       "movq %%rdi, %1       \n\t"                  \
       "movq %%rcx, %2       \n\t"                  \
       :"=m"(px),"=m"(pz),"=m"(u)                   \
       :"m"(px),"m"(pz),"m"(u),"m"(y)               \
       :"%rax","%rcx","%rdx","%rsi","%rdi"          \
    );                                              \
}

 

        ★ 总结

        算法也很简单，按照 Baseline Multiplication 的方法作就好了，注意关键的地方优化一下。下一篇讲讲平方的计算。        

 

   【回到本系列目录】 

 

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

相关文章
相关标签/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公众号
   欢迎关注本站公众号,获取更多信息
联系我们 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程