Barret 約簡以獲得 128 位數字的 64 位餘數
在 github 上有Microsoft 的 SEAL 的程式碼部分:
SEAL_ITERATE(iter(operand1, operand2, result), coeff_count, [&](auto I) { // Reduces z using base 2^64 Barrett reduction unsigned long long z[2], tmp1, tmp2[2], tmp3, carry; multiply_uint64(get<0>(I), get<1>(I), z); // Multiply input and const_ratio // Round 1 multiply_uint64_hw64(z[0], const_ratio_0, &carry); multiply_uint64(z[0], const_ratio_1, tmp2); tmp3 = tmp2[1] + add_uint64(tmp2[0], carry, &tmp1); // Round 2 multiply_uint64(z[1], const_ratio_0, tmp2); carry = tmp2[1] + add_uint64(tmp1, tmp2[0], &tmp1); // This is all we care about tmp1 = z[1] * const_ratio_1 + tmp3 + carry; // Barrett subtraction tmp3 = z[0] - tmp1 * modulus_value; // Claim: One more subtraction is enough get<2>(I) = SEAL_COND_SELECT(tmp3 >= modulus_value, tmp3 - modulus_value, tmp3); });
那應該做Barrett Reduction,一種無需除法即可計算模數的技術。
它看起來像
multiply_uint64_hw64
將兩個 64 位數字相乘並只得到 64 個最高有效位。multiply_uint64
在兩個 64 位數字中獲得一個 128 位數字。但是,我不明白正在做什麼,最重要的是,在哪裡$$ a-\left\lfloor a,s\right\rfloor,n $$
發生。這段程式碼甚至沒有 floor 函式。
問題的程式碼計算 64 位
get<2>(I)
= $ h:=f,g\bmod n $ 從輸入:
- 64 位
modulus_value
= $ n $ 和 $ n\in[2,,2^{63}] $- 64 位
get<0>(I)
= $ f $ 和 $ f\in[0,,2^{64}-1] $- 64 位
get<1>(I)
= $ g $ 和 $ g\in[0,,2^{64}-1] $ 和 $ f,g<2^{64},n $ , 滿足的條件如果 $ f,g\in[0,,n-1] $ (我猜在應用程序中總是如此)。結果 $ h $ 是歐幾里得除法的餘數 $ f,g $ 經過 $ n $ . 它在數學上定義為 $ 0\le h<m $ 和 $ \exists q\in\mathbb Z,\ f,g=q\cdot n+h $ .
程式碼首先計算 128 位
z
$ =z:=f,g $ , 那麼結果get<2>(I)
$ =h:=z\bmod n $ 通過巴雷特減少:
- 它是預先計算的(在問題程式碼的外部)
const_ratio
$ =r=\left\lfloor2^{128}/n\right\rfloor $- $ \hat q:=\left\lfloor z,r/2^{128}\right\rfloor $ , 這是正確的 $ q $ 預設情況下在一個內(注意: $ \hat q $ 是變數的最終值
tmp1
)。- $ \hat h:=z-q\cdot n $ , 這是正確的 $ h $ 在可能超過精確的範圍內 $ n $ (筆記: $ \hat h $ 是變數的最終值
tmp3
)。- $ h:=\hat h-n $ 什麼時候 $ \hat h\ge n $ , 或者 $ h:=\hat h $ 否則。
該程式碼使用小學算法執行多位數算術,從基數轉置 $ 10 $ 基地 $ 2^{64} $ (在下一節中詳細介紹了優化和變體)。數字的等價物是所謂的肢體,這裡是 64 位。
產品
z
$ =z $ 表示為兩個肢體z[0]
$ =z_0 $ (低階),z[1]
$ =z_1 $ (高階),因此與 $ z=z_0+2^{64},z_1 $ , 和 $ z_0, z_1\in[0,,2^{64}-1] $ .巴雷特乘數
const_ratio
$ =r=\left\lfloor2^{128}/n\right\rfloor $ 類似地表示為兩個肢體const_ratio_0
$ =r_0 $ 和const_ratio_1
$ =r_1 $ ,由於前提條件中區間的低端 $ n\in[2,,2^{63}] $ .中間產品 $ z,r $ 總是適合三個肢體(即使通常表示為兩個肢體的兩個數量的乘積需要四個肢體)。
預設的暫定商 $ \hat q $ 適合單肢,因為 $ \hat q\le q $ 和 $ q<2^{64} $ , 輸入前提條件的最新投保人 $ f,g<2^{64},n $ .
暫定餘數 $ \hat h $ 適合單肢,因為 $ \hat h<2n $ 和 $ 2n\le2^{64} $ , 後面要歸功於前提條件中區間的高端 $ n\in[2,,2^{63}] $ .
詳細說明程式碼的算法(如評論和賞金中所述):
我將使用十進制的插圖。在那個基礎上,由於 $ 2\le n\le10/2 $ ,
const_ratio
$ =r $ 只能是 $ \left\lfloor100/2\right\rfloor=50 $ , $ \left\lfloor100/3\right\rfloor=33 $ , $ \left\lfloor100/4\right\rfloor=25 $ , $ \left\lfloor100/5\right\rfloor=20 $ ,但我會假裝const_ratio
$ =r=29 $ 因為這是一個更有趣的例子。出於同樣的原因,我將使用z
$ =z=34 $ ,即使這不能作為兩位數的乘積獲得。產品
z
在程式碼中通過multiply_uint64(get<0>(I), get<1>(I), z)
as 兩個肢體z[0]
和z[1]
.計算的核心是 $ \hat q:=\left\lfloor z,r/2^{128}\right\rfloor $ . 這是base中的模擬 $ 2^{64} $ 的 $ 9:=\left\lfloor29\cdot34/100\right\rfloor $ 以 10 為基數。兩個參數 $ 29 $ 和 $ 34 $ 乘法是兩位數,但足夠小,以至於它們的乘積 $ 986 $ 是三位(而不是四位),我們只對右邊第三位感興趣。小學算法計算 $ 986:=29\cdot34 $ 將呈現為
2 9 const_ratio x 3 4 z ----- 1 1 6 + 8 7 ------- 9 8 6
在小學算法中有四個個位數的乘法(程式碼執行)和一些額外的操作(程式碼稍微重組):
4
次9
,36
; 寫6
,保存3
;4
次2
,8
; 加3
(保留)11
,;寫那個。3
次9
,27
; 寫7
,保存2
;3
次2
,6
; 加2
(保留)8
,;寫那個。這四個乘法中的第一個發生在程式碼片段
multiply_uint64_hw64(z[0], const_ratio_0, &carry)
中,它乘以 $ r $ 與低階肢體 $ z $ ,就像我們將 的低位數字4
與34
的低位數字9
相乘一樣29
。請注意,“write6
”在這種情況下是沒有意義的,因為它寫入的任何數字都將保持在計算的右列,沒有任何機會影響最左邊的數字,並且當我們除以 100 並向下舍入時被忽略(等效地,只保留右數第三個數字)。這就是為什麼甚至不計算 128 位產品的低位 64 位的原因,如問題中所述。3
in的等價物36
保存在 in 中carry
。第二次乘法發生在 中
multiply_uint64(z[0], const_ratio_1, tmp2)
,它乘以 的高階肢體 $ r $ 與低階肢體 $ z $ ,導致兩個肢體tmp2
;64 位tmp[0]
接收等價於8
in8
,並tmp[1]
接收等價於0
for (注意在十進制整數的傳統書寫中,前導0
被抑制)。的等價物8 plus 3
出現在 中,結果add_uint64(tmp2[0], carry, &tmp1)
的低位出現在 中,新進位出現在該函式的輸出中。這被用作 in 的右操作數(在小學算法中恰好被跳過,因為我在被抑制後採用了特定範例),從而產生了與左in的等價物。[關於輸出的注意事項:它是由1``11``tmp1``1``tmp3 = tmp2[1] + …``0``1``116``add_uint64``static_cast<unsigned char>(*result < operand1)
, 比較*result
andoperand
, 然後true
變成1
,false
to0
. Made after*result = operand1 + operand2
,它告訴我們這個加法是否產生了進位。一些編譯器辨識出這個習語,使用狀態字的 C 位,並在即將添加的內容中重用 C]。第三次乘法發生在 中
multiply_uint64(z[1], const_ratio_0, tmp2)
,它乘以 的低階肢體 $ r $ 與高階肢體 $ z $ ,結果到四肢上tmp2
,就像我們一樣3 x 9 = 27
。這次我們需要兩個肢體/數字:相當於7
去到tmp2[0]
和相當於2
去到tmp2[1]
。這裡它是小學算法的一個變體:它立即添加tmp1
(相當於中間的中間1
)116
到低階肢體中add_uint64(tmp1, tmp2[0], &tmp1)
,執行相當於1 + 7 = 8
,沒有進位。結果8
被儲存在tmp1
因為add_uint64
需要一個目的地的語義,但它真的被忽略了,因為我們不關心986
. 進位輸出 byadd_uint64
用作 中的右操作數carry = tmp2[1] + …
,相當於1 + 0 = 1
我們範例中的 。儘管名字carry
,它包含一個成熟的 64 位肢體/數字。第四次乘法發生在 中
z[1] * const_ratio_1
,它乘以 的高階肢體 $ r $ 與高階肢體 $ z $ ,就像我們一樣3 x 2 = 6
。這裡的上下文確保結果適合單個肢體,因此可以使用本地 C 運算符進行乘法運算。然後將結果用作 的左運算符… + tmp3 + carry
,執行等效的6 + 1 + 1 = 8
。上下文再次確保了這個值 $ \hat q $ ,儲存在 中tmp1
,適合單個肢體/數字。然後
tmp3 = z[0] - tmp1 * modulus_value
執行 $ \hat h:=z-q\cdot n $ . 上下文確保數學上精確的結果適合單個肢體/數字(儲存在 中tmp3
),即使 $ q\cdot n $ 才不是。這允許使用本機 C 運算符,它完全跳過計算高階肢體。然後
SEAL_COND_SELECT(tmp3 >= modulus_value, tmp3 - modulus_value, tmp3)
計算 $ h $ 從 $ \hat h $ 通過有條件地減去 $ n $ 什麼時候 $ \hat h\ge n $ . 選擇運算符隱藏在宏中。基數的兩個例子 $ 2^{64} $ (在四肢之間插入空格的大端十六進制值):
modulus_value 000076513ae0b1cd const_ratio 00000000000229e6 7f4ca82ba3a115f1 get<0>(I) 00005f0fd669f2c7 get<1>(I) 000041a1f91ef16f z 00000000185f2ae8 a455846cb7cf9b49 tmp1 000034bb854f9a8d tmp3 00000fcebfd55b60 get<2>(I) 00000fcebfd55b60 modulus_value 686f4b7702a9c775 const_ratio 0000000000000002 7387d66ffd685b82 get<0>(I) 536094611fa2b19b get<1>(I) 675ef5187093ff63 z 21aac8fcf31d6421 62e675ba16d513f1 tmp1 5287278703394bb1 tmp3 72b1d3d2b9f5e50c get<2>(I) 0a42885bb74c1d97
注:對於 $ n\in[2^{63}+1,,2^{64}-1] $ , 數量 $ \hat h $ 可能會溢出一個肢體,並且目前的程式碼會失敗。例如輸入 $ f=g=2^{32} $ ,我們得到 $ z=2^{64} $ 因此 $ \hat q=0 $ (對於任何 $ n>2^{63} $ ), 因此 $ \hat h=z=2^{64} $ 和一個輸出 $ 0 $ 而不是真實的 $ h=2^{64}-n $ . 完整原始碼有註釋“Modulus 類表示高達 61 位的非負整數模數”,因此對於大型 $ n $ 僅在呼叫程式碼出錯時發生。另外,如果我理解正確, $ n=2^{60}-2^{14}+1 $ 是主要目標。
替代方案:對於與問題程式碼執行相同功能的內容,使用 4 行短程式碼而不是 11 行程式碼,適用於所有 $ n\in[1,,2^{64}-1] $ ,不需要預先計算,可能更快,但僅與最近的 x64 編譯器 + CPU 兼容,請參閱這些程式碼片段中的第一個(第二個是沒有限制的小變體 $ f,g<2^{64},n $ )。我沒有聲明任何程式碼的恆定時間性。