2022q1 Homework5 (quiz5)

contributed by < scottxxxabc >

測驗題目

解釋上述程式碼運作原理，指出可改進之處並實作

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

isqrt

Code

/* 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 使用到查表的方式來計算平方根，clz
x << lz 的範圍就會被限定在 0x4000 0000 ~ 0xFFFF FFFF，isqrt 取 x 的最高 8 個 bit，讓表格的大小可以限定在 0x40 (

64_{10}

) ~ 0xFF (

255_{10}

) 之間。
接下來用最高 8 個 bit (x >> 56) - 64 作為 index k 查閱表格，可以查到 256 * (k + 65) - 1 的 isqrt 值。

k + 65 的值在 65~256 之間，表格的範圍落在

256 * 65 - 1

256 * 256 - 1

因為首先取了 8 個 bit，查表時再乘上 256 (

2^{8}

)，因此還需要將查表得到的 y 值左移

(64 - 8 - 8) / 2

，也就是24 bit。

利用兩次 Babylonian method 來逼近正確答案
求
$a$ 的平方根時，若是計算出的平方根
$b$ 小於實際平方根，則
$\frac{a}{b}$ 會大於實際平方根，將兩者平均可以更逼近實際數值，但根據算幾不等式

$\frac{\frac{a}{b} + b}{2} \geq \sqrt{\frac{a}{b} \times b}$
此方法的估計值始終會大於等於實際數值。

//y = (y << 7) + (x >> 41) / y;
y = ((y + x >> 48 / y) >> 1) << 8
    
//y = (y << 15) + (x >> 17) / y;
y = ((y + x >> 32 / y) >> 1) << 16

最後因為一開始將 x 左移 lz 個 bit，所以還要將 y 右移 lz / 2個 bit 。

Sieve of Eratosthenes

Code

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

generate_sieve 產生了一個 bitmask，每一個 bit 代表一個數字是否為質數。

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

max 代表回文數字的上限，因此只需要標記到 isqrt(max) 的數字為止，另外所有的偶數都不是質數，所以 bit 數減半。 bit

n

為 1 的話就代表

2 n + 1

這個數不是質數，bit

n

為 0 則是質數。

一開始先將 1 (bit 0)設為合數，接下來逐個檢查每個奇數是不是質數 :

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

因為不包含偶數，所以每一個 byte 表示的數相差 16，第2行利用 n >> 4 來取得 index，psieve[n >> 4]就是所在的 byte。
因為只有奇數，第(n >> 1) bit 代表的數字為
$2 n + 1$ 。利用 bitwise-AND mask，來取得所在 byte 的 offset，將 mask 設為 0x07，相當於除以 8 的餘數，也就是一個byte的長度。
迴圈將 n 的所有倍數的 bit 設為 1

isprime

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;
}

pquadbits 每次取出 psieve 的 64 個 bits
quad 將 pquadbits 反轉，因此 quad 的 0 代表合數，1 代表質數。& ~0b111 忽略前 3 個 bits ，從 7 開始檢查。
prev 代表前一次迴圈檢查到第 prev 個 bit，每次 for 迴圈將 prev 增加 quad 的大小 64 bits。若是 prev 大於 half_max 就結束迴圈。
i = __builtin_ctzll(quad) 找出 quad 結尾的 0 的數量。 quad 的 0 代表合數，所有的合數都可以拆解成二個或多個更小質數的乘積，所以可以跳過。
!(val % ((next << 1) + 1)) 檢查 val 是否能被
$2 * n e x t + 1$ 整除，若是可以則直接 return false，代表 val 不是質數。
quad &= ~(1ULL << i) 將已經檢查過的 bit clear 為 0。