Remove unused spec for mld_montgomery_reduce, add bounds reasoning

Context:

We previously had two specifications for mld_montgomery_reduce:

- A 'weak' one where the output upper bound matches the input upper
  bound of mld_reduce32().
- A 'strong' one where the output interval is (-Q,Q), exclusively.

Both contracts were combined into one in a somewhat ad-hoc way
(the 'weak' version being encoded as pre- and post-conditions as
usual, the 'strong' version being stated as an implication in the
post-condition).

The combined contract was used in the CBMC proofs of
- mld_poly_pointwise_montgomery
- mld_polyvecl_pointwise_acc_montgomery

There are two more functions which rely on mld_montgomery_reduce,
namely
- mld_fqmul
- mld_fqscale

Observations:
- The 'weak' contract of mld_montgomery_reduce is not needed:
  In the context of mld_poly_pointwise_montgomery and
  mld_polyvecl_pointwise_acc_montgomery, the strong version applies.
- The 'strong' version also applies in the context of mld_fqmul
  and mld_fqscale.

Accordingly, this commit removes the 'weak' version and re-formulates
the 'strong' version in traditional pre/post condition form.

The commit also removes the macro wrapper around the input bound,
which I find difficult to read. It's easier to understand if the
bound is explicit.

The 'strong' version only applies to `mld_fqscale` because of the
recent modification of the iNTT output bound from 3/4*Q to Q. Otherwise,
an even stronger spec would have been needed. At present, one can
even remove mld_fqscale and just call mld_fqmul directly, but this
commit does not do this.

The commit also takes the opportunity to add some pen-and-paper
bounds reasoning for mld_montgomery_reduce.

Fixes #602

Signed-off-by: Hanno Becker <beckphan@amazon.co.uk>

Author: hanno-becker

mldsa/mldsa_native.c

@@ -374,8 +374,6 @@
 /* mldsa/src/reduce.h */
 #undef MLD_REDUCE_H
 #undef MONT
-#undef MONTGOMERY_REDUCE_DOMAIN_MAX
-#undef MONTGOMERY_REDUCE_STRONG_DOMAIN_MAX
 #undef REDUCE32_DOMAIN_MAX
 #undef REDUCE32_RANGE_MAX
 /* mldsa/src/symmetric.h */

mldsa/src/ntt.c

@@ -27,8 +27,15 @@ __contract__(
   ensures(return_value > -MLDSA_Q && return_value < MLDSA_Q)
 )
 {
+  /* Bounds: We argue in mld_montgomery_reduce() that the reult
+   * of Montgomery reduction is < MLDSA_Q if the input is smaller
+   * than 2^31 * MLDSA_Q in absolute value. Indeed, we have:
+   *
+   *    |a * b|   = |a| * |b|
+   *              < 2^31 * MLDSA_Q_HALF
+   *              < 2^31 * MLDSA_Q
+   */
   return mld_montgomery_reduce((int64_t)a * (int64_t)b);
-  /* TODO: reason about bounds */
 }
 
 /*************************************************
@@ -50,8 +57,9 @@ __contract__(
 {
   /* check-magic: 41978 == pow(2,64-8,MLDSA_Q) */
   const int32_t f = 41978;
+  /* Bounds: MLD_INTT_BOUND is MLDSA_Q, so the bounds reasoning is just
+   * a special case of that in mld_fqmul(). */
   return mld_montgomery_reduce((int64_t)a * f);
-  /* TODO: reason about bounds */
 }
 
 #include "zetas.inc"

mldsa/src/reduce.h

@@ -21,41 +21,30 @@
 /* check-magic: 6283009 == (REDUCE32_DOMAIN_MAX - 255 * MLDSA_Q + 1) */
 #define REDUCE32_RANGE_MAX 6283009
 
-/* Absolute bound for domain of mld_montgomery_reduce() */
-#define MONTGOMERY_REDUCE_DOMAIN_MAX ((int64_t)INT32_MIN * INT32_MIN)
-
-/* Absolute bound for tight domain of mld_montgomery_reduce() */
-#define MONTGOMERY_REDUCE_STRONG_DOMAIN_MAX ((int64_t)INT32_MIN * -MLDSA_Q)
-
 /*************************************************
  * Name:        mld_montgomery_reduce
  *
- * Description:
- *      For finite field element a with
- *        -MONTGOMERY_REDUCE_DOMAIN_MAX <= a <= MONTGOMERY_REDUCE_DOMAIN_MAX,
- *      compute r == a*2^{-32} (mod MLDSA_Q) such that
- *        INT32_MIN <= r < REDUCE32_DOMAIN_MAX
+ * Description: Generic Montgomery reduction; given a 64-bit integer a, computes
+ *              32-bit integer congruent to a * R^-1 mod q, where R=2^32
  *
- *      The upper-bound on the result ensures that a result from this
- *      function can be used as an input to mld_reduce32() declared below.
+ * Arguments:   - int64_t a: input integer to be reduced, of absolute value
+ *                smaller or equal to INT64_MAX - 2^31 * MLDSA_Q.
  *
- *      Additionally, as a special case, if the input a is in range
- *        -MONTGOMERY_REDUCE_STRONG_DOMAIN_MAX < a <
- *          MONTGOMERY_REDUCE_STRONG_DOMAIN_MAX
- *      then the result satisfies -MLDSA_Q < r < MLDSA_Q.
+ * Returns:     Integer congruent to a * R^-1 modulo q, with absolute value
+ *                <= |a| / 2^32 + MLDSA_Q / 2
  *
- * Arguments:   - int64_t: finite field element a
- *
- * Returns r.
+ *              In particular, if |a| < 2^31 * MLDSA_Q, the absolute value
+ *              of the return value is < MLDSA_Q.
  **************************************************/
 static MLD_INLINE int32_t mld_montgomery_reduce(int64_t a)
 __contract__(
-  requires(a >= -MONTGOMERY_REDUCE_DOMAIN_MAX && a <= MONTGOMERY_REDUCE_DOMAIN_MAX)
-  ensures(return_value >= INT32_MIN && return_value < REDUCE32_DOMAIN_MAX)
-
-  /* Special case - for stronger input bounds, we can ensure stronger bounds on the output */
-  ensures((a >= -MONTGOMERY_REDUCE_STRONG_DOMAIN_MAX && a < MONTGOMERY_REDUCE_STRONG_DOMAIN_MAX) ==>
-          (return_value > -MLDSA_Q && return_value < MLDSA_Q))
+  /* We don't attempt to express an input-dependent output bound
+   * as the post-condition here, as all call-sites satisfy the
+   * absolute input bound 2^31 * MLDSA_Q and higher-level
+   * reasoning can be conducted using |return_value| < MLDSA_Q. */
+  requires(a > -(((int64_t)1 << 31) * MLDSA_Q) &&
+           a <  (((int64_t)1 << 31) * MLDSA_Q))
+  ensures(return_value > -MLDSA_Q && return_value < MLDSA_Q)
 )
 {
   /* check-magic: 58728449 == unsigned_mod(pow(MLDSA_Q, -1, 2^32), 2^32) */
@@ -70,14 +59,30 @@ __contract__(
 
   int64_t r;
 
-  r = a - ((int64_t)t * MLDSA_Q);
+  mld_assert(a < +(INT64_MAX - (((int64_t)1 << 31) * MLDSA_Q)) &&
+             a > -(INT64_MAX - (((int64_t)1 << 31) * MLDSA_Q)));
+
+  r = a - (int64_t)t * MLDSA_Q;
 
   /*
    * PORTABILITY: Right-shift on a signed integer is, strictly-speaking,
    * implementation-defined for negative left argument. Here,
    * we assume it's sign-preserving "arithmetic" shift right. (C99 6.5.7 (5))
    */
   r = r >> 32;
+
+  /* Bounds:
+   *
+   * By construction of the Montgomery multiplication, by the time we
+   * compute r >> 32, r is divisible by 2^32, and hence
+   *
+   *   |r >> 32|  = |r| / 2^32
+   *             <= |a| / 2^32 + MLDSA_Q / 2
+   *
+   * (In general, we would only have |x >> n| <= ceil(|x| / 2^n)).
+   *
+   * In particular, if |a| < 2^31 * MLDSA_Q, then |return_value| < MLDSA_Q.
+   */
   return (int32_t)r;
 }
 

proofs/cbmc/fqmul/Makefile

@@ -20,7 +20,7 @@ PROOF_SOURCES += $(PROOFDIR)/$(HARNESS_FILE).c
 PROJECT_SOURCES += $(SRCDIR)/mldsa/src/ntt.c
 
 CHECK_FUNCTION_CONTRACTS=mld_fqmul
-USE_FUNCTION_CONTRACTS=
+USE_FUNCTION_CONTRACTS=mld_montgomery_reduce
 APPLY_LOOP_CONTRACTS=on
 USE_DYNAMIC_FRAMES=1
 

proofs/cbmc/fqscale/Makefile

@@ -20,7 +20,7 @@ PROOF_SOURCES += $(PROOFDIR)/$(HARNESS_FILE).c
 PROJECT_SOURCES += $(SRCDIR)/mldsa/src/ntt.c
 
 CHECK_FUNCTION_CONTRACTS=mld_fqscale
-USE_FUNCTION_CONTRACTS=
+USE_FUNCTION_CONTRACTS=mld_montgomery_reduce
 APPLY_LOOP_CONTRACTS=on
 USE_DYNAMIC_FRAMES=1