2022q1 Homework5 (quiz5)

contributed by < linjohnss >

2022q1 第 5 週測驗題

測驗 1

實作利用 Sieve of Eratosthenes，嘗試找出第 1500 個回文質數，有效位數是 15

演算法: sieve of Eratosthenes

給定要找的質數範圍 n，找出 \(\sqrt n\) 內的所有質數，並依序剔除每著質數在範圍 n 內的所有倍數，留下來的數字即為和為內的所有質數。

程式碼實做原理

static ull pow10[] = {
    1,
    10,
    100,
    1000,
    10000,
    100000,
    1000000,
    10000000,
    100000000,
    1000000000,
    10000000000,
    100000000000,
    1000000000000,
    10000000000000,
    100000000000000,
    1000000000000000,
    10000000000000000,
};
static inline ull fastpow10(const int n)
{
    if (n < 0 || n > 16) {
        free(psieve);
        printf("n = %d\n", n);
        exit(1);
    }

    return pow10[n];
}

fastpow10 用查表的方式，找出 10 的冪
後續會用來快速配置 100...001

/* isqrt64_tab[k] = isqrt(256 * (k + 65) - 1) for 0 <= k < 192 */
static const uint8_t isqrt64_tab[192] = {
    128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142,
    143, 143, 144, 145, 146, 147, 148, 149, 150, 150, 151, 152, 153, 154, 155,
    155, 156, 157, 158, 159, 159, 160, 161, 162, 163, 163, 164, 165, 166, 167,
    167, 168, 169, 170, 170, 171, 172, 173, 173, 174, 175, 175, 176, 177, 178,
    178, 179, 180, 181, 181, 182, 183, 183, 184, 185, 185, 186, 187, 187, 188,
    189, 189, 190, 191, 191, 192, 193, 193, 194, 195, 195, 196, 197, 197, 198,
    199, 199, 200, 201, 201, 202, 203, 203, 204, 204, 205, 206, 206, 207, 207,
    208, 209, 209, 210, 211, 211, 212, 212, 213, 214, 214, 215, 215, 216, 217,
    217, 218, 218, 219, 219, 220, 221, 221, 222, 222, 223, 223, 224, 225, 225,
    226, 226, 227, 227, 228, 229, 229, 230, 230, 231, 231, 232, 232, 233, 234,
    234, 235, 235, 236, 236, 237, 237, 238, 238, 239, 239, 240, 241, 241, 242,
    242, 243, 243, 244, 244, 245, 245, 246, 246, 247, 247, 248, 248, 249, 249,
    250, 250, 251, 251, 252, 252, 253, 253, 254, 254, 255, 255,
};

/* integer square root of a 64-bit unsigned integer */
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);
}

isqrt 利用 lookup table (LUT) 的手法實做整數開平方根，用於找到 sieve of Eratosthenes 演算法的 \(\sqrt n\)

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 */
            }
        }
    }
}

根據位數 digits 產生對應的 bitmaps psieve

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. 首先，測試是否為偶數（代表可以被 2 整除）
2. 在 for 迴圈裡面
   • pquadbits 每次取 64 bits 的 psieve
   • quad 用於找出 psieve 的下一個質數，初始化 AND ~0b111，跳過已另外處理的 1, 3, 5 (為何要取 NOT?)
   • prev 每次加 64，用於判斷查找的進度是否小於 half_max
3. 接著在 while 迴圈裡
   • quad 中 1 表示尚未檢驗過的質數，所以我們可以用 __builtin_ctzll 來找出左邊第一格 1 之前的 0 的個數，進而推出 next
     • EXP1 = __builtin_ctzll(quad)
   • 判斷完是否可以整除後，要在 quad 標注已檢查的質數為 0，利用 i 代表位數與 quad 作 AND 運算
     • EXP2 = ~(1ULL << i)

TODO:

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

   是否有必要先將數值轉成字串？用十進位的角度處理運算是否產生額外的計算負擔？

2. Linux 核心原始程式碼 lib/math/prime_numbers.c 有類似的程式碼，請解釋其作用、裡頭 RCU 的使用情境，及針對執行時間的改進

測驗 2

程式碼行為