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==2**16 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

[cryptojedi]

Roderick Chapman ***@***.***> wrote:

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

I'm a bit worried about the t1s = t1s - (65536 * (t1s >= 32768)); This has some risk to be compiled to a branch. I assume that's what you mean by the conditional subtraction that isn't guaranteed to be constant time? Also, writing % and / operators may have a risk to compile to a DIV instruction, although I would hope that no sane compiler does that if the divisor is a power of two. I'm leaning towards the solution that Hanno proposes, i.e., add comments to document the assumptions on the compiler/environment.

[rod-chapman]

Yes... I too am concerned about the conditional subtraction being compiled to a branch. The "right" way to do it would be to generate a "mask" value using Dan Bernstein's library, but that will further erode performance.

And yes... I am assuming that % and / by a power of two will be implemented as a mask or shift reglardless of optimization level.

I am sympathetic that leaving the code alone and documenting the assumptions properly is better for now.

[hanno-becker]

There is similar implementation-defined behaviour in barrett_reduce(). We merged code with suitable PORTABILITY sections into mlkem-c-aarch64, see https://github.com/pq-code-package/mlkem-c-aarch64/blob/main/mlkem/reduce.c. Let us know if you'd like an upstream patch along those lines.