## PQC Final Project ### 2025 Spring B10902082 林秉軒 --- ## Results ---- 相乘兩個 $N$-bit 的數字可以做到比 `mpz_mul` 快大約兩倍（在 $N = 2^{25}$ 的時候）  （因為 $2^{25}$ 真的太大，此處的 `NWARMUP`, `NITERATIONS`, `NTESTS` 分別只有 1, 6, 10） ---- 在 $N = 2^{21}$ 可以略微贏過 `mpz_mul`  （此處的 `NWARMUP`, `NITERATIONS`, `NTESTS` 分別是 10, 60, 100，需要跑 8 分鐘多） ---- $N = 2^{22}$  ---- 有寫測試，所以計算結果應該(?)是對的  --- ## 整體架構 - 一個 limb 為 32-bit（$R = 2^{32}$） - 要相乘兩個 $N$-limb 的數字，先計算不做進位的捲積 - 捲積的每一項 $c_k = \sum _ {i+j=k} a_i b_j < N \cdot R^2$ - 分別模 $p_1, p_2, p_3$ 三個質數計算長度 $2N$ 的捲積 $\langle c_{k} \mod p_1 \rangle, \langle c_{k} \mod p_2 \rangle, \langle c_{k} \mod p_3 \rangle$ - 根據中國剩餘定理，只要 $p_1 \cdot p_2 \cdot p_3 \geq N \cdot R^2$ 就可以唯一還原出 $c_k$ - 用 Garner's algorithm 一邊還原一邊計算進位 ----  選定大約是 30-bit 且可以做夠長的 NTT 的三個質數（理論上這份 code 可以計算兩個 $2^{27}$-bit 的數相乘） --- ## NTT 的優化們 --- ### 採用不預先 bitrev permute 的寫法 讀取記憶體很慢(?)  每次到新的一整塊，butterfly 要乘的因子會乘上一個 $\omega ^k$。只需要預處理 $\log (N)$ 種 $\omega$ 的次方，因為 $k$ 一定是二進位 `00...01*** - 11...10***` ---- ### radix-4 的變換 手動展開兩層 butterfly $$ \begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \end{bmatrix} \mapsto \begin{bmatrix} 1 & 1 & 1 & 1 \\ 1 & -1 & 1 & -1 \\ 1 & i & -1 & -i \\ 1 & -i & -1 & i \\ \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \omega \\ a_2\omega^2 \\ a_3 \omega^3 \end{bmatrix} $$ 做的乘法、加減法次數和做兩次 radix-2 幾乎一樣。也許對 cache 以及指令重排有幫助 ---- 上面兩個優化參考自 https://github.com/atcoder/ac-library/blob/master/atcoder/convolution.hpp （基本上是以此 library 為基礎開始改的） --- ### Montgomery reduction - 事先把 $a, b$ 以及所有 NTT 要乘的那些單位根變換到 montgomery domain - 把幾乎所有乘法取模運算都改成在 montgomery domain 下的 `REDC` 操作 - 變換到 montgomery domain 只需要做 $O(N)$ 次，且都是簡單的 sequential access ---- ### vectorize - 一次做 `uint32x4_t` 的 `REDC`、加減法取模 - 一個 register 只有 128-bit，所以要做 64-bit 乘法得分成 low 和 high - radix-4 和 radix-2 都可以 vectorize ----  採用 unsigned reduction $T \mapsto (T + (-M^{-1}\cdot T\mod R) \cdot M) / R$ --- ### Incomplete NTT - 最後幾層中，每一塊的大小太小而不能 vectorize - 因此就不分解最後幾層，而是直接 schoolbook multiplication 算 $a, b$ 每一塊在 $\mathbb{Z}_p[x]/(x^{2^l}-\omega ^k)$ 逐塊相乘的捲積 - 此乘法很好 vectorize ---- ### Schoolbook multiplication 的優化 假設是模 $(x^4 - w)$  ----  ---- - `res_64` 的值域大小是 $4 \cdot M^2$，小於 $2^{64}$（我們選的 $M < 2^{30}$），所以可以在 `uint64_t` 裡面直接加，最後再一次 `REDC`（Montgomery 只要求輸入要小於 $R \cdot M$） - 此形式的乘法很容易 vectorize，且實際執行非常之快速 - 實際上我選擇的是不分解最後三層，改計算模 $(x^8 - w)$ 的逐塊相乘 - 需要在最後化約一次才能保證小於 $R \cdot M$ - 如此一來 NTT 的所有部份都被 vectorize 了(?) ---- 以上兩點優化的想法來自 https://judge.yosupo.jp/submission/270589 （這份 code 在 x86 上可以在更小的 size 就贏過 `mpz_mul`） --- ### 把乘以常數用 `vqrdmulh` 做 - 瓶頸是 NTT 的正向變換 - 每一層的變換都有大量的跟某個常數做乘法取模的操作 - 按照課堂上的簡報的兩次 mulhi 一次 mullo 的作法：`REDC(a * b) = mulhi(a, b) - mulhi(M, mullo(a, b'))` 其中 $b' = (bM^{-1} \mod R)$ ----  注意 `vqrdmulhq_s32` 有乘兩倍所以需要除回來（shift） ---- - 直接 `REDC`：`vmull_u32 * 2 + vmulq_u32 + vmlal_u32 * 2 + vshrn_n_u64 * 2` - 乘常數 `REDC`：預處理需要一個 `uint32_t` 乘法，之後只需要 `vmulq_u32 + vqrdmulhq_s32 * 2 + vshrq_n_s32` ---- ### template & constexpr & Ofast --- ### 可能的進一步優化 目前的 code 中，大部份的時候都有把數字化約到 $[0, M)$，這需要花費不少的 conditional add/sub 或 min + add/sub，也許可以在某些地方放寬到 $2M$ 或 $4M$ 之類的以減少這些化約 --- end