Bignum: Montgomery multiplication from bignum prototype

library/bignum.c

@@ -38,7 +38,6 @@
 #if defined(MBEDTLS_BIGNUM_C)
 
 #include "mbedtls/bignum.h"
-#include "bignum_internal.h"
 #include "bignum_core.h"
 #include "bn_mul.h"
 #include "mbedtls/platform_util.h"
@@ -962,40 +961,6 @@ int mbedtls_mpi_add_abs( mbedtls_mpi *X, const mbedtls_mpi *A, const mbedtls_mpi
     return( ret );
 }
 
-/**
- * Helper for mbedtls_mpi subtraction.
- *
- * Calculate l - r where l and r have the same size.
- * This function operates modulo (2^ciL)^n and returns the carry
- * (1 if there was a wraparound, i.e. if `l < r`, and 0 otherwise).
- *
- * d may be aliased to l or r.
- *
- * \param n             Number of limbs of \p d, \p l and \p r.
- * \param[out] d        The result of the subtraction.
- * \param[in] l         The left operand.
- * \param[in] r         The right operand.
- *
- * \return              1 if `l < r`.
- *                      0 if `l >= r`.
- */
-static mbedtls_mpi_uint mpi_sub_hlp( size_t n,
-                                     mbedtls_mpi_uint *d,
-                                     const mbedtls_mpi_uint *l,
-                                     const mbedtls_mpi_uint *r )
-{
-    size_t i;
-    mbedtls_mpi_uint c = 0, t, z;
-
-    for( i = 0; i < n; i++ )
-    {
-        z = ( l[i] <  c );    t = l[i] - c;
-        c = ( t < r[i] ) + z; d[i] = t - r[i];
-    }
-
-    return( c );
-}
-
 /*
  * Unsigned subtraction: X = |A| - |B|  (HAC 14.9, 14.10)
  */
@@ -1028,7 +993,7 @@ int mbedtls_mpi_sub_abs( mbedtls_mpi *X, const mbedtls_mpi *A, const mbedtls_mpi
     if( X->n > A->n )
         memset( X->p + A->n, 0, ( X->n - A->n ) * ciL );
 
-    carry = mpi_sub_hlp( n, X->p, A->p, B->p );
+    carry = mbedtls_mpi_core_sub( X->p, A->p, B->p, n );
     if( carry != 0 )
     {
         /* Propagate the carry to the first nonzero limb of X. */
@@ -1157,38 +1122,6 @@ int mbedtls_mpi_sub_int( mbedtls_mpi *X, const mbedtls_mpi *A, mbedtls_mpi_sint
     return( mbedtls_mpi_sub_mpi( X, A, &B ) );
 }
 
-mbedtls_mpi_uint mbedtls_mpi_core_mla( mbedtls_mpi_uint *d, size_t d_len,
-                                       const mbedtls_mpi_uint *s, size_t s_len,
-                                       mbedtls_mpi_uint b )
-{
-    mbedtls_mpi_uint c = 0; /* carry */
-    size_t excess_len = d_len - s_len;
-
-    size_t steps_x8 = s_len / 8;
-    size_t steps_x1 = s_len & 7;
-
-    while( steps_x8-- )
-    {
-        MULADDC_X8_INIT
-        MULADDC_X8_CORE
-        MULADDC_X8_STOP
-    }
-
-    while( steps_x1-- )
-    {
-        MULADDC_X1_INIT
-        MULADDC_X1_CORE
-        MULADDC_X1_STOP
-    }
-
-    while( excess_len-- )
-    {
-        *d += c; c = ( *d < c ); d++;
-    }
-
-    return( c );
-}
-
 /*
  * Baseline multiplication: X = A * B  (HAC 14.12)
  */
@@ -1612,21 +1545,9 @@ int mbedtls_mpi_mod_int( mbedtls_mpi_uint *r, const mbedtls_mpi *A, mbedtls_mpi_
     return( 0 );
 }
 
-/*
- * Fast Montgomery initialization (thanks to Tom St Denis)
- */
 static void mpi_montg_init( mbedtls_mpi_uint *mm, const mbedtls_mpi *N )
 {
-    mbedtls_mpi_uint x, m0 = N->p[0];
-    unsigned int i;
-
-    x  = m0;
-    x += ( ( m0 + 2 ) & 4 ) << 1;
-
-    for( i = biL; i >= 8; i /= 2 )
-        x *= ( 2 - ( m0 * x ) );
-
-    *mm = ~x + 1;
+    *mm = mbedtls_mpi_montg_init( N->p );
 }
 
 /** Montgomery multiplication: A = A * B * R^-1 mod N  (HAC 14.36)
@@ -1640,67 +1561,20 @@ static void mpi_montg_init( mbedtls_mpi_uint *mm, const mbedtls_mpi *N )
  * \param[in]       B   One of the numbers to multiply.
  *                      It must be nonzero and must not have more limbs than N
  *                      (B->n <= N->n).
- * \param[in]       N   The modulo. N must be odd.
+ * \param[in]       N   The modulus. \p N must be odd.
  * \param           mm  The value calculated by `mpi_montg_init(&mm, N)`.
  *                      This is -N^-1 mod 2^ciL.
  * \param[in,out]   T   A bignum for temporary storage.
  *                      It must be at least twice the limb size of N plus 1
  *                      (T->n >= 2 * N->n + 1).
  *                      Its initial content is unused and
  *                      its final content is indeterminate.
- *                      Note that unlike the usual convention in the library
- *                      for `const mbedtls_mpi*`, the content of T can change.
  */
-static void mpi_montmul( mbedtls_mpi *A, const mbedtls_mpi *B, const mbedtls_mpi *N, mbedtls_mpi_uint mm,
-                         const mbedtls_mpi *T )
+static void mpi_montmul( mbedtls_mpi *A, const mbedtls_mpi *B,
+                         const mbedtls_mpi *N, mbedtls_mpi_uint mm,
+                         mbedtls_mpi *T )
 {
-    size_t n, m;
-    mbedtls_mpi_uint *d;
-
-    memset( T->p, 0, T->n * ciL );
-
-    d = T->p;
-    n = N->n;
-    m = ( B->n < n ) ? B->n : n;
-
-    for( size_t i = 0; i < n; i++ )
-    {
-        mbedtls_mpi_uint u0, u1;
-
-        /*
-         * T = (T + u0*B + u1*N) / 2^biL
-         */
-        u0 = A->p[i];
-        u1 = ( d[0] + u0 * B->p[0] ) * mm;
-
-        (void) mbedtls_mpi_core_mla( d, n + 2,
-                                     B->p, m,
-                                     u0 );
-        (void) mbedtls_mpi_core_mla( d, n + 2,
-                                     N->p, n,
-                                     u1 );
-        d++;
-    }
-
-    /* At this point, d is either the desired result or the desired result
-     * plus N. We now potentially subtract N, avoiding leaking whether the
-     * subtraction is performed through side channels. */
-
-    /* Copy the n least significant limbs of d to A, so that
-     * A = d if d < N (recall that N has n limbs). */
-    memcpy( A->p, d, n * ciL );
-    /* If d >= N then we want to set A to d - N. To prevent timing attacks,
-     * do the calculation without using conditional tests. */
-    /* Set d to d0 + (2^biL)^n - N where d0 is the current value of d. */
-    d[n] += 1;
-    d[n] -= mpi_sub_hlp( n, d, d, N->p );
-    /* If d0 < N then d < (2^biL)^n
-     * so d[n] == 0 and we want to keep A as it is.
-     * If d0 >= N then d >= (2^biL)^n, and d <= (2^biL)^n + N < 2 * (2^biL)^n
-     * so d[n] == 1 and we want to set A to the result of the subtraction
-     * which is d - (2^biL)^n, i.e. the n least significant limbs of d.
-     * This exactly corresponds to a conditional assignment. */
-    mbedtls_ct_mpi_uint_cond_assign( n, A->p, d, (unsigned char) d[n] );
+    mbedtls_mpi_core_montmul( A->p, A->p, B->p, B->n, N->p, N->n, mm, T->p );
 }
 
 /*
@@ -1709,7 +1583,7 @@ static void mpi_montmul( mbedtls_mpi *A, const mbedtls_mpi *B, const mbedtls_mpi
  * See mpi_montmul() regarding constraints and guarantees on the parameters.
  */
 static void mpi_montred( mbedtls_mpi *A, const mbedtls_mpi *N,
-                         mbedtls_mpi_uint mm, const mbedtls_mpi *T )
+                         mbedtls_mpi_uint mm, mbedtls_mpi *T )
 {
     mbedtls_mpi_uint z = 1;
     mbedtls_mpi U;

library/bignum_core.c

@@ -37,6 +37,8 @@
 #endif
 
 #include "bignum_core.h"
+#include "bn_mul.h"
+#include "constant_time_internal.h"
 
 size_t mbedtls_mpi_core_clz( mbedtls_mpi_uint a )
 {
@@ -291,4 +293,139 @@ int mbedtls_mpi_core_write_be( const mbedtls_mpi_uint *X,
     return( 0 );
 }
 
+void mbedtls_mpi_core_montmul( mbedtls_mpi_uint *X,
+                               const mbedtls_mpi_uint *A,
+                               const mbedtls_mpi_uint *B,
+                               size_t B_limbs,
+                               const mbedtls_mpi_uint *N,
+                               size_t AN_limbs,
+                               mbedtls_mpi_uint mm,
+                               mbedtls_mpi_uint *T )
+{
+    memset( T, 0, ( 2 * AN_limbs + 1 ) * ciL );
+
+    for( size_t i = 0; i < AN_limbs; i++, T++ )
+    {
+        mbedtls_mpi_uint u0, u1;
+        /* T = (T + u0*B + u1*N) / 2^biL */
+        u0 = A[i];
+        u1 = ( T[0] + u0 * B[0] ) * mm;
+
+        (void) mbedtls_mpi_core_mla( T, AN_limbs + 2, B, B_limbs, u0 );
+        (void) mbedtls_mpi_core_mla( T, AN_limbs + 2, N, AN_limbs, u1 );
+    }
+
+    /* It's possible that the result in T is > N, and so we might need to subtract N */
+
+    mbedtls_mpi_uint carry  = T[AN_limbs];
+    mbedtls_mpi_uint borrow = mbedtls_mpi_core_sub( X, T, N, AN_limbs );
+
+    /*
+     * Both carry and borrow can only be 0 or 1.
+     *
+     * If carry = 1, the result in T must be > N by definition, and the subtraction
+     * using only AN_limbs limbs will create borrow, but that will have the correct
+     * final result.
+     *
+     * i.e. (carry, borrow) of (1, 1) => return X
+     *
+     * If carry = 0, then we want to use the result of the subtraction iff
+     * borrow = 0.
+     *
+     * i.e. (carry, borrow) of (0, 0) => return X
+     *                         (0, 1) => return T
+     *
+     * We've confirmed that the unit tests exercise this function with all 3 of
+     * the valid (carry, borrow) combinations (listed above), and that we don't
+     * see (carry, borrow) = (1, 0).
+     *
+     * So the correct return value is already in X if (carry ^ borrow) = 0,
+     * but is in (the lower AN_limbs limbs of) T if (carry ^ borrow) = 1.
+     */
+     mbedtls_ct_mpi_uint_cond_assign( AN_limbs, X, T, (unsigned char) ( carry ^ borrow ) );
+}
+
+/*
+ * Fast Montgomery initialization (thanks to Tom St Denis).
+ */
+mbedtls_mpi_uint mbedtls_mpi_montg_init( const mbedtls_mpi_uint *N )
+{
+    mbedtls_mpi_uint x = N[0];
+
+    x += ( ( N[0] + 2 ) & 4 ) << 1;
+
+    for( unsigned int i = biL; i >= 8; i /= 2 )
+        x *= ( 2 - ( N[0] * x ) );
+
+    return( ~x + 1 );
+}
+
+mbedtls_mpi_uint mbedtls_mpi_core_mla( mbedtls_mpi_uint *d, size_t d_len,
+                                       const mbedtls_mpi_uint *s, size_t s_len,
+                                       mbedtls_mpi_uint b )
+{
+    mbedtls_mpi_uint c = 0; /* carry */
+    if( d_len < s_len )
+        s_len = d_len;
+    size_t excess_len = d_len - s_len;
+    size_t steps_x8 = s_len / 8;
+    size_t steps_x1 = s_len & 7;
+
+    while( steps_x8-- )
+    {
+        MULADDC_X8_INIT
+        MULADDC_X8_CORE
+        MULADDC_X8_STOP
+    }
+
+    while( steps_x1-- )
+    {
+        MULADDC_X1_INIT
+        MULADDC_X1_CORE
+        MULADDC_X1_STOP
+    }
+
+    while( excess_len-- )
+    {
+        *d += c; c = ( *d < c ); d++;
+    }
+
+    return( c );
+}
+
+mbedtls_mpi_uint mbedtls_mpi_core_sub( mbedtls_mpi_uint *X,
+                                       const mbedtls_mpi_uint *A,
+                                       const mbedtls_mpi_uint *B,
+                                       size_t limbs )
+{
+    mbedtls_mpi_uint c = 0;
+
+    for( size_t i = 0; i < limbs; i++ )
+    {
+        mbedtls_mpi_uint z = ( A[i] < c );
+        mbedtls_mpi_uint t = A[i] - c;
+        c = ( t < B[i] ) + z;
+        X[i] = t - B[i];
+    }
+
+    return( c );
+}
+
+mbedtls_mpi_uint mbedtls_mpi_core_add_if( mbedtls_mpi_uint *A,
+                                          const mbedtls_mpi_uint *B,
+                                          size_t limbs,
+                                          unsigned cond )
+{
+    mbedtls_mpi_uint c = 0, t;
+    for( size_t i = 0; i < limbs; i++ )
+    {
+        mbedtls_mpi_uint add = cond * B[i];
+        t  = c;
+        t += A[i]; c  = ( t < A[i] );
+        t += add;  c += ( t < add  );
+        A[i] = t;
+    }
+    return( c );
+}
+
 #endif /* MBEDTLS_BIGNUM_C */