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. It is explained why
under the _inclusive_ input bound 2^31 * MLDSA_Q one gets an
_inclusive_ output bound by MLDSA_Q. However, a pen-and-paper
argument that the same holds for _exclusive_ input/output bounds
-- which is what the 'strong' CBMC spec states -- is still missing.

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/reduce.h

@@ -21,41 +21,33 @@
 /* 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
+ *                <= ceil(|a| / 2^32 + MLDSA_Q / 2)
  *
- * Arguments:   - int64_t: finite field element a
+ *              In particular, if |a| <= 2^31 * MLDSA_Q, the absolute value
+ *              of the return value is <= MLDSA_Q.
  *
- * Returns r.
+ *              For |a| < 2^31 * MLDSA_Q, the return value is even < MLDSA_Q
+ *              in absolute value. This is what the CBMC specification encodes.
  **************************************************/
 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)INT32_MIN * MLDSA_Q) &&
+	   a < -((int64_t)INT32_MIN * MLDSA_Q))
+  ensures(return_value > -MLDSA_Q && return_value < MLDSA_Q)
 )
 {
   /* check-magic: 58728449 == unsigned_mod(pow(MLDSA_Q, -1, 2^32), 2^32) */
@@ -68,16 +60,38 @@ __contract__(
   /* Lift to signed canonical representative mod 2^32. */
   const int32_t t = mld_cast_uint32_to_int32(a_inverted);
 
-  int64_t r;
+  int64_t r, tp;
+
+  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);
+  tp = (int64_t)t * MLDSA_Q;
+  r = a - tp;
 
   /*
    * 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: |r >> 32| <= ceil(|r| / 2^32)
+   *                   <= ceil(|a| / 2^32 + MLDSA_Q / 2)
+   *
+   * (Note that |a >> n| = ceil(|a| / 2^n) for negative a).
+   * This is the bound stated in the documentation.
+   *
+   * Note that for a = 2^31 * MLDSA_Q, we get return_value = MLDSA_Q, so the
+   * bound is strict under the assumption |a| <= 2^31 * MLDSA_Q.
+   *
+   * However, if we assume a _strict_ input inequality |a| < 2^31 * MLDSA_Q,
+   * the output inequality is _also_ strict:
+   *
+   *   |return_value| < MLDSA_Q
+   *
+   * This is what is captured by the CBMC spec.
+   *
+   * TODO: Find a pen-and-paper argument.
+   */
   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