-
Notifications
You must be signed in to change notification settings - Fork 12
/
Copy pathreduce.c
154 lines (137 loc) · 4.84 KB
/
reduce.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
/*
* Copyright (c) 2024 The mlkem-native project authors
* SPDX-License-Identifier: Apache-2.0
*/
#include "reduce.h"
#include <stdint.h>
#include "params.h"
/* QINV == -3327 converted to uint16_t == -3327 + 65536 == 62209 */
static const uint32_t QINV = 62209; /* q^-1 mod 2^16 */
/*************************************************
* Name: cast_uint16_to_int16
*
* Description: Cast uint16 value to int16
*
* Returns:
* input x in 0 .. 32767: returns value unchanged
* input x in 32768 .. 65535: returns (x - 65536)
**************************************************/
#ifdef CBMC
#pragma CPROVER check push
#pragma CPROVER check disable "conversion"
#endif
static INLINE int16_t cast_uint16_to_int16(uint16_t x)
{
/*
* PORTABILITY: This relies on uint16_t -> int16_t
* being implemented as the inverse of int16_t -> uint16_t,
* which is implementation-defined (C99 6.3.1.3 (3))
* CBMC (correctly) fails to prove this conversion is OK,
* so we have to suppress that check here
*/
return (int16_t)x;
}
#ifdef CBMC
#pragma CPROVER check pop
#endif
/*************************************************
* Name: montgomery_reduce_generic
*
* Description: Generic Montgomery reduction; given a 32-bit integer a, computes
* 16-bit integer congruent to a * R^-1 mod q, where R=2^16
*
* Arguments: - int32_t a: input integer to be reduced
*
* Returns: integer congruent to a * R^-1 modulo q
*
* Bounds: For any C such that |a| < q * C, the return value
* has absolute value < q (C/2^16 + 1/2).
*
* Notable special cases:
* - The Montgomery multiplication of a value of absolute value
* < q * C with a signed-canonical value ( < q/2 ) has
* absolute value q * (0.0254 * C + 1/2).
* - The Montgomery multiplication of a value of absolute value
* < q * C with a value t of |t| < q has absolute value
* < q * (0.0508 * C + 1/2).
* - The Montgomery multiplication of a value of absolute value
* < C with a value of abs < q has absolute value
* < q (C/2^16 + 1/2).
**************************************************/
ALWAYS_INLINE
static INLINE int16_t montgomery_reduce_generic(int32_t a)
{
/*
*Bounds on paper
* - Case |a| < q * C, for some C
* |t| <= |a|/2^16 + |t|*q/2^16
* < q * C / 2^16 + q/2
* = q (C/2^16 + 1/2)
* - Case |a| < (q/2) * C * q, for some C
* Replace C -> C * q in the above and estimate
* q / 2^17 < 0.0254.
*/
/* Compute a*q^{-1} mod 2^16 in unsigned representatives */
const uint16_t a_reduced = a & UINT16_MAX;
const uint16_t a_inverted = (a_reduced * QINV) & UINT16_MAX;
/* Lift to signed canonical representative mod 2^16. */
const int16_t t = cast_uint16_to_int16(a_inverted);
int32_t r = a - ((int32_t)t * MLKEM_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 >> 16;
return (int16_t)r;
}
int16_t montgomery_reduce(int32_t a)
{
int16_t res;
SCALAR_BOUND(a, 2 * MLKEM_Q * 32768, "montgomery_reduce input");
res = montgomery_reduce_generic(a);
SCALAR_BOUND(res, (3 * (MLKEM_Q + 1)) / 2, "montgomery_reduce output");
return res;
}
int16_t fqmul(int16_t a, int16_t b)
{
int16_t res;
SCALAR_BOUND(b, HALF_Q, "fqmul input");
res = montgomery_reduce((int32_t)a * (int32_t)b);
SCALAR_BOUND(res, MLKEM_Q, "fqmul output");
return res;
}
/*
* To divide by MLKEM_Q using Barrett multiplication, the "magic number"
* multiplier is round_to_nearest(2**26/MLKEM_Q)
*/
#define BPOWER 26
static const int32_t barrett_multiplier =
((1 << BPOWER) + MLKEM_Q / 2) / MLKEM_Q;
/*************************************************
* Name: barrett_reduce
*
* Description: Barrett reduction; given a 16-bit integer a, computes
* centered representative congruent to a mod q in
* {-(q-1)/2,...,(q-1)/2}
*
* Arguments: - int16_t a: input integer to be reduced
*
* Returns: integer in {-(q-1)/2,...,(q-1)/2} congruent to a modulo q.
**************************************************/
int16_t barrett_reduce(int16_t a)
{
/*
* Compute round_to_nearest(a/MLKEM_Q) using the multiplier
* above and shift by BPOWER places.
* 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))
*/
const int32_t t = (barrett_multiplier * a + (1 << (BPOWER - 1))) >> BPOWER;
/*
* t is in -10 .. +10, so we need 32-bit math to
* evaluate t * MLKEM_Q and the subsequent subtraction
*/
return (int16_t)(a - t * MLKEM_Q);
}