Modular-Arithmetic

Barret 約簡以獲得 128 位數字的 64 位餘數

  • September 18, 2021

在 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&lt;1&gt;(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&lt;2&gt;(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&lt;0&gt;(I), get&lt;1&gt;(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

在小學算法中有四個個位數的乘法(程式碼執行)和一些額外的操作(程式碼稍微重組):

  • 4936; 寫6,保存3
  • 428; 加3(保留)11,;寫那個。
  • 3927; 寫7,保存2
  • 326; 加2(保留) 8,;寫那個。

這四個乘法中的第一個發生在程式碼片段multiply_uint64_hw64(z[0], const_ratio_0, &carry)中,它乘以 $ r $ 與低階肢體 $ z $ ,就像我們將 的低位數字434的低位數字9相乘一樣29。請注意,“write 6”在這種情況下是沒有意義的,因為它寫入的任何數字都將保持在計算的右列,沒有任何機會影響最左邊的數字,並且當我們除以 100 並向下舍入時被忽略(等效地,只保留右數第三個數字)。這就是為什麼甚至不計算 128 位產品的低位 64 位的原因,如問題中所述。3in的等價物36保存在 in 中carry

第二次乘法發生在 中multiply_uint64(z[0], const_ratio_1, tmp2),它乘以 的高階肢體 $ r $ 與低階肢體 $ z $ ,導致兩個肢體tmp2;64 位tmp[0]接收等價於8in 8,並tmp[1]接收等價於0for (注意在十進制整數的傳統書寫中,前導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&lt;unsigned char&gt;(*result &lt; operand1), 比較*resultand operand, 然後true變成1, falseto 0. 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(相當於中間的中間1116到低階肢體中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 &gt;= modulus_value, tmp3 - modulus_value, tmp3)計算 $ h $ 從 $ \hat h $ 通過有條件地減去 $ n $ 什麼時候 $ \hat h\ge n $ . 選擇運算符隱藏在宏中。

基數的兩個例子 $ 2^{64} $ (在四肢之間插入空格的大端十六進制值):

modulus_value                      000076513ae0b1cd
const_ratio       00000000000229e6 7f4ca82ba3a115f1
get&lt;0&gt;(I)                          00005f0fd669f2c7
get&lt;1&gt;(I)                          000041a1f91ef16f
z                 00000000185f2ae8 a455846cb7cf9b49
tmp1                               000034bb854f9a8d
tmp3                               00000fcebfd55b60
get&lt;2&gt;(I)                          00000fcebfd55b60

modulus_value                      686f4b7702a9c775
const_ratio       0000000000000002 7387d66ffd685b82
get&lt;0&gt;(I)                          536094611fa2b19b
get&lt;1&gt;(I)                          675ef5187093ff63
z                 21aac8fcf31d6421 62e675ba16d513f1
tmp1                               5287278703394bb1
tmp3                               72b1d3d2b9f5e50c
get&lt;2&gt;(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 $ )。我沒有聲明任何程式碼的恆定時間性。

引用自:https://crypto.stackexchange.com/questions/94972