一、搞要
Montgomery大整数乘法 在fpga中实现资源占用高，特别是DSP资源，而且一个算法包含好多个Montgomery大整数乘法，fpga中的dsp资源总是很快就耗尽。本文讲述了如何改进算法，减少dsp使用量，并优化逻辑，提高实现的最终频率。
二、原理简述
2.1、蒙哥马利(Montgomery)模乘算法
Montgomery乘法是一种用于加速模运算的算法，其数学表达式为，其中A和B是与M同位长的大数，R是
（M为位长），
是R相对于M的模逆，即满足
。这种算法特别适用于对奇数求模的情况，因为它通过使用蒙哥马利表示法，可以在计算过程中避免直接的除法运算，从而提高计算效率。Montgomery乘法通常用于加密算法中，如椭圆曲线密码学，以加速模乘和模逆运算
下面涉汲到的Montgomery模乘算法都是指密码学中的大整数模乘算法，所谓的大整数得指位数达几百位的大数，例如256位的大整数。
2.2、在fpga中如何实现Montgomery模乘法算法
公式推导及数学原理本文不再详述。我这里就直接给出算法的实现公式。因在fpga要实现除法，耗资源，速度也慢，因而要把模乘中的除法转换成模域中的乘法，如何转换不是本文的重点。现直接给出Montgomery乘法公式。下面讲述的模乘运算基于Bn128曲线的模及模逆
模P：P=256'h73eda753299d7d483339d80809a1d80553bda402fffe5bfeffffffff00000001
模逆Q：Q=256'h3d443ab0d7bf2839181b2c170004ec0653ba5bfffffe5bfdfffffffeffffffff
Montgomery模乘算法 C=A*B mod P在FPGA中分以下4步计算(公式推导不再详述)：
1、计算中间值Ma=A*B；
2、计算中间值Mb=Ma*Q;
3、计算中间值Mc=Mb*P;
4、计算结果C=(Ma+Mc)modP;
Montgomery模乘算法完美的去除了除法，只剩下3个大整数乘法和一个大整数模加，完全可以在fpga中实现全流水线计算。
然而要实现这3个大整数乘法，单个模乘就要要消耗200多个DSP48资源，因而想办法改进算法，还是要进一步减少DSP48资源的使用量。
三、Montgomery模乘算法在fpga中优化方法
Montgomery模乘算法中的模加完全是用fpga中逻辑资源实现，在这里不再详细讲述优化方法。
三个大正数乘法的在算法方面优化改进方法如下：
1、采用karatsuba乘法减少乘法次数；
2、采用CSD编码乘法器算法把乘法改用加法减法逻辑实现；
3、当A=B时，可以进行单独优化，减少逻辑资源的使用；
4、当B=常数时，也可以进一优化，
光采用karatsuba乘法优化，要实现三个大整数乘法要消耗DSP48资源210个左右，当采用前两种方法时，即采用karatsuba乘法+CSD编码乘法，则可以降低到180个DSP48资源左右。
注意：以上数据都是基于xilinx fpga Virtex ultraScale+ VU19P上综合结果。
四、karatsuba乘法原理
Karatsuba乘法是一种快速的大数乘法算法，它的原理是利用分治法将两个大数相乘。
设要计算的两个大数为A和B，分别表示为和
，其中a、b、c、d是相应位数的小数。根据乘法的分配律，可以将A * B表示为：
在上式中，a * c、b * d和(a * d + b * c)是三个乘法运算，其中(a * d + b * c)是两个小数相乘，需要计算两次。为了减少乘法的次数，Karatsuba乘法采用如下的递归思想：
- 将A和B都分成两部分：
，
- 计算a * c，b * d和(a * d + b * c)
- 通过减少两次大数乘法的次数来计算A * B：
- 计算a * c的乘积，记为P1
- 计算b * d的乘积，记为P2
- 计算(a * d + b * c)的乘积，记为P3
- 计算
通过递归地应用这个过程，可以快速计算出大数的乘积。
Karatsuba乘法的关键思想是通过分解大数，减少乘法的次数。相较于传统的乘法算法，它具有更高的效率。然而，在实际应用中，Karatsuba乘法的效率在乘数长度较小时可能并不明显，因此通常在处理大数乘法时使用。
// 输入：n位长度的正整数x和y
// 输出：x*y的结果
函数 karatsuba (x, y) {
a和b为x的前半部分和后半部分；
c和d是y的前半部分和后半部分；
递归调用函数karatsuba解决q=a*c， 递归karatsuba解决p=b*d，递归karatsuba解决k=(a+b)(c+d)
返回10的n次方*q + 10的n/2次方 * (k - q - p) + p
}4.1、karatsuba算法实现的verilogHDL代码分为三级嵌套
256位乘法器分为256位，128位，64位等三级嵌套。如下图
y_mult_64：64位乘法器采用karatsuba算法实现的代码；
y_mult_128：128位乘法器，它包含三个y_mult_64模块；
y_mult_256：256位乘法器，它包含三个y_mult_128模块。
所以一个256位乘法器包含了3*3=9个64位乘法器。
4.2、 64位大整数乘法采用karatsuba算法实现的verilogHDL代码
64位大整数乘法，采用karatsuba算法在Xilinx fpga Virtex ultraScale+ VU19P 上的实现的最后一级64位乘法代码如下：
//
// Company: 杭州幻智星科技有限公司
// Engineer: 俞工
// email: yupc123@qq.com
// Create Date: 01/25/2021 03:45:32 PM
// Design Name: Karatsuba 64位乘法器
// Module Name: y_mult_64
// Revision 0.01
//
module y_mult_64#(
parameter AB_WIDTH = 64,
parameter DSP_STEP = 6,
parameter RESULT_WIDTH = AB_WIDTH*2
)(
input clk,
input [AB_WIDTH-1:0] a,
input [AB_WIDTH-1:0] b,
output reg [RESULT_WIDTH-1:0] result // delay 9 clk
);
wire [AB_WIDTH-1:0] z2;
wire [AB_WIDTH-1:0] z0;
wire [AB_WIDTH+1:0] z1;
reg [RESULT_WIDTH-1:0] tmp1;
reg [RESULT_WIDTH:0] tmp1_;
reg [AB_WIDTH/2:0] z1_tmp3,z1_tmp4;
reg [AB_WIDTH+1:0] z2_z0_tmp[1:0];
generate
if(AB_WIDTH==64)begin
mult32 mult32_u0 (
.CLK(clk), // input wire CLK
.A(a[AB_WIDTH-1:AB_WIDTH/2]), // input wire [31 : 0] A
.B(b[AB_WIDTH-1:AB_WIDTH/2]), // input wire [31 : 0] B
.P(z2) // output wire [63 : 0] P // delay 6 clk
);
mult32 mult32_u1 (
.CLK(clk), // input wire CLK
.A(a[(AB_WIDTH/2) -1:0]), // input wire [31 : 0] A
.B(b[(AB_WIDTH/2) -1:0]), // input wire [31 : 0] B
.P(z0) // output wire [63 : 0] P
);
mult33 mult33_u0 (
.CLK(clk), // input wire CLK
.A(z1_tmp3), // input wire [31 : 0] A
.B(z1_tmp4), // input wire [31 : 0] B
.P(z1) // output wire [63 : 0] P
);
end else begin
reg [AB_WIDTH-1:0] z2_[DSP_STEP-1:0];
reg [AB_WIDTH-1:0] z0_[DSP_STEP-1:0];
reg [AB_WIDTH+1:0] z1_[DSP_STEP-1:0];
always @ (posedge clk)begin
(*use_dsp48="yes"*) z2_[0]<=a[AB_WIDTH-1:AB_WIDTH/2]*b[AB_WIDTH-1:AB_WIDTH/2];
(*use_dsp48="yes"*) z0_[0]<=a[(AB_WIDTH/2) -1:0]*b[(AB_WIDTH/2) -1:0];
(*use_dsp48="yes"*) z1_[0]<=z1_tmp3*z1_tmp4;
if(DSP_STEP>1) begin
z2_[DSP_STEP-1:1]<=z2_[DSP_STEP-2:0];
z0_[DSP_STEP-1:1]<=z0_[DSP_STEP-2:0];
z1_[DSP_STEP-1:1]<=z1_[DSP_STEP-2:0];
end
end
assign z2=z2_[DSP_STEP-1]; // delay 6 clk
assign z1=z1_[DSP_STEP-1];
assign z0=z0_[DSP_STEP-1];
end
endgenerate
// delay 1 clk
always @ (posedge clk)begin
z1_tmp3 <= a[AB_WIDTH-1:AB_WIDTH/2] + a[(AB_WIDTH/2) -1:0]; // delay 1 clk
z1_tmp4 <= b[AB_WIDTH-1:AB_WIDTH/2] + b[(AB_WIDTH/2) -1:0];
end
// delay 2 clk
always @ (posedge clk)begin
z2_z0_tmp[0] <= z2+z0; // delay 7 clk
tmp1 <= {z2,z0};
tmp1_<={(tmp1[RESULT_WIDTH-1 :AB_WIDTH/2]+z1),tmp1[(AB_WIDTH/2) -1:0]}; // delay 8 clk
z2_z0_tmp[1]<=z2_z0_tmp[0];
end
// delay 4 clk
always @ (posedge clk)begin
result <= {tmp1_[RESULT_WIDTH :AB_WIDTH/2] - z2_z0_tmp[1],tmp1_[(AB_WIDTH/2) -1:0]}; // delay 9 clk
end
endmodule五、CSD编码乘法器原理
CSD编码（Canonic Signed Digit）是一种用于表示和计算数字的编码方法，特别适用于乘法器的设计。CSD编码将每个数字表示为一系列正负1的权重值之和，其中正数表示1，负数表示-1。
乘法器的设计可以采用CSD编码来实现高效的乘法运算。具体步骤如下：
- 将输入的乘数和被乘数转换为CSD编码。这涉及到将二进制表示的数字转换为CSD形式。
- 实现部分乘法器来计算每个位上的部分乘积。这可以通过位与运算和位移运算来实现。
- 对于每个部分乘积，根据CSD编码的规则，将其转换为CSD形式。
- 将所有部分乘积相加，得到最终结果。在这一步中，可以使用加法器和位移运算来实现。
通过使用CSD编码，可以减少乘法器中的部分乘积和加法器的数量，从而提高乘法器的效率。然而，CSD编码的实现相对复杂，需要额外的逻辑电路。因此，在实际中，需要综合考虑CSD编码的优点和实现的复杂性来选择乘法器的设计方法。
5.1、CSD编码乘法器的verilogHDL代码例子
此代码是Montgomery的乘法器的第二个乘法即：计算中间值Mb=Ma*Q;
//
// Company: 杭州幻智星科技有限公司
// Engineer: 俞工
// email: yupc123@qq.com
// Create Date: 01/25/2021 03:45:32 PM
// Design Name: CSD编码256位乘法器
// Module Name: y_CSD_mul2_256
// Revision 0.01
// Description: 计算中间值Mb=Ma*Q
// Q=256'h3d443ab0d7bf2839181b2c170004ec0653ba5bfffffe5bfdfffffffeffffffff
//
module y_CSD_mul2_256(
input clk,
input in_val,
input[255:0] a,
output reg out_val=0,
output reg[255:0] result=0 //delay 4 clk
);
wire[255:0] SL_254;
wire[255:0] SL_247;
wire[255:0] SL_245;
wire[255:0] SL_241;
wire[255:0] SL_237;
wire[255:0] SL_223; //30
wire[255:0] SL_204;
wire[255:0] SL_202;
wire[255:0] SL_197;
wire[255:0] SL_191;
wire[255:0] SL_187;
wire[255:0] SL_186;
wire[255:0] SL_180;
wire[255:0] SL_172;
wire[255:0] SL_170;
wire[255:0] SL_169;
wire[255:0] SL_163;
wire[255:0] SL_162;
wire[255:0] SL_145;
wire[255:0] SL_143;
wire[255:0] SL_129;
wire[255:0] SL_128;
wire[255:0] SL_125;
wire[255:0] SL_123;
wire[255:0] SL_121;
wire[255:0] SL_112;
wire[255:0] SL_110;
wire[255:0] SL_78;
//-1 383
wire[255:0] SL_249;
wire[255:0] SL_233;
wire[255:0] SL_231;
wire[255:0] SL_229;
wire[255:0] SL_227;
wire[255:0] SL_220;
wire[255:0] SL_218;
wire[255:0] SL_213;
wire[255:0] SL_207;
wire[255:0] SL_194;
wire[255:0] SL_177;
wire[255:0] SL_175;
wire[255:0] SL_159;
wire[255:0] SL_139;
wire[255:0] SL_137;
wire[255:0] SL_117;
wire[255:0] SL_114;
wire[255:0] SL_108;
wire[255:0] SL_105;
wire[255:0] SL_80;
wire[255:0] SL_76;
wire[255:0] SL_73;
wire[255:0] SL_64;
wire[255:0] SL_31;
wire[255:0] SL_0;
reg [255:0] r_tmp1,r_tmp2,r_tmp3,r_tmp4,r_tmp5,r_tmp6,r_tmp7,r_tmp8,r_tmp9;
reg [255:0] r1_tmp1,r1_tmp2,r1_tmp3,r1_tmp4,r1_tmp5,r1_tmp6,r1_tmp7,r1_tmp8,r1_tmp9;
reg [2:0] shift_csdval=0;
//1
//assign SL_384 = a<<384;
assign SL_254= a<<254;
assign SL_247= a<<248;
assign SL_245= a<<246;
assign SL_241= a<<242;
assign SL_237= a<<238;
assign SL_223= a<<224;
assign SL_204= a<<205;
assign SL_202= a<<203;
assign SL_197= a<<198;
assign SL_191= a<<192;
assign SL_187= a<<188;
assign SL_186= a<<187;
assign SL_180= a<<181;
assign SL_172= a<<173;
assign SL_170= a<<171;
assign SL_169= a<<170;
assign SL_163= a<<164;
assign SL_162= a<<163;
assign SL_145= a<<146;
assign SL_143= a<<144;
assign SL_129= a<<130;
assign SL_128= a<<129;
assign SL_125= a<<126;
assign SL_123= a<<124;
assign SL_121= a<<122;
assign SL_112= a<<113;
assign SL_110= a<<111;
assign SL_78 = a<<79;
//-1 = a<<
assign SL_249 = a<<250;
assign SL_233 = a<<234;
assign SL_231 = a<<232;
assign SL_229 = a<<230;
assign SL_227 = a<<228;
assign SL_220 = a<<221;
assign SL_218 = a<<219;
assign SL_213 = a<<214;
assign SL_207 = a<<208;
assign SL_194 = a<<195;
assign SL_177 = a<<178;
assign SL_175 = a<<176;
assign SL_159 = a<<160;
assign SL_139 = a<<140;
assign SL_137 = a<<138;
assign SL_117 = a<<118;
assign SL_114 = a<<115;
assign SL_108 = a<<109;
assign SL_105 = a<<106;
assign SL_80 = a<<81;
assign SL_76 = a<<77;
assign SL_73 = a<<74;
assign SL_64 = a<<65;
assign SL_31 = a<<32;
assign SL_0 = a;
always @ (posedge clk)begin
//add tree level-1
//r_tmp0 <= SL_384 + SL_380 + SL_367;
r_tmp1 <= SL_125 + SL_163 + SL_187 + SL_237 + SL_241;
r_tmp2 <= SL_123 + SL_162 + SL_186 + SL_223 + SL_245;
r_tmp3 <= SL_121 + SL_145 + SL_180 + SL_204 + SL_247;
r_tmp4 <= SL_112 + SL_143 + SL_172 + SL_202 + SL_254;
r_tmp5 <= SL_110 + SL_129 + SL_170 + SL_197;
r_tmp6 <= SL_78 + SL_128 + SL_169 + SL_191;
//r1_tmp0 <= SL_382 + SL_376 + SL_374;
r1_tmp1 <= SL_80 + SL_139 + SL_213 + SL_233 + SL_249;
r1_tmp2 <= SL_76 + SL_137 + SL_207 + SL_231;
r1_tmp3 <= SL_73 + SL_117 + SL_194 + SL_229;
r1_tmp4 <= SL_64 + SL_114 + SL_177 + SL_227;
r1_tmp5 <= SL_31 + SL_108 + SL_175 + SL_220;
r1_tmp6 <= SL_0 + SL_105 + SL_159 + SL_218;
//add tree level-2
r_tmp7 <= r_tmp1 + r_tmp2 + r_tmp3;
r_tmp8 <= r_tmp4 + r_tmp5 + r_tmp6;
r1_tmp7 <= r1_tmp1 + r1_tmp3 + r1_tmp5;
r1_tmp8 <= r1_tmp2 + r1_tmp4 + r1_tmp6;
//add tree level-3
r_tmp9<=r_tmp7+r_tmp8;
r1_tmp9<=r1_tmp7+r1_tmp8;
//add tree level-4
result<=(r_tmp9)-(r1_tmp9);
end
always @ (posedge clk)begin
{out_val,shift_csdval}<={shift_csdval,in_val};
end
endmoduleMontgomery的乘法器的第三步计算中间值Mc=Mb*P，也适合用CSD编码实现，但常数模P中连续的1太少，使得用CSD编码实现的话，消耗过多的逻辑资源，最终得不偿失，最终放弃。
六、 karatsuba算法同CSD编码其同实现256位乘法
上一节所述的 Montgomery的乘法器的第三步计算中间值Mc=Mb*P不适合单独用CSD编码实现，但是可以用karatsuba算法结合CSD编码一起实现，即达到了不消耗过多的逻辑资源，也能减少DSP48的资源。
6.1、karatsuba算法的三级嵌套
6.2、karatsuba算法和CSD编码混合编程代码
//
// Company: 杭州幻智星科技有限公司
// Engineer: 俞工
// email: yupc123@qq.com
// Create Date: 01/25/2021 03:45:32 PM
// Design Name: Karatsuba算法同CSD编码一起实现的64位乘法器
// Module Name: y_mult_c64
// Revision 0.01
// Description: 计算中间值Mc=Mb*P
//
module y_mult_c64#(
parameter AB_WIDTH = 96,
parameter DSP_STEP = 6,
parameter DSP_EN = 0,
parameter RESULT_WIDTH = AB_WIDTH*2,
parameter [AB_WIDTH-1:0] CONST_B = -1,
parameter [AB_WIDTH/2:0] z1_tmp4 = CONST_B[AB_WIDTH-1:AB_WIDTH/2] + CONST_B[(AB_WIDTH/2) -1:0]
)(
input clk,
input [AB_WIDTH-1:0] a,
output reg [RESULT_WIDTH-1:0] result //delay 4clk
);
reg [RESULT_WIDTH-1:0] result_delay[DSP_STEP-2:0];
generate
if(AB_WIDTH==64&&CONST_B==64'b1111111111111111111111111111111100000000000000000000000000000001)begin
reg [AB_WIDTH-1:0] tmp1,tmp2;
reg [RESULT_WIDTH-1:0] result_;
always @ (posedge clk)begin
tmp1<=a;
tmp2<=tmp1;
result_<={tmp2,tmp2}-{tmp2,32'd0};
if(DSP_STEP==1) begin
result<=result_;
end else begin
if(DSP_STEP>2) begin
result_delay[DSP_STEP-2:0]<={result_delay[DSP_STEP-3:0],result_};
end
result<=result_delay[DSP_STEP-2];
end
end
end else if(DSP_EN==0&&AB_WIDTH==64&&CONST_B==64'b0101001110111101101001000000001011111111111111100101101111111110)begin
reg [AB_WIDTH+48:0] tmp1;
reg [AB_WIDTH+44:0] tmp2;
reg [AB_WIDTH+45:0] tmp3;
reg [AB_WIDTH+20:0] tmp4;
reg [AB_WIDTH-1:0] tmp;
reg [RESULT_WIDTH-1:0] result_;
always @ (posedge clk)begin
tmp<=a;
tmp1<={tmp}+{tmp,12'd0}+{tmp,31'd0}+{tmp,48'd0};
tmp2<={tmp}+{tmp,7'd0}+{tmp,37'd0}+{tmp,44'd0};
tmp3<={tmp}+{tmp,19'd0}+{tmp,30'd0}+{tmp,45'd0};
tmp4<={tmp}+{tmp,16'd0}+{tmp,20'd0};
result_<=({tmp4,42'd0}+{tmp3,15'd0})-({tmp1,1'd0}+{tmp2,10'd0});
if(DSP_STEP==1) begin
result<=result_;
end else begin
if(DSP_STEP>2) begin
result_delay[DSP_STEP-2:0]<={result_delay[DSP_STEP-3:0],result_};
end
result<=result_delay[DSP_STEP-2];
end
end
end else if(DSP_EN==0&&AB_WIDTH==64&&CONST_B==64'b1010100111011110110100100000000011111111111111110010110111111111)begin
// 1010101000h0000h0h01001000000001000000000000000h010h00h00000000h
reg [AB_WIDTH+48:0] tmp1;
reg [AB_WIDTH+53:0] tmp2;
reg [AB_WIDTH+61:0] tmp3;
reg [AB_WIDTH+63:0] tmp4;
reg [AB_WIDTH-1:0] tmp;
reg [RESULT_WIDTH-1:0] result_;
always @ (posedge clk)begin
tmp<=a;
tmp1<={tmp}+{tmp,9'd0}+{tmp,16'd0}+{tmp,48'd0};
tmp2<={tmp,12'd0}+{tmp,46'd0}+{tmp,53'd0};
tmp3<={tmp,14'd0}+{tmp,41'd0}+{tmp,57'd0}+{tmp,61'd0};
tmp4<={tmp,32'd0}+{tmp,44'd0}+{tmp,59'd0}+{tmp,63'd0};
result_<=(tmp4+tmp3)-(tmp1+tmp2);
if(DSP_STEP==1) begin
result<=result_;
end else begin
if(DSP_STEP>2) begin
result_delay[DSP_STEP-2:0]<={result_delay[DSP_STEP-3:0],result_};
end
result<=result_delay[DSP_STEP-2];
end
end
end else
begin
wire [AB_WIDTH-1:0] z2;
wire [AB_WIDTH-1:0] z0;
wire [AB_WIDTH+1:0] z1;
reg [AB_WIDTH+1:0] z1_tmp;
reg [RESULT_WIDTH-1:0] tmp1;
reg [RESULT_WIDTH:0] tmp1_;
reg [AB_WIDTH/2:0] z1_tmp3;
reg [AB_WIDTH+1:0] z2_z0_tmp[1:0];
if(AB_WIDTH==64)begin
mult32 mult32_u0 (
.CLK(clk), // input wire CLK
.A(a[AB_WIDTH-1:AB_WIDTH/2]), // input wire [31 : 0] A
.B(CONST_B[AB_WIDTH-1:AB_WIDTH/2]), // input wire [31 : 0] B
.P(z2) // output wire [63 : 0] P
);
mult32 mult32_u1 (
.CLK(clk), // input wire CLK
.A(a[(AB_WIDTH/2) -1:0]), // input wire [31 : 0] A
.B(CONST_B[(AB_WIDTH/2) -1:0]), // input wire [31 : 0] B
.P(z0) // output wire [63 : 0] P
);
mult33 mult33_u0 (
.CLK(clk), // input wire CLK
.A(z1_tmp3), // input wire [31 : 0] A
.B(z1_tmp4), // input wire [31 : 0] B
.P(z1) // output wire [63 : 0] P
);
end else begin
reg [AB_WIDTH-1:0] z2_[DSP_STEP-1:0];
reg [AB_WIDTH-1:0] z0_[DSP_STEP-1:0];
reg [AB_WIDTH+1:0] z1_[DSP_STEP-1:0];
always @ (posedge clk)begin
(*use_dsp48="yes"*) z2_[0]<=a[AB_WIDTH-1:AB_WIDTH/2]*CONST_B[AB_WIDTH-1:AB_WIDTH/2];
(*use_dsp48="yes"*) z0_[0]<=a[(AB_WIDTH/2) -1:0]*CONST_B[(AB_WIDTH/2) -1:0];
(*use_dsp48="yes"*) z1_[0]<=z1_tmp3*z1_tmp4;
if(DSP_STEP>1) begin
z2_[DSP_STEP-1:1]<=z2_[DSP_STEP-2:0];
z0_[DSP_STEP-1:1]<=z0_[DSP_STEP-2:0];
z1_[DSP_STEP-1:1]<=z1_[DSP_STEP-2:0];
end
end
assign z2=z2_[DSP_STEP-1]; // delay 6 clk
assign z1=z1_[DSP_STEP-1];
assign z0=z0_[DSP_STEP-1];
end
always @ (posedge clk)begin
z1_tmp3 <= a[AB_WIDTH-1:AB_WIDTH/2] + a[(AB_WIDTH/2) -1:0];
end
// delay 2 clk
always @ (posedge clk)begin
z2_z0_tmp[0] <= z2+z0;
tmp1 <= {z2,z0};
tmp1_<={(tmp1[RESULT_WIDTH-1 :AB_WIDTH/2]+z1),tmp1[(AB_WIDTH/2) -1:0]};
z2_z0_tmp[1]<=z2_z0_tmp[0];// delay 3 clk
end
// delay 3 clk
always @ (posedge clk)begin
result <= {tmp1_[RESULT_WIDTH :AB_WIDTH/2] - z2_z0_tmp[1],tmp1_[(AB_WIDTH/2) -1:0]}; // delay 4 clk
end
end
endgenerate
endmodule七、Montgomery模乘算法的设计与实现
下图为Montgomery模乘算法层级结构，y_mul_mont2、y_CSD_mul2_256、y_mult_c256是三个256位的大整数乘法器，用上面所述的三种方法实现。
7.1、Montgomery模乘算法顶层代码
全流水算法设计实现，流水级数为40级。
//
// Company: 杭州幻智星科技有限公司
// Engineer: 俞工
// email: yupc123@qq.com
// Create Date: 01/25/2021 03:45:32 PM
// Design Name: Montgomery算法实现的256位乘法器top层代码
// Module Name: y_mul_mont
// Revision 0.01
//
module y_mul_mont#(
parameter P =256'h73eda753299d7d483339d80809a1d80553bda402fffe5bfeffffffff00000001,
//parameter P = 256'd21888242871839275222246405745257275088696311157297823662689037894645226208583,
parameter INV = bn128_pkg::MONT_FACTOR,//256'd2613481120475770576337551280382610583496597584173309950404576975523838984063,
parameter CTL_WIDTH = 1,
parameter DSP_EN=0,
parameter DSP_STEP = 6,
parameter DAT_BITS = 256
)
(
input clk,
input in_val,
input wire [CTL_WIDTH-1:0] in_ctl,
input [DAT_BITS-1:0] a,
input [DAT_BITS-1:0] b,
output reg out_val=0,
output wire [CTL_WIDTH-1:0] out_ctl,
output reg [DAT_BITS-1:0] result=0 //delay DSP_STEP*2+28 clk
);
wire [DAT_BITS*2-1:0] m;
reg [DAT_BITS-1:0] r_a=0;
reg [DAT_BITS-1:0] r_b=0;
reg in_valid=0;
wire m_val;
logic [CTL_WIDTH-1:0] shift_ctl[DSP_STEP*2+28:0];
assign out_ctl=shift_ctl[DSP_STEP*2+28];
always @ (posedge clk)begin
r_a<=a;
r_b<=b;
end
initial begin
for(int i=0;i<(DSP_STEP*2+29);i=i+1)
shift_ctl[i]=0;
end
always @ (posedge clk)begin
in_valid<=in_val;
shift_ctl<={shift_ctl[DSP_STEP*2+27:0],in_ctl};
end
y_mult_256 #(
.DSP_STEP(DSP_STEP),
.AB_WIDTH(DAT_BITS)
) U0_m(
.clk(clk),
.in_val(in_valid),
.a(r_a),
.b(r_b),
.out_val(m_val),
.result(m) //delay 16 clk
);
wire [DAT_BITS-1:0] m_f;
wire m_f_val;
y_CSD_mul2_256 U1_m_f(
.clk(clk),
.in_val(m_val),
.a(m[DAT_BITS-1:0]),
.out_val(m_f_val),
.result(m_f) //delay 20 clk
);
wire[DAT_BITS*2-1:0] m_f_P;
wire m_f_p_val;
y_mult_c256 #( .DSP_EN(DSP_EN),
.DSP_STEP(DSP_STEP),
.AB_WIDTH(DAT_BITS),
.CONST_B(P)
)U2_m_f_P(
.clk(clk),
.in_val(m_f_val),
.a(m_f[DAT_BITS-1:0]),
.out_val(m_f_p_val),
.result(m_f_P) //delay 36 clk
);
wire[DAT_BITS*2-1:0] result_tmp1;
fifo_shift512 fifo_shift_u (
.clk(clk), // input wire clk
.din(m), // input wire [767 : 0] din
.wr_en(m_val), // input wire wr_en
.rd_en(m_f_p_val), // input wire rd_en
.dout(result_tmp1), // output wire [767 : 0] dout
.full(), // output wire full
.empty() // output wire empty
);
//assign w_result_tmp=(m_f_P+result_tmp1); //513=512+512
reg [256:0]add_result_low=0;
reg [256:0]add_result_high=0;
always@(posedge clk)begin
add_result_low<=m_f_P[255:0]+result_tmp1[255:0]; //delay 37 clk
add_result_high<=m_f_P[511:256]+result_tmp1[511:256];
end
reg [256:0]add_result_high_full=0;
always@(posedge clk)add_result_high_full<=add_result_high+add_result_low[256]; //delay 38 clk
reg[DAT_BITS:0] result_tmp=0;
always @ (posedge clk) result_tmp<=add_result_high_full; //delay 39 clk
always @ (posedge clk)begin
if(result_tmp>=P)begin
result <= result_tmp - P; //delay 40 clk
end
else begin
result <= result_tmp;
end
end
reg [2:0] valid_out_reg=0;
always@(posedge clk) begin
{out_val,valid_out_reg} <= {valid_out_reg,m_f_p_val};
end
endmodule八、在区块链filecoin P2阶段中fpga加速中的应用
上述的Montgomery模乘算法算法用在filecon矿机的 P2阶段，用Xilinx Alveo U200加速卡实现。本项目难点在于如何简化降低单个256位模乘算子的DSP资源使用量，如何将算法逻辑布局约束在不同区域，优化矩阵运算并减少模乘算子的使用量等，最终经过6个月的不断努力，把单个模乘算子的DSP使用量减少到180个，模乘算子使用量减少到25个。在Alveo U200上主要资源使用情况如下。
1、逻辑资源使用率:70%;
2、DSP使用率：91%；
3、运行频率：324M;
工程布线综合结果如下图：
九、bln128模乘源代码链接
下载地址：https://download.csdn.net/download/weixin_46897017/89584817
4554
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
