---
tags: linux2022
---

# 2022q1 Homework5 (quiz5)
contributed by < `scottxxxabc` >
> [測驗題目](https://hackmd.io/@sysprog/linux2022-quiz5)
:::success
1. 解釋上述程式碼運作原理，指出可改進之處並實作
> 是否有必要先將數值轉成字串？用十進位的角度處理運算是否產生額外的計算負擔？
:::
### isqrt
:::spoiler 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](https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Babylonian_method) 來逼近正確答案
求 $a$ 的平方根時，若是計算出的平方根 $b$ 小於實際平方根，則 $\dfrac{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

:::spoiler 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 的話就代表 $2n+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 代表的數字為 $2n+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*next+1$ 整除，若是可以則直接 return false，代表 `val` 不是質數。
* `quad &= ~(1ULL << i)` 將已經檢查過的 bit clear 為 0。