# 2025q1 Homework4 (quiz3+4)
## quiz3
### 測驗一
#### 1. 解釋上述程式碼運作原理
**核心資料結構與型別設計**
* 每個 `mpi_t` 都是一個可以動態擴充的「大整數」。
* `31 bits per word`：你用 `32-bit int` 但只存 `31 bits`，有效避開加法溢位問題（詳見 `INTMAX = 0x7fffffff`）。
大數資料實際儲存在 `data` ，當你需要更大（乘法、進位、I/O），就 realloc 把 data 擴大。
每一格 `data[n]` 都只存 31 bits 有效資料，如 `data[0]` 是最低 31 bit，`data[1]` 是再高 31 bit，以此類推。
```c
typedef struct {
uint32_t *data;
size_t capacity;
} mpi_t[1];
```
這裡的 `mpi_t` 不是單純「一個陣列」，而是「一個只有一格的 `struct` 陣列」。這種寫法是C 語言一種經典技巧，多用於實現物件導向感覺的 `pass-by-reference`，讓你用 `mpi_t x`; 就可以像指標一樣用 `x` 傳參數。
以 capacity = 3 為例，資料內容：
```graphviz
digraph mpi_struct {
rankdir=TB;
node [shape=record, fontname="monospace", fontsize=13];
// stack struct node
struct [label="{ struct (on stack) | capacity = 3 | data (pointer) }", width=2];
// heap array node
heap [label="{ <d0> 0x6F23ABCD | <d1> 0x00327894 | <d2> 0x00000001 } | heap (malloc array)", shape=record, width=5];
struct:data -> heap:d0 [label=" data pointer", fontsize=16];
}
```
**初始化和記憶體配置**
```c
void mpi_init(mpi_t rop)
{
rop->capacity = 0;
rop->data = NULL;
}
```
設定 rop->data = NULL, capacity = 0
預備之後分配空間。
```c
void mpi_enlarge(mpi_t rop, size_t capacity)
{
if (capacity > rop->capacity) {
size_t min = rop->capacity;
rop->capacity = capacity;
rop->data = realloc(rop->data, capacity * 4);
if (!rop->data && capacity != 0) {
fprintf(stderr, "Out of memory (%zu words requested)\n", capacity);
abort();
}
for (size_t n = min; n < capacity; ++n)
rop->data[n] = 0;
}
}
```
若新容量比原本大，就 realloc 一塊更大的空間，而新增的那幾格清成 0。
```c
void mpi_clear(mpi_t rop)
{
free(rop->data);
}
```
釋放 rop->data，避免 memory leak。
```c
void mpi_compact(mpi_t rop)
{
size_t capacity = rop->capacity;
if (rop->capacity == 0)
return;
for (size_t i = rop->capacity - 1; i != (size_t) -1; --i) {
if (rop->data[i])
break;
capacity--;
}
assert(capacity != (size_t) -1);
rop->data = realloc(rop->data, capacity * 4);
rop->capacity = capacity;
if (!rop->data && capacity != 0) {
fprintf(stderr, "Out of memory (%zu words requested)\n", capacity);
abort();
}
}
```
去除高位都是 0 的部分，節省記憶體（壓縮成「正確位數」）。
**基本設定與數字轉換**
```c
void mpi_set_u64(mpi_t rop, uint64_t op)
{
size_t capacity = ceil_div(64, 31);
mpi_enlarge(rop, capacity);
for (size_t n = 0; n < capacity; ++n) {
rop->data[n] = op & INTMAX;
op >>= 31;
}
for (size_t n = capacity; n < rop->capacity; ++n)
rop->data[n] = 0;
}
```
把 64-bit 標準整數，拆成兩個（或三個）31-bit word 寫進 `rop->data[]`。
每次只取 31 bits，`op >>= 31` 進到下一輪。
**歐基里德演算法**
```c
void mpi_gcd(mpi_t rop, const mpi_t op1, const mpi_t op2)
{
if (mpi_cmp_u32(op2, 0) == 0) {
mpi_set(rop, op1);
return;
}
mpi_t q, r;
mpi_init(q);
mpi_init(r);
mpi_fdiv_qr(q, r, op1, op2);
mpi_gcd(rop, op2, r);
mpi_clear(q);
mpi_clear(r);
}
```
**定理敘述**
對任意兩個非負整數 $a, b$（其中 $b \neq 0$），
$$
\gcd(a, b) = \gcd(b,\, a \bmod b)
$$其中 $\gcd$ 表示最大公因數，$\bmod$ 為取餘數運算。
**證明：**
設 $d$ 是 $a$ 與 $b$ 的公因數，則 $d$ 也會是 $b$ 與 $r$ 的公因數，反之亦然。
其中 $r = a \bmod b$，依整數除法定理，存在唯一整數 $q, r$ 使得
$$
a = q \cdot b + r,\quad 0 \leq r < b
$$
**若 $d \mid a$ 且 $d \mid b$，**
因為 $a = q \cdot b + r$， 所以 $r = a - q \cdot b$。
由於 $d \mid a$ 且 $d \mid b$， 則 $d \mid (a - q \cdot b) = r$，因此，$d$ 也整除 $r$。
**若 $d \mid b$ 且 $d \mid r，$**
同理，$a = q \cdot b + r$。
因為 $d \mid b$，則 $d \mid (q \cdot b)$。
又 $d \mid r$，所以 $d \mid (q \cdot b + r) = a$。
因此，$d$ 也整除 $a$。
**歸納結論**
$a, b$ 的所有公因數等同於 $b, r$ 的所有公因數。
故最大公因數也相等：
$$
\gcd(a, b) = \gcd(b,\, a \bmod b)
$$
**參考**
- [Wikipedia: Euclidean algorithm](https://en.wikipedia.org/wiki/Euclidean_algorithm#Proof_of_correctness)
- 《Introduction to Algorithms》（CLRS）
#### 降低記憶體開銷
目前的 `mpi_gcd` 實作是每次遞迴都 `mpi_init/mpi_clear q、r`：
```c
void mpi_gcd(mpi_t rop, const mpi_t op1, const mpi_t op2)
{
if (mpi_cmp_u32(op2, 0) == 0) {
mpi_set(rop, op1);
return;
}
mpi_t q, r;
mpi_init(q); mpi_init(r);
mpi_fdiv_qr(q, r, op1, op2);
mpi_gcd(rop, op2, r);
mpi_clear(q); mpi_clear(r);
}
```
但這會出現一個問題，就是遞迴的深度月多，`heap` 裡面的暫存量就會越來越大，很明顯的這樣分配及釋放資源開銷很大，且每次遞迴都需要 `copy` 跟分配新的 `q&r` 這樣明顯浪費空間。
所以我想直接接用 `while-loop` 的寫法，這樣不單只要存一組暫存，且多次 `reuse` ，不用遞迴。
```c
void mpi_gcd(mpi_t rop, const mpi_t op1, const mpi_t op2)
{
mpi_t a, b, q, r;
mpi_init(a); mpi_init(b); mpi_init(q); mpi_init(r);
mpi_set(a, op1);
mpi_set(b, op2);
while (mpi_cmp_u32(b, 0) != 0) {
mpi_fdiv_qr(q, r, a, b);
mpi_set(a, b);
mpi_set(b, r);
}
mpi_set(rop, a);
mpi_clear(a); mpi_clear(b); mpi_clear(q); mpi_clear(r);
}
```
目前的重點就是，不管 `gcd` 做幾層，記憶體只需要分配一次，共四組 `buffer` ，不再有隨遞迴深度而改變 `heap` 的開銷， `a b q r` 這四個變數在每次迴圈都可以 `reuse` 。
而因為我目前的想法只有改掉遞迴，所以我有使用 `Chatgpt` 提供給我更多的想法 。
* **`In-place/Buffer-pool` 減少動態分配**
可以將多數 `MPI` 運算改成 `in-place buffer`（如 `mpi_fdiv_qr` 支援原地更新）。
或維護一個 `static buffer pool`。
* **`Reduce realloc/shrink` 頻率**
`mpi_compact` 壓縮記憶體，只在需要時（運算結束/最終輸出）才呼叫。
* **用 `stack buffer` 暫存小數**
對於小範圍運算，可以考慮用固定長度的 `local buffer`（如 128-bit/256-bit array 作為臨時變數），再需要大數時才用 `heap`。
**改進最大公因數的運算效率**
目前是使用 Euclidean Algorithm ，利用重複做 mod 運算，並使用遞迴或迴圈，這方法算出來的解答會是正確的，但在多精度下會有幾個效能的缺點
1. 大規模運算成本高
多精度的除法與 mod 運算，遠比加法、減法、位元位移還要慢很多，所以在位數很大的時候 mod 與除法運算可佔全部 gcd 時間的 80% 或更多。
> GNU Multiple Precision Arithmetic Library 官方手冊（Section 15.3）提到：
"Euclidean GCD involves many divisions and modulus operations, which are much more expensive than additions or subtractions, especially for very large numbers. That's why the binary GCD, which only uses shifts and subtractions, is often preferred for bignums."
> 論文 Fast Algorithms for Computing Large Integers GCD 中提到：
"The cost of a division between large integers dominates the overall cost of the Euclidean algorithm. In practice, 80% or more of the computation time is spent in division (modulo) operations."
Sign in with Wallet
Connect another wallet