# [2019q3](http://wiki.csie.ncku.edu.tw/sysprog/schedule) 第 8 週測驗題
:::info
目的: 檢驗學員對 CS:APP 第 5 章的認知
:::

---

### 測驗 `1`

[bc](https://linux.die.net/man/1/bc) 工具很好用，可讓我們快速計算特定數值表示，如:

$ echo "2^64" | bc -l

預期輸出為:

18446744073709551616

這個十進位的輸出需要幾個字元來展現呢？

$ echo "2^64" | bc -l | wc -c

其中 [wc](https://linux.die.net/man/1/wc) 可計算字元數量，於是我們知道 2^64^ 用十進位表示，需要 21 個位數。倘若我們要計算 48 個位數的十進位運算，顯然不能用 long 或 long long (在 [LP64 data model](https://en.wikipedia.org/wiki/64-bit_computing#64-bit_data_models) 中，這兩者都是 64 位元的寬度) 型態來處理，這時候我們需要 [大數運算](https://openhome.cc/Gossip/AlgorithmGossip/BigNumber.htm) 的實作。

考慮到 64 位元的有效數值上限是 `0xFFFFFFFFFFFFFFF`，而十進位 $99999999×99999999$ 的結果在這有效範圍以內，被乘數 (前者) 和乘數 (後者) 都是 8 個位數 (有利於 x86-64 架構)。於是下方程式碼嘗試依照每 8 個十進位數字，進行分割運算: (`big-mul`):

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

static inline void add_with_carry(char *input, char num, int idx) {
    char input_num = input[idx] + num;
    if (input_num <= 9) {
        input[idx] = input_num;
        return;
    }
    char carry = input_num / 10;
    input[idx] = input_num % 10;
    add_with_carry(input, carry, idx - 1);
}

#define max(a, b)               \
    ({                          \
        __typeof__(a) _a = (a); \
        __typeof__(b) _b = (b); \
        _a > _b ? _a : _b;      \
    })

static int add(const char *a1, const char *a2, char *ret, int *ppos) {
    int len1 = strlen(a1), len2 = strlen(a2);
    int len = max(len1, len2);

    int i;
    for (i = len - 1; i >= 0; i--) {
        unsigned char tmp;
        int len_diff = len1 - len2;
        if (len_diff >= 0) {
            tmp = ARR1[i] - '0';
            if (i > len_diff - 1)
                tmp += ARR2[i - len_diff] - '0';
        } else {
            tmp = ARR3[i] - '0';
            if (i > -len_diff - 1)
                tmp += ARR4[i + len_diff] - '0';
        }
        add_with_carry(ret, tmp, i + 1);
    }

    *ppos = 1;
    if (ret[BBB]) {
        len++;
        *ppos = 0;
    }

    for (i = 0; i < len + 1; i++)
        ret[i] += '0';

    return len;
}

char *mul(const char *m1, const char *m2) {
    int len = strlen(m1);
    int ret_len = len + strlen(m2);
    int n = 0, pos = 0, nblocks = len / 8;
    char str1[10], str2[10], resstr[20];

    char *tmp_ret = calloc(ret_len, 1);
    char *pow1 = calloc(len * 2, 1), *pow2 = calloc(len * 2, 1);
    char *ret = calloc(ret_len, 1);

    /* split into several blocks, 8 digits per block */
    for (int i = 0; i < nblocks; i++) {
        memcpy(str1, m1 + len - i * 8 - 8, 8);
        unsigned long im1 = atoi(str1);

        for (int j = 0; j < nblocks; j++) {
            memcpy(str2, m2 + len - j * 8 - 8, 8);
            unsigned long im2 = atoi(str2);
            unsigned long tmp_res = im1 * im2;

            /* multiplying pow(10,n) equals to preserving trailing '0' */
            int pow1size = i * 8, pow2size = j * 8;

            int reslen = sprintf(resstr, "%lu", tmp_res);
            int totlen = reslen + pow2size;

            memset(pow2, '0', totlen);
            memcpy(pow2, resstr, reslen);

            reslen = totlen;
            totlen += pow1size;
            memset(pow1, '0', totlen);
            pow1[totlen] = '\0';
            memcpy(pow1, pow2, reslen);

            memset(ret, 0, n + pos);

            n = add(pow1, tmp_ret, ret, &pos);
            memcpy(tmp_ret, ret + pos, n);
        }
    }
    memset(ret, 0, n + pos);
    memcpy(ret, tmp_ret, n);

    free(tmp_ret);
    free(pow1);
    free(pow2);

    return ret;
}

int main(int argc, char **argv) {
    if (argc < 3)
        return -1;
    char *m1 = argv[1], *m2 = argv[2];

    char *ret = mul(m1, m2);
    printf("%s\n", ret);
    free(ret);

    return 0;
}

參考執行方式和輸出:

./mul 119334567890334449388883313579158334567098134455 \
      667908995633221198765432134678040000123411113456

預期可得:

79704631383957730438879843848804741889926116047138197998269353980447530720116354515912427726480

其中被乘數和乘數都用 48 個字元 (再加 1 個 null terminator) 的字串來表示 (可執行 `$ echo -n "119334567890334449388883313579158334567098134455" | wc -c`)，最終乘法結果用 95 個字元 (再加 1 個 null terminator) 的字串來表示。

請補完程式碼。

==作答區==

ARR1 = ?
* `(a)` a1
* `(b)` a2
* `(c)` ret

ARR2 = ?
* `(a)` a1
* `(b)` a2
* `(c)` ret

ARR3 = ?
* `(a)` a1
* `(b)` a2
* `(c)` ret

ARR4 = ?
* `(a)` a1
* `(b)` a2
* `(c)` ret

BBB = ?
* `(a)` len
* `(b)` `len - 1`
* `(c)` i
* `(d)` 1
* `(e)` 0

:::success
延伸問題:
1. 解釋上述程式碼運作原理，並指出影響運算正確的缺陷;
2. 修正上述程式碼的缺陷，嘗試依據 CS:APP 第 5 章的手法予以分析執行時期的成本，從而改善效能;
3. 研究 [Karatsuba algorithm](https://en.wikipedia.org/wiki/Karatsuba_algorithm) 快速相乘演算法，嘗試和本實作比較，被乘數和乘數至少該是十進位 48 個位數
- 解說影片: [Fast Multiplication: From Grade-School Multiplication To Karatsuba's Algorithm](https://youtu.be/-dfsxsiGoC8)
- 善用 [GMP: GNU Multiple Precision Arithmetic Library](https://gmplib.org/)
4. 研讀 [隱藏的傅立葉分析](https://hackmd.io/@sysprog/S1jR0mHfN) 提到透過傅立葉轉換來進行乘法運算，嘗試引入到大數運算並比較正確性及效能表現;
- 參照 [FFT based multiplication of large numbers](http://numbers.computation.free.fr/Constants/Algorithms/fft.html) 及 [Multiplying Huge Integers using Fourier Transforms](http://www.cs.rug.nl/~ando/pdfs/Ando_Emerencia_multiplying_huge_integers_using_fourier_transforms_paper.pdf)
> A plot showing the average time spent per multiplication operation as a function of the sum of the lengths of the
integers to be multiplied. Both axes are dawn on a logarithmic
scale (base 10)
:::