Commit 1ca3a48b873dc8cc00ad8dfb74794856ecd58764 | zig

author	Frank Denis <124872+jedisct1@users.noreply.github.com> 2024-03-10 15:30:13+01:00
committer	Frank Denis <124872+jedisct1@users.noreply.github.com> 2024-03-10 15:48:38+01:00
log	1ca3a48b873dc8cc00ad8dfb74794856ecd58764
tree	d5678ee365fe01c777799f59b70289d419bd0a17
parent	4ba4f94c93d5eb1945f1b2c8c53a45cbee609d3b

std.crypto: add support for ML-KEM

ML-KEM is the Kyber post-quantum secure key encapsulation mechanism,
as being standardized by NIST.

Too bad, they decided to rename it; the "Kyber" name was so much
better!

This implements the current draft (NIST FIPS-203), which is already
being deployed even though the specification is not finalized.

3 files changed, 1832 insertions(+), 1784 deletions(-)

lib
- std
  - crypto.zig
  - crypto
    kyber_d00.zig
    ml_kem.zig

lib/std/crypto.zig+2-1

...	...	@@ -70,7 +70,8 @@ pub const dh = struct {
70	70
71	71	/// Key Encapsulation Mechanisms.
72	72	pub const kem = struct {
73		pub const kyber_d00 = @import("crypto/kyber_d00.zig");
	73	pub const kyber_d00 = @import("crypto/ml_kem.zig").kyber_d00;
	74	pub const ml_kem_01 = @import("crypto/ml_kem.zig").ml_kem_01;
74	75	};
75	76
76	77	/// Elliptic-curve arithmetic.

lib/std/crypto/kyber_d00.zig deleted-1783

...	...	@@ -1,1783 +0,0 @@
1		//! Implementation of the IND-CCA2 post-quantum secure key encapsulation
2		//! mechanism (KEM) CRYSTALS-Kyber, as submitted to the third round of the NIST
3		//! Post-Quantum Cryptography (v3.02/"draft00"), and selected for standardisation.
4		//!
5		//! Kyber will likely change before final standardisation.
6		//!
7		//! The namespace suffix (currently `_d00`) refers to the version currently
8		//! implemented, in accordance with the draft. It may not be updated if new
9		//! versions of the draft only include editorial changes.
10		//!
11		//! The suffix will eventually be removed once Kyber is finalized.
12		//!
13		//! Quoting from the CFRG I-D:
14		//!
15		//! Kyber is not a Diffie-Hellman (DH) style non-interactive key
16		//! agreement, but instead, Kyber is a Key Encapsulation Method (KEM).
17		//! In essence, a KEM is a Public-Key Encryption (PKE) scheme where the
18		//! plaintext cannot be specified, but is generated as a random key as
19		//! part of the encryption. A KEM can be transformed into an unrestricted
20		//! PKE using HPKE (RFC9180). On its own, a KEM can be used as a key
21		//! agreement method in TLS.
22		//!
23		//! Kyber is an IND-CCA2 secure KEM. It is constructed by applying a
24		//! Fujisaki--Okamato style transformation on InnerPKE, which is the
25		//! underlying IND-CPA secure Public Key Encryption scheme. We cannot
26		//! use InnerPKE directly, as its ciphertexts are malleable.
27		//!
28		//! ```
29		//! F.O. transform
30		//! InnerPKE ----------------------> Kyber
31		//! IND-CPA IND-CCA2
32		//! ```
33		//!
34		//! Kyber is a lattice-based scheme. More precisely, its security is
35		//! based on the learning-with-errors-and-rounding problem in module
36		//! lattices (MLWER). The underlying polynomial ring R (defined in
37		//! Section 5) is chosen such that multiplication is very fast using the
38		//! number theoretic transform (NTT, see Section 5.1.3).
39		//!
40		//! An InnerPKE private key is a vector _s_ over R of length k which is
41		//! _small_ in a particular way. Here k is a security parameter akin to
42		//! the size of a prime modulus. For Kyber512, which targets AES-128's
43		//! security level, the value of k is 2.
44		//!
45		//! The public key consists of two values:
46		//!
47		//! * _A_ a uniformly sampled k by k matrix over R _and_
48		//!
49		//! * _t = A s + e_, where e is a suitably small masking vector.
50		//!
51		//! Distinguishing between such A s + e and a uniformly sampled t is the
52		//! module learning-with-errors (MLWE) problem. If that is hard, then it
53		//! is also hard to recover the private key from the public key as that
54		//! would allow you to distinguish between those two.
55		//!
56		//! To save space in the public key, A is recomputed deterministically
57		//! from a seed _rho_.
58		//!
59		//! A ciphertext for a message m under this public key is a pair (c_1,
60		//! c_2) computed roughly as follows:
61		//!
62		//! c_1 = Compress(A^T r + e_1, d_u)
63		//! c_2 = Compress(t^T r + e_2 + Decompress(m, 1), d_v)
64		//!
65		//! where
66		//!
67		//! * e_1, e_2 and r are small blinds;
68		//!
69		//! * Compress(-, d) removes some information, leaving d bits per
70		//! coefficient and Decompress is such that Compress after Decompress
71		//! does nothing and
72		//!
73		//! * d_u, d_v are scheme parameters.
74		//!
75		//! Distinguishing such a ciphertext and uniformly sampled (c_1, c_2) is
76		//! an example of the full MLWER problem, see section 4.4 of [KyberV302].
77		//!
78		//! To decrypt the ciphertext, one computes
79		//!
80		//! m = Compress(Decompress(c_2, d_v) - s^T Decompress(c_1, d_u), 1).
81		//!
82		//! It it not straight-forward to see that this formula is correct. In
83		//! fact, there is negligible but non-zero probability that a ciphertext
84		//! does not decrypt correctly given by the DFP column in Table 4. This
85		//! failure probability can be computed by a careful automated analysis
86		//! of the probabilities involved, see kyber_failure.py of [SecEst].
87		//!
88		//! [KyberV302](https://pq-crystals.org/kyber/data/kyber-specification-round3-20210804.pdf)
89		//! [I-D](https://github.com/bwesterb/draft-schwabe-cfrg-kyber)
90		//! [SecEst](https://github.com/pq-crystals/security-estimates)
91
92		// TODO
93		//
94		// - The bottleneck in Kyber are the various hash/xof calls:
95		// - Optimize Zig's keccak implementation.
96		// - Use SIMD to compute keccak in parallel.
97		// - Can we track bounds of coefficients using comptime types without
98		// duplicating code?
99		// - Would be neater to have tests closer to the thing under test.
100		// - When generating a keypair, we have a copy of the inner public key with
101		// its large matrix A in both the public key and the private key. In Go we
102		// can just have a pointer in the private key to the public key, but
103		// how do we do this elegantly in Zig?
104
105		const std = @import("std");
106		const builtin = @import("builtin");
107
108		const testing = std.testing;
109		const assert = std.debug.assert;
110		const crypto = std.crypto;
111		const math = std.math;
112		const mem = std.mem;
113		const RndGen = std.Random.DefaultPrng;
114		const sha3 = crypto.hash.sha3;
115
116		// Q is the parameter q ≡ 3329 = 2¹¹ + 2¹⁰ + 2⁸ + 1.
117		const Q: i16 = 3329;
118
119		// Montgomery R
120		const R: i32 = 1 << 16;
121
122		// Parameter n, degree of polynomials.
123		const N: usize = 256;
124
125		// Size of "small" vectors used in encryption blinds.
126		const eta2: u8 = 2;
127
128		const Params = struct {
129		name: []const u8,
130
131		// Width and height of the matrix A.
132		k: u8,
133
134		// Size of "small" vectors used in private key and encryption blinds.
135		eta1: u8,
136
137		// How many bits to retain of u, the private-key independent part
138		// of the ciphertext.
139		du: u8,
140
141		// How many bits to retain of v, the private-key dependent part
142		// of the ciphertext.
143		dv: u8,
144		};
145
146		pub const Kyber512 = Kyber(.{
147		.name = "Kyber512",
148		.k = 2,
149		.eta1 = 3,
150		.du = 10,
151		.dv = 4,
152		});
153
154		pub const Kyber768 = Kyber(.{
155		.name = "Kyber768",
156		.k = 3,
157		.eta1 = 2,
158		.du = 10,
159		.dv = 4,
160		});
161
162		pub const Kyber1024 = Kyber(.{
163		.name = "Kyber1024",
164		.k = 4,
165		.eta1 = 2,
166		.du = 11,
167		.dv = 5,
168		});
169
170		const modes = [_]type{ Kyber512, Kyber768, Kyber1024 };
171		const h_length: usize = 32;
172		const inner_seed_length: usize = 32;
173		const common_encaps_seed_length: usize = 32;
174		const common_shared_key_size: usize = 32;
175
176		fn Kyber(comptime p: Params) type {
177		return struct {
178		// Size of a ciphertext, in bytes.
179		pub const ciphertext_length = Poly.compressedSize(p.du) * p.k + Poly.compressedSize(p.dv);
180
181		const Self = @This();
182		const V = Vec(p.k);
183		const M = Mat(p.k);
184
185		/// Length (in bytes) of a shared secret.
186		pub const shared_length = common_shared_key_size;
187		/// Length (in bytes) of a seed for deterministic encapsulation.
188		pub const encaps_seed_length = common_encaps_seed_length;
189		/// Length (in bytes) of a seed for key generation.
190		pub const seed_length: usize = inner_seed_length + shared_length;
191		/// Algorithm name.
192		pub const name = p.name;
193
194		/// A shared secret, and an encapsulated (encrypted) representation of it.
195		pub const EncapsulatedSecret = struct {
196		shared_secret: [shared_length]u8,
197		ciphertext: [ciphertext_length]u8,
198		};
199
200		/// A Kyber public key.
201		pub const PublicKey = struct {
202		pk: InnerPk,
203
204		// Cached
205		hpk: [h_length]u8, // H(pk)
206
207		/// Size of a serialized representation of the key, in bytes.
208		pub const bytes_length = InnerPk.bytes_length;
209
210		/// Generates a shared secret, and encapsulates it for the public key.
211		/// If `seed` is `null`, a random seed is used. This is recommended.
212		/// If `seed` is set, encapsulation is deterministic.
213		pub fn encaps(pk: PublicKey, seed_: ?[encaps_seed_length]u8) EncapsulatedSecret {
214		const seed = seed_ orelse seed: {
215		var random_seed: [encaps_seed_length]u8 = undefined;
216		crypto.random.bytes(&random_seed);
217		break :seed random_seed;
218		};
219
220		var m: [inner_plaintext_length]u8 = undefined;
221
222		// m = H(seed)
223		var h = sha3.Sha3_256.init(.{});
224		h.update(&seed);
225		h.final(&m);
226
227		// (K', r) = G(m ‖ H(pk))
228		var kr: [inner_plaintext_length + h_length]u8 = undefined;
229		var g = sha3.Sha3_512.init(.{});
230		g.update(&m);
231		g.update(&pk.hpk);
232		g.final(&kr);
233
234		// c = innerEncrypy(pk, m, r)
235		const ct = pk.pk.encrypt(&m, kr[32..64]);
236
237		// Compute H(c) and put in second slot of kr, which will be (K', H(c)).
238		h = sha3.Sha3_256.init(.{});
239		h.update(&ct);
240		h.final(kr[32..64]);
241
242		// K = KDF(K' ‖ H(c))
243		var kdf = sha3.Shake256.init(.{});
244		kdf.update(&kr);
245		var ss: [shared_length]u8 = undefined;
246		kdf.squeeze(&ss);
247
248		return EncapsulatedSecret{
249		.shared_secret = ss,
250		.ciphertext = ct,
251		};
252		}
253
254		/// Serializes the key into a byte array.
255		pub fn toBytes(pk: PublicKey) [bytes_length]u8 {
256		return pk.pk.toBytes();
257		}
258
259		/// Deserializes the key from a byte array.
260		pub fn fromBytes(buf: *const [bytes_length]u8) !PublicKey {
261		var ret: PublicKey = undefined;
262		ret.pk = InnerPk.fromBytes(buf[0..InnerPk.bytes_length]);
263
264		var h = sha3.Sha3_256.init(.{});
265		h.update(buf);
266		h.final(&ret.hpk);
267		return ret;
268		}
269		};
270
271		/// A Kyber secret key.
272		pub const SecretKey = struct {
273		sk: InnerSk,
274		pk: InnerPk,
275		hpk: [h_length]u8, // H(pk)
276		z: [shared_length]u8,
277
278		/// Size of a serialized representation of the key, in bytes.
279		pub const bytes_length: usize =
280		InnerSk.bytes_length + InnerPk.bytes_length + h_length + shared_length;
281
282		/// Decapsulates the shared secret within ct using the private key.
283		pub fn decaps(sk: SecretKey, ct: *const [ciphertext_length]u8) ![shared_length]u8 {
284		// m' = innerDec(ct)
285		const m2 = sk.sk.decrypt(ct);
286
287		// (K'', r') = G(m' ‖ H(pk))
288		var kr2: [64]u8 = undefined;
289		var g = sha3.Sha3_512.init(.{});
290		g.update(&m2);
291		g.update(&sk.hpk);
292		g.final(&kr2);
293
294		// ct' = innerEnc(pk, m', r')
295		const ct2 = sk.pk.encrypt(&m2, kr2[32..64]);
296
297		// Compute H(ct) and put in the second slot of kr2 which will be (K'', H(ct)).
298		var h = sha3.Sha3_256.init(.{});
299		h.update(ct);
300		h.final(kr2[32..64]);
301
302		// Replace K'' by z in the first slot of kr2 if ct ≠ ct'.
303		cmov(32, kr2[0..32], sk.z, ctneq(ciphertext_length, ct.*, ct2));
304
305		// K = KDF(K''/z, H(c))
306		var kdf = sha3.Shake256.init(.{});
307		var ss: [shared_length]u8 = undefined;
308		kdf.update(&kr2);
309		kdf.squeeze(&ss);
310		return ss;
311		}
312
313		/// Serializes the key into a byte array.
314		pub fn toBytes(sk: SecretKey) [bytes_length]u8 {
315		return sk.sk.toBytes() ++ sk.pk.toBytes() ++ sk.hpk ++ sk.z;
316		}
317
318		/// Deserializes the key from a byte array.
319		pub fn fromBytes(buf: *const [bytes_length]u8) !SecretKey {
320		var ret: SecretKey = undefined;
321		comptime var s: usize = 0;
322		ret.sk = InnerSk.fromBytes(buf[s .. s + InnerSk.bytes_length]);
323		s += InnerSk.bytes_length;
324		ret.pk = InnerPk.fromBytes(buf[s .. s + InnerPk.bytes_length]);
325		s += InnerPk.bytes_length;
326		ret.hpk = buf[s..][0..h_length].*;
327		s += h_length;
328		ret.z = buf[s..][0..shared_length].*;
329		return ret;
330		}
331		};
332
333		/// A Kyber key pair.
334		pub const KeyPair = struct {
335		secret_key: SecretKey,
336		public_key: PublicKey,
337
338		/// Create a new key pair.
339		/// If seed is null, a random seed will be generated.
340		/// If a seed is provided, the key pair will be determinsitic.
341		pub fn create(seed_: ?[seed_length]u8) !KeyPair {
342		const seed = seed_ orelse sk: {
343		var random_seed: [seed_length]u8 = undefined;
344		crypto.random.bytes(&random_seed);
345		break :sk random_seed;
346		};
347		var ret: KeyPair = undefined;
348		ret.secret_key.z = seed[inner_seed_length..seed_length].*;
349
350		// Generate inner key
351		innerKeyFromSeed(
352		seed[0..inner_seed_length].*,
353		&ret.public_key.pk,
354		&ret.secret_key.sk,
355		);
356		ret.secret_key.pk = ret.public_key.pk;
357
358		// Copy over z from seed.
359		ret.secret_key.z = seed[inner_seed_length..seed_length].*;
360
361		// Compute H(pk)
362		var h = sha3.Sha3_256.init(.{});
363		h.update(&ret.public_key.pk.toBytes());
364		h.final(&ret.secret_key.hpk);
365		ret.public_key.hpk = ret.secret_key.hpk;
366
367		return ret;
368		}
369		};
370
371		// Size of plaintexts of the in
372		const inner_plaintext_length: usize = Poly.compressedSize(1);
373
374		const InnerPk = struct {
375		rho: [32]u8, // ρ, the seed for the matrix A
376		th: V, // NTT(t), normalized
377
378		// Cached values
379		aT: M,
380
381		const bytes_length = V.bytes_length + 32;
382
383		fn encrypt(
384		pk: InnerPk,
385		pt: *const [inner_plaintext_length]u8,
386		seed: *const [32]u8,
387		) [ciphertext_length]u8 {
388		// Sample r, e₁ and e₂ appropriately
389		const rh = V.noise(p.eta1, 0, seed).ntt().barrettReduce();
390		const e1 = V.noise(eta2, p.k, seed);
391		const e2 = Poly.noise(eta2, 2 * p.k, seed);
392
393		// Next we compute u = Aᵀ r + e₁. First Aᵀ.
394		var u: V = undefined;
395		for (0..p.k) \|i\| {
396		// Note that coefficients of r are bounded by q and those of Aᵀ
397		// are bounded by 4.5q and so their product is bounded by 2¹⁵q
398		// as required for multiplication.
399		u.ps[i] = pk.aT.vs[i].dotHat(rh);
400		}
401
402		// Aᵀ and r were not in Montgomery form, so the Montgomery
403		// multiplications in the inner product added a factor R⁻¹ which
404		// the InvNTT cancels out.
405		u = u.barrettReduce().invNTT().add(e1).normalize();
406
407		// Next, compute v = <t, r> + e₂ + Decompress_q(m, 1)
408		const v = pk.th.dotHat(rh).barrettReduce().invNTT()
409		.add(Poly.decompress(1, pt)).add(e2).normalize();
410
411		return u.compress(p.du) ++ v.compress(p.dv);
412		}
413
414		fn toBytes(pk: InnerPk) [bytes_length]u8 {
415		return pk.th.toBytes() ++ pk.rho;
416		}
417
418		fn fromBytes(buf: *const [bytes_length]u8) InnerPk {
419		var ret: InnerPk = undefined;
420		ret.th = V.fromBytes(buf[0..V.bytes_length]).normalize();
421		ret.rho = buf[V.bytes_length..bytes_length].*;
422		ret.aT = M.uniform(ret.rho, true);
423		return ret;
424		}
425		};
426
427		// Private key of the inner PKE
428		const InnerSk = struct {
429		sh: V, // NTT(s), normalized
430		const bytes_length = V.bytes_length;
431
432		fn decrypt(sk: InnerSk, ct: *const [ciphertext_length]u8) [inner_plaintext_length]u8 {
433		const u = V.decompress(p.du, ct[0..comptime V.compressedSize(p.du)]);
434		const v = Poly.decompress(
435		p.dv,
436		ct[comptime V.compressedSize(p.du)..ciphertext_length],
437		);
438
439		// Compute m = v - <s, u>
440		return v.sub(sk.sh.dotHat(u.ntt()).barrettReduce().invNTT())
441		.normalize().compress(1);
442		}
443
444		fn toBytes(sk: InnerSk) [bytes_length]u8 {
445		return sk.sh.toBytes();
446		}
447
448		fn fromBytes(buf: *const [bytes_length]u8) InnerSk {
449		var ret: InnerSk = undefined;
450		ret.sh = V.fromBytes(buf).normalize();
451		return ret;
452		}
453		};
454
455		// Derives inner PKE keypair from given seed.
456		fn innerKeyFromSeed(seed: [inner_seed_length]u8, pk: InnerPk, sk: InnerSk) void {
457		var expanded_seed: [64]u8 = undefined;
458
459		var h = sha3.Sha3_512.init(.{});
460		h.update(&seed);
461		h.final(&expanded_seed);
462		pk.rho = expanded_seed[0..32].*;
463		const sigma = expanded_seed[32..64];
464		pk.aT = M.uniform(pk.rho, false); // Expand ρ to A; we'll transpose later on
465
466		// Sample secret vector s.
467		sk.sh = V.noise(p.eta1, 0, sigma).ntt().normalize();
468
469		const eh = Vec(p.k).noise(p.eta1, p.k, sigma).ntt(); // sample blind e.
470		var th: V = undefined;
471
472		// Next, we compute t = A s + e.
473		for (0..p.k) \|i\| {
474		// Note that coefficients of s are bounded by q and those of A
475		// are bounded by 4.5q and so their product is bounded by 2¹⁵q
476		// as required for multiplication.
477		// A and s were not in Montgomery form, so the Montgomery
478		// multiplications in the inner product added a factor R⁻¹ which
479		// we'll cancel out with toMont(). This will also ensure the
480		// coefficients of th are bounded in absolute value by q.
481		th.ps[i] = pk.aT.vs[i].dotHat(sk.sh).toMont();
482		}
483
484		pk.th = th.add(eh).normalize(); // bounded by 8q
485		pk.aT = pk.aT.transpose();
486		}
487		};
488		}
489
490		// R mod q
491		const r_mod_q: i32 = @rem(@as(i32, R), Q);
492
493		// R² mod q
494		const r2_mod_q: i32 = @rem(r_mod_q * r_mod_q, Q);
495
496		// ζ is the degree 256 primitive root of unity used for the NTT.
497		const zeta: i16 = 17;
498
499		// (128)⁻¹ R². Used in inverse NTT.
500		const r2_over_128: i32 = @mod(invertMod(128, Q) * r2_mod_q, Q);
501
502		// zetas lists precomputed powers of the primitive root of unity in
503		// Montgomery representation used for the NTT:
504		//
505		// zetas[i] = ζᵇʳᵛ⁽ⁱ⁾ R mod q
506		//
507		// where ζ = 17, brv(i) is the bitreversal of a 7-bit number and R=2¹⁶ mod q.
508		const zetas = computeZetas();
509
510		// invNTTReductions keeps track of which coefficients to apply Barrett
511		// reduction to in Poly.invNTT().
512		//
513		// Generated lazily: once a butterfly is computed which is about to
514		// overflow the i16, the largest coefficient is reduced. If that is
515		// not enough, the other coefficient is reduced as well.
516		//
517		// This is actually optimal, as proven in https://eprint.iacr.org/2020/1377.pdf
518		// TODO generate comptime?
519		const inv_ntt_reductions = [_]i16{
520		-1, // after layer 1
521		-1, // after layer 2
522		16,
523		17,
524		48,
525		49,
526		80,
527		81,
528		112,
529		113,
530		144,
531		145,
532		176,
533		177,
534		208,
535		209,
536		240, 241, -1, // after layer 3
537		0, 1, 32,
538		33, 34, 35,
539		64, 65, 96,
540		97, 98, 99,
541		128, 129,
542		160, 161, 162, 163, 192, 193, 224, 225, 226, 227, -1, // after layer 4
543		2, 3, 66, 67, 68, 69, 70, 71, 130, 131, 194,
544		195, 196, 197,
545		198, 199, -1, // after layer 5
546		4, 5, 6,
547		7, 132, 133,
548		134, 135, 136,
549		137, 138, 139,
550		140, 141,
551		142, 143, -1, // after layer 6
552		-1, // after layer 7
553		};
554
555		test "invNTTReductions bounds" {
556		// Checks whether the reductions proposed by invNTTReductions
557		// don't overflow during invNTT().
558		var xs = [_]i32{1} ** 256; // start at \|x\| ≤ q
559
560		var r: usize = 0;
561		var layer: math.Log2Int(usize) = 1;
562		while (layer < 8) : (layer += 1) {
563		const w = @as(usize, 1) << layer;
564		var i: usize = 0;
565
566		while (i + w < 256) {
567		xs[i] = xs[i] + xs[i + w];
568		try testing.expect(xs[i] <= 9); // we can't exceed 9q
569		xs[i + w] = 1;
570		i += 1;
571		if (@mod(i, w) == 0) {
572		i += w;
573		}
574		}
575
576		while (true) {
577		const j = inv_ntt_reductions[r];
578		r += 1;
579		if (j < 0) {
580		break;
581		}
582		xs[@as(usize, @intCast(j))] = 1;
583		}
584		}
585		}
586
587		// Extended euclidean algorithm.
588		//
589		// For a, b finds x, y such that x a + y b = gcd(a, b). Used to compute
590		// modular inverse.
591		fn eea(a: anytype, b: @TypeOf(a)) EeaResult(@TypeOf(a)) {
592		if (a == 0) {
593		return .{ .gcd = b, .x = 0, .y = 1 };
594		}
595		const r = eea(@rem(b, a), a);
596		return .{ .gcd = r.gcd, .x = r.y - @divTrunc(b, a) * r.x, .y = r.x };
597		}
598
599		fn EeaResult(comptime T: type) type {
600		return struct { gcd: T, x: T, y: T };
601		}
602
603		// Returns least common multiple of a and b.
604		fn lcm(a: anytype, b: @TypeOf(a)) @TypeOf(a) {
605		const r = eea(a, b);
606		return a * b / r.gcd;
607		}
608
609		// Invert modulo p.
610		fn invertMod(a: anytype, p: @TypeOf(a)) @TypeOf(a) {
611		const r = eea(a, p);
612		assert(r.gcd == 1);
613		return r.x;
614		}
615
616		// Reduce mod q for testing.
617		fn modQ32(x: i32) i16 {
618		var y = @as(i16, @intCast(@rem(x, @as(i32, Q))));
619		if (y < 0) {
620		y += Q;
621		}
622		return y;
623		}
624
625		// Given -2¹⁵ q ≤ x < 2¹⁵ q, returns -q < y < q with x 2⁻¹⁶ = y (mod q).
626		fn montReduce(x: i32) i16 {
627		const qInv = comptime invertMod(@as(i32, Q), R);
628		// This is Montgomery reduction with R=2¹⁶.
629		//
630		// Note gcd(2¹⁶, q) = 1 as q is prime. Write q' := 62209 = q⁻¹ mod R.
631		// First we compute
632		//
633		// m := ((x mod R) q') mod R
634		// = x q' mod R
635		// = int16(x q')
636		// = int16(int32(x) * int32(q'))
637		//
638		// Note that x q' might be as big as 2³² and could overflow the int32
639		// multiplication in the last line. However for any int32s a and b,
640		// we have int32(int64(a)int64(b)) = int32(ab) and so the result is ok.
641		const m: i16 = @truncate(@as(i32, @truncate(x *% qInv)));
642
643		// Note that x - m q is divisible by R; indeed modulo R we have
644		//
645		// x - m q ≡ x - x q' q ≡ x - x q⁻¹ q ≡ x - x = 0.
646		//
647		// We return y := (x - m q) / R. Note that y is indeed correct as
648		// modulo q we have
649		//
650		// y ≡ x R⁻¹ - m q R⁻¹ = x R⁻¹
651		//
652		// and as both 2¹⁵ q ≤ m q, x < 2¹⁵ q, we have
653		// 2¹⁶ q ≤ x - m q < 2¹⁶ and so q ≤ (x - m q) / R < q as desired.
654		const yR = x - @as(i32, m) * @as(i32, Q);
655		return @bitCast(@as(u16, @truncate(@as(u32, @bitCast(yR)) >> 16)));
656		}
657
658		test "Test montReduce" {
659		var rnd = RndGen.init(0);
660		for (0..1000) \|_\| {
661		const bound = comptime @as(i32, Q) * (1 << 15);
662		const x = rnd.random().intRangeLessThan(i32, -bound, bound);
663		const y = montReduce(x);
664		try testing.expect(-Q < y and y < Q);
665		try testing.expectEqual(modQ32(x), modQ32(@as(i32, y) * R));
666		}
667		}
668
669		// Given any x, return x R mod q where R=2¹⁶.
670		fn feToMont(x: i16) i16 {
671		// Note \|1353 x\| ≤ 1353 2¹⁵ ≤ 13318 q ≤ 2¹⁵ q and so we're within
672		// the bounds of montReduce.
673		return montReduce(@as(i32, x) * r2_mod_q);
674		}
675
676		test "Test feToMont" {
677		var x: i32 = -(1 << 15);
678		while (x < 1 << 15) : (x += 1) {
679		const y = feToMont(@as(i16, @intCast(x)));
680		try testing.expectEqual(modQ32(@as(i32, y)), modQ32(x * r_mod_q));
681		}
682		}
683
684		// Given any x, compute 0 ≤ y ≤ q with x = y (mod q).
685		//
686		// Beware: we might have feBarrettReduce(x) = q ≠ 0 for some x. In fact,
687		// this happens if and only if x = -nq for some positive integer n.
688		fn feBarrettReduce(x: i16) i16 {
689		// This is standard Barrett reduction.
690		//
691		// For any x we have x mod q = x - ⌊x/q⌋ q. We will use 20159/2²⁶ as
692		// an approximation of 1/q. Note that 0 ≤ 20159/2²⁶ - 1/q ≤ 0.135/2²⁶
693		// and so \| x 20156/2²⁶ - x/q \| ≤ 2⁻¹⁰ for \|x\| ≤ 2¹⁶. For all x
694		// not a multiple of q, the number x/q is further than 1/q from any integer
695		// and so ⌊x 20156/2²⁶⌋ = ⌊x/q⌋. If x is a multiple of q and x is positive,
696		// then x 20156/2²⁶ is larger than x/q so ⌊x 20156/2²⁶⌋ = ⌊x/q⌋ as well.
697		// Finally, if x is negative multiple of q, then ⌊x 20156/2²⁶⌋ = ⌊x/q⌋-1.
698		// Thus
699		// [ q if x=-nq for pos. integer n
700		// x - ⌊x 20156/2²⁶⌋ q = [
701		// [ x mod q otherwise
702		//
703		// To actually compute this, note that
704		//
705		// ⌊x 20156/2²⁶⌋ = (20159 x) >> 26.
706		return x -% @as(i16, @intCast((@as(i32, x) * 20159) >> 26)) *% Q;
707		}
708
709		test "Test Barrett reduction" {
710		var x: i32 = -(1 << 15);
711		while (x < 1 << 15) : (x += 1) {
712		var y1 = feBarrettReduce(@as(i16, @intCast(x)));
713		const y2 = @mod(@as(i16, @intCast(x)), Q);
714		if (x < 0 and @rem(-x, Q) == 0) {
715		y1 -= Q;
716		}
717		try testing.expectEqual(y1, y2);
718		}
719		}
720
721		// Returns x if x < q and x - q otherwise. Assumes x ≥ -29439.
722		fn csubq(x: i16) i16 {
723		var r = x;
724		r -= Q;
725		r += (r >> 15) & Q;
726		return r;
727		}
728
729		test "Test csubq" {
730		var x: i32 = -29439;
731		while (x < 1 << 15) : (x += 1) {
732		const y1 = csubq(@as(i16, @intCast(x)));
733		var y2 = @as(i16, @intCast(x));
734		if (@as(i16, @intCast(x)) >= Q) {
735		y2 -= Q;
736		}
737		try testing.expectEqual(y1, y2);
738		}
739		}
740
741		// Compute a^s mod p.
742		fn mpow(a: anytype, s: @TypeOf(a), p: @TypeOf(a)) @TypeOf(a) {
743		var ret: @TypeOf(a) = 1;
744		var s2 = s;
745		var a2 = a;
746
747		while (true) {
748		if (s2 & 1 == 1) {
749		ret = @mod(ret * a2, p);
750		}
751		s2 >>= 1;
752		if (s2 == 0) {
753		break;
754		}
755		a2 = @mod(a2 * a2, p);
756		}
757		return ret;
758		}
759
760		// Computes zetas table used by ntt and invNTT.
761		fn computeZetas() [128]i16 {
762		@setEvalBranchQuota(10000);
763		var ret: [128]i16 = undefined;
764		for (&ret, 0..) \|*r, i\| {
765		const t = @as(i16, @intCast(mpow(@as(i32, zeta), @bitReverse(@as(u7, @intCast(i))), Q)));
766		r.* = csubq(feBarrettReduce(feToMont(t)));
767		}
768		return ret;
769		}
770
771		// An element of our base ring R which are polynomials over ℤ_q
772		// modulo the equation Xᴺ = -1, where q=3329 and N=256.
773		//
774		// This type is also used to store NTT-transformed polynomials,
775		// see Poly.NTT().
776		//
777		// Coefficients aren't always reduced. See Normalize().
778		const Poly = struct {
779		cs: [N]i16,
780
781		const bytes_length = N / 2 * 3;
782		const zero: Poly = .{ .cs = .{0} ** N };
783
784		fn add(a: Poly, b: Poly) Poly {
785		var ret: Poly = undefined;
786		for (0..N) \|i\| {
787		ret.cs[i] = a.cs[i] + b.cs[i];
788		}
789		return ret;
790		}
791
792		fn sub(a: Poly, b: Poly) Poly {
793		var ret: Poly = undefined;
794		for (0..N) \|i\| {
795		ret.cs[i] = a.cs[i] - b.cs[i];
796		}
797		return ret;
798		}
799
800		// For testing, generates a random polynomial with for each
801		// coefficient \|x\| ≤ q.
802		fn randAbsLeqQ(rnd: anytype) Poly {
803		var ret: Poly = undefined;
804		for (0..N) \|i\| {
805		ret.cs[i] = rnd.random().intRangeAtMost(i16, -Q, Q);
806		}
807		return ret;
808		}
809
810		// For testing, generates a random normalized polynomial.
811		fn randNormalized(rnd: anytype) Poly {
812		var ret: Poly = undefined;
813		for (0..N) \|i\| {
814		ret.cs[i] = rnd.random().intRangeLessThan(i16, 0, Q);
815		}
816		return ret;
817		}
818
819		// Executes a forward "NTT" on p.
820		//
821		// Assumes the coefficients are in absolute value ≤q. The resulting
822		// coefficients are in absolute value ≤7q. If the input is in Montgomery
823		// form, then the result is in Montgomery form and so (by linearity of the NTT)
824		// if the input is in regular form, then the result is also in regular form.
825		fn ntt(a: Poly) Poly {
826		// Note that ℤ_q does not have a primitive 512ᵗʰ root of unity (as 512
827		// does not divide into q-1) and so we cannot do a regular NTT. ℤ_q
828		// does have a primitive 256ᵗʰ root of unity, the smallest of which
829		// is ζ := 17.
830		//
831		// Recall that our base ring R := ℤ_q[x] / (x²⁵⁶ + 1). The polynomial
832		// x²⁵⁶+1 will not split completely (as its roots would be 512ᵗʰ roots
833		// of unity.) However, it does split almost (using ζ¹²⁸ = -1):
834		//
835		// x²⁵⁶ + 1 = (x²)¹²⁸ - ζ¹²⁸
836		// = ((x²)⁶⁴ - ζ⁶⁴)((x²)⁶⁴ + ζ⁶⁴)
837		// = ((x²)³² - ζ³²)((x²)³² + ζ³²)((x²)³² - ζ⁹⁶)((x²)³² + ζ⁹⁶)
838		// ⋮
839		// = (x² - ζ)(x² + ζ)(x² - ζ⁶⁵)(x² + ζ⁶⁵) … (x² + ζ¹²⁷)
840		//
841		// Note that the powers of ζ that appear (from the second line down) are
842		// in binary
843		//
844		// 0100000 1100000
845		// 0010000 1010000 0110000 1110000
846		// 0001000 1001000 0101000 1101000 0011000 1011000 0111000 1111000
847		// …
848		//
849		// That is: brv(2), brv(3), brv(4), …, where brv(x) denotes the 7-bit
850		// bitreversal of x. These powers of ζ are given by the Zetas array.
851		//
852		// The polynomials x² ± ζⁱ are irreducible and coprime, hence by
853		// the Chinese Remainder Theorem we know
854		//
855		// ℤ_q[x]/(x²⁵⁶+1) → ℤ_q[x]/(x²-ζ) x … x ℤ_q[x]/(x²+ζ¹²⁷)
856		//
857		// given by a ↦ ( a mod x²-ζ, …, a mod x²+ζ¹²⁷ )
858		// is an isomorphism, which is the "NTT". It can be efficiently computed by
859		//
860		//
861		// a ↦ ( a mod (x²)⁶⁴ - ζ⁶⁴, a mod (x²)⁶⁴ + ζ⁶⁴ )
862		// ↦ ( a mod (x²)³² - ζ³², a mod (x²)³² + ζ³²,
863		// a mod (x²)⁹⁶ - ζ⁹⁶, a mod (x²)⁹⁶ + ζ⁹⁶ )
864		//
865		// et cetera
866		// If N was 8 then this can be pictured in the following diagram:
867		//
868		// https://cnx.org/resources/17ee4dfe517a6adda05377b25a00bf6e6c93c334/File0026.png
869		//
870		// Each cross is a Cooley-Tukey butterfly: it's the map
871		//
872		// (a, b) ↦ (a + ζb, a - ζb)
873		//
874		// for the appropriate power ζ for that column and row group.
875		var p = a;
876		var k: usize = 0; // index into zetas
877
878		var l = N >> 1;
879		while (l > 1) : (l >>= 1) {
880		// On the nᵗʰ iteration of the l-loop, the absolute value of the
881		// coefficients are bounded by nq.
882
883		// offset effectively loops over the row groups in this column; it is
884		// the first row in the row group.
885		var offset: usize = 0;
886		while (offset < N - l) : (offset += 2 * l) {
887		k += 1;
888		const z = @as(i32, zetas[k]);
889
890		// j loops over each butterfly in the row group.
891		for (offset..offset + l) \|j\| {
892		const t = montReduce(z * @as(i32, p.cs[j + l]));
893		p.cs[j + l] = p.cs[j] - t;
894		p.cs[j] += t;
895		}
896		}
897		}
898
899		return p;
900		}
901
902		// Executes an inverse "NTT" on p and multiply by the Montgomery factor R.
903		//
904		// Assumes the coefficients are in absolute value ≤q. The resulting
905		// coefficients are in absolute value ≤q. If the input is in Montgomery
906		// form, then the result is in Montgomery form and so (by linearity)
907		// if the input is in regular form, then the result is also in regular form.
908		fn invNTT(a: Poly) Poly {
909		var k: usize = 127; // index into zetas
910		var r: usize = 0; // index into invNTTReductions
911		var p = a;
912
913		// We basically do the oppposite of NTT, but postpone dividing by 2 in the
914		// inverse of the Cooley-Tukey butterfly and accumulate that into a big
915		// division by 2⁷ at the end. See the comments in the ntt() function.
916
917		var l: usize = 2;
918		while (l < N) : (l <<= 1) {
919		var offset: usize = 0;
920		while (offset < N - l) : (offset += 2 * l) {
921		// As we're inverting, we need powers of ζ⁻¹ (instead of ζ).
922		// To be precise, we need ζᵇʳᵛ⁽ᵏ⁾⁻¹²⁸. However, as ζ⁻¹²⁸ = -1,
923		// we can use the existing zetas table instead of
924		// keeping a separate invZetas table as in Dilithium.
925
926		const minZeta = @as(i32, zetas[k]);
927		k -= 1;
928
929		for (offset..offset + l) \|j\| {
930		// Gentleman-Sande butterfly: (a, b) ↦ (a + b, ζ(a-b))
931		const t = p.cs[j + l] - p.cs[j];
932		p.cs[j] += p.cs[j + l];
933		p.cs[j + l] = montReduce(minZeta * @as(i32, t));
934
935		// Note that if we had \|a\| < αq and \|b\| < βq before the
936		// butterfly, then now we have \|a\| < (α+β)q and \|b\| < q.
937		}
938		}
939
940		// We let the invNTTReductions instruct us which coefficients to
941		// Barrett reduce.
942		while (true) {
943		const i = inv_ntt_reductions[r];
944		r += 1;
945		if (i < 0) {
946		break;
947		}
948		p.cs[@as(usize, @intCast(i))] = feBarrettReduce(p.cs[@as(usize, @intCast(i))]);
949		}
950		}
951
952		for (0..N) \|j\| {
953		// Note 1441 = (128)⁻¹ R². The coefficients are bounded by 9q, so
954		// as 1441 * 9 ≈ 2¹⁴ < 2¹⁵, we're within the required bounds
955		// for montReduce().
956		p.cs[j] = montReduce(r2_over_128 * @as(i32, p.cs[j]));
957		}
958
959		return p;
960		}
961
962		// Normalizes coefficients.
963		//
964		// Ensures each coefficient is in {0, …, q-1}.
965		fn normalize(a: Poly) Poly {
966		var ret: Poly = undefined;
967		for (0..N) \|i\| {
968		ret.cs[i] = csubq(feBarrettReduce(a.cs[i]));
969		}
970		return ret;
971		}
972
973		// Put p in Montgomery form.
974		fn toMont(a: Poly) Poly {
975		var ret: Poly = undefined;
976		for (0..N) \|i\| {
977		ret.cs[i] = feToMont(a.cs[i]);
978		}
979		return ret;
980		}
981
982		// Barret reduce coefficients.
983		//
984		// Beware, this does not fully normalize coefficients.
985		fn barrettReduce(a: Poly) Poly {
986		var ret: Poly = undefined;
987		for (0..N) \|i\| {
988		ret.cs[i] = feBarrettReduce(a.cs[i]);
989		}
990		return ret;
991		}
992
993		fn compressedSize(comptime d: u8) usize {
994		return @divTrunc(N * d, 8);
995		}
996
997		// Returns packed Compress_q(p, d).
998		//
999		// Assumes p is normalized.
1000		fn compress(p: Poly, comptime d: u8) [compressedSize(d)]u8 {
1001		@setEvalBranchQuota(10000);
1002		const q_over_2: u32 = comptime @divTrunc(Q, 2); // (q-1)/2
1003		const two_d_min_1: u32 = comptime (1 << d) - 1; // 2ᵈ-1
1004		var in_off: usize = 0;
1005		var out_off: usize = 0;
1006
1007		const batch_size: usize = comptime lcm(@as(i16, d), 8);
1008		const in_batch_size: usize = comptime batch_size / d;
1009		const out_batch_size: usize = comptime batch_size / 8;
1010
1011		const out_length: usize = comptime @divTrunc(N * d, 8);
1012		comptime assert(out_length * 8 == d * N);
1013		var out = [_]u8{0} ** out_length;
1014
1015		while (in_off < N) {
1016		// First we compress into in.
1017		var in: [in_batch_size]u16 = undefined;
1018		inline for (0..in_batch_size) \|i\| {
1019		// Compress_q(x, d) = ⌈(2ᵈ/q)x⌋ mod⁺ 2ᵈ
1020		// = ⌊(2ᵈ/q)x+½⌋ mod⁺ 2ᵈ
1021		// = ⌊((x << d) + q/2) / q⌋ mod⁺ 2ᵈ
1022		// = DIV((x << d) + q/2, q) & ((1<<d) - 1)
1023		const t = @as(u24, @intCast(p.cs[in_off + i])) << d;
1024		// Division by invariant multiplication, equivalent to DIV(t + q/2, q).
1025		// A division may not be a constant-time operation, even with a constant denominator.
1026		// Here, side channels would leak information about the shared secret, see https://kyberslash.cr.yp.to
1027		// Multiplication, on the other hand, is a constant-time operation on the CPUs we currently support.
1028		comptime assert(d <= 11);
1029		comptime assert(((20642679 * @as(u64, Q)) >> 36) == 1);
1030		const u: u32 = @intCast((@as(u64, t + q_over_2) * 20642679) >> 36);
1031		in[i] = @intCast(u & two_d_min_1);
1032		}
1033
1034		// Now we pack the d-bit integers from `in' into out as bytes.
1035		comptime var in_shift: usize = 0;
1036		comptime var j: usize = 0;
1037		comptime var i: usize = 0;
1038		inline while (i < in_batch_size) : (j += 1) {
1039		comptime var todo: usize = 8;
1040		inline while (todo > 0) {
1041		const out_shift = comptime 8 - todo;
1042		out[out_off + j] \|= @as(u8, @truncate((in[i] >> in_shift) << out_shift));
1043
1044		const done = comptime @min(@min(d, todo), d - in_shift);
1045		todo -= done;
1046		in_shift += done;
1047
1048		if (in_shift == d) {
1049		in_shift = 0;
1050		i += 1;
1051		}
1052		}
1053		}
1054
1055		in_off += in_batch_size;
1056		out_off += out_batch_size;
1057		}
1058
1059		return out;
1060		}
1061
1062		// Set p to Decompress_q(m, d).
1063		fn decompress(comptime d: u8, in: *const [compressedSize(d)]u8) Poly {
1064		@setEvalBranchQuota(10000);
1065		const inLen = comptime @divTrunc(N * d, 8);
1066		comptime assert(inLen * 8 == d * N);
1067		var ret: Poly = undefined;
1068		var in_off: usize = 0;
1069		var out_off: usize = 0;
1070
1071		const batch_size: usize = comptime lcm(@as(i16, d), 8);
1072		const in_batch_size: usize = comptime batch_size / 8;
1073		const out_batch_size: usize = comptime batch_size / d;
1074
1075		while (out_off < N) {
1076		comptime var in_shift: usize = 0;
1077		comptime var j: usize = 0;
1078		comptime var i: usize = 0;
1079		inline while (i < out_batch_size) : (i += 1) {
1080		// First, unpack next coefficient.
1081		comptime var todo = d;
1082		var out: u16 = 0;
1083
1084		inline while (todo > 0) {
1085		const out_shift = comptime d - todo;
1086		const m = comptime (1 << d) - 1;
1087		out \|= (@as(u16, in[in_off + j] >> in_shift) << out_shift) & m;
1088
1089		const done = comptime @min(@min(8, todo), 8 - in_shift);
1090		todo -= done;
1091		in_shift += done;
1092
1093		if (in_shift == 8) {
1094		in_shift = 0;
1095		j += 1;
1096		}
1097		}
1098
1099		// Decompress_q(x, d) = ⌈(q/2ᵈ)x⌋
1100		// = ⌊(q/2ᵈ)x+½⌋
1101		// = ⌊(qx + 2ᵈ⁻¹)/2ᵈ⌋
1102		// = (qx + (1<<(d-1))) >> d
1103		const qx = @as(u32, out) * @as(u32, Q);
1104		ret.cs[out_off + i] = @as(i16, @intCast((qx + (1 << (d - 1))) >> d));
1105		}
1106
1107		in_off += in_batch_size;
1108		out_off += out_batch_size;
1109		}
1110
1111		return ret;
1112		}
1113
1114		// Returns the "pointwise" multiplication a o b.
1115		//
1116		// That is: invNTT(a o b) = invNTT(a) * invNTT(b). Assumes a and b are in
1117		// Montgomery form. Products between coefficients of a and b must be strictly
1118		// bounded in absolute value by 2¹⁵q. a o b will be in Montgomery form and
1119		// bounded in absolute value by 2q.
1120		fn mulHat(a: Poly, b: Poly) Poly {
1121		// Recall from the discussion in ntt(), that a transformed polynomial is
1122		// an element of ℤ_q[x]/(x²-ζ) x … x ℤ_q[x]/(x²+ζ¹²⁷);
1123		// that is: 128 degree-one polynomials instead of simply 256 elements
1124		// from ℤ_q as in the regular NTT. So instead of pointwise multiplication,
1125		// we multiply the 128 pairs of degree-one polynomials modulo the
1126		// right equation:
1127		//
1128		// (a₁ + a₂x)(b₁ + b₂x) = a₁b₁ + a₂b₂ζ' + (a₁b₂ + a₂b₁)x,
1129		//
1130		// where ζ' is the appropriate power of ζ.
1131
1132		var p: Poly = undefined;
1133		var k: usize = 64;
1134		var i: usize = 0;
1135		while (i < N) : (i += 4) {
1136		const z = @as(i32, zetas[k]);
1137		k += 1;
1138
1139		const a1b1 = montReduce(@as(i32, a.cs[i + 1]) * @as(i32, b.cs[i + 1]));
1140		const a0b0 = montReduce(@as(i32, a.cs[i]) * @as(i32, b.cs[i]));
1141		const a1b0 = montReduce(@as(i32, a.cs[i + 1]) * @as(i32, b.cs[i]));
1142		const a0b1 = montReduce(@as(i32, a.cs[i]) * @as(i32, b.cs[i + 1]));
1143
1144		p.cs[i] = montReduce(a1b1 * z) + a0b0;
1145		p.cs[i + 1] = a0b1 + a1b0;
1146
1147		const a3b3 = montReduce(@as(i32, a.cs[i + 3]) * @as(i32, b.cs[i + 3]));
1148		const a2b2 = montReduce(@as(i32, a.cs[i + 2]) * @as(i32, b.cs[i + 2]));
1149		const a3b2 = montReduce(@as(i32, a.cs[i + 3]) * @as(i32, b.cs[i + 2]));
1150		const a2b3 = montReduce(@as(i32, a.cs[i + 2]) * @as(i32, b.cs[i + 3]));
1151
1152		p.cs[i + 2] = a2b2 - montReduce(a3b3 * z);
1153		p.cs[i + 3] = a2b3 + a3b2;
1154		}
1155
1156		return p;
1157		}
1158
1159		// Sample p from a centered binomial distribution with n=2η and p=½ - viz:
1160		// coefficients are in {-η, …, η} with probabilities
1161		//
1162		// {ncr(0, 2η)/2^2η, ncr(1, 2η)/2^2η, …, ncr(2η,2η)/2^2η}
1163		fn noise(comptime eta: u8, nonce: u8, seed: *const [32]u8) Poly {
1164		var h = sha3.Shake256.init(.{});
1165		const suffix: [1]u8 = .{nonce};
1166		h.update(seed);
1167		h.update(&suffix);
1168
1169		// The distribution at hand is exactly the same as that
1170		// of (a₁ + a₂ + … + a_η) - (b₁ + … + b_η) where a_i,b_i~U(1).
1171		// Thus we need 2η bits per coefficient.
1172		const buf_len = comptime 2 * eta * N / 8;
1173		var buf: [buf_len]u8 = undefined;
1174		h.squeeze(&buf);
1175
1176		// buf is interpreted as a₁…a_ηb₁…b_ηa₁…a_ηb₁…b_η…. We process
1177		// multiple coefficients in one batch.
1178
1179		const T = switch (builtin.target.cpu.arch) {
1180		.x86_64, .x86 => u32, // Generates better code on Intel CPUs
1181		else => u64, // u128 might be faster on some other CPUs.
1182		};
1183
1184		comptime var batch_count: usize = undefined;
1185		comptime var batch_bytes: usize = undefined;
1186		comptime var mask: T = 0;
1187		comptime {
1188		batch_count = @bitSizeOf(T) / @as(usize, 2 * eta);
1189		while (@rem(N, batch_count) != 0 and batch_count > 0) : (batch_count -= 1) {}
1190		assert(batch_count > 0);
1191		assert(@rem(2 * eta * batch_count, 8) == 0);
1192		batch_bytes = 2 * eta * batch_count / 8;
1193
1194		for (0..2 * eta * batch_count) \|_\| {
1195		mask <<= eta;
1196		mask \|= 1;
1197		}
1198		}
1199
1200		var ret: Poly = undefined;
1201		for (0..comptime N / batch_count) \|i\| {
1202		// Read coefficients into t. In the case of η=3,
1203		// we have t = a₁ + 2a₂ + 4a₃ + 8b₁ + 16b₂ + …
1204		var t: T = 0;
1205		inline for (0..batch_bytes) \|j\| {
1206		t \|= @as(T, buf[batch_bytes * i + j]) << (8 * j);
1207		}
1208
1209		// Accumelate `a's and `b's together by masking them out, shifting
1210		// and adding. For η=3, we have d = a₁ + a₂ + a₃ + 8(b₁ + b₂ + b₃) + …
1211		var d: T = 0;
1212		inline for (0..eta) \|j\| {
1213		d += (t >> j) & mask;
1214		}
1215
1216		// Extract each a and b separately and set coefficient in polynomial.
1217		inline for (0..batch_count) \|j\| {
1218		const mask2 = comptime (1 << eta) - 1;
1219		const a = @as(i16, @intCast((d >> (comptime (2 * j * eta))) & mask2));
1220		const b = @as(i16, @intCast((d >> (comptime ((2 * j + 1) * eta))) & mask2));
1221		ret.cs[batch_count * i + j] = a - b;
1222		}
1223		}
1224
1225		return ret;
1226		}
1227
1228		// Sample p uniformly from the given seed and x and y coordinates.
1229		fn uniform(seed: [32]u8, x: u8, y: u8) Poly {
1230		var h = sha3.Shake128.init(.{});
1231		const suffix: [2]u8 = .{ x, y };
1232		h.update(&seed);
1233		h.update(&suffix);
1234
1235		const buf_len = sha3.Shake128.block_length; // rate SHAKE-128
1236		var buf: [buf_len]u8 = undefined;
1237
1238		var ret: Poly = undefined;
1239		var i: usize = 0; // index into ret.cs
1240		outer: while (true) {
1241		h.squeeze(&buf);
1242
1243		var j: usize = 0; // index into buf
1244		while (j < buf_len) : (j += 3) {
1245		const b0 = @as(u16, buf[j]);
1246		const b1 = @as(u16, buf[j + 1]);
1247		const b2 = @as(u16, buf[j + 2]);
1248
1249		const ts: [2]u16 = .{
1250		b0 \| ((b1 & 0xf) << 8),
1251		(b1 >> 4) \| (b2 << 4),
1252		};
1253
1254		inline for (ts) \|t\| {
1255		if (t < Q) {
1256		ret.cs[i] = @as(i16, @intCast(t));
1257		i += 1;
1258
1259		if (i == N) {
1260		break :outer;
1261		}
1262		}
1263		}
1264		}
1265		}
1266
1267		return ret;
1268		}
1269
1270		// Packs p.
1271		//
1272		// Assumes p is normalized (and not just Barrett reduced).
1273		fn toBytes(p: Poly) [bytes_length]u8 {
1274		var ret: [bytes_length]u8 = undefined;
1275		for (0..comptime N / 2) \|i\| {
1276		const t0 = @as(u16, @intCast(p.cs[2 * i]));
1277		const t1 = @as(u16, @intCast(p.cs[2 * i + 1]));
1278		ret[3 * i] = @as(u8, @truncate(t0));
1279		ret[3 * i + 1] = @as(u8, @truncate((t0 >> 8) \| (t1 << 4)));
1280		ret[3 * i + 2] = @as(u8, @truncate(t1 >> 4));
1281		}
1282		return ret;
1283		}
1284
1285		// Unpacks a Poly from buf.
1286		//
1287		// p will not be normalized; instead 0 ≤ p[i] < 4096.
1288		fn fromBytes(buf: *const [bytes_length]u8) Poly {
1289		var ret: Poly = undefined;
1290		for (0..comptime N / 2) \|i\| {
1291		const b0 = @as(i16, buf[3 * i]);
1292		const b1 = @as(i16, buf[3 * i + 1]);
1293		const b2 = @as(i16, buf[3 * i + 2]);
1294		ret.cs[2 * i] = b0 \| ((b1 & 0xf) << 8);
1295		ret.cs[2 * i + 1] = (b1 >> 4) \| b2 << 4;
1296		}
1297		return ret;
1298		}
1299		};
1300
1301		// A vector of K polynomials.
1302		fn Vec(comptime K: u8) type {
1303		return struct {
1304		ps: [K]Poly,
1305
1306		const Self = @This();
1307		const bytes_length = K * Poly.bytes_length;
1308
1309		fn compressedSize(comptime d: u8) usize {
1310		return Poly.compressedSize(d) * K;
1311		}
1312
1313		fn ntt(a: Self) Self {
1314		var ret: Self = undefined;
1315		for (0..K) \|i\| {
1316		ret.ps[i] = a.ps[i].ntt();
1317		}
1318		return ret;
1319		}
1320
1321		fn invNTT(a: Self) Self {
1322		var ret: Self = undefined;
1323		for (0..K) \|i\| {
1324		ret.ps[i] = a.ps[i].invNTT();
1325		}
1326		return ret;
1327		}
1328
1329		fn normalize(a: Self) Self {
1330		var ret: Self = undefined;
1331		for (0..K) \|i\| {
1332		ret.ps[i] = a.ps[i].normalize();
1333		}
1334		return ret;
1335		}
1336
1337		fn barrettReduce(a: Self) Self {
1338		var ret: Self = undefined;
1339		for (0..K) \|i\| {
1340		ret.ps[i] = a.ps[i].barrettReduce();
1341		}
1342		return ret;
1343		}
1344
1345		fn add(a: Self, b: Self) Self {
1346		var ret: Self = undefined;
1347		for (0..K) \|i\| {
1348		ret.ps[i] = a.ps[i].add(b.ps[i]);
1349		}
1350		return ret;
1351		}
1352
1353		fn sub(a: Self, b: Self) Self {
1354		var ret: Self = undefined;
1355		for (0..K) \|i\| {
1356		ret.ps[i] = a.ps[i].sub(b.ps[i]);
1357		}
1358		return ret;
1359		}
1360
1361		// Samples v[i] from centered binomial distribution with the given η,
1362		// seed and nonce+i.
1363		fn noise(comptime eta: u8, nonce: u8, seed: *const [32]u8) Self {
1364		var ret: Self = undefined;
1365		for (0..K) \|i\| {
1366		ret.ps[i] = Poly.noise(eta, nonce + @as(u8, @intCast(i)), seed);
1367		}
1368		return ret;
1369		}
1370
1371		// Sets p to the inner product of a and b using "pointwise" multiplication.
1372		//
1373		// See MulHat() and NTT() for a description of the multiplication.
1374		// Assumes a and b are in Montgomery form. p will be in Montgomery form,
1375		// and its coefficients will be bounded in absolute value by 2kq.
1376		// If a and b are not in Montgomery form, then the action is the same
1377		// as "pointwise" multiplication followed by multiplying by R⁻¹, the inverse
1378		// of the Montgomery factor.
1379		fn dotHat(a: Self, b: Self) Poly {
1380		var ret: Poly = Poly.zero;
1381		for (0..K) \|i\| {
1382		ret = ret.add(a.ps[i].mulHat(b.ps[i]));
1383		}
1384		return ret;
1385		}
1386
1387		fn compress(v: Self, comptime d: u8) [compressedSize(d)]u8 {
1388		const cs = comptime Poly.compressedSize(d);
1389		var ret: [compressedSize(d)]u8 = undefined;
1390		inline for (0..K) \|i\| {
1391		ret[i * cs .. (i + 1) * cs].* = v.ps[i].compress(d);
1392		}
1393		return ret;
1394		}
1395
1396		fn decompress(comptime d: u8, buf: *const [compressedSize(d)]u8) Self {
1397		const cs = comptime Poly.compressedSize(d);
1398		var ret: Self = undefined;
1399		inline for (0..K) \|i\| {
1400		ret.ps[i] = Poly.decompress(d, buf[i * cs .. (i + 1) * cs]);
1401		}
1402		return ret;
1403		}
1404
1405		/// Serializes the key into a byte array.
1406		fn toBytes(v: Self) [bytes_length]u8 {
1407		var ret: [bytes_length]u8 = undefined;
1408		inline for (0..K) \|i\| {
1409		ret[i * Poly.bytes_length .. (i + 1) * Poly.bytes_length].* = v.ps[i].toBytes();
1410		}
1411		return ret;
1412		}
1413
1414		/// Deserializes the key from a byte array.
1415		fn fromBytes(buf: *const [bytes_length]u8) Self {
1416		var ret: Self = undefined;
1417		inline for (0..K) \|i\| {
1418		ret.ps[i] = Poly.fromBytes(
1419		buf[i * Poly.bytes_length .. (i + 1) * Poly.bytes_length],
1420		);
1421		}
1422		return ret;
1423		}
1424		};
1425		}
1426
1427		// A matrix of K vectors
1428		fn Mat(comptime K: u8) type {
1429		return struct {
1430		const Self = @This();
1431		vs: [K]Vec(K),
1432
1433		fn uniform(seed: [32]u8, comptime transposed: bool) Self {
1434		var ret: Self = undefined;
1435		var i: u8 = 0;
1436		while (i < K) : (i += 1) {
1437		var j: u8 = 0;
1438		while (j < K) : (j += 1) {
1439		ret.vs[i].ps[j] = Poly.uniform(
1440		seed,
1441		if (transposed) i else j,
1442		if (transposed) j else i,
1443		);
1444		}
1445		}
1446		return ret;
1447		}
1448
1449		// Returns transpose of A
1450		fn transpose(m: Self) Self {
1451		var ret: Self = undefined;
1452		for (0..K) \|i\| {
1453		for (0..K) \|j\| {
1454		ret.vs[i].ps[j] = m.vs[j].ps[i];
1455		}
1456		}
1457		return ret;
1458		}
1459		};
1460		}
1461
1462		// Returns `true` if a ≠ b.
1463		fn ctneq(comptime len: usize, a: [len]u8, b: [len]u8) u1 {
1464		return 1 - @intFromBool(crypto.utils.timingSafeEql([len]u8, a, b));
1465		}
1466
1467		// Copy src into dst given b = 1.
1468		fn cmov(comptime len: usize, dst: *[len]u8, src: [len]u8, b: u1) void {
1469		const mask = @as(u8, 0) -% b;
1470		for (0..len) \|i\| {
1471		dst[i] ^= mask & (dst[i] ^ src[i]);
1472		}
1473		}
1474
1475		test "MulHat" {
1476		var rnd = RndGen.init(0);
1477
1478		for (0..100) \|_\| {
1479		const a = Poly.randAbsLeqQ(&rnd);
1480		const b = Poly.randAbsLeqQ(&rnd);
1481
1482		const p2 = a.ntt().mulHat(b.ntt()).barrettReduce().invNTT().normalize();
1483		var p: Poly = undefined;
1484
1485		@memset(&p.cs, 0);
1486
1487		for (0..N) \|i\| {
1488		for (0..N) \|j\| {
1489		var v = montReduce(@as(i32, a.cs[i]) * @as(i32, b.cs[j]));
1490		var k = i + j;
1491		if (k >= N) {
1492		// Recall Xᴺ = -1.
1493		k -= N;
1494		v = -v;
1495		}
1496		p.cs[k] = feBarrettReduce(v + p.cs[k]);
1497		}
1498		}
1499
1500		p = p.toMont().normalize();
1501
1502		try testing.expectEqual(p, p2);
1503		}
1504		}
1505
1506		test "NTT" {
1507		var rnd = RndGen.init(0);
1508
1509		for (0..1000) \|_\| {
1510		var p = Poly.randAbsLeqQ(&rnd);
1511		const q = p.toMont().normalize();
1512		p = p.ntt();
1513
1514		for (0..N) \|i\| {
1515		try testing.expect(p.cs[i] <= 7 * Q and -7 * Q <= p.cs[i]);
1516		}
1517
1518		p = p.normalize().invNTT();
1519		for (0..N) \|i\| {
1520		try testing.expect(p.cs[i] <= Q and -Q <= p.cs[i]);
1521		}
1522
1523		p = p.normalize();
1524
1525		try testing.expectEqual(p, q);
1526		}
1527		}
1528
1529		test "Compression" {
1530		var rnd = RndGen.init(0);
1531		inline for (.{ 1, 4, 5, 10, 11 }) \|d\| {
1532		for (0..1000) \|_\| {
1533		const p = Poly.randNormalized(&rnd);
1534		const pp = p.compress(d);
1535		const pq = Poly.decompress(d, &pp).compress(d);
1536		try testing.expectEqual(pp, pq);
1537		}
1538		}
1539		}
1540
1541		test "noise" {
1542		var seed: [32]u8 = undefined;
1543		for (&seed, 0..) \|*s, i\| {
1544		s.* = @as(u8, @intCast(i));
1545		}
1546		try testing.expectEqual(Poly.noise(3, 37, &seed).cs, .{
1547		0, 0, 1, -1, 0, 2, 0, -1, -1, 3, 0, 1, -2, -2, 0, 1, -2,
1548		1, 0, -2, 3, 0, 0, 0, 1, 3, 1, 1, 2, 1, -1, -1, -1, 0,
1549		1, 0, 1, 0, 2, 0, 1, -2, 0, -1, -1, -2, 1, -1, -1, 2, -1,
1550		1, 1, 2, -3, -1, -1, 0, 0, 0, 0, 1, -1, -2, -2, 0, -2, 0,
1551		0, 0, 1, 0, -1, -1, 1, -2, 2, 0, 0, 2, -2, 0, 1, 0, 1,
1552		1, 1, 0, 1, -2, -1, -2, -1, 1, 0, 0, 0, 0, 0, 1, 0, -1,
1553		-1, 0, -1, 1, 0, 1, 0, -1, -1, 0, -2, 2, 0, -2, 1, -1, 0,
1554		1, -1, -1, 2, 1, 0, 0, -2, -1, 2, 0, 0, 0, -1, -1, 3, 1,
1555		0, 1, 0, 1, 0, 2, 1, 0, 0, 1, 0, 1, 0, 0, -1, -1, -1,
1556		0, 1, 3, 1, 0, 1, 0, 1, -1, -1, -1, -1, 0, 0, -2, -1, -1,
1557		2, 0, 1, 0, 1, 0, 2, -2, 0, 1, 1, -3, -1, -2, -1, 0, 1,
1558		0, 1, -2, 2, 2, 1, 1, 0, -1, 0, -1, -1, 1, 0, -1, 2, 1,
1559		-1, 1, 2, -2, 1, 2, 0, 1, 2, 1, 0, 0, 2, 1, 2, 1, 0,
1560		2, 1, 0, 0, -1, -1, 1, -1, 0, 1, -1, 2, 2, 0, 0, -1, 1,
1561		1, 1, 1, 0, 0, -2, 0, -1, 1, 2, 0, 0, 1, 1, -1, 1, 0,
1562		1,
1563		});
1564		try testing.expectEqual(Poly.noise(2, 37, &seed).cs, .{
1565		1, 0, 1, -1, -1, -2, -1, -1, 2, 0, -1, 0, 0, -1,
1566		1, 1, -1, 1, 0, 2, -2, 0, 1, 2, 0, 0, -1, 1,
1567		0, -1, 1, -1, 1, 2, 1, 1, 0, -1, 1, -1, -2, -1,
1568		1, -1, -1, -1, 2, -1, -1, 0, 0, 1, 1, -1, 1, 1,
1569		1, 1, -1, -2, 0, 1, 0, 0, 2, 1, -1, 2, 0, 0,
1570		1, 1, 0, -1, 0, 0, -1, -1, 2, 0, 1, -1, 2, -1,
1571		-1, -1, -1, 0, -2, 0, 2, 1, 0, 0, 0, -1, 0, 0,
1572		0, -1, -1, 0, -1, -1, 0, -1, 0, 0, -2, 1, 1, 0,
1573		1, 0, 1, 0, 1, 1, -1, 2, 0, 1, -1, 1, 2, 0,
1574		0, 0, 0, -1, -1, -1, 0, 1, 0, -1, 2, 0, 0, 1,
1575		1, 1, 0, 1, -1, 1, 2, 1, 0, 2, -1, 1, -1, -2,
1576		-1, -2, -1, 1, 0, -2, -2, -1, 1, 0, 0, 0, 0, 1,
1577		0, 0, 0, 2, 2, 0, 1, 0, -1, -1, 0, 2, 0, 0,
1578		-2, 1, 0, 2, 1, -1, -2, 0, 0, -1, 1, 1, 0, 0,
1579		2, 0, 1, 1, -2, 1, -2, 1, 1, 0, 2, 0, -1, 0,
1580		-1, 0, 1, 2, 0, 1, 0, -2, 1, -2, -2, 1, -1, 0,
1581		-1, 1, 1, 0, 0, 0, 1, 0, -1, 1, 1, 0, 0, 0,
1582		0, 1, 0, 1, -1, 0, 1, -1, -1, 2, 0, 0, 1, -1,
1583		0, 1, -1, 0,
1584		});
1585		}
1586
1587		test "uniform sampling" {
1588		var seed: [32]u8 = undefined;
1589		for (&seed, 0..) \|*s, i\| {
1590		s.* = @as(u8, @intCast(i));
1591		}
1592		try testing.expectEqual(Poly.uniform(seed, 1, 0).cs, .{
1593		797, 993, 161, 6, 2608, 2385, 2096, 2661, 1676, 247, 2440,
1594		342, 634, 194, 1570, 2848, 986, 684, 3148, 3208, 2018, 351,
1595		2288, 612, 1394, 170, 1521, 3119, 58, 596, 2093, 1549, 409,
1596		2156, 1934, 1730, 1324, 388, 446, 418, 1719, 2202, 1812, 98,
1597		1019, 2369, 214, 2699, 28, 1523, 2824, 273, 402, 2899, 246,
1598		210, 1288, 863, 2708, 177, 3076, 349, 44, 949, 854, 1371,
1599		957, 292, 2502, 1617, 1501, 254, 7, 1761, 2581, 2206, 2655,
1600		1211, 629, 1274, 2358, 816, 2766, 2115, 2985, 1006, 2433, 856,
1601		2596, 3192, 1, 1378, 2345, 707, 1891, 1669, 536, 1221, 710,
1602		2511, 120, 1176, 322, 1897, 2309, 595, 2950, 1171, 801, 1848,
1603		695, 2912, 1396, 1931, 1775, 2904, 893, 2507, 1810, 2873, 253,
1604		1529, 1047, 2615, 1687, 831, 1414, 965, 3169, 1887, 753, 3246,
1605		1937, 115, 2953, 586, 545, 1621, 1667, 3187, 1654, 1988, 1857,
1606		512, 1239, 1219, 898, 3106, 391, 1331, 2228, 3169, 586, 2412,
1607		845, 768, 156, 662, 478, 1693, 2632, 573, 2434, 1671, 173,
1608		969, 364, 1663, 2701, 2169, 813, 1000, 1471, 720, 2431, 2530,
1609		3161, 733, 1691, 527, 2634, 335, 26, 2377, 1707, 767, 3020,
1610		950, 502, 426, 1138, 3208, 2607, 2389, 44, 1358, 1392, 2334,
1611		875, 2097, 173, 1697, 2578, 942, 1817, 974, 1165, 2853, 1958,
1612		2973, 3282, 271, 1236, 1677, 2230, 673, 1554, 96, 242, 1729,
1613		2518, 1884, 2272, 71, 1382, 924, 1807, 1610, 456, 1148, 2479,
1614		2152, 238, 2208, 2329, 713, 1175, 1196, 757, 1078, 3190, 3169,
1615		708, 3117, 154, 1751, 3225, 1364, 154, 23, 2842, 1105, 1419,
1616		79, 5, 2013,
1617		});
1618		}
1619
1620		test "Polynomial packing" {
1621		var rnd = RndGen.init(0);
1622
1623		for (0..1000) \|_\| {
1624		const p = Poly.randNormalized(&rnd);
1625		try testing.expectEqual(Poly.fromBytes(&p.toBytes()), p);
1626		}
1627		}
1628
1629		test "Test inner PKE" {
1630		var seed: [32]u8 = undefined;
1631		var pt: [32]u8 = undefined;
1632		for (&seed, &pt, 0..) \|s, p, i\| {
1633		s.* = @as(u8, @intCast(i));
1634		p.* = @as(u8, @intCast(i + 32));
1635		}
1636		inline for (modes) \|mode\| {
1637		for (0..100) \|i\| {
1638		var pk: mode.InnerPk = undefined;
1639		var sk: mode.InnerSk = undefined;
1640		seed[0] = @as(u8, @intCast(i));
1641		mode.innerKeyFromSeed(seed, &pk, &sk);
1642		for (0..10) \|j\| {
1643		seed[1] = @as(u8, @intCast(j));
1644		try testing.expectEqual(sk.decrypt(&pk.encrypt(&pt, &seed)), pt);
1645		}
1646		}
1647		}
1648		}
1649
1650		test "Test happy flow" {
1651		var seed: [64]u8 = undefined;
1652		for (&seed, 0..) \|*s, i\| {
1653		s.* = @as(u8, @intCast(i));
1654		}
1655		inline for (modes) \|mode\| {
1656		for (0..100) \|i\| {
1657		seed[0] = @as(u8, @intCast(i));
1658		const kp = try mode.KeyPair.create(seed);
1659		const sk = try mode.SecretKey.fromBytes(&kp.secret_key.toBytes());
1660		try testing.expectEqual(sk, kp.secret_key);
1661		const pk = try mode.PublicKey.fromBytes(&kp.public_key.toBytes());
1662		try testing.expectEqual(pk, kp.public_key);
1663		for (0..10) \|j\| {
1664		seed[1] = @as(u8, @intCast(j));
1665		const e = pk.encaps(seed[0..32].*);
1666		try testing.expectEqual(e.shared_secret, try sk.decaps(&e.ciphertext));
1667		}
1668		}
1669		}
1670		}
1671
1672		// Code to test NIST Known Answer Tests (KAT), see PQCgenKAT.c.
1673
1674		const sha2 = crypto.hash.sha2;
1675
1676		test "NIST KAT test" {
1677		inline for (.{
1678		.{ Kyber512, "e9c2bd37133fcb40772f81559f14b1f58dccd1c816701be9ba6214d43baf4547" },
1679		.{ Kyber1024, "89248f2f33f7f4f7051729111f3049c409a933ec904aedadf035f30fa5646cd5" },
1680		.{ Kyber768, "a1e122cad3c24bc51622e4c242d8b8acbcd3f618fee4220400605ca8f9ea02c2" },
1681		}) \|modeHash\| {
1682		const mode = modeHash[0];
1683		var seed: [48]u8 = undefined;
1684		for (&seed, 0..) \|*s, i\| {
1685		s.* = @as(u8, @intCast(i));
1686		}
1687		var f = sha2.Sha256.init(.{});
1688		const fw = f.writer();
1689		var g = NistDRBG.init(seed);
1690		try std.fmt.format(fw, "# {s}\n\n", .{mode.name});
1691		for (0..100) \|i\| {
1692		g.fill(&seed);
1693		try std.fmt.format(fw, "count = {}\n", .{i});
1694		try std.fmt.format(fw, "seed = {s}\n", .{std.fmt.fmtSliceHexUpper(&seed)});
1695		var g2 = NistDRBG.init(seed);
1696
1697		// This is not equivalent to g2.fill(kseed[:]). As the reference
1698		// implementation calls randombytes twice generating the keypair,
1699		// we have to do that as well.
1700		var kseed: [64]u8 = undefined;
1701		var eseed: [32]u8 = undefined;
1702		g2.fill(kseed[0..32]);
1703		g2.fill(kseed[32..64]);
1704		g2.fill(&eseed);
1705		const kp = try mode.KeyPair.create(kseed);
1706		const e = kp.public_key.encaps(eseed);
1707		const ss2 = try kp.secret_key.decaps(&e.ciphertext);
1708		try testing.expectEqual(ss2, e.shared_secret);
1709		try std.fmt.format(fw, "pk = {s}\n", .{std.fmt.fmtSliceHexUpper(&kp.public_key.toBytes())});
1710		try std.fmt.format(fw, "sk = {s}\n", .{std.fmt.fmtSliceHexUpper(&kp.secret_key.toBytes())});
1711		try std.fmt.format(fw, "ct = {s}\n", .{std.fmt.fmtSliceHexUpper(&e.ciphertext)});
1712		try std.fmt.format(fw, "ss = {s}\n\n", .{std.fmt.fmtSliceHexUpper(&e.shared_secret)});
1713		}
1714
1715		var out: [32]u8 = undefined;
1716		f.final(&out);
1717		var outHex: [64]u8 = undefined;
1718		_ = try std.fmt.bufPrint(&outHex, "{s}", .{std.fmt.fmtSliceHexLower(&out)});
1719		try testing.expectEqual(outHex, modeHash[1].*);
1720		}
1721		}
1722
1723		const NistDRBG = struct {
1724		key: [32]u8,
1725		v: [16]u8,
1726
1727		fn incV(g: *NistDRBG) void {
1728		var j: usize = 15;
1729		while (j >= 0) : (j -= 1) {
1730		if (g.v[j] == 255) {
1731		g.v[j] = 0;
1732		} else {
1733		g.v[j] += 1;
1734		break;
1735		}
1736		}
1737		}
1738
1739		// AES256_CTR_DRBG_Update(pd, &g.key, &g.v).
1740		fn update(g: *NistDRBG, pd: ?[48]u8) void {
1741		var buf: [48]u8 = undefined;
1742		const ctx = crypto.core.aes.Aes256.initEnc(g.key);
1743		var i: usize = 0;
1744		while (i < 3) : (i += 1) {
1745		g.incV();
1746		var block: [16]u8 = undefined;
1747		ctx.encrypt(&block, &g.v);
1748		buf[i * 16 ..][0..16].* = block;
1749		}
1750		if (pd) \|p\| {
1751		for (&buf, p) \|*b, x\| {
1752		b.* ^= x;
1753		}
1754		}
1755		g.key = buf[0..32].*;
1756		g.v = buf[32..48].*;
1757		}
1758
1759		// randombytes.
1760		fn fill(g: *NistDRBG, out: []u8) void {
1761		var block: [16]u8 = undefined;
1762		var dst = out;
1763
1764		const ctx = crypto.core.aes.Aes256.initEnc(g.key);
1765		while (dst.len > 0) {
1766		g.incV();
1767		ctx.encrypt(&block, &g.v);
1768		if (dst.len < 16) {
1769		@memcpy(dst, block[0..dst.len]);
1770		break;
1771		}
1772		dst[0..block.len].* = block;
1773		dst = dst[16..dst.len];
1774		}
1775		g.update(null);
1776		}
1777
1778		fn init(seed: [48]u8) NistDRBG {
1779		var ret: NistDRBG = .{ .key = .{0} 32, .v = .{0} 16 };
1780		ret.update(seed);
1781		return ret;
1782		}
1783		};

lib/std/crypto/ml_kem.zig created+1830

...	...	@@ -0,0 +1,1830 @@
	1	//! Implementation of the IND-CCA2 post-quantum secure key encapsulation mechanism (KEM)
	2	//! ML-KEM (NIST FIPS-203 publication) and CRYSTALS-Kyber (v3.02/"draft00" CFRG draft).
	3	//!
	4	//! The schemes are not finalized yet, and are still subject to breaking changes.
	5	//!
	6	//! The Kyber namespace suffix (currently `_d00`) refers to the version currently
	7	//! implemented, in accordance with the draft.
	8	//! The ML-KEM namespace suffix (currently `_01`) refers to the NIST FIPS-203 draft
	9	//! published on August 24, 2023, with the unintentional transposition of Â having been reverted.
	10	//!
	11	//! Suffixes may not be updated if new versions of the documents only include editorial changes.
	12	//! The suffixes will be removed once the schemes are finalized.
	13	//!
	14	//! Quoting from the CFRG I-D:
	15	//!
	16	//! Kyber is not a Diffie-Hellman (DH) style non-interactive key
	17	//! agreement, but instead, Kyber is a Key Encapsulation Method (KEM).
	18	//! In essence, a KEM is a Public-Key Encryption (PKE) scheme where the
	19	//! plaintext cannot be specified, but is generated as a random key as
	20	//! part of the encryption. A KEM can be transformed into an unrestricted
	21	//! PKE using HPKE (RFC9180). On its own, a KEM can be used as a key
	22	//! agreement method in TLS.
	23	//!
	24	//! Kyber is an IND-CCA2 secure KEM. It is constructed by applying a
	25	//! Fujisaki--Okamato style transformation on InnerPKE, which is the
	26	//! underlying IND-CPA secure Public Key Encryption scheme. We cannot
	27	//! use InnerPKE directly, as its ciphertexts are malleable.
	28	//!
	29	//! ```
	30	//! F.O. transform
	31	//! InnerPKE ----------------------> Kyber
	32	//! IND-CPA IND-CCA2
	33	//! ```
	34	//!
	35	//! Kyber is a lattice-based scheme. More precisely, its security is
	36	//! based on the learning-with-errors-and-rounding problem in module
	37	//! lattices (MLWER). The underlying polynomial ring R (defined in
	38	//! Section 5) is chosen such that multiplication is very fast using the
	39	//! number theoretic transform (NTT, see Section 5.1.3).
	40	//!
	41	//! An InnerPKE private key is a vector _s_ over R of length k which is
	42	//! _small_ in a particular way. Here k is a security parameter akin to
	43	//! the size of a prime modulus. For Kyber512, which targets AES-128's
	44	//! security level, the value of k is 2.
	45	//!
	46	//! The public key consists of two values:
	47	//!
	48	//! * _A_ a uniformly sampled k by k matrix over R _and_
	49	//!
	50	//! * _t = A s + e_, where e is a suitably small masking vector.
	51	//!
	52	//! Distinguishing between such A s + e and a uniformly sampled t is the
	53	//! module learning-with-errors (MLWE) problem. If that is hard, then it
	54	//! is also hard to recover the private key from the public key as that
	55	//! would allow you to distinguish between those two.
	56	//!
	57	//! To save space in the public key, A is recomputed deterministically
	58	//! from a seed _rho_.
	59	//!
	60	//! A ciphertext for a message m under this public key is a pair (c_1,
	61	//! c_2) computed roughly as follows:
	62	//!
	63	//! c_1 = Compress(A^T r + e_1, d_u)
	64	//! c_2 = Compress(t^T r + e_2 + Decompress(m, 1), d_v)
	65	//!
	66	//! where
	67	//!
	68	//! * e_1, e_2 and r are small blinds;
	69	//!
	70	//! * Compress(-, d) removes some information, leaving d bits per
	71	//! coefficient and Decompress is such that Compress after Decompress
	72	//! does nothing and
	73	//!
	74	//! * d_u, d_v are scheme parameters.
	75	//!
	76	//! Distinguishing such a ciphertext and uniformly sampled (c_1, c_2) is
	77	//! an example of the full MLWER problem, see section 4.4 of [KyberV302].
	78	//!
	79	//! To decrypt the ciphertext, one computes
	80	//!
	81	//! m = Compress(Decompress(c_2, d_v) - s^T Decompress(c_1, d_u), 1).
	82	//!
	83	//! It it not straight-forward to see that this formula is correct. In
	84	//! fact, there is negligible but non-zero probability that a ciphertext
	85	//! does not decrypt correctly given by the DFP column in Table 4. This
	86	//! failure probability can be computed by a careful automated analysis
	87	//! of the probabilities involved, see kyber_failure.py of [SecEst].
	88	//!
	89	//! [KyberV302](https://pq-crystals.org/kyber/data/kyber-specification-round3-20210804.pdf)
	90	//! [I-D](https://github.com/bwesterb/draft-schwabe-cfrg-kyber)
	91	//! [SecEst](https://github.com/pq-crystals/security-estimates)
	92
	93	// TODO
	94	//
	95	// - The bottleneck in Kyber are the various hash/xof calls:
	96	// - Optimize Zig's keccak implementation.
	97	// - Use SIMD to compute keccak in parallel.
	98	// - Can we track bounds of coefficients using comptime types without
	99	// duplicating code?
	100	// - Would be neater to have tests closer to the thing under test.
	101	// - When generating a keypair, we have a copy of the inner public key with
	102	// its large matrix A in both the public key and the private key. In Go we
	103	// can just have a pointer in the private key to the public key, but
	104	// how do we do this elegantly in Zig?
	105
	106	const std = @import("std");
	107	const builtin = @import("builtin");
	108
	109	const testing = std.testing;
	110	const assert = std.debug.assert;
	111	const crypto = std.crypto;
	112	const errors = std.crypto.errors;
	113	const math = std.math;
	114	const mem = std.mem;
	115	const RndGen = std.Random.DefaultPrng;
	116	const sha3 = crypto.hash.sha3;
	117
	118	// Q is the parameter q ≡ 3329 = 2¹¹ + 2¹⁰ + 2⁸ + 1.
	119	const Q: i16 = 3329;
	120
	121	// Montgomery R
	122	const R: i32 = 1 << 16;
	123
	124	// Parameter n, degree of polynomials.
	125	const N: usize = 256;
	126
	127	// Size of "small" vectors used in encryption blinds.
	128	const eta2: u8 = 2;
	129
	130	const Params = struct {
	131	name: []const u8,
	132
	133	// NIST ML-KEM variant instead of Kyber as originally submitted.
	134	ml_kem: bool = false,
	135
	136	// Width and height of the matrix A.
	137	k: u8,
	138
	139	// Size of "small" vectors used in private key and encryption blinds.
	140	eta1: u8,
	141
	142	// How many bits to retain of u, the private-key independent part
	143	// of the ciphertext.
	144	du: u8,
	145
	146	// How many bits to retain of v, the private-key dependent part
	147	// of the ciphertext.
	148	dv: u8,
	149	};
	150
	151	pub const kyber_d00 = struct {
	152	pub const Kyber512 = Kyber(.{
	153	.name = "Kyber512",
	154	.k = 2,
	155	.eta1 = 3,
	156	.du = 10,
	157	.dv = 4,
	158	});
	159
	160	pub const Kyber768 = Kyber(.{
	161	.name = "Kyber768",
	162	.k = 3,
	163	.eta1 = 2,
	164	.du = 10,
	165	.dv = 4,
	166	});
	167
	168	pub const Kyber1024 = Kyber(.{
	169	.name = "Kyber1024",
	170	.k = 4,
	171	.eta1 = 2,
	172	.du = 11,
	173	.dv = 5,
	174	});
	175	};
	176
	177	pub const ml_kem_01 = struct {
	178	pub const MLKem512 = Kyber(.{
	179	.name = "ML-KEM-512",
	180	.ml_kem = true,
	181	.k = 2,
	182	.eta1 = 3,
	183	.du = 10,
	184	.dv = 4,
	185	});
	186
	187	pub const MLKem768 = Kyber(.{
	188	.name = "ML-KEM-768",
	189	.ml_kem = true,
	190	.k = 3,
	191	.eta1 = 2,
	192	.du = 10,
	193	.dv = 4,
	194	});
	195
	196	pub const MLKem1024 = Kyber(.{
	197	.name = "ML-KEM-1024",
	198	.ml_kem = true,
	199	.k = 4,
	200	.eta1 = 2,
	201	.du = 11,
	202	.dv = 5,
	203	});
	204	};
	205
	206	const modes = [_]type{
	207	kyber_d00.Kyber512,
	208	kyber_d00.Kyber768,
	209	kyber_d00.Kyber1024,
	210	ml_kem_01.MLKem512,
	211	ml_kem_01.MLKem768,
	212	ml_kem_01.MLKem1024,
	213	};
	214	const h_length: usize = 32;
	215	const inner_seed_length: usize = 32;
	216	const common_encaps_seed_length: usize = 32;
	217	const common_shared_key_size: usize = 32;
	218
	219	fn Kyber(comptime p: Params) type {
	220	return struct {
	221	// Size of a ciphertext, in bytes.
	222	pub const ciphertext_length = Poly.compressedSize(p.du) * p.k + Poly.compressedSize(p.dv);
	223
	224	const Self = @This();
	225	const V = Vec(p.k);
	226	const M = Mat(p.k);
	227
	228	/// Length (in bytes) of a shared secret.
	229	pub const shared_length = common_shared_key_size;
	230	/// Length (in bytes) of a seed for deterministic encapsulation.
	231	pub const encaps_seed_length = common_encaps_seed_length;
	232	/// Length (in bytes) of a seed for key generation.
	233	pub const seed_length: usize = inner_seed_length + shared_length;
	234	/// Algorithm name.
	235	pub const name = p.name;
	236
	237	/// A shared secret, and an encapsulated (encrypted) representation of it.
	238	pub const EncapsulatedSecret = struct {
	239	shared_secret: [shared_length]u8,
	240	ciphertext: [ciphertext_length]u8,
	241	};
	242
	243	/// A Kyber public key.
	244	pub const PublicKey = struct {
	245	pk: InnerPk,
	246
	247	// Cached
	248	hpk: [h_length]u8, // H(pk)
	249
	250	/// Size of a serialized representation of the key, in bytes.
	251	pub const bytes_length = InnerPk.bytes_length;
	252
	253	/// Generates a shared secret, and encapsulates it for the public key.
	254	/// If `seed` is `null`, a random seed is used. This is recommended.
	255	/// If `seed` is set, encapsulation is deterministic.
	256	pub fn encaps(pk: PublicKey, seed_: ?[encaps_seed_length]u8) EncapsulatedSecret {
	257	var m: [inner_plaintext_length]u8 = undefined;
	258
	259	if (seed_) \|seed\| {
	260	if (p.ml_kem) {
	261	@memcpy(&m, &seed);
	262	} else {
	263	// m = H(seed)
	264	sha3.Sha3_256.hash(&seed, &m, .{});
	265	}
	266	} else {
	267	crypto.random.bytes(&m);
	268	}
	269
	270	// (K', r) = G(m ‖ H(pk))
	271	var kr: [inner_plaintext_length + h_length]u8 = undefined;
	272	var g = sha3.Sha3_512.init(.{});
	273	g.update(&m);
	274	g.update(&pk.hpk);
	275	g.final(&kr);
	276
	277	// c = innerEncrypy(pk, m, r)
	278	const ct = pk.pk.encrypt(&m, kr[32..64]);
	279
	280	// Compute H(c) and put in second slot of kr, which will be (K', H(c)).
	281	sha3.Sha3_256.hash(&ct, kr[32..], .{});
	282
	283	if (p.ml_kem) {
	284	return EncapsulatedSecret{
	285	.shared_secret = kr[0..shared_length].*, // ML-KEM: K = K'
	286	.ciphertext = ct,
	287	};
	288	} else {
	289	var ss: [shared_length]u8 = undefined;
	290	sha3.Shake256.hash(&kr, &ss, .{});
	291	return EncapsulatedSecret{
	292	.shared_secret = ss, // Kyber: K = KDF(K' ‖ H(c))
	293	.ciphertext = ct,
	294	};
	295	}
	296	}
	297
	298	/// Serializes the key into a byte array.
	299	pub fn toBytes(pk: PublicKey) [bytes_length]u8 {
	300	return pk.pk.toBytes();
	301	}
	302
	303	/// Deserializes the key from a byte array.
	304	pub fn fromBytes(buf: *const [bytes_length]u8) errors.NonCanonicalError!PublicKey {
	305	var ret: PublicKey = undefined;
	306	ret.pk = try InnerPk.fromBytes(buf[0..InnerPk.bytes_length]);
	307	sha3.Sha3_256.hash(buf, &ret.hpk, .{});
	308	return ret;
	309	}
	310	};
	311
	312	/// A Kyber secret key.
	313	pub const SecretKey = struct {
	314	sk: InnerSk,
	315	pk: InnerPk,
	316	hpk: [h_length]u8, // H(pk)
	317	z: [shared_length]u8,
	318
	319	/// Size of a serialized representation of the key, in bytes.
	320	pub const bytes_length: usize =
	321	InnerSk.bytes_length + InnerPk.bytes_length + h_length + shared_length;
	322
	323	/// Decapsulates the shared secret within ct using the private key.
	324	pub fn decaps(sk: SecretKey, ct: *const [ciphertext_length]u8) ![shared_length]u8 {
	325	// m' = innerDec(ct)
	326	const m2 = sk.sk.decrypt(ct);
	327
	328	// (K'', r') = G(m' ‖ H(pk))
	329	var kr2: [64]u8 = undefined;
	330	var g = sha3.Sha3_512.init(.{});
	331	g.update(&m2);
	332	g.update(&sk.hpk);
	333	g.final(&kr2);
	334
	335	// ct' = innerEnc(pk, m', r')
	336	const ct2 = sk.pk.encrypt(&m2, kr2[32..64]);
	337
	338	// Compute H(ct) and put in the second slot of kr2 which will be (K'', H(ct)).
	339	sha3.Sha3_256.hash(ct, kr2[32..], .{});
	340
	341	// Replace K'' by z in the first slot of kr2 if ct ≠ ct'.
	342	cmov(32, kr2[0..32], sk.z, ctneq(ciphertext_length, ct.*, ct2));
	343
	344	if (p.ml_kem) {
	345	// ML-KEM: K = K''/z
	346	return kr2[0..shared_length].*;
	347	} else {
	348	// Kyber: K = KDF(K''/z ‖ H(c))
	349	var ss: [shared_length]u8 = undefined;
	350	sha3.Shake256.hash(&kr2, &ss, .{});
	351	return ss;
	352	}
	353	}
	354
	355	/// Serializes the key into a byte array.
	356	pub fn toBytes(sk: SecretKey) [bytes_length]u8 {
	357	return sk.sk.toBytes() ++ sk.pk.toBytes() ++ sk.hpk ++ sk.z;
	358	}
	359
	360	/// Deserializes the key from a byte array.
	361	pub fn fromBytes(buf: *const [bytes_length]u8) errors.NonCanonicalError!SecretKey {
	362	var ret: SecretKey = undefined;
	363	comptime var s: usize = 0;
	364	ret.sk = InnerSk.fromBytes(buf[s .. s + InnerSk.bytes_length]);
	365	s += InnerSk.bytes_length;
	366	ret.pk = try InnerPk.fromBytes(buf[s .. s + InnerPk.bytes_length]);
	367	s += InnerPk.bytes_length;
	368	ret.hpk = buf[s..][0..h_length].*;
	369	s += h_length;
	370	ret.z = buf[s..][0..shared_length].*;
	371	return ret;
	372	}
	373	};
	374
	375	/// A Kyber key pair.
	376	pub const KeyPair = struct {
	377	secret_key: SecretKey,
	378	public_key: PublicKey,
	379
	380	/// Create a new key pair.
	381	/// If seed is null, a random seed will be generated.
	382	/// If a seed is provided, the key pair will be determinsitic.
	383	pub fn create(seed_: ?[seed_length]u8) !KeyPair {
	384	const seed = seed_ orelse sk: {
	385	var random_seed: [seed_length]u8 = undefined;
	386	crypto.random.bytes(&random_seed);
	387	break :sk random_seed;
	388	};
	389	var ret: KeyPair = undefined;
	390	ret.secret_key.z = seed[inner_seed_length..seed_length].*;
	391
	392	// Generate inner key
	393	innerKeyFromSeed(
	394	seed[0..inner_seed_length].*,
	395	&ret.public_key.pk,
	396	&ret.secret_key.sk,
	397	);
	398	ret.secret_key.pk = ret.public_key.pk;
	399
	400	// Copy over z from seed.
	401	ret.secret_key.z = seed[inner_seed_length..seed_length].*;
	402
	403	// Compute H(pk)
	404	sha3.Sha3_256.hash(&ret.public_key.pk.toBytes(), &ret.secret_key.hpk, .{});
	405	ret.public_key.hpk = ret.secret_key.hpk;
	406
	407	return ret;
	408	}
	409	};
	410
	411	// Size of plaintexts of the in
	412	const inner_plaintext_length: usize = Poly.compressedSize(1);
	413
	414	const InnerPk = struct {
	415	rho: [32]u8, // ρ, the seed for the matrix A
	416	th: V, // NTT(t), normalized
	417
	418	// Cached values
	419	aT: M,
	420
	421	const bytes_length = V.bytes_length + 32;
	422
	423	fn encrypt(
	424	pk: InnerPk,
	425	pt: *const [inner_plaintext_length]u8,
	426	seed: *const [32]u8,
	427	) [ciphertext_length]u8 {
	428	// Sample r, e₁ and e₂ appropriately
	429	const rh = V.noise(p.eta1, 0, seed).ntt().barrettReduce();
	430	const e1 = V.noise(eta2, p.k, seed);
	431	const e2 = Poly.noise(eta2, 2 * p.k, seed);
	432
	433	// Next we compute u = Aᵀ r + e₁. First Aᵀ.
	434	var u: V = undefined;
	435	for (0..p.k) \|i\| {
	436	// Note that coefficients of r are bounded by q and those of Aᵀ
	437	// are bounded by 4.5q and so their product is bounded by 2¹⁵q
	438	// as required for multiplication.
	439	u.ps[i] = pk.aT.vs[i].dotHat(rh);
	440	}
	441
	442	// Aᵀ and r were not in Montgomery form, so the Montgomery
	443	// multiplications in the inner product added a factor R⁻¹ which
	444	// the InvNTT cancels out.
	445	u = u.barrettReduce().invNTT().add(e1).normalize();
	446
	447	// Next, compute v = <t, r> + e₂ + Decompress_q(m, 1)
	448	const v = pk.th.dotHat(rh).barrettReduce().invNTT()
	449	.add(Poly.decompress(1, pt)).add(e2).normalize();
	450
	451	return u.compress(p.du) ++ v.compress(p.dv);
	452	}
	453
	454	fn toBytes(pk: InnerPk) [bytes_length]u8 {
	455	return pk.th.toBytes() ++ pk.rho;
	456	}
	457
	458	fn fromBytes(buf: *const [bytes_length]u8) errors.NonCanonicalError!InnerPk {
	459	var ret: InnerPk = undefined;
	460
	461	const th_bytes = buf[0..V.bytes_length];
	462	ret.th = V.fromBytes(th_bytes).normalize();
	463
	464	if (p.ml_kem) {
	465	// Verify that the coefficients used a canonical representation.
	466	if (!mem.eql(u8, &ret.th.toBytes(), th_bytes)) {
	467	return error.NonCanonical;
	468	}
	469	}
	470
	471	ret.rho = buf[V.bytes_length..bytes_length].*;
	472	ret.aT = M.uniform(ret.rho, true);
	473	return ret;
	474	}
	475	};
	476
	477	// Private key of the inner PKE
	478	const InnerSk = struct {
	479	sh: V, // NTT(s), normalized
	480	const bytes_length = V.bytes_length;
	481
	482	fn decrypt(sk: InnerSk, ct: *const [ciphertext_length]u8) [inner_plaintext_length]u8 {
	483	const u = V.decompress(p.du, ct[0..comptime V.compressedSize(p.du)]);
	484	const v = Poly.decompress(
	485	p.dv,
	486	ct[comptime V.compressedSize(p.du)..ciphertext_length],
	487	);
	488
	489	// Compute m = v - <s, u>
	490	return v.sub(sk.sh.dotHat(u.ntt()).barrettReduce().invNTT())
	491	.normalize().compress(1);
	492	}
	493
	494	fn toBytes(sk: InnerSk) [bytes_length]u8 {
	495	return sk.sh.toBytes();
	496	}
	497
	498	fn fromBytes(buf: *const [bytes_length]u8) InnerSk {
	499	var ret: InnerSk = undefined;
	500	ret.sh = V.fromBytes(buf).normalize();
	501	return ret;
	502	}
	503	};
	504
	505	// Derives inner PKE keypair from given seed.
	506	fn innerKeyFromSeed(seed: [inner_seed_length]u8, pk: InnerPk, sk: InnerSk) void {
	507	var expanded_seed: [64]u8 = undefined;
	508	sha3.Sha3_512.hash(&seed, &expanded_seed, .{});
	509	pk.rho = expanded_seed[0..32].*;
	510	const sigma = expanded_seed[32..64];
	511	pk.aT = M.uniform(pk.rho, false); // Expand ρ to A; we'll transpose later on
	512
	513	// Sample secret vector s.
	514	sk.sh = V.noise(p.eta1, 0, sigma).ntt().normalize();
	515
	516	const eh = Vec(p.k).noise(p.eta1, p.k, sigma).ntt(); // sample blind e.
	517	var th: V = undefined;
	518
	519	// Next, we compute t = A s + e.
	520	for (0..p.k) \|i\| {
	521	// Note that coefficients of s are bounded by q and those of A
	522	// are bounded by 4.5q and so their product is bounded by 2¹⁵q
	523	// as required for multiplication.
	524	// A and s were not in Montgomery form, so the Montgomery
	525	// multiplications in the inner product added a factor R⁻¹ which
	526	// we'll cancel out with toMont(). This will also ensure the
	527	// coefficients of th are bounded in absolute value by q.
	528	th.ps[i] = pk.aT.vs[i].dotHat(sk.sh).toMont();
	529	}
	530
	531	pk.th = th.add(eh).normalize(); // bounded by 8q
	532	pk.aT = pk.aT.transpose();
	533	}
	534	};
	535	}
	536
	537	// R mod q
	538	const r_mod_q: i32 = @rem(@as(i32, R), Q);
	539
	540	// R² mod q
	541	const r2_mod_q: i32 = @rem(r_mod_q * r_mod_q, Q);
	542
	543	// ζ is the degree 256 primitive root of unity used for the NTT.
	544	const zeta: i16 = 17;
	545
	546	// (128)⁻¹ R². Used in inverse NTT.
	547	const r2_over_128: i32 = @mod(invertMod(128, Q) * r2_mod_q, Q);
	548
	549	// zetas lists precomputed powers of the primitive root of unity in
	550	// Montgomery representation used for the NTT:
	551	//
	552	// zetas[i] = ζᵇʳᵛ⁽ⁱ⁾ R mod q
	553	//
	554	// where ζ = 17, brv(i) is the bitreversal of a 7-bit number and R=2¹⁶ mod q.
	555	const zetas = computeZetas();
	556
	557	// invNTTReductions keeps track of which coefficients to apply Barrett
	558	// reduction to in Poly.invNTT().
	559	//
	560	// Generated lazily: once a butterfly is computed which is about to
	561	// overflow the i16, the largest coefficient is reduced. If that is
	562	// not enough, the other coefficient is reduced as well.
	563	//
	564	// This is actually optimal, as proven in https://eprint.iacr.org/2020/1377.pdf
	565	// TODO generate comptime?
	566	const inv_ntt_reductions = [_]i16{
	567	-1, // after layer 1
	568	-1, // after layer 2
	569	16,
	570	17,
	571	48,
	572	49,
	573	80,
	574	81,
	575	112,
	576	113,
	577	144,
	578	145,
	579	176,
	580	177,
	581	208,
	582	209,
	583	240, 241, -1, // after layer 3
	584	0, 1, 32,
	585	33, 34, 35,
	586	64, 65, 96,
	587	97, 98, 99,
	588	128, 129,
	589	160, 161, 162, 163, 192, 193, 224, 225, 226, 227, -1, // after layer 4
	590	2, 3, 66, 67, 68, 69, 70, 71, 130, 131, 194,
	591	195, 196, 197,
	592	198, 199, -1, // after layer 5
	593	4, 5, 6,
	594	7, 132, 133,
	595	134, 135, 136,
	596	137, 138, 139,
	597	140, 141,
	598	142, 143, -1, // after layer 6
	599	-1, // after layer 7
	600	};
	601
	602	test "invNTTReductions bounds" {
	603	// Checks whether the reductions proposed by invNTTReductions
	604	// don't overflow during invNTT().
	605	var xs = [_]i32{1} ** 256; // start at \|x\| ≤ q
	606
	607	var r: usize = 0;
	608	var layer: math.Log2Int(usize) = 1;
	609	while (layer < 8) : (layer += 1) {
	610	const w = @as(usize, 1) << layer;
	611	var i: usize = 0;
	612
	613	while (i + w < 256) {
	614	xs[i] = xs[i] + xs[i + w];
	615	try testing.expect(xs[i] <= 9); // we can't exceed 9q
	616	xs[i + w] = 1;
	617	i += 1;
	618	if (@mod(i, w) == 0) {
	619	i += w;
	620	}
	621	}
	622
	623	while (true) {
	624	const j = inv_ntt_reductions[r];
	625	r += 1;
	626	if (j < 0) {
	627	break;
	628	}
	629	xs[@as(usize, @intCast(j))] = 1;
	630	}
	631	}
	632	}
	633
	634	// Extended euclidean algorithm.
	635	//
	636	// For a, b finds x, y such that x a + y b = gcd(a, b). Used to compute
	637	// modular inverse.
	638	fn eea(a: anytype, b: @TypeOf(a)) EeaResult(@TypeOf(a)) {
	639	if (a == 0) {
	640	return .{ .gcd = b, .x = 0, .y = 1 };
	641	}
	642	const r = eea(@rem(b, a), a);
	643	return .{ .gcd = r.gcd, .x = r.y - @divTrunc(b, a) * r.x, .y = r.x };
	644	}
	645
	646	fn EeaResult(comptime T: type) type {
	647	return struct { gcd: T, x: T, y: T };
	648	}
	649
	650	// Returns least common multiple of a and b.
	651	fn lcm(a: anytype, b: @TypeOf(a)) @TypeOf(a) {
	652	const r = eea(a, b);
	653	return a * b / r.gcd;
	654	}
	655
	656	// Invert modulo p.
	657	fn invertMod(a: anytype, p: @TypeOf(a)) @TypeOf(a) {
	658	const r = eea(a, p);
	659	assert(r.gcd == 1);
	660	return r.x;
	661	}
	662
	663	// Reduce mod q for testing.
	664	fn modQ32(x: i32) i16 {
	665	var y = @as(i16, @intCast(@rem(x, @as(i32, Q))));
	666	if (y < 0) {
	667	y += Q;
	668	}
	669	return y;
	670	}
	671
	672	// Given -2¹⁵ q ≤ x < 2¹⁵ q, returns -q < y < q with x 2⁻¹⁶ = y (mod q).
	673	fn montReduce(x: i32) i16 {
	674	const qInv = comptime invertMod(@as(i32, Q), R);
	675	// This is Montgomery reduction with R=2¹⁶.
	676	//
	677	// Note gcd(2¹⁶, q) = 1 as q is prime. Write q' := 62209 = q⁻¹ mod R.
	678	// First we compute
	679	//
	680	// m := ((x mod R) q') mod R
	681	// = x q' mod R
	682	// = int16(x q')
	683	// = int16(int32(x) * int32(q'))
	684	//
	685	// Note that x q' might be as big as 2³² and could overflow the int32
	686	// multiplication in the last line. However for any int32s a and b,
	687	// we have int32(int64(a)int64(b)) = int32(ab) and so the result is ok.
	688	const m: i16 = @truncate(@as(i32, @truncate(x *% qInv)));
	689
	690	// Note that x - m q is divisible by R; indeed modulo R we have
	691	//
	692	// x - m q ≡ x - x q' q ≡ x - x q⁻¹ q ≡ x - x = 0.
	693	//
	694	// We return y := (x - m q) / R. Note that y is indeed correct as
	695	// modulo q we have
	696	//
	697	// y ≡ x R⁻¹ - m q R⁻¹ = x R⁻¹
	698	//
	699	// and as both 2¹⁵ q ≤ m q, x < 2¹⁵ q, we have
	700	// 2¹⁶ q ≤ x - m q < 2¹⁶ and so q ≤ (x - m q) / R < q as desired.
	701	const yR = x - @as(i32, m) * @as(i32, Q);
	702	return @bitCast(@as(u16, @truncate(@as(u32, @bitCast(yR)) >> 16)));
	703	}
	704
	705	test "Test montReduce" {
	706	var rnd = RndGen.init(0);
	707	for (0..1000) \|_\| {
	708	const bound = comptime @as(i32, Q) * (1 << 15);
	709	const x = rnd.random().intRangeLessThan(i32, -bound, bound);
	710	const y = montReduce(x);
	711	try testing.expect(-Q < y and y < Q);
	712	try testing.expectEqual(modQ32(x), modQ32(@as(i32, y) * R));
	713	}
	714	}
	715
	716	// Given any x, return x R mod q where R=2¹⁶.
	717	fn feToMont(x: i16) i16 {
	718	// Note \|1353 x\| ≤ 1353 2¹⁵ ≤ 13318 q ≤ 2¹⁵ q and so we're within
	719	// the bounds of montReduce.
	720	return montReduce(@as(i32, x) * r2_mod_q);
	721	}
	722
	723	test "Test feToMont" {
	724	var x: i32 = -(1 << 15);
	725	while (x < 1 << 15) : (x += 1) {
	726	const y = feToMont(@as(i16, @intCast(x)));
	727	try testing.expectEqual(modQ32(@as(i32, y)), modQ32(x * r_mod_q));
	728	}
	729	}
	730
	731	// Given any x, compute 0 ≤ y ≤ q with x = y (mod q).
	732	//
	733	// Beware: we might have feBarrettReduce(x) = q ≠ 0 for some x. In fact,
	734	// this happens if and only if x = -nq for some positive integer n.
	735	fn feBarrettReduce(x: i16) i16 {
	736	// This is standard Barrett reduction.
	737	//
	738	// For any x we have x mod q = x - ⌊x/q⌋ q. We will use 20159/2²⁶ as
	739	// an approximation of 1/q. Note that 0 ≤ 20159/2²⁶ - 1/q ≤ 0.135/2²⁶
	740	// and so \| x 20156/2²⁶ - x/q \| ≤ 2⁻¹⁰ for \|x\| ≤ 2¹⁶. For all x
	741	// not a multiple of q, the number x/q is further than 1/q from any integer
	742	// and so ⌊x 20156/2²⁶⌋ = ⌊x/q⌋. If x is a multiple of q and x is positive,
	743	// then x 20156/2²⁶ is larger than x/q so ⌊x 20156/2²⁶⌋ = ⌊x/q⌋ as well.
	744	// Finally, if x is negative multiple of q, then ⌊x 20156/2²⁶⌋ = ⌊x/q⌋-1.
	745	// Thus
	746	// [ q if x=-nq for pos. integer n
	747	// x - ⌊x 20156/2²⁶⌋ q = [
	748	// [ x mod q otherwise
	749	//
	750	// To actually compute this, note that
	751	//
	752	// ⌊x 20156/2²⁶⌋ = (20159 x) >> 26.
	753	return x -% @as(i16, @intCast((@as(i32, x) * 20159) >> 26)) *% Q;
	754	}
	755
	756	test "Test Barrett reduction" {
	757	var x: i32 = -(1 << 15);
	758	while (x < 1 << 15) : (x += 1) {
	759	var y1 = feBarrettReduce(@as(i16, @intCast(x)));
	760	const y2 = @mod(@as(i16, @intCast(x)), Q);
	761	if (x < 0 and @rem(-x, Q) == 0) {
	762	y1 -= Q;
	763	}
	764	try testing.expectEqual(y1, y2);
	765	}
	766	}
	767
	768	// Returns x if x < q and x - q otherwise. Assumes x ≥ -29439.
	769	fn csubq(x: i16) i16 {
	770	var r = x;
	771	r -= Q;
	772	r += (r >> 15) & Q;
	773	return r;
	774	}
	775
	776	test "Test csubq" {
	777	var x: i32 = -29439;
	778	while (x < 1 << 15) : (x += 1) {
	779	const y1 = csubq(@as(i16, @intCast(x)));
	780	var y2 = @as(i16, @intCast(x));
	781	if (@as(i16, @intCast(x)) >= Q) {
	782	y2 -= Q;
	783	}
	784	try testing.expectEqual(y1, y2);
	785	}
	786	}
	787
	788	// Compute a^s mod p.
	789	fn mpow(a: anytype, s: @TypeOf(a), p: @TypeOf(a)) @TypeOf(a) {
	790	var ret: @TypeOf(a) = 1;
	791	var s2 = s;
	792	var a2 = a;
	793
	794	while (true) {
	795	if (s2 & 1 == 1) {
	796	ret = @mod(ret * a2, p);
	797	}
	798	s2 >>= 1;
	799	if (s2 == 0) {
	800	break;
	801	}
	802	a2 = @mod(a2 * a2, p);
	803	}
	804	return ret;
	805	}
	806
	807	// Computes zetas table used by ntt and invNTT.
	808	fn computeZetas() [128]i16 {
	809	@setEvalBranchQuota(10000);
	810	var ret: [128]i16 = undefined;
	811	for (&ret, 0..) \|*r, i\| {
	812	const t = @as(i16, @intCast(mpow(@as(i32, zeta), @bitReverse(@as(u7, @intCast(i))), Q)));
	813	r.* = csubq(feBarrettReduce(feToMont(t)));
	814	}
	815	return ret;
	816	}
	817
	818	// An element of our base ring R which are polynomials over ℤ_q
	819	// modulo the equation Xᴺ = -1, where q=3329 and N=256.
	820	//
	821	// This type is also used to store NTT-transformed polynomials,
	822	// see Poly.NTT().
	823	//
	824	// Coefficients aren't always reduced. See Normalize().
	825	const Poly = struct {
	826	cs: [N]i16,
	827
	828	const bytes_length = N / 2 * 3;
	829	const zero: Poly = .{ .cs = .{0} ** N };
	830
	831	fn add(a: Poly, b: Poly) Poly {
	832	var ret: Poly = undefined;
	833	for (0..N) \|i\| {
	834	ret.cs[i] = a.cs[i] + b.cs[i];
	835	}
	836	return ret;
	837	}
	838
	839	fn sub(a: Poly, b: Poly) Poly {
	840	var ret: Poly = undefined;
	841	for (0..N) \|i\| {
	842	ret.cs[i] = a.cs[i] - b.cs[i];
	843	}
	844	return ret;
	845	}
	846
	847	// For testing, generates a random polynomial with for each
	848	// coefficient \|x\| ≤ q.
	849	fn randAbsLeqQ(rnd: anytype) Poly {
	850	var ret: Poly = undefined;
	851	for (0..N) \|i\| {
	852	ret.cs[i] = rnd.random().intRangeAtMost(i16, -Q, Q);
	853	}
	854	return ret;
	855	}
	856
	857	// For testing, generates a random normalized polynomial.
	858	fn randNormalized(rnd: anytype) Poly {
	859	var ret: Poly = undefined;
	860	for (0..N) \|i\| {
	861	ret.cs[i] = rnd.random().intRangeLessThan(i16, 0, Q);
	862	}
	863	return ret;
	864	}
	865
	866	// Executes a forward "NTT" on p.
	867	//
	868	// Assumes the coefficients are in absolute value ≤q. The resulting
	869	// coefficients are in absolute value ≤7q. If the input is in Montgomery
	870	// form, then the result is in Montgomery form and so (by linearity of the NTT)
	871	// if the input is in regular form, then the result is also in regular form.
	872	fn ntt(a: Poly) Poly {
	873	// Note that ℤ_q does not have a primitive 512ᵗʰ root of unity (as 512
	874	// does not divide into q-1) and so we cannot do a regular NTT. ℤ_q
	875	// does have a primitive 256ᵗʰ root of unity, the smallest of which
	876	// is ζ := 17.
	877	//
	878	// Recall that our base ring R := ℤ_q[x] / (x²⁵⁶ + 1). The polynomial
	879	// x²⁵⁶+1 will not split completely (as its roots would be 512ᵗʰ roots
	880	// of unity.) However, it does split almost (using ζ¹²⁸ = -1):
	881	//
	882	// x²⁵⁶ + 1 = (x²)¹²⁸ - ζ¹²⁸
	883	// = ((x²)⁶⁴ - ζ⁶⁴)((x²)⁶⁴ + ζ⁶⁴)
	884	// = ((x²)³² - ζ³²)((x²)³² + ζ³²)((x²)³² - ζ⁹⁶)((x²)³² + ζ⁹⁶)
	885	// ⋮
	886	// = (x² - ζ)(x² + ζ)(x² - ζ⁶⁵)(x² + ζ⁶⁵) … (x² + ζ¹²⁷)
	887	//
	888	// Note that the powers of ζ that appear (from the second line down) are
	889	// in binary
	890	//
	891	// 0100000 1100000
	892	// 0010000 1010000 0110000 1110000
	893	// 0001000 1001000 0101000 1101000 0011000 1011000 0111000 1111000
	894	// …
	895	//
	896	// That is: brv(2), brv(3), brv(4), …, where brv(x) denotes the 7-bit
	897	// bitreversal of x. These powers of ζ are given by the Zetas array.
	898	//
	899	// The polynomials x² ± ζⁱ are irreducible and coprime, hence by
	900	// the Chinese Remainder Theorem we know
	901	//
	902	// ℤ_q[x]/(x²⁵⁶+1) → ℤ_q[x]/(x²-ζ) x … x ℤ_q[x]/(x²+ζ¹²⁷)
	903	//
	904	// given by a ↦ ( a mod x²-ζ, …, a mod x²+ζ¹²⁷ )
	905	// is an isomorphism, which is the "NTT". It can be efficiently computed by
	906	//
	907	//
	908	// a ↦ ( a mod (x²)⁶⁴ - ζ⁶⁴, a mod (x²)⁶⁴ + ζ⁶⁴ )
	909	// ↦ ( a mod (x²)³² - ζ³², a mod (x²)³² + ζ³²,
	910	// a mod (x²)⁹⁶ - ζ⁹⁶, a mod (x²)⁹⁶ + ζ⁹⁶ )
	911	//
	912	// et cetera
	913	// If N was 8 then this can be pictured in the following diagram:
	914	//
	915	// https://cnx.org/resources/17ee4dfe517a6adda05377b25a00bf6e6c93c334/File0026.png
	916	//
	917	// Each cross is a Cooley-Tukey butterfly: it's the map
	918	//
	919	// (a, b) ↦ (a + ζb, a - ζb)
	920	//
	921	// for the appropriate power ζ for that column and row group.
	922	var p = a;
	923	var k: usize = 0; // index into zetas
	924
	925	var l = N >> 1;
	926	while (l > 1) : (l >>= 1) {
	927	// On the nᵗʰ iteration of the l-loop, the absolute value of the
	928	// coefficients are bounded by nq.
	929
	930	// offset effectively loops over the row groups in this column; it is
	931	// the first row in the row group.
	932	var offset: usize = 0;
	933	while (offset < N - l) : (offset += 2 * l) {
	934	k += 1;
	935	const z = @as(i32, zetas[k]);
	936
	937	// j loops over each butterfly in the row group.
	938	for (offset..offset + l) \|j\| {
	939	const t = montReduce(z * @as(i32, p.cs[j + l]));
	940	p.cs[j + l] = p.cs[j] - t;
	941	p.cs[j] += t;
	942	}
	943	}
	944	}
	945
	946	return p;
	947	}
	948
	949	// Executes an inverse "NTT" on p and multiply by the Montgomery factor R.
	950	//
	951	// Assumes the coefficients are in absolute value ≤q. The resulting
	952	// coefficients are in absolute value ≤q. If the input is in Montgomery
	953	// form, then the result is in Montgomery form and so (by linearity)
	954	// if the input is in regular form, then the result is also in regular form.
	955	fn invNTT(a: Poly) Poly {
	956	var k: usize = 127; // index into zetas
	957	var r: usize = 0; // index into invNTTReductions
	958	var p = a;
	959
	960	// We basically do the oppposite of NTT, but postpone dividing by 2 in the
	961	// inverse of the Cooley-Tukey butterfly and accumulate that into a big
	962	// division by 2⁷ at the end. See the comments in the ntt() function.
	963
	964	var l: usize = 2;
	965	while (l < N) : (l <<= 1) {
	966	var offset: usize = 0;
	967	while (offset < N - l) : (offset += 2 * l) {
	968	// As we're inverting, we need powers of ζ⁻¹ (instead of ζ).
	969	// To be precise, we need ζᵇʳᵛ⁽ᵏ⁾⁻¹²⁸. However, as ζ⁻¹²⁸ = -1,
	970	// we can use the existing zetas table instead of
	971	// keeping a separate invZetas table as in Dilithium.
	972
	973	const minZeta = @as(i32, zetas[k]);
	974	k -= 1;
	975
	976	for (offset..offset + l) \|j\| {
	977	// Gentleman-Sande butterfly: (a, b) ↦ (a + b, ζ(a-b))
	978	const t = p.cs[j + l] - p.cs[j];
	979	p.cs[j] += p.cs[j + l];
	980	p.cs[j + l] = montReduce(minZeta * @as(i32, t));
	981
	982	// Note that if we had \|a\| < αq and \|b\| < βq before the
	983	// butterfly, then now we have \|a\| < (α+β)q and \|b\| < q.
	984	}
	985	}
	986
	987	// We let the invNTTReductions instruct us which coefficients to
	988	// Barrett reduce.
	989	while (true) {
	990	const i = inv_ntt_reductions[r];
	991	r += 1;
	992	if (i < 0) {
	993	break;
	994	}
	995	p.cs[@as(usize, @intCast(i))] = feBarrettReduce(p.cs[@as(usize, @intCast(i))]);
	996	}
	997	}
	998
	999	for (0..N) \|j\| {
	1000	// Note 1441 = (128)⁻¹ R². The coefficients are bounded by 9q, so
	1001	// as 1441 * 9 ≈ 2¹⁴ < 2¹⁵, we're within the required bounds
	1002	// for montReduce().
	1003	p.cs[j] = montReduce(r2_over_128 * @as(i32, p.cs[j]));
	1004	}
	1005
	1006	return p;
	1007	}
	1008
	1009	// Normalizes coefficients.
	1010	//
	1011	// Ensures each coefficient is in {0, …, q-1}.
	1012	fn normalize(a: Poly) Poly {
	1013	var ret: Poly = undefined;
	1014	for (0..N) \|i\| {
	1015	ret.cs[i] = csubq(feBarrettReduce(a.cs[i]));
	1016	}
	1017	return ret;
	1018	}
	1019
	1020	// Put p in Montgomery form.
	1021	fn toMont(a: Poly) Poly {
	1022	var ret: Poly = undefined;
	1023	for (0..N) \|i\| {
	1024	ret.cs[i] = feToMont(a.cs[i]);
	1025	}
	1026	return ret;
	1027	}
	1028
	1029	// Barret reduce coefficients.
	1030	//
	1031	// Beware, this does not fully normalize coefficients.
	1032	fn barrettReduce(a: Poly) Poly {
	1033	var ret: Poly = undefined;
	1034	for (0..N) \|i\| {
	1035	ret.cs[i] = feBarrettReduce(a.cs[i]);
	1036	}
	1037	return ret;
	1038	}
	1039
	1040	fn compressedSize(comptime d: u8) usize {
	1041	return @divTrunc(N * d, 8);
	1042	}
	1043
	1044	// Returns packed Compress_q(p, d).
	1045	//
	1046	// Assumes p is normalized.
	1047	fn compress(p: Poly, comptime d: u8) [compressedSize(d)]u8 {
	1048	@setEvalBranchQuota(10000);
	1049	const q_over_2: u32 = comptime @divTrunc(Q, 2); // (q-1)/2
	1050	const two_d_min_1: u32 = comptime (1 << d) - 1; // 2ᵈ-1
	1051	var in_off: usize = 0;
	1052	var out_off: usize = 0;
	1053
	1054	const batch_size: usize = comptime lcm(@as(i16, d), 8);
	1055	const in_batch_size: usize = comptime batch_size / d;
	1056	const out_batch_size: usize = comptime batch_size / 8;
	1057
	1058	const out_length: usize = comptime @divTrunc(N * d, 8);
	1059	comptime assert(out_length * 8 == d * N);
	1060	var out = [_]u8{0} ** out_length;
	1061
	1062	while (in_off < N) {
	1063	// First we compress into in.
	1064	var in: [in_batch_size]u16 = undefined;
	1065	inline for (0..in_batch_size) \|i\| {
	1066	// Compress_q(x, d) = ⌈(2ᵈ/q)x⌋ mod⁺ 2ᵈ
	1067	// = ⌊(2ᵈ/q)x+½⌋ mod⁺ 2ᵈ
	1068	// = ⌊((x << d) + q/2) / q⌋ mod⁺ 2ᵈ
	1069	// = DIV((x << d) + q/2, q) & ((1<<d) - 1)
	1070	const t = @as(u24, @intCast(p.cs[in_off + i])) << d;
	1071	// Division by invariant multiplication, equivalent to DIV(t + q/2, q).
	1072	// A division may not be a constant-time operation, even with a constant denominator.
	1073	// Here, side channels would leak information about the shared secret, see https://kyberslash.cr.yp.to
	1074	// Multiplication, on the other hand, is a constant-time operation on the CPUs we currently support.
	1075	comptime assert(d <= 11);
	1076	comptime assert(((20642679 * @as(u64, Q)) >> 36) == 1);
	1077	const u: u32 = @intCast((@as(u64, t + q_over_2) * 20642679) >> 36);
	1078	in[i] = @intCast(u & two_d_min_1);
	1079	}
	1080
	1081	// Now we pack the d-bit integers from `in' into out as bytes.
	1082	comptime var in_shift: usize = 0;
	1083	comptime var j: usize = 0;
	1084	comptime var i: usize = 0;
	1085	inline while (i < in_batch_size) : (j += 1) {
	1086	comptime var todo: usize = 8;
	1087	inline while (todo > 0) {
	1088	const out_shift = comptime 8 - todo;
	1089	out[out_off + j] \|= @as(u8, @truncate((in[i] >> in_shift) << out_shift));
	1090
	1091	const done = comptime @min(@min(d, todo), d - in_shift);
	1092	todo -= done;
	1093	in_shift += done;
	1094
	1095	if (in_shift == d) {
	1096	in_shift = 0;
	1097	i += 1;
	1098	}
	1099	}
	1100	}
	1101
	1102	in_off += in_batch_size;
	1103	out_off += out_batch_size;
	1104	}
	1105
	1106	return out;
	1107	}
	1108
	1109	// Set p to Decompress_q(m, d).
	1110	fn decompress(comptime d: u8, in: *const [compressedSize(d)]u8) Poly {
	1111	@setEvalBranchQuota(10000);
	1112	const inLen = comptime @divTrunc(N * d, 8);
	1113	comptime assert(inLen * 8 == d * N);
	1114	var ret: Poly = undefined;
	1115	var in_off: usize = 0;
	1116	var out_off: usize = 0;
	1117
	1118	const batch_size: usize = comptime lcm(@as(i16, d), 8);
	1119	const in_batch_size: usize = comptime batch_size / 8;
	1120	const out_batch_size: usize = comptime batch_size / d;
	1121
	1122	while (out_off < N) {
	1123	comptime var in_shift: usize = 0;
	1124	comptime var j: usize = 0;
	1125	comptime var i: usize = 0;
	1126	inline while (i < out_batch_size) : (i += 1) {
	1127	// First, unpack next coefficient.
	1128	comptime var todo = d;
	1129	var out: u16 = 0;
	1130
	1131	inline while (todo > 0) {
	1132	const out_shift = comptime d - todo;
	1133	const m = comptime (1 << d) - 1;
	1134	out \|= (@as(u16, in[in_off + j] >> in_shift) << out_shift) & m;
	1135
	1136	const done = comptime @min(@min(8, todo), 8 - in_shift);
	1137	todo -= done;
	1138	in_shift += done;
	1139
	1140	if (in_shift == 8) {
	1141	in_shift = 0;
	1142	j += 1;
	1143	}
	1144	}
	1145
	1146	// Decompress_q(x, d) = ⌈(q/2ᵈ)x⌋
	1147	// = ⌊(q/2ᵈ)x+½⌋
	1148	// = ⌊(qx + 2ᵈ⁻¹)/2ᵈ⌋
	1149	// = (qx + (1<<(d-1))) >> d
	1150	const qx = @as(u32, out) * @as(u32, Q);
	1151	ret.cs[out_off + i] = @as(i16, @intCast((qx + (1 << (d - 1))) >> d));
	1152	}
	1153
	1154	in_off += in_batch_size;
	1155	out_off += out_batch_size;
	1156	}
	1157
	1158	return ret;
	1159	}
	1160
	1161	// Returns the "pointwise" multiplication a o b.
	1162	//
	1163	// That is: invNTT(a o b) = invNTT(a) * invNTT(b). Assumes a and b are in
	1164	// Montgomery form. Products between coefficients of a and b must be strictly
	1165	// bounded in absolute value by 2¹⁵q. a o b will be in Montgomery form and
	1166	// bounded in absolute value by 2q.
	1167	fn mulHat(a: Poly, b: Poly) Poly {
	1168	// Recall from the discussion in ntt(), that a transformed polynomial is
	1169	// an element of ℤ_q[x]/(x²-ζ) x … x ℤ_q[x]/(x²+ζ¹²⁷);
	1170	// that is: 128 degree-one polynomials instead of simply 256 elements
	1171	// from ℤ_q as in the regular NTT. So instead of pointwise multiplication,
	1172	// we multiply the 128 pairs of degree-one polynomials modulo the
	1173	// right equation:
	1174	//
	1175	// (a₁ + a₂x)(b₁ + b₂x) = a₁b₁ + a₂b₂ζ' + (a₁b₂ + a₂b₁)x,
	1176	//
	1177	// where ζ' is the appropriate power of ζ.
	1178
	1179	var p: Poly = undefined;
	1180	var k: usize = 64;
	1181	var i: usize = 0;
	1182	while (i < N) : (i += 4) {
	1183	const z = @as(i32, zetas[k]);
	1184	k += 1;
	1185
	1186	const a1b1 = montReduce(@as(i32, a.cs[i + 1]) * @as(i32, b.cs[i + 1]));
	1187	const a0b0 = montReduce(@as(i32, a.cs[i]) * @as(i32, b.cs[i]));
	1188	const a1b0 = montReduce(@as(i32, a.cs[i + 1]) * @as(i32, b.cs[i]));
	1189	const a0b1 = montReduce(@as(i32, a.cs[i]) * @as(i32, b.cs[i + 1]));
	1190
	1191	p.cs[i] = montReduce(a1b1 * z) + a0b0;
	1192	p.cs[i + 1] = a0b1 + a1b0;
	1193
	1194	const a3b3 = montReduce(@as(i32, a.cs[i + 3]) * @as(i32, b.cs[i + 3]));
	1195	const a2b2 = montReduce(@as(i32, a.cs[i + 2]) * @as(i32, b.cs[i + 2]));
	1196	const a3b2 = montReduce(@as(i32, a.cs[i + 3]) * @as(i32, b.cs[i + 2]));
	1197	const a2b3 = montReduce(@as(i32, a.cs[i + 2]) * @as(i32, b.cs[i + 3]));
	1198
	1199	p.cs[i + 2] = a2b2 - montReduce(a3b3 * z);
	1200	p.cs[i + 3] = a2b3 + a3b2;
	1201	}
	1202
	1203	return p;
	1204	}
	1205
	1206	// Sample p from a centered binomial distribution with n=2η and p=½ - viz:
	1207	// coefficients are in {-η, …, η} with probabilities
	1208	//
	1209	// {ncr(0, 2η)/2^2η, ncr(1, 2η)/2^2η, …, ncr(2η,2η)/2^2η}
	1210	fn noise(comptime eta: u8, nonce: u8, seed: *const [32]u8) Poly {
	1211	var h = sha3.Shake256.init(.{});
	1212	const suffix: [1]u8 = .{nonce};
	1213	h.update(seed);
	1214	h.update(&suffix);
	1215
	1216	// The distribution at hand is exactly the same as that
	1217	// of (a₁ + a₂ + … + a_η) - (b₁ + … + b_η) where a_i,b_i~U(1).
	1218	// Thus we need 2η bits per coefficient.
	1219	const buf_len = comptime 2 * eta * N / 8;
	1220	var buf: [buf_len]u8 = undefined;
	1221	h.squeeze(&buf);
	1222
	1223	// buf is interpreted as a₁…a_ηb₁…b_ηa₁…a_ηb₁…b_η…. We process
	1224	// multiple coefficients in one batch.
	1225
	1226	const T = switch (builtin.target.cpu.arch) {
	1227	.x86_64, .x86 => u32, // Generates better code on Intel CPUs
	1228	else => u64, // u128 might be faster on some other CPUs.
	1229	};
	1230
	1231	comptime var batch_count: usize = undefined;
	1232	comptime var batch_bytes: usize = undefined;
	1233	comptime var mask: T = 0;
	1234	comptime {
	1235	batch_count = @bitSizeOf(T) / @as(usize, 2 * eta);
	1236	while (@rem(N, batch_count) != 0 and batch_count > 0) : (batch_count -= 1) {}
	1237	assert(batch_count > 0);
	1238	assert(@rem(2 * eta * batch_count, 8) == 0);
	1239	batch_bytes = 2 * eta * batch_count / 8;
	1240
	1241	for (0..2 * eta * batch_count) \|_\| {
	1242	mask <<= eta;
	1243	mask \|= 1;
	1244	}
	1245	}
	1246
	1247	var ret: Poly = undefined;
	1248	for (0..comptime N / batch_count) \|i\| {
	1249	// Read coefficients into t. In the case of η=3,
	1250	// we have t = a₁ + 2a₂ + 4a₃ + 8b₁ + 16b₂ + …
	1251	var t: T = 0;
	1252	inline for (0..batch_bytes) \|j\| {
	1253	t \|= @as(T, buf[batch_bytes * i + j]) << (8 * j);
	1254	}
	1255
	1256	// Accumelate `a's and `b's together by masking them out, shifting
	1257	// and adding. For η=3, we have d = a₁ + a₂ + a₃ + 8(b₁ + b₂ + b₃) + …
	1258	var d: T = 0;
	1259	inline for (0..eta) \|j\| {
	1260	d += (t >> j) & mask;
	1261	}
	1262
	1263	// Extract each a and b separately and set coefficient in polynomial.
	1264	inline for (0..batch_count) \|j\| {
	1265	const mask2 = comptime (1 << eta) - 1;
	1266	const a = @as(i16, @intCast((d >> (comptime (2 * j * eta))) & mask2));
	1267	const b = @as(i16, @intCast((d >> (comptime ((2 * j + 1) * eta))) & mask2));
	1268	ret.cs[batch_count * i + j] = a - b;
	1269	}
	1270	}
	1271
	1272	return ret;
	1273	}
	1274
	1275	// Sample p uniformly from the given seed and x and y coordinates.
	1276	fn uniform(seed: [32]u8, x: u8, y: u8) Poly {
	1277	var h = sha3.Shake128.init(.{});
	1278	const suffix: [2]u8 = .{ x, y };
	1279	h.update(&seed);
	1280	h.update(&suffix);
	1281
	1282	const buf_len = sha3.Shake128.block_length; // rate SHAKE-128
	1283	var buf: [buf_len]u8 = undefined;
	1284
	1285	var ret: Poly = undefined;
	1286	var i: usize = 0; // index into ret.cs
	1287	outer: while (true) {
	1288	h.squeeze(&buf);
	1289
	1290	var j: usize = 0; // index into buf
	1291	while (j < buf_len) : (j += 3) {
	1292	const b0 = @as(u16, buf[j]);
	1293	const b1 = @as(u16, buf[j + 1]);
	1294	const b2 = @as(u16, buf[j + 2]);
	1295
	1296	const ts: [2]u16 = .{
	1297	b0 \| ((b1 & 0xf) << 8),
	1298	(b1 >> 4) \| (b2 << 4),
	1299	};
	1300
	1301	inline for (ts) \|t\| {
	1302	if (t < Q) {
	1303	ret.cs[i] = @as(i16, @intCast(t));
	1304	i += 1;
	1305
	1306	if (i == N) {
	1307	break :outer;
	1308	}
	1309	}
	1310	}
	1311	}
	1312	}
	1313
	1314	return ret;
	1315	}
	1316
	1317	// Packs p.
	1318	//
	1319	// Assumes p is normalized (and not just Barrett reduced).
	1320	fn toBytes(p: Poly) [bytes_length]u8 {
	1321	var ret: [bytes_length]u8 = undefined;
	1322	for (0..comptime N / 2) \|i\| {
	1323	const t0 = @as(u16, @intCast(p.cs[2 * i]));
	1324	const t1 = @as(u16, @intCast(p.cs[2 * i + 1]));
	1325	ret[3 * i] = @as(u8, @truncate(t0));
	1326	ret[3 * i + 1] = @as(u8, @truncate((t0 >> 8) \| (t1 << 4)));
	1327	ret[3 * i + 2] = @as(u8, @truncate(t1 >> 4));
	1328	}
	1329	return ret;
	1330	}
	1331
	1332	// Unpacks a Poly from buf.
	1333	//
	1334	// p will not be normalized; instead 0 ≤ p[i] < 4096.
	1335	fn fromBytes(buf: *const [bytes_length]u8) Poly {
	1336	var ret: Poly = undefined;
	1337	for (0..comptime N / 2) \|i\| {
	1338	const b0 = @as(i16, buf[3 * i]);
	1339	const b1 = @as(i16, buf[3 * i + 1]);
	1340	const b2 = @as(i16, buf[3 * i + 2]);
	1341	ret.cs[2 * i] = b0 \| ((b1 & 0xf) << 8);
	1342	ret.cs[2 * i + 1] = (b1 >> 4) \| b2 << 4;
	1343	}
	1344	return ret;
	1345	}
	1346	};
	1347
	1348	// A vector of K polynomials.
	1349	fn Vec(comptime K: u8) type {
	1350	return struct {
	1351	ps: [K]Poly,
	1352
	1353	const Self = @This();
	1354	const bytes_length = K * Poly.bytes_length;
	1355
	1356	fn compressedSize(comptime d: u8) usize {
	1357	return Poly.compressedSize(d) * K;
	1358	}
	1359
	1360	fn ntt(a: Self) Self {
	1361	var ret: Self = undefined;
	1362	for (0..K) \|i\| {
	1363	ret.ps[i] = a.ps[i].ntt();
	1364	}
	1365	return ret;
	1366	}
	1367
	1368	fn invNTT(a: Self) Self {
	1369	var ret: Self = undefined;
	1370	for (0..K) \|i\| {
	1371	ret.ps[i] = a.ps[i].invNTT();
	1372	}
	1373	return ret;
	1374	}
	1375
	1376	fn normalize(a: Self) Self {
	1377	var ret: Self = undefined;
	1378	for (0..K) \|i\| {
	1379	ret.ps[i] = a.ps[i].normalize();
	1380	}
	1381	return ret;
	1382	}
	1383
	1384	fn barrettReduce(a: Self) Self {
	1385	var ret: Self = undefined;
	1386	for (0..K) \|i\| {
	1387	ret.ps[i] = a.ps[i].barrettReduce();
	1388	}
	1389	return ret;
	1390	}
	1391
	1392	fn add(a: Self, b: Self) Self {
	1393	var ret: Self = undefined;
	1394	for (0..K) \|i\| {
	1395	ret.ps[i] = a.ps[i].add(b.ps[i]);
	1396	}
	1397	return ret;
	1398	}
	1399
	1400	fn sub(a: Self, b: Self) Self {
	1401	var ret: Self = undefined;
	1402	for (0..K) \|i\| {
	1403	ret.ps[i] = a.ps[i].sub(b.ps[i]);
	1404	}
	1405	return ret;
	1406	}
	1407
	1408	// Samples v[i] from centered binomial distribution with the given η,
	1409	// seed and nonce+i.
	1410	fn noise(comptime eta: u8, nonce: u8, seed: *const [32]u8) Self {
	1411	var ret: Self = undefined;
	1412	for (0..K) \|i\| {
	1413	ret.ps[i] = Poly.noise(eta, nonce + @as(u8, @intCast(i)), seed);
	1414	}
	1415	return ret;
	1416	}
	1417
	1418	// Sets p to the inner product of a and b using "pointwise" multiplication.
	1419	//
	1420	// See MulHat() and NTT() for a description of the multiplication.
	1421	// Assumes a and b are in Montgomery form. p will be in Montgomery form,
	1422	// and its coefficients will be bounded in absolute value by 2kq.
	1423	// If a and b are not in Montgomery form, then the action is the same
	1424	// as "pointwise" multiplication followed by multiplying by R⁻¹, the inverse
	1425	// of the Montgomery factor.
	1426	fn dotHat(a: Self, b: Self) Poly {
	1427	var ret: Poly = Poly.zero;
	1428	for (0..K) \|i\| {
	1429	ret = ret.add(a.ps[i].mulHat(b.ps[i]));
	1430	}
	1431	return ret;
	1432	}
	1433
	1434	fn compress(v: Self, comptime d: u8) [compressedSize(d)]u8 {
	1435	const cs = comptime Poly.compressedSize(d);
	1436	var ret: [compressedSize(d)]u8 = undefined;
	1437	inline for (0..K) \|i\| {
	1438	ret[i * cs .. (i + 1) * cs].* = v.ps[i].compress(d);
	1439	}
	1440	return ret;
	1441	}
	1442
	1443	fn decompress(comptime d: u8, buf: *const [compressedSize(d)]u8) Self {
	1444	const cs = comptime Poly.compressedSize(d);
	1445	var ret: Self = undefined;
	1446	inline for (0..K) \|i\| {
	1447	ret.ps[i] = Poly.decompress(d, buf[i * cs .. (i + 1) * cs]);
	1448	}
	1449	return ret;
	1450	}
	1451
	1452	/// Serializes the key into a byte array.
	1453	fn toBytes(v: Self) [bytes_length]u8 {
	1454	var ret: [bytes_length]u8 = undefined;
	1455	inline for (0..K) \|i\| {
	1456	ret[i * Poly.bytes_length .. (i + 1) * Poly.bytes_length].* = v.ps[i].toBytes();
	1457	}
	1458	return ret;
	1459	}
	1460
	1461	/// Deserializes the key from a byte array.
	1462	fn fromBytes(buf: *const [bytes_length]u8) Self {
	1463	var ret: Self = undefined;
	1464	inline for (0..K) \|i\| {
	1465	ret.ps[i] = Poly.fromBytes(
	1466	buf[i * Poly.bytes_length .. (i + 1) * Poly.bytes_length],
	1467	);
	1468	}
	1469	return ret;
	1470	}
	1471	};
	1472	}
	1473
	1474	// A matrix of K vectors
	1475	fn Mat(comptime K: u8) type {
	1476	return struct {
	1477	const Self = @This();
	1478	vs: [K]Vec(K),
	1479
	1480	fn uniform(seed: [32]u8, comptime transposed: bool) Self {
	1481	var ret: Self = undefined;
	1482	var i: u8 = 0;
	1483	while (i < K) : (i += 1) {
	1484	var j: u8 = 0;
	1485	while (j < K) : (j += 1) {
	1486	ret.vs[i].ps[j] = Poly.uniform(
	1487	seed,
	1488	if (transposed) i else j,
	1489	if (transposed) j else i,
	1490	);
	1491	}
	1492	}
	1493	return ret;
	1494	}
	1495
	1496	// Returns transpose of A
	1497	fn transpose(m: Self) Self {
	1498	var ret: Self = undefined;
	1499	for (0..K) \|i\| {
	1500	for (0..K) \|j\| {
	1501	ret.vs[i].ps[j] = m.vs[j].ps[i];
	1502	}
	1503	}
	1504	return ret;
	1505	}
	1506	};
	1507	}
	1508
	1509	// Returns `true` if a ≠ b.
	1510	fn ctneq(comptime len: usize, a: [len]u8, b: [len]u8) u1 {
	1511	return 1 - @intFromBool(crypto.utils.timingSafeEql([len]u8, a, b));
	1512	}
	1513
	1514	// Copy src into dst given b = 1.
	1515	fn cmov(comptime len: usize, dst: *[len]u8, src: [len]u8, b: u1) void {
	1516	const mask = @as(u8, 0) -% b;
	1517	for (0..len) \|i\| {
	1518	dst[i] ^= mask & (dst[i] ^ src[i]);
	1519	}
	1520	}
	1521
	1522	test "MulHat" {
	1523	var rnd = RndGen.init(0);
	1524
	1525	for (0..100) \|_\| {
	1526	const a = Poly.randAbsLeqQ(&rnd);
	1527	const b = Poly.randAbsLeqQ(&rnd);
	1528
	1529	const p2 = a.ntt().mulHat(b.ntt()).barrettReduce().invNTT().normalize();
	1530	var p: Poly = undefined;
	1531
	1532	@memset(&p.cs, 0);
	1533
	1534	for (0..N) \|i\| {
	1535	for (0..N) \|j\| {
	1536	var v = montReduce(@as(i32, a.cs[i]) * @as(i32, b.cs[j]));
	1537	var k = i + j;
	1538	if (k >= N) {
	1539	// Recall Xᴺ = -1.
	1540	k -= N;
	1541	v = -v;
	1542	}
	1543	p.cs[k] = feBarrettReduce(v + p.cs[k]);
	1544	}
	1545	}
	1546
	1547	p = p.toMont().normalize();
	1548
	1549	try testing.expectEqual(p, p2);
	1550	}
	1551	}
	1552
	1553	test "NTT" {
	1554	var rnd = RndGen.init(0);
	1555
	1556	for (0..1000) \|_\| {
	1557	var p = Poly.randAbsLeqQ(&rnd);
	1558	const q = p.toMont().normalize();
	1559	p = p.ntt();
	1560
	1561	for (0..N) \|i\| {
	1562	try testing.expect(p.cs[i] <= 7 * Q and -7 * Q <= p.cs[i]);
	1563	}
	1564
	1565	p = p.normalize().invNTT();
	1566	for (0..N) \|i\| {
	1567	try testing.expect(p.cs[i] <= Q and -Q <= p.cs[i]);
	1568	}
	1569
	1570	p = p.normalize();
	1571
	1572	try testing.expectEqual(p, q);
	1573	}
	1574	}
	1575
	1576	test "Compression" {
	1577	var rnd = RndGen.init(0);
	1578	inline for (.{ 1, 4, 5, 10, 11 }) \|d\| {
	1579	for (0..1000) \|_\| {
	1580	const p = Poly.randNormalized(&rnd);
	1581	const pp = p.compress(d);
	1582	const pq = Poly.decompress(d, &pp).compress(d);
	1583	try testing.expectEqual(pp, pq);
	1584	}
	1585	}
	1586	}
	1587
	1588	test "noise" {
	1589	var seed: [32]u8 = undefined;
	1590	for (&seed, 0..) \|*s, i\| {
	1591	s.* = @as(u8, @intCast(i));
	1592	}
	1593	try testing.expectEqual(Poly.noise(3, 37, &seed).cs, .{
	1594	0, 0, 1, -1, 0, 2, 0, -1, -1, 3, 0, 1, -2, -2, 0, 1, -2,
	1595	1, 0, -2, 3, 0, 0, 0, 1, 3, 1, 1, 2, 1, -1, -1, -1, 0,
	1596	1, 0, 1, 0, 2, 0, 1, -2, 0, -1, -1, -2, 1, -1, -1, 2, -1,
	1597	1, 1, 2, -3, -1, -1, 0, 0, 0, 0, 1, -1, -2, -2, 0, -2, 0,
	1598	0, 0, 1, 0, -1, -1, 1, -2, 2, 0, 0, 2, -2, 0, 1, 0, 1,
	1599	1, 1, 0, 1, -2, -1, -2, -1, 1, 0, 0, 0, 0, 0, 1, 0, -1,
	1600	-1, 0, -1, 1, 0, 1, 0, -1, -1, 0, -2, 2, 0, -2, 1, -1, 0,
	1601	1, -1, -1, 2, 1, 0, 0, -2, -1, 2, 0, 0, 0, -1, -1, 3, 1,
	1602	0, 1, 0, 1, 0, 2, 1, 0, 0, 1, 0, 1, 0, 0, -1, -1, -1,
	1603	0, 1, 3, 1, 0, 1, 0, 1, -1, -1, -1, -1, 0, 0, -2, -1, -1,
	1604	2, 0, 1, 0, 1, 0, 2, -2, 0, 1, 1, -3, -1, -2, -1, 0, 1,
	1605	0, 1, -2, 2, 2, 1, 1, 0, -1, 0, -1, -1, 1, 0, -1, 2, 1,
	1606	-1, 1, 2, -2, 1, 2, 0, 1, 2, 1, 0, 0, 2, 1, 2, 1, 0,
	1607	2, 1, 0, 0, -1, -1, 1, -1, 0, 1, -1, 2, 2, 0, 0, -1, 1,
	1608	1, 1, 1, 0, 0, -2, 0, -1, 1, 2, 0, 0, 1, 1, -1, 1, 0,
	1609	1,
	1610	});
	1611	try testing.expectEqual(Poly.noise(2, 37, &seed).cs, .{
	1612	1, 0, 1, -1, -1, -2, -1, -1, 2, 0, -1, 0, 0, -1,
	1613	1, 1, -1, 1, 0, 2, -2, 0, 1, 2, 0, 0, -1, 1,
	1614	0, -1, 1, -1, 1, 2, 1, 1, 0, -1, 1, -1, -2, -1,
	1615	1, -1, -1, -1, 2, -1, -1, 0, 0, 1, 1, -1, 1, 1,
	1616	1, 1, -1, -2, 0, 1, 0, 0, 2, 1, -1, 2, 0, 0,
	1617	1, 1, 0, -1, 0, 0, -1, -1, 2, 0, 1, -1, 2, -1,
	1618	-1, -1, -1, 0, -2, 0, 2, 1, 0, 0, 0, -1, 0, 0,
	1619	0, -1, -1, 0, -1, -1, 0, -1, 0, 0, -2, 1, 1, 0,
	1620	1, 0, 1, 0, 1, 1, -1, 2, 0, 1, -1, 1, 2, 0,
	1621	0, 0, 0, -1, -1, -1, 0, 1, 0, -1, 2, 0, 0, 1,
	1622	1, 1, 0, 1, -1, 1, 2, 1, 0, 2, -1, 1, -1, -2,
	1623	-1, -2, -1, 1, 0, -2, -2, -1, 1, 0, 0, 0, 0, 1,
	1624	0, 0, 0, 2, 2, 0, 1, 0, -1, -1, 0, 2, 0, 0,
	1625	-2, 1, 0, 2, 1, -1, -2, 0, 0, -1, 1, 1, 0, 0,
	1626	2, 0, 1, 1, -2, 1, -2, 1, 1, 0, 2, 0, -1, 0,
	1627	-1, 0, 1, 2, 0, 1, 0, -2, 1, -2, -2, 1, -1, 0,
	1628	-1, 1, 1, 0, 0, 0, 1, 0, -1, 1, 1, 0, 0, 0,
	1629	0, 1, 0, 1, -1, 0, 1, -1, -1, 2, 0, 0, 1, -1,
	1630	0, 1, -1, 0,
	1631	});
	1632	}
	1633
	1634	test "uniform sampling" {
	1635	var seed: [32]u8 = undefined;
	1636	for (&seed, 0..) \|*s, i\| {
	1637	s.* = @as(u8, @intCast(i));
	1638	}
	1639	try testing.expectEqual(Poly.uniform(seed, 1, 0).cs, .{
	1640	797, 993, 161, 6, 2608, 2385, 2096, 2661, 1676, 247, 2440,
	1641	342, 634, 194, 1570, 2848, 986, 684, 3148, 3208, 2018, 351,
	1642	2288, 612, 1394, 170, 1521, 3119, 58, 596, 2093, 1549, 409,
	1643	2156, 1934, 1730, 1324, 388, 446, 418, 1719, 2202, 1812, 98,
	1644	1019, 2369, 214, 2699, 28, 1523, 2824, 273, 402, 2899, 246,
	1645	210, 1288, 863, 2708, 177, 3076, 349, 44, 949, 854, 1371,
	1646	957, 292, 2502, 1617, 1501, 254, 7, 1761, 2581, 2206, 2655,
	1647	1211, 629, 1274, 2358, 816, 2766, 2115, 2985, 1006, 2433, 856,
	1648	2596, 3192, 1, 1378, 2345, 707, 1891, 1669, 536, 1221, 710,
	1649	2511, 120, 1176, 322, 1897, 2309, 595, 2950, 1171, 801, 1848,
	1650	695, 2912, 1396, 1931, 1775, 2904, 893, 2507, 1810, 2873, 253,
	1651	1529, 1047, 2615, 1687, 831, 1414, 965, 3169, 1887, 753, 3246,
	1652	1937, 115, 2953, 586, 545, 1621, 1667, 3187, 1654, 1988, 1857,
	1653	512, 1239, 1219, 898, 3106, 391, 1331, 2228, 3169, 586, 2412,
	1654	845, 768, 156, 662, 478, 1693, 2632, 573, 2434, 1671, 173,
	1655	969, 364, 1663, 2701, 2169, 813, 1000, 1471, 720, 2431, 2530,
	1656	3161, 733, 1691, 527, 2634, 335, 26, 2377, 1707, 767, 3020,
	1657	950, 502, 426, 1138, 3208, 2607, 2389, 44, 1358, 1392, 2334,
	1658	875, 2097, 173, 1697, 2578, 942, 1817, 974, 1165, 2853, 1958,
	1659	2973, 3282, 271, 1236, 1677, 2230, 673, 1554, 96, 242, 1729,
	1660	2518, 1884, 2272, 71, 1382, 924, 1807, 1610, 456, 1148, 2479,
	1661	2152, 238, 2208, 2329, 713, 1175, 1196, 757, 1078, 3190, 3169,
	1662	708, 3117, 154, 1751, 3225, 1364, 154, 23, 2842, 1105, 1419,
	1663	79, 5, 2013,
	1664	});
	1665	}
	1666
	1667	test "Polynomial packing" {
	1668	var rnd = RndGen.init(0);
	1669
	1670	for (0..1000) \|_\| {
	1671	const p = Poly.randNormalized(&rnd);
	1672	try testing.expectEqual(Poly.fromBytes(&p.toBytes()), p);
	1673	}
	1674	}
	1675
	1676	test "Test inner PKE" {
	1677	var seed: [32]u8 = undefined;
	1678	var pt: [32]u8 = undefined;
	1679	for (&seed, &pt, 0..) \|s, p, i\| {
	1680	s.* = @as(u8, @intCast(i));
	1681	p.* = @as(u8, @intCast(i + 32));
	1682	}
	1683	inline for (modes) \|mode\| {
	1684	for (0..100) \|i\| {
	1685	var pk: mode.InnerPk = undefined;
	1686	var sk: mode.InnerSk = undefined;
	1687	seed[0] = @as(u8, @intCast(i));
	1688	mode.innerKeyFromSeed(seed, &pk, &sk);
	1689	for (0..10) \|j\| {
	1690	seed[1] = @as(u8, @intCast(j));
	1691	try testing.expectEqual(sk.decrypt(&pk.encrypt(&pt, &seed)), pt);
	1692	}
	1693	}
	1694	}
	1695	}
	1696
	1697	test "Test happy flow" {
	1698	var seed: [64]u8 = undefined;
	1699	for (&seed, 0..) \|*s, i\| {
	1700	s.* = @as(u8, @intCast(i));
	1701	}
	1702	inline for (modes) \|mode\| {
	1703	for (0..100) \|i\| {
	1704	seed[0] = @as(u8, @intCast(i));
	1705	const kp = try mode.KeyPair.create(seed);
	1706	const sk = try mode.SecretKey.fromBytes(&kp.secret_key.toBytes());
	1707	try testing.expectEqual(sk, kp.secret_key);
	1708	const pk = try mode.PublicKey.fromBytes(&kp.public_key.toBytes());
	1709	try testing.expectEqual(pk, kp.public_key);
	1710	for (0..10) \|j\| {
	1711	seed[1] = @as(u8, @intCast(j));
	1712	const e = pk.encaps(seed[0..32].*);
	1713	try testing.expectEqual(e.shared_secret, try sk.decaps(&e.ciphertext));
	1714	}
	1715	}
	1716	}
	1717	}
	1718
	1719	// Code to test NIST Known Answer Tests (KAT), see PQCgenKAT.c.
	1720
	1721	const sha2 = crypto.hash.sha2;
	1722
	1723	test "NIST KAT test" {
	1724	inline for (.{
	1725	.{ kyber_d00.Kyber512, "e9c2bd37133fcb40772f81559f14b1f58dccd1c816701be9ba6214d43baf4547" },
	1726	.{ kyber_d00.Kyber1024, "89248f2f33f7f4f7051729111f3049c409a933ec904aedadf035f30fa5646cd5" },
	1727	.{ kyber_d00.Kyber768, "a1e122cad3c24bc51622e4c242d8b8acbcd3f618fee4220400605ca8f9ea02c2" },
	1728	}) \|modeHash\| {
	1729	const mode = modeHash[0];
	1730	var seed: [48]u8 = undefined;
	1731	for (&seed, 0..) \|*s, i\| {
	1732	s.* = @as(u8, @intCast(i));
	1733	}
	1734	var f = sha2.Sha256.init(.{});
	1735	const fw = f.writer();
	1736	var g = NistDRBG.init(seed);
	1737	try std.fmt.format(fw, "# {s}\n\n", .{mode.name});
	1738	for (0..100) \|i\| {
	1739	g.fill(&seed);
	1740	try std.fmt.format(fw, "count = {}\n", .{i});
	1741	try std.fmt.format(fw, "seed = {s}\n", .{std.fmt.fmtSliceHexUpper(&seed)});
	1742	var g2 = NistDRBG.init(seed);
	1743
	1744	// This is not equivalent to g2.fill(kseed[:]). As the reference
	1745	// implementation calls randombytes twice generating the keypair,
	1746	// we have to do that as well.
	1747	var kseed: [64]u8 = undefined;
	1748	var eseed: [32]u8 = undefined;
	1749	g2.fill(kseed[0..32]);
	1750	g2.fill(kseed[32..64]);
	1751	g2.fill(&eseed);
	1752	const kp = try mode.KeyPair.create(kseed);
	1753	const e = kp.public_key.encaps(eseed);
	1754	const ss2 = try kp.secret_key.decaps(&e.ciphertext);
	1755	try testing.expectEqual(ss2, e.shared_secret);
	1756	try std.fmt.format(fw, "pk = {s}\n", .{std.fmt.fmtSliceHexUpper(&kp.public_key.toBytes())});
	1757	try std.fmt.format(fw, "sk = {s}\n", .{std.fmt.fmtSliceHexUpper(&kp.secret_key.toBytes())});
	1758	try std.fmt.format(fw, "ct = {s}\n", .{std.fmt.fmtSliceHexUpper(&e.ciphertext)});
	1759	try std.fmt.format(fw, "ss = {s}\n\n", .{std.fmt.fmtSliceHexUpper(&e.shared_secret)});
	1760	}
	1761
	1762	var out: [32]u8 = undefined;
	1763	f.final(&out);
	1764	var outHex: [64]u8 = undefined;
	1765	_ = try std.fmt.bufPrint(&outHex, "{s}", .{std.fmt.fmtSliceHexLower(&out)});
	1766	try testing.expectEqual(outHex, modeHash[1].*);
	1767	}
	1768	}
	1769
	1770	const NistDRBG = struct {
	1771	key: [32]u8,
	1772	v: [16]u8,
	1773
	1774	fn incV(g: *NistDRBG) void {
	1775	var j: usize = 15;
	1776	while (j >= 0) : (j -= 1) {
	1777	if (g.v[j] == 255) {
	1778	g.v[j] = 0;
	1779	} else {
	1780	g.v[j] += 1;
	1781	break;
	1782	}
	1783	}
	1784	}
	1785
	1786	// AES256_CTR_DRBG_Update(pd, &g.key, &g.v).
	1787	fn update(g: *NistDRBG, pd: ?[48]u8) void {
	1788	var buf: [48]u8 = undefined;
	1789	const ctx = crypto.core.aes.Aes256.initEnc(g.key);
	1790	var i: usize = 0;
	1791	while (i < 3) : (i += 1) {
	1792	g.incV();
	1793	var block: [16]u8 = undefined;
	1794	ctx.encrypt(&block, &g.v);
	1795	buf[i * 16 ..][0..16].* = block;
	1796	}
	1797	if (pd) \|p\| {
	1798	for (&buf, p) \|*b, x\| {
	1799	b.* ^= x;
	1800	}
	1801	}
	1802	g.key = buf[0..32].*;
	1803	g.v = buf[32..48].*;
	1804	}
	1805
	1806	// randombytes.
	1807	fn fill(g: *NistDRBG, out: []u8) void {
	1808	var block: [16]u8 = undefined;
	1809	var dst = out;
	1810
	1811	const ctx = crypto.core.aes.Aes256.initEnc(g.key);
	1812	while (dst.len > 0) {
	1813	g.incV();
	1814	ctx.encrypt(&block, &g.v);
	1815	if (dst.len < 16) {
	1816	@memcpy(dst, block[0..dst.len]);
	1817	break;
	1818	}
	1819	dst[0..block.len].* = block;
	1820	dst = dst[16..dst.len];
	1821	}
	1822	g.update(null);
	1823	}
	1824
	1825	fn init(seed: [48]u8) NistDRBG {
	1826	var ret: NistDRBG = .{ .key = .{0} 32, .v = .{0} 16 };
	1827	ret.update(seed);
	1828	return ret;
	1829	}
	1830	};