# 2022q1 Homework5 (quiz5)
contributed by < [Nomad1230](https://github.com/Nomad1230) >
## 測驗 `1`
### 解釋程式碼原理
#### 建立質數表
:::spoiler generate_sieve
```c=
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 */
}
}
}
}
```
:::
此段程式碼基於 Sieve of Eratosthenes 演算法來建立質數表，參數 `n` 代表欲搜尋之有效位數區間。
第 3 ~ 16 行決定質數表的大小，首先根據有效位數取出此區間最大的數值作為 `max` ，例如本題要計算有效位數為 15 的區間，則 `max` 為 10 ^ 15 - 1 ，再來將 `max` 取 isqrt ，因為在進行試除法時只需判斷到 $\sqrt{N}$ 。
:::info
在進行試除法時會將目標數 $N$ 對 1 ~ $\sqrt{N}$ 一一進行試除，若整除則代表 ${N}$ 有一個小於 $\sqrt{N}$ 及一個大於 $\sqrt{N}$ 的因數，因此不用再對大於 $\sqrt{N}$ 的數進行試除。
:::
此時 `max` 為預期的質數表大小，但由於所有的偶數中只有 2 為質數，故可以把所有偶數忽略不計，因此再將 `max` 除以 2 。
最後將 `max` 除以 8 就是要分配給記憶體的 bytes 數。
在分配完空間之後便是要將質數填入表中以方便查詢，第 27 ~ 43 行的程式從 3 ~ `max` 進行試除，若試除成功則將其所有的倍數設成 COMPOSITE。
值得注意的是質數表的儲存方式，一個 `psieve` 指標所指向的值可以儲存 16 個數的質數表，因為偶數忽略不計，剩下 8 個奇數的 state 以 bitwise 的方式儲存在 `uint8_t` 的變數中，如圖所示:
```graphviz
graph graph1{
label = "psieve[0]"
node [shape = record]
num1 [label = "1|3|5|7|9|11|13|15"]
}
```
```graphviz
graph graph2{
label = "psieve[1]"
node [shape = record]
num1 [label = "17|19|21|23|25|27|29|31"]
}
```
故 `psieve` 的索引值都需先 >> 4 ，儲存質數表，且在儲存質數表時都要以 `|=` 運算子來儲存特定的 bit 。
#### 用查表輔助 64 bit 開方根
:::spoiler isqrt
```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);
}
```
:::
首先可以看到定義了 `isqrt64_tab` 來儲存 `isqrt(256 * (k + 65) - 1)` for 0 <= k < 192 的結果，看到第 25 、 26 行中在查表之前已經先把 `x` left shift 掉 clz & 62 的個數，若 `x` 的 clz 為奇數，則 shift 完最左邊兩 bits 為 01 ，為偶數則 shift 完最左邊的 bit 為 1 ，為甚麼要做這樣的分別呢?
:::info
在查表時是取 `x` 最左邊的 8 bits 作為 index ，此 index 必為 64 以上，故 table 只需儲存 256 - 64 = 192 個元素。
:::
原因就出在回傳值的部分:
```c
return y >> (lz >> 1);
```
一開始取出 clz 是為了將數值都對齊左側方便接下來直式開方法的運算，所以最後計算出的結果 `y` 要右移 (clz / 2) 個位數(因為開方所以除2)，假設 clz 的值為 27 ，由於開方前先把數值乘上了 $2^{27}$ ，所以最後結果應該要右移掉 $2^{13.5}$ ，這樣會造成計算上的問題，故多用了一個變數 `lz` 且使其必為偶數，來確保運算中的正確性與便利性。
第 26 ~ 29 行做的是直式開方法，值得注意的是這邊一次取 8 bits 為單位在做開方法。
先看到一般十進位的直式開方法流程:
[圖源](https://highscope.ch.ntu.edu.tw/wordpress/?p=51628)
以上圖為例，若要做 144 的開方，先把面積減掉邊長為 10 的正方形，剩下 44 的面積為兩個矩形及一個小正方形，以數學公式來看就是 $2ab + b^2 = b(2a + b) = 44$ 。
這時再選定小正方形的長度使得小正方形 + 兩個矩形的面積逼近 44 ，若面積還有剩餘則再選定一個更小的正方形一直迭代下去，最後所有正方形加起來的長度就會逼近於開方根。
從以上流程可以寫出直式開方法的演算法流程:
(1).先求最高位數的平方根 a
(2).選定一數 b 使得 b(2a + b) 接近但不超過餘數
(3).若餘數不為 0 則重複第二步
再來看到程式碼，可以看出其是將 64 位元的數切成 8 等份，每一份都是一個 byte 的大小，因此可以看成是 256 進位的 8 位數的直式開方法。
第 26 行為步驟 1 ， 第 27 、 28 行為步驟 2 的迭代。
以第 27 行為例:
```c
y = (y << 7) + (x >> 41) / y;
```
`y` 在第 26 行的結果是最高位的 byte 的開方根，
#### 質數判斷
:::spoiler isprime
```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); //EXP1;
next = prev + i;
if (!(val % ((next << 1) + 1)))
return false;
quad &= ~(1ULL << i); //EXP2;
}
prev += 64;
}
return true;
}
```
:::
Sign in with Wallet
Connect another wallet