# 2022q1 Homework5 (quiz5)
contributed by < `Kevin-Shih` >
> [第 5 週測驗題](https://hackmd.io/@sysprog/linux2022-quiz5#測驗-1)
## 測驗 `1`
**題目**
利用 [Sieve of Eratosthenes](https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes)，嘗試找出第 1500 個回文質數，有效位數是 15。 預期執行結果為 `100191191191001`。
完整程式碼: [prime_palindrome.c](https://gist.github.com/Kevin-Shih/5dc5ac521875f56958721ba14f80f3a8)
:::warning
:warning: 不用複製原題目，專注於探討程式原理，包含你的推測及驗證，只要列出關鍵程式碼。善用 GitHub 或 gist 保存完整程式碼及相關測試。
:notes: jserv
:::
:::info
收到，另外放在 gist 中。
:ballot_box_with_check: Kevin Shih
:::
#### 運作原理與解題
部分較為重要的片段及本題答案
**以查表方式求 64-bit uint 的整數平方根**
> 從 [第 3 週測驗題 `11`](https://hackmd.io/@sysprog/linux2022-quiz3#%E6%B8%AC%E9%A9%97-11) 改寫裡頭的整數開平方根實作，藉由 lookup table (LUT) 的手法降低分支的數量 --quiz5 題目敘述
```c
/* 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);
}
```
:::warning
TODO:
解釋上述程式碼的原理
:notes: jserv
:::
**建構 bitmaps 用於進行 Sieve of Eratosthenes algorithm**
```c
/* 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 的因數)
:::warning
TODO:
補上質數測試時，為何只要檢驗到 $\sqrt n$ 的數學證明。
:notes: jserv
:::
關於 bit offset
```graphviz
digraph "reverse-bits" {
node [shape=record];
rankdir=LR;
a [label="{<0> 1 |<1> 3 |<2> 5 |<3> 7 |<4> 9 |<5> 11 |<6> 13 |<7> 15 }"];
b [label="{<0> 1 |<1> 1 |<2> 1 |<3> 0 |<4> 0 |<5> 0 |<6> 0 |<7> 0 }"]
node [shape="plain"]
"uint8_t psieve[0]"->a:0:w
"mask" -> b:0:w
}
```
`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` 後，指派到該元素。
**透過篩選表檢驗給定的回文是否為質數**
```c
/* 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;
}
```
:::warning
不用列出 `EXP1` 和 `EXP2`
:notes: jserv
:::
:::info
留下答案已移到 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。
**產生下一個回文**
```c
/* 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 的部分。
**產生回文質數**
```c
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`)。
Sign in with Wallet
Connect another wallet