MbedTLS中的Montgomery算法实现解析（三）

MbedTLS中的Montgomery算法实现解析（三）

在前面的两篇博客中，我们详细讨论了Montgomery算法的数学原理以及MbedTLS中是如何快速求得-N^-1(mod n)。在这一篇博客中，我们将探讨MbedTLS中的蒙哥马利约减（Montgomery reduction）和蒙哥马利模乘（Montgomery multiplication）是如何实现的。

在阅读本博客之前，建议先掌握MbedTLS的大数表示方式

蒙哥马利约减Montgomery reduction

在我们的数学推导中，蒙哥马利约减是蒙哥马利模乘的基础，也就是说我们首先有：

X ≡ a × R^-1 (mod n)

然后才有：

X = a × R × b × R × R^-1 (mod n)

但是MbedTLS将这两者结合了起来，在MbedTLS中定义的蒙哥马利模乘式为：

X = A × B × R^-1 (mod n)

既包含了模乘（Montgomery multiplication），又包含了约减（Montgomery reduction），这样的好处就是约减可以写成模乘的一种特殊情况，即当B=1时：

X = A × B × R^-1 (mod n) = A × R^-1 (mod n)

因此在MbenTLS中，约减就是模乘的一种特殊情况。

下面我们来看代码。

// Montgomery reduction: A = A * R^-1 mod N
static int mpi_montred( mbedtls_mpi *A, const mbedtls_mpi *N, mbedtls_mpi_uint mm, const mbedtls_mpi *T )
{
    mbedtls_mpi_uint z = 1;
    mbedtls_mpi U;

    U.n = U.s = (int) z;
    U.p = &z;

    return( mpi_montmul( A, &U, N, mm, T ) );
}

可以看到这里初始化了一个值为1的大数变量U，随后调用了蒙哥马利模乘函数mpi_montmul

蒙哥马利模乘Montgomery multiplication

下面我们来看看这个蒙哥马利模乘函数mpi_montmul：

// Montgomery multiplication: A = A * B * R^-1 mod N
static int mpi_montmul( mbedtls_mpi *A, const mbedtls_mpi *B, const mbedtls_mpi *N, mbedtls_mpi_uint mm,const mbedtls_mpi *T )
{
    size_t i, n, m;
    mbedtls_mpi_uint u0, u1, *d;

    if( T->n < N->n + 1 || T->p == NULL )
        return( MBEDTLS_ERR_MPI_BAD_INPUT_DATA );

    memset( T->p, 0, T->n * ciL );

    d = T->p;
    n = N->n;
    m = ( B->n < n ) ? B->n : n;

    for( i = 0; i < n; i++ )
    {
        /*
         * T = (T + u0*B + u1*N) / 2^biL
         */
        u0 = A->p[i];
        u1 = ( d[0] + u0 * B->p[0] ) * mm;

        mpi_mul_hlp( m, B->p, d, u0 );
        mpi_mul_hlp( n, N->p, d, u1 );

        *d++ = u0; d[n + 1] = 0;
    }

    memcpy( A->p, d, ( n + 1 ) * ciL );

    if( mbedtls_mpi_cmp_abs( A, N ) >= 0 )
        mpi_sub_hlp( n, N->p, A->p );
    else
        /* prevent timing attacks */
        mpi_sub_hlp( n, A->p, T->p );

    return( 0 );
}

函数有五个形参：
mbedtls_mpi*A ：第一个乘数
const mbedtls_mpi*B：第二个乘数
const mbedtls_mpi*N：模除的大数N
mbedtls_mpi_uint mm：预计算结果-N^-1(mod n)
const mbedtls_mpi*T：用来保存中间计算结果的T

函数没有提到R的取值，在后面的实现中函数采用2^32×n作为R。而在每一轮的运算中，我们针对的是2³²，之所以这样实现是因为MbedTLS中的大整数采用的字长是32位，这样在涉及到模除运算的时候，我们只需要关注大整数的最后一个字即可。

我们一步一步来看代码，首先是：

 if( T->n < N->n + 1 || T->p == NULL )
        return( MBEDTLS_ERR_MPI_BAD_INPUT_DATA );

首先T是用来保存两个大数相乘的中间变量，T的limbs数必须足够大，当不够大或者为空的时候要返回错误。

// 将T的内容全部置为0
memset( T->p, 0, T->n * ciL );

d = T->p; //d指向T->p
n = N->n; //n记录N->n
m = ( B->n < n ) ? B->n : n; //m记录N->n和B->n中较小的一个

接下来是Montgomery multiplication的核心，一个循环：

for( i = 0; i < n; i++ )
    {
        /*
         * T = (T + u0*B + u1*N) / 2^biL
         */
        u0 = A->p[i];
        u1 = ( d[0] + u0 * B->p[0] ) * mm;

        mpi_mul_hlp( m, B->p, d, u0 );
        mpi_mul_hlp( n, N->p, d, u1 );

        *d++ = u0; d[n + 1] = 0;
    }

首先，有一点我们需要明确的是，函数接收的A和B其实是：

A = a × R (mod n)
B = b × R (mod n)

这也就是说明：

0 < A , B < n - 1

这意味着A和B只是低n×Bil上的内容是有意义的，在高于n×Bil位上的内容全部是0.

因此我们在处理的时候只需要考虑A和B的低n×Bil上的内容即可。

在这个循环中，u₀=A->p[i]，经过循环后u₀会遍历A->p中的内容。

u1计算的内容是我们在博客一中提到的：s ≡ - a₀ × n^-1 (mod R)

d[0] + u₀ × B->p[0] = A->p[i] × B->p[0] + T->p[0]

相当于每一轮计算的低32位再加上之前计算结果中的低32位。为什么要只计算低32位呢？这就是蒙哥马利算法很巧妙的地方，因为R的取值是2³²，这说明一个数在模2³²之后只用看他的低32位内容即可，也就是说我们只用关注第一个字即可（MbedTLS中低位放在前面的字中）。

再乘上mm，我们在博客二中提到了mm ≡ -n^-1 (mod R)

于是我们现在就得到了u₁ ≡ - T₀ × n^-1 (mod R)

读者可能会注意到，为什么u₁是两个字相乘的结果，应该占2个字的存储空间，为什么我们只用了一个字来存放？这是因为高位字的内容模R等于0，我们只需要关注低位字就行了。

mpi_mul_hlp( m, B->p, d, u0 );
mpi_mul_hlp( n, N->p, d, u1 );

这里计算了：

d + = u₀ × B + u₁ × N

在这里我们讨论一下d现在有多少个字长？
首先u₀和u₁都是一个字长，而B的长度小于等于n个字，N的长度等于n个字。那么d的长度应该不超过n+1个字。

图解如下：

到这里我们就只差最后一步了，这个时候的d[0]的内容已经全是0（原因可以见博客一，这里就不解释了），下面我们需要进行右移操作，将整个d向右移动一个字，在这里MbedTLS采用了一种非常聪明的办法，并不是右移数字，而是左移指针。

*d++ = u0;

最后d[n+1]=0，是为了防止内存中出现其它的数据干扰计算结果。
因为我们通过图解可以看到d的数据不可能超过d[n]，d[n+1]用来保存的是
d^’[n]的内容，在下一个循环开始之前，我们将这个内容清空，方便保存d^’[n]。

每左移一次指针，就相当于除以2³²，一共左移了n次指针，一共就相当于除以了2^32×n，就相当于除以R。

在这个循环结束后，我们就得到了一个长度为n+1个字的d，此时的
d_终=d[n]。

此时的d[0]到d[n-1]也记录了A->p的n个字的内容。

在这个循环结束之后，我们就已经计算到了：
X = A × B × R^-1 (mod n)

下面我们将这(n+1)个字的值copy给A：

 memcpy( A->p, d, ( n + 1 ) * ciL );

就相当于完成了：

A = A × B × R^-1 (mod n)

最后我们需要检查一下A是否大于N，如果大于N的话就直接 A=A-N就行，在数学上有严格的证明可以证明 A < 2×N 。

if( mbedtls_mpi_cmp_abs( A, N ) >= 0 )
        mpi_sub_hlp( n, N->p, A->p );
    else
        /* prevent timing attacks */
        mpi_sub_hlp( n, A->p, T->p );

最后将T->p中的低n个字（也就是A->p）的内容清除是为了防止时序攻击，与算法的实现原理关系不大。

写在最后

从刚开始研究蒙哥马利到写完这篇博客，已经断断续续地过了一个星期的时间，在这个时间里我慢慢了解到了蒙哥马利算法的精巧，同时在实现的过程中也有很多小技巧需要去琢磨。在师兄们的无私帮助之下，我总算写完了这三篇博客，博客内也可能有一些错误，欢迎大家的批评指正。