Merge pull request #1408 from mpg/improve-gcd-3.6

[3.6] Make GCD (a lot) less leaky

ChangeLog.d/gcd-sign.txt

@@ -0,0 +1,5 @@
+Changes
+   * The function mbedtls_mpi_gcd() now always gives a non-negative output.
+     Previously the output was negative when B = 0 and A < 0, which was not
+     documented, and inconsistent as all other inputs resulted in a non-negative
+     output.

include/mbedtls/bignum.h

@@ -974,8 +974,7 @@ int mbedtls_mpi_random(mbedtls_mpi *X,
  * \brief          Compute the greatest common divisor: G = gcd(A, B)
  *
  * \param G        The destination MPI. This must point to an initialized MPI.
- *                 This will be positive unless \p B is 0, in which case \p A
+ *                 This will always be positive or 0.
- *                 will be returned, where \p A could be negative.
  * \param A        The first operand. This must point to an initialized MPI.
  * \param B        The second operand. This must point to an initialized MPI.
  *

library/bignum.c

@@ -430,13 +430,6 @@ cleanup:
     return ret;
 }
 
-/*
- * Return the number of less significant zero-bits
- */
-size_t mbedtls_mpi_lsb(const mbedtls_mpi *X)
-{
-    size_t i;
-
 #if defined(__has_builtin)
 #if (MBEDTLS_MPI_UINT_MAX == UINT_MAX) && __has_builtin(__builtin_ctz)
     #define mbedtls_mpi_uint_ctz __builtin_ctz
@@ -447,22 +440,34 @@ size_t mbedtls_mpi_lsb(const mbedtls_mpi *X)
 #endif
 #endif
 
-#if defined(mbedtls_mpi_uint_ctz)
+#if !defined(mbedtls_mpi_uint_ctz)
+static size_t mbedtls_mpi_uint_ctz(mbedtls_mpi_uint x)
+{
+    size_t count = 0;
+    mbedtls_ct_condition_t done = MBEDTLS_CT_FALSE;
+
+    for (size_t i = 0; i < biL; i++) {
+        mbedtls_ct_condition_t non_zero = mbedtls_ct_bool((x >> i) & 1);
+        done = mbedtls_ct_bool_or(done, non_zero);
+        count = mbedtls_ct_size_if(done, count, i + 1);
+    }
+
+    return count;
+}
+#endif
+
+/*
+ * Return the number of less significant zero-bits
+ */
+size_t mbedtls_mpi_lsb(const mbedtls_mpi *X)
+{
+    size_t i;
+
     for (i = 0; i < X->n; i++) {
         if (X->p[i] != 0) {
             return i * biL + mbedtls_mpi_uint_ctz(X->p[i]);
         }
     }
-#else
-    size_t count = 0;
-    for (i = 0; i < X->n; i++) {
-        for (size_t j = 0; j < biL; j++, count++) {
-            if (((X->p[i] >> j) & 1) != 0) {
-                return count;
-            }
-        }
-    }
-#endif
 
     return 0;
 }
@@ -1814,103 +1819,51 @@ cleanup:
 }
 
 /*
- * Greatest common divisor: G = gcd(A, B)  (HAC 14.54)
+ * Greatest common divisor: G = gcd(A, B)
+ * Wrapper around mbedtls_mpi_gcd_modinv() that removes its restrictions.
  */
 int mbedtls_mpi_gcd(mbedtls_mpi *G, const mbedtls_mpi *A, const mbedtls_mpi *B)
 {
     int ret = MBEDTLS_ERR_ERROR_CORRUPTION_DETECTED;
-    size_t lz, lzt;
     mbedtls_mpi TA, TB;
 
     mbedtls_mpi_init(&TA); mbedtls_mpi_init(&TB);
 
+    /* Make copies and take absolute values */
     MBEDTLS_MPI_CHK(mbedtls_mpi_copy(&TA, A));
     MBEDTLS_MPI_CHK(mbedtls_mpi_copy(&TB, B));
+    TA.s = TB.s = 1;
 
-    lz = mbedtls_mpi_lsb(&TA);
+    /* Make the two values the same (non-zero) number of limbs.
-    lzt = mbedtls_mpi_lsb(&TB);
+     * This is needed to use mbedtls_mpi_core functions below. */
+    MBEDTLS_MPI_CHK(mbedtls_mpi_grow(&TA, TB.n != 0 ? TB.n : 1));
+    MBEDTLS_MPI_CHK(mbedtls_mpi_grow(&TB, TA.n)); // non-zero from above
 
-    /* The loop below gives the correct result when A==0 but not when B==0.
+    /* Handle special cases (that don't happen in crypto usage) */
-     * So have a special case for B==0. Leverage the fact that we just
+    if (mbedtls_mpi_core_check_zero_ct(TA.p, TA.n) == MBEDTLS_CT_FALSE) {
-     * calculated the lsb and lsb(B)==0 iff B is odd or 0 to make the test
+        MBEDTLS_MPI_CHK(mbedtls_mpi_copy(G, &TB)); // GCD(0, B) = abs(B)
-     * slightly more efficient than cmp_int(). */
+        goto cleanup;
-    if (lzt == 0 && mbedtls_mpi_get_bit(&TB, 0) == 0) {
+    }
-        ret = mbedtls_mpi_copy(G, A);
+    if (mbedtls_mpi_core_check_zero_ct(TB.p, TB.n) == MBEDTLS_CT_FALSE) {
+        MBEDTLS_MPI_CHK(mbedtls_mpi_copy(G, &TA)); // GCD(A, 0) = abs(A)
         goto cleanup;
     }
 
-    if (lzt < lz) {
+    /* Make boths inputs odd by putting powers of 2 on the side */
-        lz = lzt;
+    const size_t za = mbedtls_mpi_lsb(&TA);
-    }
+    const size_t zb = mbedtls_mpi_lsb(&TB);
+    MBEDTLS_MPI_CHK(mbedtls_mpi_shift_r(&TA, za));
+    MBEDTLS_MPI_CHK(mbedtls_mpi_shift_r(&TB, zb));
 
-    TA.s = TB.s = 1;
+    /* Ensure A <= B: if B < A, swap them */
+    mbedtls_ct_condition_t swap = mbedtls_mpi_core_lt_ct(TB.p, TA.p, TA.n);
+    mbedtls_mpi_core_cond_swap(TA.p, TB.p, TA.n, swap);
 
-    /* We mostly follow the procedure described in HAC 14.54, but with some
+    MBEDTLS_MPI_CHK(mbedtls_mpi_gcd_modinv_odd(G, NULL, &TA, &TB));
-     * minor differences:
-     * - Sequences of multiplications or divisions by 2 are grouped into a
-     *   single shift operation.
-     * - The procedure in HAC assumes that 0 < TB <= TA.
-     *     - The condition TB <= TA is not actually necessary for correctness.
-     *       TA and TB have symmetric roles except for the loop termination
-     *       condition, and the shifts at the beginning of the loop body
-     *       remove any significance from the ordering of TA vs TB before
-     *       the shifts.
-     *     - If TA = 0, the loop goes through 0 iterations and the result is
-     *       correctly TB.
-     *     - The case TB = 0 was short-circuited above.
-     *
-     * For the correctness proof below, decompose the original values of
-     * A and B as
-     *   A = sa * 2^a * A' with A'=0 or A' odd, and sa = +-1
-     *   B = sb * 2^b * B' with B'=0 or B' odd, and sb = +-1
-     * Then gcd(A, B) = 2^{min(a,b)} * gcd(A',B'),
-     * and gcd(A',B') is odd or 0.
-     *
-     * At the beginning, we have TA = |A| and TB = |B| so gcd(A,B) = gcd(TA,TB).
-     * The code maintains the following invariant:
-     *     gcd(A,B) = 2^k * gcd(TA,TB) for some k   (I)
-     */
 
-    /* Proof that the loop terminates:
+    /* Re-inject the power of 2 we had previously put aside */
-     * At each iteration, either the right-shift by 1 is made on a nonzero
+    size_t zg = za > zb ? zb : za; // zg = min(za, zb)
-     * value and the nonnegative integer bitlen(TA) + bitlen(TB) decreases
+    MBEDTLS_MPI_CHK(mbedtls_mpi_shift_l(G, zg));
-     * by at least 1, or the right-shift by 1 is made on zero and then
-     * TA becomes 0 which ends the loop (TB cannot be 0 if it is right-shifted
-     * since in that case TB is calculated from TB-TA with the condition TB>TA).
-     */
-    while (mbedtls_mpi_cmp_int(&TA, 0) != 0) {
-        /* Divisions by 2 preserve the invariant (I). */
-        MBEDTLS_MPI_CHK(mbedtls_mpi_shift_r(&TA, mbedtls_mpi_lsb(&TA)));
-        MBEDTLS_MPI_CHK(mbedtls_mpi_shift_r(&TB, mbedtls_mpi_lsb(&TB)));
-
-        /* Set either TA or TB to |TA-TB|/2. Since TA and TB are both odd,
-         * TA-TB is even so the division by 2 has an integer result.
-         * Invariant (I) is preserved since any odd divisor of both TA and TB
-         * also divides |TA-TB|/2, and any odd divisor of both TA and |TA-TB|/2
-         * also divides TB, and any odd divisor of both TB and |TA-TB|/2 also
-         * divides TA.
-         */
-        if (mbedtls_mpi_cmp_mpi(&TA, &TB) >= 0) {
-            MBEDTLS_MPI_CHK(mbedtls_mpi_sub_abs(&TA, &TA, &TB));
-            MBEDTLS_MPI_CHK(mbedtls_mpi_shift_r(&TA, 1));
-        } else {
-            MBEDTLS_MPI_CHK(mbedtls_mpi_sub_abs(&TB, &TB, &TA));
-            MBEDTLS_MPI_CHK(mbedtls_mpi_shift_r(&TB, 1));
-        }
-        /* Note that one of TA or TB is still odd. */
-    }
-
-    /* By invariant (I), gcd(A,B) = 2^k * gcd(TA,TB) for some k.
-     * At the loop exit, TA = 0, so gcd(TA,TB) = TB.
-     * - If there was at least one loop iteration, then one of TA or TB is odd,
-     *   and TA = 0, so TB is odd and gcd(TA,TB) = gcd(A',B'). In this case,
-     *   lz = min(a,b) so gcd(A,B) = 2^lz * TB.
-     * - If there was no loop iteration, then A was 0, and gcd(A,B) = B.
-     *   In this case, lz = 0 and B = TB so gcd(A,B) = B = 2^lz * TB as well.
-     */
-
-    MBEDTLS_MPI_CHK(mbedtls_mpi_shift_l(&TB, lz));
-    MBEDTLS_MPI_CHK(mbedtls_mpi_copy(G, &TB));
 
 cleanup:
 

tests/suites/test_suite_bignum.misc.data

@@ -1466,10 +1466,10 @@ Test GCD: 6, 0 (1 limb)
 mpi_gcd:"06":"00":"6"
 
 Test GCD: negative, 0 (null)
-mpi_gcd:"-50000":"":"-50000"
+mpi_gcd:"-50000":"":"50000"
 
 Test GCD: negative, 0 (1 limb)
-mpi_gcd:"-a782374b2ee927df28802745833a":"00":"-a782374b2ee927df28802745833a"
+mpi_gcd:"-a782374b2ee927df28802745833a":"00":"a782374b2ee927df28802745833a"
 
 Test GCD: 0 (null), negative
 mpi_gcd:"":"-50000":"50000"