一. Barrett模乘

1.1.前言

Barrett算法是利用移位代替除法。
举个例子：已知X M，求Z(三者均为正整数，Z<M)
Z=X mod M

1. 方法一：

一个很简单的方法就是用X=X-M，一直减下去，直到最后结果<M。

2. 方法二

我们令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.证明

需要注意的是，最后一张图片的式子有误，第５行应为：

Z<——Z-k'M

没有链接，只好选择原创了，侵删，论文：NTT处理器的研究与实现_宋鹏飞

二. Montgomery模乘

2.1 算法原理

算法原理可以参见：https://blog.csdn.net/a675115471/article/details/107553091?utm_medium=distribute.pc_relevant.none-task-blog-2_defaultbaidujs_baidulandingword~default-4.opensearchhbase&spm=1001.2101.3001.4242.3
这里，给出算法原文献：

以a*b mod q为例

1. 计算 $T_0=ab\cdot q^{-1}\mathrm{m}\mathrm{o}\mathrm{d}R$

2. 计算$T_1=T_0\cdot q$

3. 计算$T_2=ab-T_1$

4. 计算$T_3=\frac{T_2}{R}$
   注意：第1步必须mod R !
   上述四步得到结果为$ab\cdot R^{-1}\mathrm{m}\mathrm{o}\mathrm{d}q$，下面将$R^2$乘上后再进行计算

5. 计算 $T_4=T_3\cdot (R^2\mathrm{m}\mathrm{o}\mathrm{d}q)\cdot q^{-1}\mathrm{m}\mathrm{o}\mathrm{d}R$

6. 计算$T_5=T_4\cdot q$

7. 计算$T_6=T_3\cdot (R^2\mathrm{m}\mathrm{o}\mathrm{d}q)-T_5$

8. 计算$T_7=\frac{T_6}{R}$
   注意：第5步必须mod R !

下面举个实例：a = 123, b = 456, q = 677, $q^{-1}=613$, R=1000, $R^2\mathrm{m}\mathrm{o}\mathrm{d}q=71$(注：因为这里的数是十进制表示，所以取R是一千，而不是$2^{something}$，原理一样)

1. 计算 $T_0=123\ast 456\cdot 613\mathrm{m}\mathrm{o}\mathrm{d}1000=944$
2. 计算$T_1=944\cdot 677=639088$
3. 计算$T_2=123\ast 456-639088=-583000$
4. 计算$T_3=\frac{-583000}{1000}=-583$
5. 计算 $T_4=-583\ast 71\cdot 613\mathrm{m}\mathrm{o}\mathrm{d}1000=-909$
6. 计算$T_5=-909\cdot 677=-615393$
7. 计算$T_6=-583\ast 71-(-615393)=574000$
8. 计算$T_7=\frac{574000}{1000}=574$

print(123*456%677)

结果与python直接计算一致

2.2 算法描述与硬件实现

摘自：Verbauwhede, Ingrid MR, ed. Secure integrated circuits and systems. Berlin: Springer, 2010.
蒙哥马利模乘算法是最常用的一种模乘算法。计算在蒙哥马利域中进行，蒙哥马利域定义为：对于元素$a\in F$和$R(p<2^R)$，有映射$a\to a\cdot 2^R\mathrm{m}\mathrm{o}\mathrm{d}p$。蒙哥马利域允许仅用乘法进行有效约简。但在计算之前，所有的输入值必须转换为蒙哥马利域（并在计算出结果后再次转换回去），增加了计算前后的额外复杂度。该方法的优点是可以减少成本高昂的约简运算，将其替换为除以2（位移）操作。假定X、Y为蒙哥马来利坐标中的两个因子，即$\hat{X}=X\cdot 2^R\mathrm{m}\mathrm{o}\mathrm{d}p$和$\hat{Y}=Y\cdot 2^R\mathrm{m}\mathrm{o}\mathrm{d}p$，可使用标准乘法计算$\hat{Z}=\hat{X}\cdot \hat{Y}=XY\cdot 2^{2R}$。需要注意的是，该计算结果既不是蒙哥马利域，也不是标准域需要进行修正。在使用特殊蒙哥马利模乘算法，用$X\cdot Y\cdot 2^{-R}\mathrm{m}\mathrm{o}\mathrm{d}p$替代$X\cdot Y\mathrm{m}\mathrm{o}\mathrm{d}p$体时需要考虑到该情况。

需要指出的是，当涉及多个重复的模乘运算时，可以忽略蒙哥马利计算的额外转换步骤，如进行模幂运算的情况。这些局部乘法运算结果会按照最低有效位到最高有效位的顺序连续相加。在每次迭代时判定中间结果是奇数还是偶数。因此，可以检查中间结果中的最低有效位，如果此位等于"1"，则模数会加到中间和上，这就确保了和总是偶数。在每次迭代结束时，中间结果会除以2，这样就可以避免增加中间结果规模的复杂度。

该算法要求每次循环进行两次相加。通过引入冗余表示，改进算法，构建包含CSA加法器的蒙哥马利模乘的高效架构