大整數算法[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教程