Skip to content

Commit

Permalink
feat(reductions): add barrett and montgomery reductions
Browse files Browse the repository at this point in the history
  • Loading branch information
dynm committed Oct 30, 2024
1 parent ee4ed6e commit 45d08d3
Show file tree
Hide file tree
Showing 2 changed files with 178 additions and 0 deletions.
92 changes: 92 additions & 0 deletions reductions/barrett.zh.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,92 @@
# Barrett Reduction

Barrett Reduction是一种在计算模运算(求余数)时使用的优化算法,可以避免昂贵的除法运算。它的基本思想相当直接。

当要计算 $a \bmod N$(a 模 N)的余数 r 时,常规做法是首先计算商 $q = a \div N$,然后计算余数 $r = a - q \times N$。

然而,计算 q 时需要一次除法,而除法对 CPU 的运算开销比较大。Barrett Reduction 通过引入一个预计算的常数来避免这个除法运算。

## 算法核心思想

1. 选择一个足够大的 $R = 2^k$,使得 $R > N$。
2. 预计算 $\mu = \lfloor R / N \rfloor$。($\lfloor \rfloor$ 表示向下取整)
3. 当需要计算 $a \bmod N$ 时,用以下步骤代替:
- 计算 $q = \lfloor(a \times \mu) / R\rfloor$
- 计算 $r = a - q \times N$
- 如果 $r \geq N$,则 $r = r - N$

这里的关键优化是,由于 R 是 2 的幂,$(a \times \mu) / R$ 可以通过简单的位移操作实现:$q = (a \times \mu) \gg k$。

对于一个固定的有限域,N 通常是固定的,那么 $\mu$ 也是固定的,可以预先计算并存储。

## 误差分析

这种方法可能会导致 q 略小于实际的商,因为 $\mu$ 是 $R/N$ 的下界(由于向下取整)。这就是为什么我们需要在计算 $r = a - q \times N$ 后,判断 r 是否小于 N,如果大于或等于 N 则需要再减去 N。

关于最差情况下需要减去 N 的次数:

- 误差的上界可以被证明是 2。也就是说,计算出的 r 最多比实际的余数大 2N。
- 因此,最多需要减去 N 两次就可以得到正确的余数。

## 算法优势

1. 避免了昂贵的除法运算,替换为乘法和位移操作。
2. 对于固定的模数 N,$\mu$ 可以预计算,进一步提高效率。
3. 即使在最坏情况下,也只需要进行有限次(最多两次)的额外减法。

## 计算μ的误差范围

我们设 $q'$ 为实际的 $a/N$ 值,q 为经过 Barrett Reduction 的值
$q' = \lfloor a / N \rfloor$,$q = \lfloor(a \times \mu) / R\rfloor$

1. $R = 2^k$,其中 k 是一个整数,且 $R > N$。
2. $\mu = \lfloor R / N \rfloor$
3. 假设 $a < R$(这是一个重要的假设,在实际应用中通常成立)

### 上界:$q \leq q'$

我们知道 $\mu = \lfloor R / N \rfloor$,所以 $\mu \leq R / N$。

$$
\begin{aligned}
q &= \lfloor(a \times \mu) / R\rfloor \\
&\leq \lfloor(a \times (R/N)) / R\rfloor \quad \text{(因为} \mu \leq R/N\text{)} \\
&= \lfloor a / N \rfloor \\
&= q'
\end{aligned}
$$

所以 $q \leq q'$

### 下界:$q > q' - 2$

首先,我们知道 $R/N - 1 < \mu \leq R/N$ (因为 $\mu$ 是 $R/N$ 的向下取整)
所以:

$$
\begin{aligned}
q &= \lfloor(a \times \mu) / R\rfloor \\
&> \lfloor(a \times (R/N - 1)) / R\rfloor \quad \text{(因为} \mu > R/N - 1\text{)} \\
&= \lfloor a/N - a/R \rfloor
\end{aligned}
$$

接下来,我们可以写出:
$q > a/N - a/R - 1$ (因为对任何 x,$\lfloor x \rfloor > x - 1$)

现在,利用假设 $a < R$,我们有 $a/R < 1$,所以:

$$
\begin{aligned}
q &> a/N - 1 - 1 \\
&= a/N - 2
\end{aligned}
$$

但是 $q' = \lfloor a/N \rfloor$,所以 $a/N - 1 < q' \leq a/N$
将这个不等式代入上面的式子:
$q > (q' - 1) - 2 = q' - 3$

因为 q 是整数,所以 $q \geq q' - 2$

结合之前得到的 $q \leq q'$,我们可以得出 $q' - 2 < q \leq q'$
86 changes: 86 additions & 0 deletions reductions/montgomery.zh.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,86 @@
# Montgomery Reduction

在有限域运算中,我们需要频繁使用mod N运算。对于涉及连乘的指数运算,Montgomery Reduction能提供更好的性能优化。

Montgomery Reduction的核心思想是将运算数转换到一个特殊的Montgomery Form,其所有的除法运算都可以被显著简化。

让我们通过一个具体例子来理解:

设模数 N = 97,我们需要计算 $(a \cdot b) \bmod N$。传统方法需要对$a \cdot b$的结果进行求余运算,这涉及到开销较大的除法运算。

### Montgomery Form转换
1. 选择R = 100(通常选择为2的幂,这里为便于演示选用100,我们仅需在10进制下取尾数和抹除尾部的00即可)
2. 对于输入 a = 15, b = 32:
- $T_a = aR \bmod N = 45$
- $T_b = bR \bmod N = 96$

### Montgomery乘法过程
1. 计算Montgomery Form的乘积:
$T_r = T_a \cdot T_b = 45 \cdot 96 = 4320$

2. 结果转换:
- 由于经过了两次Montgomery Form转换,结果包含 $R^2$ 项,需要消除一个 $R$
- 我们需要在不改变 $mod N$ 结果的情况下调整 $T_r$ ,使其能被 $R$ 整除
- 我们可以给T加上若干次 $N$ ,这样$mod N$的结果仍然不会被影响,于是我们有了 $T + xN$
- 我们又需要合适的x,使得 $T + xN$ 能被 $R$ 整除
- 通过同余等式 $T + xN \equiv 0 \pmod{R}$,即 $x = T \cdot (-\frac{1}{N}) \bmod R$
- 我们把 $-\frac{1}{N}$ 记为 $N'$,$N' = -\frac{1}{N} \bmod R = 67$(提前计算)
- 上面求得的 $x$ 即为Montgomery Reduction中的 $m$, $t = T + mN$

3. 具体计算:
第一次Montgomery Reduction ($aR \cdot bR/R$):

$m = T_r \cdot N' \bmod R = 4320 \cdot 67 \bmod 100 = 289440 \bmod 100 = 40$
(注:mod 100在10进制下就是取2位尾数)

$t = (T_r + m \cdot N) / R = (4320 + 40 \cdot 97) / 100 = 8200 / 100 = 82 = abR$
(注:除以100在10进制下就是抹掉尾部的00)

第二次Montgomery Reduction ($abR/R$):

$m' = t \cdot N' \bmod R = 82 \cdot 67 \bmod 100 = 94$

$t' = (t + m' \cdot N) / R = (82 + 94 \cdot 97) / 100 = 92 = ab$

### 优势分析
虽然单次运算时Montgomery Reduction似乎并未节省太多计算资源(转换到Montgomery Form时仍需mod N运算),但在需要连续模乘运算的场景下(如模幂运算),其优势就非常明显:

1. 模幂运算示例:计算 $a^5 \bmod N$
- 传统方法需要在每次乘法后进行mod N运算
- 使用Montgomery Form后:
- 初始转换:$T_a = aR \bmod N$
- 中间运算无需除法:$T_a \cdot T_a / R = (aR \cdot aR) / R = a^2R$
- 只需在最后转换回原始形式

2. 主要优势:
- 避免中间步骤的mod N运算
- 消除大部分除法运算
- 适合硬件实现(R选择2的幂时,除法和取模都可以用移位和位与运算完成)
- 特别适合需要连续模乘运算的场景

### 工程实现
在实际工程中,我们通常通过将元素a乘以$R^2$来将其编码到Montgomery Form。这种方式的优点是:

1. Montgomery Reduction和乘法可以合并为统一的mul函数:
- $\text{mul}(a, b) = ab/R$

2. 这样encode和decode都可以用同一个mul函数:
- $\text{encode}(a) = \text{mul}(a, R^2) = aR^2/R = aR$
- $\text{decode}(T_a) = \text{mul}(T_a, 1) = aR/R = a$

这种统一的实现方式使代码更简洁,同时当R选择为2的幂时,所有的除法和取模运算都可以通过移位和位与运算高效完成。

对于模幂运算,我们可以看到Mont乘的优势:
1. 传统算法:
$a^5 = ((((a \cdot a) \bmod N) \cdot a) \bmod N) \cdot a) \bmod N) \cdot a) \bmod N$
每一步都需要进行昂贵的除法运算来求余

2. Montgomery算法:
- $T_a = aR \bmod N$
- $T_a \cdot T_a / R = a^2R$
- $(a^2R \cdot aR) / R = a^3R$
- $(a^3R \cdot aR) / R = a^4R$
- $(a^4R \cdot aR) / R = a^5R$
- $\text{decode}(a^5R) = a^5$

所有中间步骤都只需要简单的移位操作,大大提升了性能。

0 comments on commit 45d08d3

Please sign in to comment.