# 2022q1 Homework5 (quiz5)

contributed by < [tommy2234](https://github.com/tommy2234) >

題目連結 [quiz5](https://hackmd.io/@sysprog/linux2022-quiz5)

## 測驗 1
本題利用 [Sieve of Eratosthenes](https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes)，嘗試找出第 1500 個回文質數，有效位數是 15

- 改寫[第 3 週測驗題](https://hackmd.io/@sysprog/linux2022-quiz3)裡頭的整數開平方根實作，藉由 [lookup table](https://en.wikipedia.org/wiki/Lookup_table) (LUT) 的手法降低分支的數量: (但這不必然較快)

static ull isqrt(ull x)
{
    if (x == 0) return 0;

    int lz = __builtin_clzll(x) & 62;
    x <<= lz;
    uint32_t y = isqrt64_tab[(x >> 56) - 64];
    y = (y << 7) + (x >> 41) / y;
    y = (y << 15) + (x >> 17) / y;
    y -= x < (uint64_t) y * y;
    return y >> (lz >> 1);
}

### 產生一個能用來檢驗質數的 bitmap
此 bitmap 將會用在 **Sieve of Eratosthenes algorithm** 之中。
其原理是是從2開始，將每個質數的各個倍數標記成合數。

psieve 只需要檢驗到最大 palindrome (max) 的開根號 ，並且所有 2 的倍數都被忽略。

/* Generates a table of prime numbers up to the maximum having digits.
 * All even numbers (including 2) are ignored. Bit n in the bitmap
 * represents the odd number (2n + 1). This enables the table to be
 * half the size of a naive Sieve of Eratosthenes implementation.
 */
static void generate_sieve(int digits)
{
    ull max = 0;
    for (int count = 0; count < digits; ++count)
        max = max * 10 + 9;

    max = isqrt(max);
    half_max = max >> 1;

    /* We need half the space as multiples of 2 can be omitted */
    int bytes = (max >> 1) + (max & 0x1);

    /* Calculate the actual number of bytes required */
    bytes = (bytes >> 3) + (bytes & 0x1);

    bytes = ALIGN(bytes, 16); /* Align-up to 16-byte */
    psieve = realloc(psieve, bytes);
    if (!psieve) {
        printf("realloc() failed!\n");
        exit(1);
    }
    memset(psieve, 0, bytes);

    /* In psieve bit 0 -> 1, 1 -> 3, 2 -> 5, 3 -> 7 and so on... */
    /* Set the 0th bit representing 1 to COMPOSITE
     */
    psieve[0] |= COMPOSITE << (1 >> 1);

    unsigned char mask = 0x7;
    for (ull n = 3; n <= max; n += 2) {
        if (((psieve[n >> 4] >> ((n >> 1) & mask)) & 0x1) == PRIME) {
            for (ull mul = (n << 1); mul < max; mul += n) {
                /* Skip the evens: there is no representation in psieve */
                if (!(mul & 0x1))
                    continue;

                /* Set offset of mul in psieve */
                psieve[mul >> 4] |= COMPOSITE
                                    << ((mul >> 1) & mask); /* bit offset */
            }
        }
    }
}

psieve 中的 bit n 代表的是奇數 (2n - 1)。
一開始將 bit 0 設為 1 代表和數

psieve[0] |= COMPOSITE << (1 >> 1);

然後逐個數字檢驗質數

if (((psieve[n >> 4] >> ((n >> 1) & mask)) & 0x1) == PRIME) {      

- 由於一個 psieve 的元素有 8 bits，可以覆蓋 16 個數字 ( 2 的倍數都忽略 )，因此 `n >> 4` 得到 n 所對應的 bit 位在 psieve 中的 index，`psieve[n >> 4]` 就得到對應的元素。

- 前面已經得到 n 所在的 index ， `((n >> 1) & mask))` 則是計算 n 在此元素中對應到的 bit (offset)。
psieve 已經將 2 的倍數都忽略不紀錄，因此 n >> 1（除以 2）可以得到 n 對應 bitmap 中第幾個 bit。然後我們可以透過 and 一個 0x7 的 mask 來達到 mod 8 的效果，以取得 offset 。

- 組合起來變成 `(psieve[n >> 4] >> ((n >> 1) & mask)) & 0x1)` 。將 psieve[index] 右移 offset 再 and 0x1 得到 LSB ，終於取得了 n 對應的 bit 。

- 最後我們檢查 n 是否為質數，如果是就將其所有倍數都標記為合數，這就形成了一個 Eratosthenes 篩表。

/* Set offset of mul in psieve */
psieve[mul >> 4] |= COMPOSITE 
            << ((mul >> 1) & mask); /* bit offset */

### 運用篩表檢驗質數
透過篩檢驗傳入的 val 是否被篩表中某個質數整除，如果表中沒有任何質數能整除 val ，那 val 就是質數。

/* Check if a number is prime. */
static bool isprime(const ull val)
{
    if (!(val & 0x1)) /* Test for divisibility by 2 */
        return false;

    ull *pquadbits = (ull *) psieve;
    ull next = 3; /* start at 7 (i.e. 3 * 2 + 1) */

    for (ull quad = ~*pquadbits & ~0b111, prev = 0; prev <= half_max;
         quad = ~*++pquadbits) {
        if (!quad) {
            prev += 64;
            continue;
        }

        while (quad) {
            ull i = EXP1;
            next = prev + i;
            if (!(val % ((next << 1) + 1)))
                return false;
            quad &= EXP2;
        }

        prev += 64;
    }

    return true;
}

函式中有幾個重要變數
1. `pquadbits`

    ull *pquadbits = (ull *) psieve;

ull 型態的 pquadbits 一次從 psieve 中拿取 64 個 bit 做檢驗，也就是一部份的篩表。

2. `next`
在下一個要檢驗的質數在 psieve 中對應的 bit ，所以 bit next 對應的數字為 next * 2 + 1 。
呼叫 isprime 前對 2、3、5 已做過檢驗，因此從 7 開始檢驗，next 起始為 3 (3 * 2 + 1 = 7)。

3. `prev`
前一個檢查的 bit。

4. `quad`
對 \*pquadbits 取 complement ， 目的是為了套用 `ctz` 快速定位下一個要檢驗的質數，檢驗過後透過 `quad &= ~(1ULL << i);` 將前方的 bits 都清成 0 ， 下次再用 `ctz` 又可以找到下個待檢驗質數。

:::warning
TODO:
1. 指出可改進之處並實作

> 是否有必要先將數值轉成字串？用十進位的角度處理運算是否產生額外的計算負擔？
2. Linux 核心原始程式碼 [lib/math/prime_numbers.c](https://github.com/torvalds/linux/blob/master/lib/math/prime_numbers.c) 有類似的程式碼，請解釋其作用、裡頭 [RCU](https://hackmd.io/@sysprog/linux-rcu) 的使用情境，及針對執行時間的改進
:::