RSA 硬件实现 : radix-2 montgomery multiplication issues

Question

我正在硬件(xilinx ZYNQ FPGA)中实现 RSA 1024，但我无法解决一些奇怪的问题。最值得注意的是，我发现我的实现仅适用于某些基数/指数/模数组合，但没有找到任何原因。

注意:我正在使用 Xilinx HLS(本质上是合成到硬件中的 C 代码)实现算法。为了这篇文章，把它当作一个标准的 C 实现来对待，除了我可以有高达 4096 位宽的变量。我还没有将它并行化，所以它应该像标准 C 代码一样运行。

---

问题

我的问题是，我能够得到某些模幂测试问题的正确答案，但前提是底数、指数和模数的值可以用比实际 1024 位少得多的位数编写操作数宽度(即它们被零填充)。

当我使用从 SSH-keygen 生成的实际 1024 位值时，我不再得到正确的结果。

例如，如果我的输入参数是

uint1024_t base     = 1570
uint1024_t exponent = 1019
uint1024_t modulus  = 3337

我正确地得到了 1570^1029 mod(3337) = 688 的结果

但是，当我实际使用占据所有(或几乎所有)1024 位的值作为输入时...

uint1024_t base     = 0x00be5416af9696937b7234421f7256f78dba8001c80a5fdecdb4ed761f2b7f955946ec920399f23ce9627f66286239d3f20e7a46df185946c6c8482e227b9ce172dd518202381706ed0f91b53c5436f233dec27e8cb46c4478f0398d2c254021a7c21596b30f77e9886e2fd2a081cadd3faf83c86bfdd6e9daad12559f8d2747
uint1024_t exponent = 0x6f1e6ab386677cdc86a18f24f42073b328847724fbbd293eee9cdec29ac4dfe953a4256d7e6b9abee426db3b4ddc367a9fcf68ff168a7000d3a7fa8b9d9064ef4f271865045925660fab620fad0aeb58f946e33bdff6968f4c29ac62bd08cf53cb8be2116f2c339465a64fd02517f2bafca72c9f3ca5bbf96b24c1345eb936d1
uint1024_t modulus  = 0xb4d92132b03210f62e52129ae31ef25e03c2dd734a7235efd36bad80c28885f3a9ee1ab626c30072bb3fd9906bf89a259ffd9d5fd75f87a30d75178b9579b257b5dca13ca7546866ad9f2db0072d59335fb128b7295412dd5c43df2c4f2d2f9c1d59d2bb444e6dac1d9cef27190a97aae7030c5c004c5aea3cf99afe89b86d6d

我错误地得到了一个巨大的数字，而不是 29 (0x1D) 的正确答案

我已经对这两种算法进行了一百万次检查，并尝试了不同的初始值和循环边界，但似乎没有任何效果。

---

我的实现

我使用标准的平方和乘法方法进行模幂运算，我选择使用 Tenca-Koc radix-2 算法进行蒙哥马利乘法，下面的伪代码中有详细说明...

/* Tenca-Koc radix2 montgomery multiplication */
Z = 0
for i = 0 to n-1
    Z = Z + X[i]*Y
    if Z is odd then Z = Z + M
    Z = Z/2  // left shift in radix2
if (S >= M) then S = S - M

我的蒙哥马利乘法实现如下:

void montMult(uint1024_t X, uint1024_t Y, uint1024_t M, uint1024_t* outData)
{
    ap_uint<2*NUM_BITS> S = 0; 

    for (int i=0; i<NUM_BITS; i++)
    {
        // add product of X.get_bit(i) and Y to partial sum
        S += X[i]*Y; 

        // if S is even, add modulus to partial sum
        if (S.test(0))  
            S += M;     

        // rightshift 1 bit (divide by 2)
        S = S >> 1;
    }

    // bring back to under 1024 bits by subtracting modulus
    if (S >= M)
        S -= M;

    // write output data
    *outData = S.range(NUM_BITS-1,0); 

我的顶级模幂运算如下，其中(转换符号!)...

// k: number of bits
// r = 2^k (radix)
// M: base
// e: exponent
// n: modulus
// Mbar: (precomputed residue) M*r mod(n)
// xbar: (precomputed initial residue) 1*r mod(n)

void ModExp(uint1024_t M, uint1024_t e, uint1024_t n, 
            uint1024_t Mbar, uint1024_t xbar, uint1024_t* out)
{
    for (int i=NUM_BITS-1; i>=0; i--)
    {
        // square
        montMult(xbar,xbar,n,&xbar);

        // multiply   
        if (e.test(i)) // if (e.bit(i) == 1)
            montMult(Mbar,xbar,n,&xbar);
    }
        // undo montgomery residue transformation
        montMult(xbar,1,n,out);
}

我这辈子都弄不明白为什么这适用于除实际 1024 位值以外的所有内容。任何帮助将不胜感激

Answer 1

我已经替换了我的答案，因为我错了。您的原始代码完全正确。我已经使用我自己的 BigInteger 库对其进行了测试，其中包括蒙哥马利算术，一切都非常有效。这是我的代码:

const
  base1     =
 '0x00be5416af9696937b7234421f7256f78dba8001c80a5fdecdb4ed761f2b7f955946ec9203'+
 '99f23ce9627f66286239d3f20e7a46df185946c6c8482e227b9ce172dd518202381706ed0f91'+
 'b53c5436f233dec27e8cb46c4478f0398d2c254021a7c21596b30f77e9886e2fd2a081cadd3f'+
 'af83c86bfdd6e9daad12559f8d2747';
  exponent1 =
 '0x6f1e6ab386677cdc86a18f24f42073b328847724fbbd293eee9cdec29ac4dfe953a4256d7e'+
 '6b9abee426db3b4ddc367a9fcf68ff168a7000d3a7fa8b9d9064ef4f271865045925660fab62'+
 '0fad0aeb58f946e33bdff6968f4c29ac62bd08cf53cb8be2116f2c339465a64fd02517f2bafc'+
 'a72c9f3ca5bbf96b24c1345eb936d1';
  modulus1  =
 '0xb4d92132b03210f62e52129ae31ef25e03c2dd734a7235efd36bad80c28885f3a9ee1ab626'+
 'c30072bb3fd9906bf89a259ffd9d5fd75f87a30d75178b9579b257b5dca13ca7546866ad9f2d'+
 'b0072d59335fb128b7295412dd5c43df2c4f2d2f9c1d59d2bb444e6dac1d9cef27190a97aae7'+
 '030c5c004c5aea3cf99afe89b86d6d';

function MontMult(X, Y, N: BigInteger): BigInteger;
var
  I: Integer;
begin
  Result:= 0;
  for I:= 0 to 1023 do begin
    if not X.IsEven then Result:= Result + Y;
    if not Result.IsEven then Result:= Result + N;
    Result:= Result shr 1;
    X:= X shr 1;
  end;
  if Result >= N then Result:= Result - N;
end;

function ModExp(B, E, N: BigInteger): BigInteger;
var
  R, MontB: BigInteger;
  I: Integer;

begin
  R:= BigInteger.PowerOfTwo(1024) mod N;
  MontB:= (B * R) mod N;
  for I:= 1023 downto 0 do begin
    R:= MontMult(R, R, N);
    if not (E shr I).IsEven then
      R:= MontMult(MontB, R, N);
  end;
  Result:= MontMult(R, 1, N);
end;

procedure TestMontMult;
var
  Base, Expo, Modulus: BigInteger;
  MontBase, MontExpo: BigInteger;
  X, Y, R: BigInteger;
  Mont: TMont;

begin
// convert to BigInteger
  Base:= BigInteger.Parse(base1);
  Expo:= BigInteger.Parse(exponent1);
  Modulus:= BigInteger.Parse(modulus1);

  R:= BigInteger.PowerOfTwo(1024) mod Modulus;
// Convert into Montgomery form
  MontBase:= (Base * R) mod Modulus;
  MontExpo:= (Expo * R) mod Modulus;
  Writeln;

// MontMult test, all 3 versions output
//  '0x146005377258684F3FFD8D9A70D723BDD3A2E3A160E11B7AD35A7106D4D903AB9D14A9201'+
//  'D0907CE2FC2E04A69656C38CE64AA0BADF2376AEFB19D8732CE2B3650466E31BB78CF24F4E3'+
//  '774A78575738B668DA0E40C8DDDA972CE101E0CADC5D4CCFF6EF2E4E97AF02F34E3AB7258A7'+
//  '323E472FC051825FFC72ADC53B0DAF3C4';
  Writeln('Using MontMult');
  Writeln(MontMult(MontMult(MontBase, MontExpo, Modulus), 1, Modulus).ToHexString);
// same using TMont instance
  Writeln('Using TMont.Multiply');
  Mont:= TMont.GetInstance(Modulus);
  Writeln(Mont.Reduce(Mont.Multiply(MontBase, MontExpo)).ToHexString);
  Writeln('Using TMont.ModMul');
  Writeln(Mont.ModMul(Base,Expo).ToHexString);

// ModExp test, all 3 versions output 29
  Writeln('Using ModExp');
  Writeln(ModExp(Base, Expo, Modulus).ToString);
  Writeln('Using BigInteger.ModPow');
  Writeln(BigInteger.ModPow(Base, Expo, Modulus).ToString);
  Writeln('Using TMont.ModPow');
  Writeln(Mont.ModPow(Base, Expo).ToString);
end;