L04: fibdrv

主講人: jserv / 課程討論區: 2023 年系統軟體課程

返回「Linux 核心設計/實作」課程進度表

sysprog21/bignum 程式碼分析

sysprog21/bignum 在 API 的設計想法，區分

• 高階 bn API
• 低階 apm API

apm 表示 Arbitrary-Precision arithmetic (高精度運算)，利用數字陣列來進行大數運算，uint8_t 的最大值是 255，該如何表示 300 呢？

1. 令一個 uint8_t *digit 的數字陣列來儲存(同時也可以指定此陣列的 size，本例中假設 size = 3)
2. 已知
   $300=0\times {256}^2+1\times {256}^1+44\times {256}^0$
3. 則 300 可以用 digit[2] = 0,digit[1] = 1,digit[0] = 44 表示

視覺化表示就是

$$\begin{array}{rrrr}300:& 00000000& 00000001& 00101100\\ digit:& digit[2]& digit[1]& digit[0]\\ value:& 0& 1& 44\end{array}$$

最後要輸出要經過格式轉換，才能成為十進位，可見 format.c。

由於除以 2 可以看作是進行 binary code 右移，結合 clz / ctz 的指令搭配 fastdoubling (先去除以數字 MSB 起算的開頭 0 位元，因為真正的數字不包含 leading 0s)，最後呈現出 sysprog21/bignum 計算費氏數列的實作

    bn *a1 = fib; /* Use output param fib as a1 */

    bn_t a0, tmp, a;
    bn_init_u32(a0, 0); /*  a0 = 0 */
    bn_set_u32(a1, 1);  /*  a1 = 1 */
    bn_init(tmp);       /* tmp = 0 */
    bn_init(a);

    /* Start at second-highest bit set. */
    for (uint64_t k = ((uint64_t) 1) << (62 - __builtin_clzll(n)); k; k >>= 1) {
        /* Both ways use two squares, two adds, one multipy and one shift. */
        bn_lshift(a0, 1, a); /* a03 = a0 * 2 */
        bn_add(a, a1, a);    /*   ... + a1 */
        bn_sqr(a1, tmp);     /* tmp = a1^2 */
        bn_sqr(a0, a0);      /* a0 = a0 * a0 */
        bn_add(a0, tmp, a0); /*  ... + a1 * a1 */
        bn_mul(a1, a, a1);   /*  a1 = a1 * a */
        if (k & n) {
            bn_swap(a1, a0);    /*  a1 <-> a0 */
            bn_add(a0, a1, a1); /*  a1 += a0 */
        }
    }
    /* Now a1 (alias of output parameter fib) = F[n] */

Karatsuba 演算法

Karatsuba 的概念是將

$x$、

$y$ 以第

$n$ 位數為界，拆成兩半

$x_1$、

$x_0$、

$y_1$、

$y_0$，把這他們視為較小的數相乘，然後再透過左移補回

$x_1$、

$y_1$ 損失的位數，以十進位為例：

$$\begin{gathered} x=x_1\ast {10}^n+x_0 \\ y=y_1\ast {10}^n+y_0 \end{gathered}$$

則

$x\ast y$ 可以化為：

$\underset{z_2}{\underbrace{x_1{y}_1}}\ast {10}^{2n}+\underset{z_1}{\underbrace{(x_1{y}_0+y_1{x}_0)}}\ast {10}^n+\underset{z_0}{\underbrace{x_0{y}_0}}$

上述算法計算

$z_2$、

$z_1$、

$z_0$ 需要四次乘法，我們還可以透過以下技巧優化成三次乘法：

觀察

$(x_1+x_0)(y_1+y_0)$ 展開的結果

$(x_1+x_0)(y_1+y_0)=\underset{z_2}{\underbrace{x_1{y}_1}}+\underset{z_1}{\underbrace{x_1{y}_0+x_0{y}_1}}+\underset{z_0}{\underbrace{x_0{y}_0}}$

移項之後，我們就能利用

$(x_1+x_0)(y_1+y_0)$、

$z_0$、

$z_2$ 來計算

$z_1$

$z_2=x_1{y}_1$

$z_0=x_0{y}_0$

$z_1=(x_1+x_0)(y_1+y_0)-z_0-z_2$

最後計算

$z_2\ast {10}^{2n}+z_1\ast {10}^n+z_0$ 便能得到

$x$

$y$ 相乘的結果。

接著以 2 個 8 位元數值相乘變成 16 bit 為例，解說 Karatsuba 演算法（假設所採用的處理器只支援 8 位元乘法）：

令

$x={01001001}_2={73}_{10},y={10000011}_2={131}_{10}$

$x=x_1\ast 2^4+x_0=0100\ast 2^4+1001$

$y=y_1\ast 2^4+y_0=1000\ast 2^4+0011$

$z_2=x_1{y}_1=0100\ast 1000=00100000$

$z_0=x_0{y}_0=1001\ast 0011=00011011$

$\begin{array}{rl}{z}_1& =(x_1+x_0)(y_1+y_0)-z_0-z_2\\ & =(0100+1001)(1000+0011)-00100000-00011011\\ & =10001111-00100000-00011011\\ & =01010100\end{array}$

$\begin{array}{rl}x\ast y& =z_2\ast 2^8+z_1\ast 2^4+z_0\\ & =0010010101011011\\ & ={9563}_{10}\end{array}$

其中

$\times 2^n$ 可以用左移運算代替。

接續上一個例子，當

$x$、

$y$ 超過 8 bit 時，可以透過分治法實作 Karatsuba。

$x_1$、

$x_0$、

$y_1$、

$y_0$ 的位元數超出處理器的乘法的位數時，就把他們再切為

$x_{11}$、

$x_{10}$、

$x_{01}$、

$x_{00}$、…，再使用 Karatsuba 計算。以下以兩個 16 bit 數值相乘變成 32 bit 來演示（假設所採用的處理器只支援 8 位元乘法）：

由上圖可以看出計算

$z_2$、

$z_1$、

$z_0$ 時，透過分治法將

$x_1$、

$x_0$、

$y_1$、

$y_0$ 切成更小的數字執行乘法運算。最後再用左移與加法計算

$z_2\ast 2^{16}+z_1\ast 2^8+z_0$ 即可求得結果。

至此可透過分治法，運用 Karatsuba 演算法計算任意位數的大數。

使用 Karatsuba algorithm 來加速乘法與平方運算，乘法與平方由於運算較複雜，占用程式執行時間的比例也最高，因此sysprog21/bignum 使用 Karatsuba algorithm 來加速乘法與平方運算。對 sysprog21/bignum 使用 perf 進行分析如下 (只擷取部分)

$ sudo perf report --stdio -g graph,0.5,caller
# To display the perf.data header info, please use --header/--header-only options.
#
#
# Total Lost Samples: 0
#
# Samples: 5K of event 'cycles'
# Event count (approx.): 7673070456
#
# Children      Self  Command    Shared Object      Symbol                        >
# ........  ........  .........  .................  ..............................>
#
    75.85%     2.07%  fibonacci  fibonacci          [.] main
            |          
            |--73.78%--main
            |          |          
            |          |--26.61%--bn_sqr
            |          |          |          
            |          |          |--9.77%--_apm_mul_base
            |          |          |          |          
            |          |          |          |--2.89%--apm_dmul
            |          |          |          |          
            |          |          |           --2.13%--apm_dmul_add
            |          |          |          
            |          |          |--3.40%--__GI___libc_realloc (inlined)
            |          |          |          |          
            |          |          |           --2.85%--_int_realloc
            |          |          |                     |          
            |          |          |                     |--0.93%--_int_free
            |          |          |                     |          
            |          |          |                      --0.85%--_int_malloc
            |          |          |          
            |          |          |--2.27%--__GI___libc_malloc (inlined)
            |          |          |          |          
            |          |          |           --0.52%--tcache_get (inlined)
            |          |          |          
            |          |          |--2.04%--apm_sqr
            |          |          |          
            |          |          |--1.97%--_int_free
            |          |          |          
            |          |          |--1.37%--__memmove_avx_unaligned_erms
            |          |          |          
            |          |           --0.52%--__GI___libc_free (inlined)
            |          |          
            |          |--19.35%--bn_mul
            |          |          |          
            |          |          |--4.56%--_apm_mul_base
            |          |          |          |          
            |          |          |          |--1.31%--apm_dmul
            |          |          |          |          
            |          |          |           --0.97%--apm_dmul_add
            |          |          |          
            |          |          |--3.90%--__GI___libc_malloc (inlined)
            |          |          |          |          
            |          |          |          |--1.20%--tcache_get (inlined)
            |          |          |          |          
            |          |          |           --0.56%--checked_request2size (inlin>
            |          |          |          
            |          |          |--2.56%--_int_free
            |          |          |          
            |          |          |--1.80%--apm_mul
            |          |          |          
            |          |          |--1.25%--__GI___libc_realloc (inlined)
            |          |          |          |          
            |          |          |           --0.96%--_int_realloc
            |          |          |          
            |          |          |--1.12%--__memmove_avx_unaligned_erms
            |          |          |          
            |          |           --0.68%--__GI___libc_free (inlined)
            |          |          
            |          |--16.17%--bn_add
            |          |          |          
            |          |          |--3.82%--apm_add_n
            |          |          |          
            |          |          |--3.45%--__GI___libc_realloc (inlined)
            |          |          |          |          
            |          |          |           --2.77%--_int_realloc
            |          |          |                     |          
            |          |          |                     |--0.90%--_int_malloc
            |          |          |                     |          
            |          |          |                      --0.60%--_int_free
            |          |          |          
            |          |           --3.08%--apm_add
            |          |          
            |          |--3.87%--bn_lshift
            |          |          |          
            |          |          |--1.31%--apm_lshift
            |          |          |          
            |          |           --0.77%--__memset_avx2_unaligned_erms
            |          |          
            |          |--1.99%--bn_init
            |          |          |          
            |          |           --1.54%--__libc_calloc
            |          |                     |          
            |          |                      --0.71%--_int_malloc
            |          |          
            |          |--1.83%--_int_free
            |          |          
            |          |--1.62%--bn_swap
            |          |          
            |           --1.21%--bn_init_u32
            |                     |          
            |                      --1.08%--bn_init
            |                                |          
            |                                 --0.91%--__libc_calloc
            |          
            |--1.30%--_start
            |          __libc_start_main
            |          main
            |          
             --0.77%--0xffffffffffffffff
                       main

可見使用 Karatsuba algorithm 的 sysprog21/bignum，其乘法與平方是程式執行時間佔比最高的運算，分別為 19.35% 與 26.61%，合計近 46% 的時間占比。

Schönhage–Strassen Algorithm

Schönhage–Strassen 的概念是將

$x$、

$y$ 以

$n$ 個位數為單位，將大數拆成若干個較小的數

$x_0$、

$x_1$…、

$y_0$、

$y_1$、…，對

$(x_0,x_1,...)$ 和

$(y_0,y_1,...)$ 線性摺積，執行完線性摺積後將所有元素依據自己的位數左移並相加就可以完成乘法運算。舉例：

令

$x=26699,y=188$

將

$x$、

$y$ 以 2 個位數為界分割會得到以下兩個數列：

$(2,66,99),(1,88)$

對數列線性摺積：

         2    66    99
x              1    88
-----------------------
       176  5808  8712
   2    66    99
-----------------------
   2   242  5907  8712

$(2,242,5907,8712)$ 為

$(2,66,99)$ 與

$(1,88)$ 的線性摺積，接下來透過左移和加法將線性摺積的結果相加即可完成大數乘法：

        8712
      5907
     242
+    2
--------------
     5019412