一. Barrett模乘
1.1.前言
Barrett算法是利用移位代替除法。
举个例子:已知X M,求Z(三者均为正整数,Z<M)
Z=X mod M
- 方法一:
一个很简单的方法就是用X=X-M,一直减下去,直到最后结果<M。
- 方法二
我们令X=kM+Z ,我们只需求出k,然后用X减kM即可。利用高斯取整函数k=[X/M],但是其中用了除了法。我们要避免除法这种长时间运算的计算,所以出现了一个叫做Barrett算法,就是利用移位来代替除法来求得一个k’,(这里的移位是任意进制下的)。经Barrett计算得到的k’要么等于k,要么等于k-1
1.2.验证
python3:
import random
def Barrett(x,n,mod):
r=10#按论文,这里应取2的次幂,但python里不好进行二进制处理,就只好取10了(\悲哀)
z=0
u=int(r**(2*n+1)/mod)
k1=int(x/(r**(n-2)))
k2=int( k1*(u/r**(n+3)))
if(k2!=int(x/mod)):
print(k2,int(x/mod))
z=x-k2*mod
while z>=mod:
z=z-mod
return z
if __name__=="__main__":
mod=16891
for i in range(1000):
x=random.randint(30000,99999)
right=x%mod
z=Barrett(x,5,mod)
if(right!=z):
print("fuck")
verilog:
`include "param.v"
module MOD_X(
input wire[`SIZE:0] x,
output wire [`SIZE-1:0] mod_x
);
wire [`SIZE-1:0] z1;
wire [`SIZE-1:0] z2;
wire [`SIZE-1:0] z3;
wire [`SIZE-1:0] k1;
wire [`SIZE-1:0] k2;
assign k1=x >>(`SIZE);
assign k2= k1*(`BARRETT_U>>(`SIZE));
assign z1=x-k2*`P;
assign z2=z1-`P;
assign z3=z2-`P;
assign mod_x=(z1>=`P)? (z2>=`P? ((z3>=`P)?z3-`P:z3):z2 ) :z1;endmodule
param.v:
`define BARRETT_R 2'b10
`define BARRETT_U 8736
有中间变量z1 z2 z3是因为参数的选择,代码中参数的选择,根据式子,到少要判断三次,当然可以选择其他参数达到只判断一次的效果。但是乘法位数将会更大,所以自己折中吧。
1.3.证明
需要注意的是,最后一张图片的式子有误,第5行应为:
Z<——Z-k'M
没有链接,只好选择原创了,侵删,论文:NTT处理器的研究与实现_宋鹏飞
二. Montgomery模乘
2.1 算法原理
以a*b mod q为例
计算 T 0 = a b ⋅ q − 1 m o d R T_0=ab\cdot q^{-1}\mod R T0=ab⋅q−1modR
计算 T 1 = T 0 ⋅ q T_1=T_0\cdot q T1=T0⋅q
计算 T 2 = a b − T 1 T_2=ab-T_1 T2=ab−T1
计算 T 3 = T 2 R T_3=\frac{T_2}{R} T3=RT2
注意:第1步必须mod R !
上述四步得到结果为 a b ⋅ R − 1 m o d q ab\cdot R^{-1}\mod q ab⋅R−1modq,下面将 R 2 R^2 R2乘上后再进行计算计算 T 4 = T 3 ⋅ ( R 2 m o d q ) ⋅ q − 1 m o d R T_4=T_3\cdot (R^2\mod q)\cdot q^{-1}\mod R T4=T3⋅(R2modq)⋅q−1modR
计算 T 5 = T 4 ⋅ q T_5=T_4\cdot q T5=T4⋅q
计算 T 6 = T 3 ⋅ ( R 2 m o d q ) − T 5 T_6=T_3\cdot (R^2\mod q)-T_5 T6=T3⋅(R2modq)−T5
计算 T 7 = T 6 R T_7=\frac{T_6}{R} T7=RT6
注意:第5步必须mod R !
下面举个实例:a = 123, b = 456, q = 677, q − 1 = 613 q^{-1}=613 q−1=613, R=1000, R 2 m o d q = 71 R^2\mod q=71 R2modq=71(注:因为这里的数是十进制表示,所以取R是一千,而不是 2 s o m e t h i n g 2^{something} 2something,原理一样)
- 计算 T 0 = 123 ∗ 456 ⋅ 613 m o d 1000 = 944 T_0=123*456\cdot 613\mod 1000=944 T0=123∗456⋅613mod1000=944
- 计算 T 1 = 944 ⋅ 677 = 639088 T_1=944\cdot 677=639088 T1=944⋅677=639088
- 计算 T 2 = 123 ∗ 456 − 639088 = − 583000 T_2=123*456-639088=-583000 T2=123∗456−639088=−583000
- 计算 T 3 = − 583000 1000 = − 583 T_3=\frac{-583000}{1000}=-583 T3=1000−583000=−583
- 计算 T 4 = − 583 ∗ 71 ⋅ 613 m o d 1000 = − 909 T_4=-583*71\cdot 613\mod 1000=-909 T4=−583∗71⋅613mod1000=−909
- 计算 T 5 = − 909 ⋅ 677 = − 615393 T_5=-909\cdot 677=-615393 T5=−909⋅677=−615393
- 计算 T 6 = − 583 ∗ 71 − ( − 615393 ) = 574000 T_6=-583*71-(-615393)=574000 T6=−583∗71−(−615393)=574000
- 计算 T 7 = 574000 1000 = 574 T_7=\frac{574000}{1000}=574 T7=1000574000=574
print(123*456%677)
结果与python直接计算一致
2.2 算法描述与硬件实现
摘自:Verbauwhede, Ingrid MR, ed. Secure integrated circuits and systems. Berlin: Springer, 2010.
蒙哥马利模乘算法是最常用的一种模乘算法。计算在蒙哥马利域中进行,蒙哥马利域定义为:对于元素
a
∈
F
a\in F
a∈F和
R
(
p
<
2
R
)
R(p<2^R)
R(p<2R),有映射
a
→
a
⋅
2
R
m
o
d
p
a\rightarrow a\cdot 2^R\mod p
a→a⋅2Rmodp。蒙哥马利域允许仅用乘法进行有效约简。但在计算之前,所有的输入值必须转换为蒙哥马利域(并在计算出结果后再次转换回去),增加了计算前后的额外复杂度。该方法的优点是可以减少成本高昂的约简运算,将其替换为除以2(位移)操作。假定X、Y为蒙哥马来利坐标中的两个因子,即
X
^
=
X
⋅
2
R
m
o
d
p
\hat{X}=X\cdot 2^R\mod p
X^=X⋅2Rmodp和
Y
^
=
Y
⋅
2
R
m
o
d
p
\hat{Y}=Y\cdot 2^R\mod p
Y^=Y⋅2Rmodp,可使用标准乘法计算
Z
^
=
X
^
⋅
Y
^
=
X
Y
⋅
2
2
R
\hat{Z}=\hat{X}\cdot\hat{Y}=XY\cdot 2^{2R}
Z^=X^⋅Y^=XY⋅22R。需要注意的是,该计算结果既不是蒙哥马利域,也不是标准域需要进行修正。在使用特殊蒙哥马利模乘算法,用
X
⋅
Y
⋅
2
−
R
m
o
d
p
X\cdot Y\cdot 2^{-R}\mod p
X⋅Y⋅2−Rmodp替代
X
⋅
Y
m
o
d
p
X\cdot Y\mod p
X⋅Ymodp体时需要考虑到该情况。
需要指出的是,当涉及多个重复的模乘运算时,可以忽略蒙哥马利计算的额外转换步骤,如进行模幂运算的情况。这些局部乘法运算结果会按照最低有效位到最高有效位的顺序连续相加。在每次迭代时判定中间结果是奇数还是偶数。因此,可以检查中间结果中的最低有效位,如果此位等于"1",则模数会加到中间和上,这就确保了和总是偶数。在每次迭代结束时,中间结果会除以2,这样就可以避免增加中间结果规模的复杂度。
该算法要求每次循环进行两次相加。通过引入冗余表示,改进算法,构建包含CSA加法器的蒙哥马利模乘的高效架构