2022q1 Homework5 (quiz5)

contributed by < Kevin-Shih >

第 5 週測驗題

測驗 1

題目
利用 Sieve of Eratosthenes，嘗試找出第 1500 個回文質數，有效位數是 15。 預期執行結果為 100191191191001。

完整程式碼: prime_palindrome.c

不用複製原題目，專注於探討程式原理，包含你的推測及驗證，只要列出關鍵程式碼。善用 GitHub 或 gist 保存完整程式碼及相關測試。
jserv

收到，另外放在 gist 中。

Kevin Shih

運作原理與解題

部分較為重要的片段及本題答案

以查表方式求 64-bit uint 的整數平方根

從 第 3 週測驗題 11 改寫裡頭的整數開平方根實作，藉由 lookup table (LUT) 的手法降低分支的數量 –quiz5 題目敘述

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

TODO:
解釋上述程式碼的原理

jserv

建構 bitmaps 用於進行 Sieve of Eratosthenes algorithm

/* 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); // why not psieve[0] |= COMPOSITE

    unsigned char mask = 0x7;
    for (ull n = 3; n <= max; n += 2) {
        /* 若 n 是質數 */
        if (((psieve[n >> 4] >> ((n >> 1) & mask)) & 0x1) == PRIME) { 
            /* 所有篩選表中 n 的倍數所代表的 bit 均設為 COMPOSITE (從 2n 到 max 之間) */
            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] |= psieve[0] |= COMPOSITE
                                    << ((mul >> 1) & mask); /* bit offset */
            }
        }
    }
}

Sieve of Eratosthenes 是透過依序刪去表中質數的倍數，最後剩下的就是質數的一種搜尋質數的方式，由於 2 以外的偶數均不是質數，在建立篩選表時只保留奇數可以節省一半的空間。 這個 bitmap 的 max 是最大的回文的 isqrt，因為後續檢驗質數不需查大於 sqrt 值 (若有大於 sqrt 的因數必然有對應的小於 sqrt 的因數)

TODO:
補上質數測試時，為何只要檢驗到

$\sqrt{n}$ 的數學證明。
jserv

關於 bit offset

psieve 陣列中的每個元素都可表示 16 個數 (只存奇數)，因此 n >> 4 取得對應 n 所在的元素；第 (n - 1) / 2 個 bit 代表 n (即第 n >> 1 bit)，前面已經處理過 n 所在的元素，因此只取 (n >> 1) % 8 即可，這邊透過 AND 一個 0x7 的 mask 達成，最後 AND 0x1 取得最低為的 bit 即為代表 n 的 bit，而若是指派數值給 n 則找到 n 所在的元素，並將要給定的值右移 (n >> 1) & mask 後，指派到該元素。

透過篩選表檢驗給定的回文是否為質數

/* 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 = __builtin_ctzll(quad);
            next = prev + i;
            if (!(val % ((next << 1) + 1)))
                return false;
            quad &= ~(1ULL << i);
        }

        prev += 64;
    }

    return true;
}

不用列出 EXP1 和 EXP2

jserv

留下答案已移到 gist 了。 (其實寫在這單純為了方便查閱、複習)
Kevin Shih

這個函式用來檢驗找出的回文是否為質數，對 2, 3, 5 是否可整除已經處理過了，後續迴圈內只要處理是否被其他質數整除即可。

• prev: 先前檢驗到這個 bit
• next: 下一個要檢驗的質數對應 psieve 中的第幾個 bit (即檢驗 val 是否被 next * 2 + 1 整除)
• quad: 先前建立的篩選表，每次取 64-bit 取 NOT，用於查找下一個要檢驗的質數，起始值再與 ~0b111 做 AND 以忽略以另外處理的前三個奇數(1, 3, 5)

quad = ~*pquadbits 得到的(部分)篩選表中 0 表示非質數或檢驗過的質數 1 表示尚未檢驗過的質數，如此每次只要取 ctz 就可以得到要檢驗的質數對應 psieve 中的第幾個 bit，檢驗後不論是否能整除都將該 bit 設為 1 以表示檢驗過了即可。 當 quad 為 0 時表示這 64-bit 中所有質數都已被檢驗過了，故將 prev 加 64 並進入下一次迴圈。 因 quad 為 ull 型態故 EXP1 = __builtin_ctzll(quad)，而用於標註檢驗過的質數的 EXP2 = ~(1ULL << i)，即尾隨 i 個 0 就左移 i bit。

產生下一個回文

/* Generate the next palindrome after a palindrome
 *   123454321 -> 123464321
 *   123494321 -> 123505321
 *   123999321 -> 124000421
 */
static ull getnextpalin(char *buf, const int *len)
{
    midindex = (*len >> 1) - 1;

    /* Handle the case of odd number of digits.
     * If the central digit is 9, reset it to 0 and process adjacent digit,
     * otherwise increment it and return.
     */
    if ((*len & 0x1) == 0x1) {
        if (buf[midindex + 1] - '0' != 9) {
            buf[midindex + 1] += 1; 
            return atol(buf);
        }
        buf[midindex + 1] = '0';
    }

    /* Adjust other digits */
    for (; midindex >= 0; --midindex) {
        if (buf[midindex] - '0' != 9) {
            buf[midindex] += 1;
            buf[*len - midindex - 1] += 1;
            return atol(buf);
        }
        buf[midindex] = '0';
        buf[*len - midindex - 1] = '0';
    }

    /* If we have exhausted numbers in *len digits, increase the number of
     * digits and return the first palindrome of the form 10..0..01
     */
    *len += 2;
    /* The last palin will always be of the form 99...99 */
    return atol(buf) + 2;
}

當位數不夠時應修改 *len 以觸發 main 中 refresh buffer 的部分。

產生回文質數

int main()
{
    int count = LIMIT;
    ull i = 100000000000001;
    int len = 0;
    char *buf = ltoa(i, &len);

    if (len < N_DIGITS) {
        printf("len: %d\n", len);
        exit(1);
    }

    generate_sieve(N_DIGITS);

    /* We started at 1000000000001. */
    while (1) {
        /* If number of digits are even, all palindromes are divisible by 11.
         * Get first palindrome of next odd number of digits.
         */
        if (!(len & 0x1)) {
            i = 1ULL | fastpow10(len);
            buf = ltoa(i, &len);
            continue;
        }

        /* Check if number is divisible by 5 or 3 */
        if ((buf[len - 1] != '5') && !is_divisible_by3(buf, len)) {
            if (isprime(i) && (--count == 0)) {
                printf("%s\n", buf);
                return 0;
            }
        }

        int oldlen = len;
        i = getnextpalin(buf, &len);
        /* Refresh buffer if number of digits has increased */
        if (oldlen != len)
            buf = ltoa(i, &len);
    }

    free(psieve);

    return 0;
}

發現並不是找出第 1500 個回文質數，而是第 1500 個 15 位數的回文質數。 因為當 LIMIT 設為 1 時其輸出為 100000323000001 而非 2。 因為 len 一開始就是 15 位數，而非從個位數慢慢往上找(第一個回文 i 設為 100000000000001)。