Write a proof of correctness for mbedtls_mpi_gcd
Signed-off-by: Gilles Peskine <Gilles.Peskine@arm.com>
diff --git a/library/bignum.c b/library/bignum.c
index 3b126e2..bae0392 100644
--- a/library/bignum.c
+++ b/library/bignum.c
@@ -2312,18 +2312,49 @@
TA.s = TB.s = 1;
- /* We follow the procedure described in HAC 14.54, except that sequences
- * of divisions by 2 are grouped into a single shift. The procedure in HAC
- * assumes that the numbers are initially positive. The case B=0 was
- * short-circuited above. If A=0, the loop goes through 0 iterations
- * and the result is correctly B.
+ /* We mostly follow the procedure described in HAC 14.54, but with some
+ * minor differences:
+ * - Sequences of multiplications or divisions by 2 are grouped into a
+ * single shift operation.
+ * - The procedure in HAC assumes that 0 < A <= B.
+ * - The condition A <= B is not actually necessary for correctness;
+ * the first round through the loop results in TA < TB.
+ * - If A = 0, the loop goes through 0 iterations and the result is
+ * correctly B.
+ * - The case B=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|/2^a and TB = |B|/2^b.
+ * The code maintains the following invariant:
+ * gcd(A,B) = 2^k * gcd(TA,TB) for some k (I)
*/
+ /* Proof that the loop terminates:
+ * At each iteration, either the right-shift by 1 is made on a nonzero
+ * value and the nonnegative integer bitlen(TA) + bitlen(TB) decreases
+ * 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 divisior 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 ) );
@@ -2334,8 +2365,18 @@
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, B = 2^lz * TB so gcd(A,B) = 2^lz * TB as well.
+ */
+
MBEDTLS_MPI_CHK( mbedtls_mpi_shift_l( &TB, lz ) );
MBEDTLS_MPI_CHK( mbedtls_mpi_copy( G, &TB ) );