Implementation-defined behaviour in montgomery_reduce()

[rod-chapman]

Our attempts at formal verification of the montgomery_reduce() function in ref/reduce.c reveal that this function depends on a subtle implementation-defined behaviour - namely a cast from int32_t to int16_t when the valued being converted lies outside the range of int16_t. (See C17 or C23 standards, para 6.3.1.3 for details.)

This is unfortunate from a portability point of view. It also confounds understanding of the code, since the actual semantics typically implemented by gcc and clang is rather non-obvious.

For example, line 20 of reduce.c: t = (int16_t)a * QINV; the cast of a to int16_t exhibits implementation-defined behaviour. There is also a second (invisible, implicit) cast there to int16_t that follows the evaluation of the multiplication, just before the assignment to t takes place.

I have developed an alternative implementation that is free from such non-portable behaviour.

[hanno-becker]

Conceptually, a * QINV is an unsigned computation mod 2^16, so I think the conversion int32_t -> int16_t should be replaced by int32_t -> uint16_t, which is well-defined. However, the result must then be lifted to a signed canonical representative, and I have not come across a C implementation not implementing uint16_t -> int16_t as exactly that (anyone else?). Still, @rod-chapman's right that this is strictly speaking implementation-defined.

I'd be inclined to merely document this limitation (in the code and at the top-level stating the C language requirements for this repository) and change the code to something like this:

int16_t montgomery_reduce(int32_t a)
{
  uint16_t u;
  int16_t t;
  // Compute a*q^{-1} mod 2^16 in unsigned representatives
  u = (uint16_t)a * QINV;
  // Lift to signed canonical representative mod 2^16.
  // PORTABILITY: This relies on uint16_t -> int16_t
  // being implemented as the inverse of int16_t -> uint16_t,
  // which is not mandated by the standard.
  t = (int16_t) u;
  // By construction, the LHS is divisible by 2^16
  t = (a - (int32_t)t*KYBER_Q) >> 16;
  return t;
}

(NB: IIUC, semantically the arithmetic (uint16_t) a * QINV would still be signed (if sizeof(int) > 2) due to integer promotion, but that doesn't matter, again since cast to an unsigned type is always well-defined)

[cryptojedi]

@rod-chapman Sorry that I didn't follow up on this for so long. Can you share the version you wrote that avoids the implementation-defined behavior?

[rod-chapman]

I'll have a look...

[rod-chapman]

Let's start at the beginning - can we arrive at a post-condition for montgomery_reduce(), expressed as a legal and well-defined C expression? I'm struggling a bit with this. R-1 (mod Q) where R==216 is 169 as far as I can tell, so for r = montgomery_reduce(x); we require r % Q == (x * 169) % Q but that doesn't seem to work because C's % operator doesn't work that way for negative left operand.

Any ideas?

[cryptojedi]

The two post-conditions are (r*(1<<16)-x)%KYBER_Q == 0 and r >= -(KYBER_Q+1) && r <= KYBER_Q-1

[rod-chapman]

The comment in kyber/ref/reduce.c says Returns: integer in {-q+1,...,q-1} congruent to a * R^-1 modulo q. so is it really supposed to be r >= -(KYBER_Q+1) // as you have above or is it r >= (-KYBER_Q)+1

??

[rod-chapman]

OK... I have an implementation which I think is well-defined, portable, and proves correct with CBMC. BUT... it's slow - 14 instructions on AArch64 compared to 7 for the original, using gcc 13.1.0 -O3 It also depends on a conditional subtraction that is not guaranteed 100% to be constant-time.

[rod-chapman]

Here's the declaration with contracts: https://github.com/rod-chapman/cbmc-examples/blob/426f81b8cc7220a355cfa28efe688ab2a2d1b2cd/pqc/crypto.h#L30

and the implementation: https://github.com/rod-chapman/cbmc-examples/blob/426f81b8cc7220a355cfa28efe688ab2a2d1b2cd/pqc/crypto.c#L88