Implementación de FPGA RSA: problemas de multiplicación de Montgomery

Estoy implementando RSA 1024 en hardware (xilinx ZYNQ FPGA), y no puedo resolver algunos problemas curiosos. En particular, estoy descubriendo que mi implementación solo funciona para ciertas combinaciones de base / exponente / módulo, pero no he encontrado ninguna razón por la que este sea el caso.

Nota: estoy implementando el algoritmo usando Xilinx HLS (esencialmente código C que se sintetiza en hardware). Por el bien de esta publicación, trátela como una implementación estándar de C, excepto que puedo tener variables de hasta 4096 bits de ancho. Todavía no lo he paralizado, por lo que debería comportarse como el código C estándar.

El problema

Mi problema es que puedo obtener la respuesta correcta para ciertos problemas de prueba de exponenciación modular, pero solo si los valores para la base, el exponente y el módulo se pueden escribir en mucho menos bits que los 1024 bits reales ancho del operando (es decir, están rellenados con cero).

Cuando uso valores reales de 1024 bits generados desde SSH-keygen, ya no obtengo los resultados correctos.

Por ejemplo, si mis argumentos de entrada son

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

Obtengo correctamente un resultado de 1570 ^ 1029 mod (3337) = 688

Sin embargo, cuando realmente uso valores que ocupan todos (o aproximadamente todos) 1024 bits para las entradas ...

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

Obtengo incorrectamente un número masivo, en lugar de la respuesta correcta de 29 (0x1D)

He comprobado ambos algoritmos un millón de veces y he experimentado con valores iniciales y límites de bucle diferentes, pero nada parece funcionar.

Mi implementación

Estoy usando el método estándar de cuadrado y multiplicación para la exponenciación modular, y elegí usar el algoritmo Tenca-Koc radix-2 para la multiplicación de Montgomery, que se detalla en el pseudocódigo a continuación ...

/* 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

La implementación de la multiplicación de My Montgomery es la siguiente:

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); 

}

y mi exponenciación modular de nivel superior es la siguiente, donde (cambio de notación) ...

// 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);
}

Por mi vida no puedo entender por qué esto funciona para todo, excepto un valor real de 1024 bits. Cualquier ayuda sería muy apreciada