| ... | @@ -0,0 +1,1780 @@ |
| 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 negligable 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.rand.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 | mem.copy(u8, &ret.hpk, buf[s .. s + h_length]); |
| 327 | s += h_length; |
| 328 | mem.copy(u8, &ret.z, buf[s .. s + 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 | mem.copy(u8, &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 | mem.copy(u8, &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 | mem.copy(u8, &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 | mem.copy(u8, &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[@intCast(usize, 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 = @intCast(i16, @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(a*b) and so the result is ok. |
| 641 | const m = @truncate(i16, @truncate(i32, x *% qInv)); |
| 642 | |
| 643 | // Note that x - m q is divisable 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(i16, @truncate(u16, @bitCast(u32, 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(@intCast(i16, 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 -% @intCast(i16, (@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(@intCast(i16, x)); |
| 713 | const y2 = @mod(@intCast(i16, 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(@intCast(i16, x)); |
| 733 | var y2 = @intCast(i16, x); |
| 734 | if (@intCast(i16, 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 = @intCast(i16, mpow(@as(i32, zeta), @bitReverse(@intCast(u7, 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[@intCast(usize, i)] = feBarrettReduce(p.cs[@intCast(usize, 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 = @intCast(u32, p.cs[in_off + i]) << d; |
| 1024 | in[i] = @intCast(u16, @divFloor(t + q_over_2, Q) & two_d_min_1); |
| 1025 | } |
| 1026 | |
| 1027 | // Now we pack the d-bit integers from `in' into out as bytes. |
| 1028 | comptime var in_shift: usize = 0; |
| 1029 | comptime var j: usize = 0; |
| 1030 | comptime var i: usize = 0; |
| 1031 | inline while (i < in_batch_size) : (j += 1) { |
| 1032 | comptime var todo: usize = 8; |
| 1033 | inline while (todo > 0) { |
| 1034 | const out_shift = comptime 8 - todo; |
| 1035 | out[out_off + j] |= @truncate(u8, (in[i] >> in_shift) << out_shift); |
| 1036 | |
| 1037 | const done = comptime @min(@min(d, todo), d - in_shift); |
| 1038 | todo -= done; |
| 1039 | in_shift += done; |
| 1040 | |
| 1041 | if (in_shift == d) { |
| 1042 | in_shift = 0; |
| 1043 | i += 1; |
| 1044 | } |
| 1045 | } |
| 1046 | } |
| 1047 | |
| 1048 | in_off += in_batch_size; |
| 1049 | out_off += out_batch_size; |
| 1050 | } |
| 1051 | |
| 1052 | return out; |
| 1053 | } |
| 1054 | |
| 1055 | // Set p to Decompress_q(m, d). |
| 1056 | fn decompress(comptime d: u8, in: *const [compressedSize(d)]u8) Poly { |
| 1057 | @setEvalBranchQuota(10000); |
| 1058 | const inLen = comptime @divTrunc(N * d, 8); |
| 1059 | comptime assert(inLen * 8 == d * N); |
| 1060 | var ret: Poly = undefined; |
| 1061 | var in_off: usize = 0; |
| 1062 | var out_off: usize = 0; |
| 1063 | |
| 1064 | const batch_size: usize = comptime lcm(@as(i16, d), 8); |
| 1065 | const in_batch_size: usize = comptime batch_size / 8; |
| 1066 | const out_batch_size: usize = comptime batch_size / d; |
| 1067 | |
| 1068 | while (out_off < N) { |
| 1069 | comptime var in_shift: usize = 0; |
| 1070 | comptime var j: usize = 0; |
| 1071 | comptime var i: usize = 0; |
| 1072 | inline while (i < out_batch_size) : (i += 1) { |
| 1073 | // First, unpack next coefficient. |
| 1074 | comptime var todo = d; |
| 1075 | var out: u16 = 0; |
| 1076 | |
| 1077 | inline while (todo > 0) { |
| 1078 | const out_shift = comptime d - todo; |
| 1079 | const m = comptime (1 << d) - 1; |
| 1080 | out |= (@as(u16, in[in_off + j] >> in_shift) << out_shift) & m; |
| 1081 | |
| 1082 | const done = comptime @min(@min(8, todo), 8 - in_shift); |
| 1083 | todo -= done; |
| 1084 | in_shift += done; |
| 1085 | |
| 1086 | if (in_shift == 8) { |
| 1087 | in_shift = 0; |
| 1088 | j += 1; |
| 1089 | } |
| 1090 | } |
| 1091 | |
| 1092 | // Decompress_q(x, d) = ⌈(q/2ᵈ)x⌋ |
| 1093 | // = ⌊(q/2ᵈ)x+½⌋ |
| 1094 | // = ⌊(qx + 2ᵈ⁻¹)/2ᵈ⌋ |
| 1095 | // = (qx + (1<<(d-1))) >> d |
| 1096 | const qx = @as(u32, out) * @as(u32, Q); |
| 1097 | ret.cs[out_off + i] = @intCast(i16, (qx + (1 << (d - 1))) >> d); |
| 1098 | } |
| 1099 | |
| 1100 | in_off += in_batch_size; |
| 1101 | out_off += out_batch_size; |
| 1102 | } |
| 1103 | |
| 1104 | return ret; |
| 1105 | } |
| 1106 | |
| 1107 | // Returns the "pointwise" multiplication a o b. |
| 1108 | // |
| 1109 | // That is: invNTT(a o b) = invNTT(a) * invNTT(b). Assumes a and b are in |
| 1110 | // Montgomery form. Products between coefficients of a and b must be strictly |
| 1111 | // bounded in absolute value by 2¹⁵q. a o b will be in Montgomery form and |
| 1112 | // bounded in absolute value by 2q. |
| 1113 | fn mulHat(a: Poly, b: Poly) Poly { |
| 1114 | // Recall from the discussion in ntt(), that a transformed polynomial is |
| 1115 | // an element of ℤ_q[x]/(x²-ζ) x … x ℤ_q[x]/(x²+ζ¹²⁷); |
| 1116 | // that is: 128 degree-one polynomials instead of simply 256 elements |
| 1117 | // from ℤ_q as in the regular NTT. So instead of pointwise multiplication, |
| 1118 | // we multiply the 128 pairs of degree-one polynomials modulo the |
| 1119 | // right equation: |
| 1120 | // |
| 1121 | // (a₁ + a₂x)(b₁ + b₂x) = a₁b₁ + a₂b₂ζ' + (a₁b₂ + a₂b₁)x, |
| 1122 | // |
| 1123 | // where ζ' is the appropriate power of ζ. |
| 1124 | |
| 1125 | var p: Poly = undefined; |
| 1126 | var k: usize = 64; |
| 1127 | var i: usize = 0; |
| 1128 | while (i < N) : (i += 4) { |
| 1129 | const z = @as(i32, zetas[k]); |
| 1130 | k += 1; |
| 1131 | |
| 1132 | const a1b1 = montReduce(@as(i32, a.cs[i + 1]) * @as(i32, b.cs[i + 1])); |
| 1133 | const a0b0 = montReduce(@as(i32, a.cs[i]) * @as(i32, b.cs[i])); |
| 1134 | const a1b0 = montReduce(@as(i32, a.cs[i + 1]) * @as(i32, b.cs[i])); |
| 1135 | const a0b1 = montReduce(@as(i32, a.cs[i]) * @as(i32, b.cs[i + 1])); |
| 1136 | |
| 1137 | p.cs[i] = montReduce(a1b1 * z) + a0b0; |
| 1138 | p.cs[i + 1] = a0b1 + a1b0; |
| 1139 | |
| 1140 | const a3b3 = montReduce(@as(i32, a.cs[i + 3]) * @as(i32, b.cs[i + 3])); |
| 1141 | const a2b2 = montReduce(@as(i32, a.cs[i + 2]) * @as(i32, b.cs[i + 2])); |
| 1142 | const a3b2 = montReduce(@as(i32, a.cs[i + 3]) * @as(i32, b.cs[i + 2])); |
| 1143 | const a2b3 = montReduce(@as(i32, a.cs[i + 2]) * @as(i32, b.cs[i + 3])); |
| 1144 | |
| 1145 | p.cs[i + 2] = a2b2 - montReduce(a3b3 * z); |
| 1146 | p.cs[i + 3] = a2b3 + a3b2; |
| 1147 | } |
| 1148 | |
| 1149 | return p; |
| 1150 | } |
| 1151 | |
| 1152 | // Sample p from a centered binomial distribution with n=2η and p=½ - viz: |
| 1153 | // coefficients are in {-η, …, η} with probabilities |
| 1154 | // |
| 1155 | // {ncr(0, 2η)/2^2η, ncr(1, 2η)/2^2η, …, ncr(2η,2η)/2^2η} |
| 1156 | fn noise(comptime eta: u8, nonce: u8, seed: *const [32]u8) Poly { |
| 1157 | var h = sha3.Shake256.init(.{}); |
| 1158 | const suffix: [1]u8 = .{nonce}; |
| 1159 | h.update(seed); |
| 1160 | h.update(&suffix); |
| 1161 | |
| 1162 | // The distribution at hand is exactly the same as that |
| 1163 | // of (a₁ + a₂ + … + a_η) - (b₁ + … + b_η) where a_i,b_i~U(1). |
| 1164 | // Thus we need 2η bits per coefficient. |
| 1165 | const buf_len = comptime 2 * eta * N / 8; |
| 1166 | var buf: [buf_len]u8 = undefined; |
| 1167 | h.squeeze(&buf); |
| 1168 | |
| 1169 | // buf is interpreted as a₁…a_ηb₁…b_ηa₁…a_ηb₁…b_η…. We process |
| 1170 | // multiple coefficients in one batch. |
| 1171 | |
| 1172 | const T = switch (builtin.target.cpu.arch) { |
| 1173 | .x86_64, .x86 => u32, // Generates better code on Intel CPUs |
| 1174 | else => u64, // u128 might be faster on some other CPUs. |
| 1175 | }; |
| 1176 | |
| 1177 | comptime var batch_count: usize = undefined; |
| 1178 | comptime var batch_bytes: usize = undefined; |
| 1179 | comptime var mask: T = 0; |
| 1180 | comptime { |
| 1181 | batch_count = @bitSizeOf(T) / @as(usize, 2 * eta); |
| 1182 | while (@rem(N, batch_count) != 0 and batch_count > 0) : (batch_count -= 1) {} |
| 1183 | assert(batch_count > 0); |
| 1184 | assert(@rem(2 * eta * batch_count, 8) == 0); |
| 1185 | batch_bytes = 2 * eta * batch_count / 8; |
| 1186 | |
| 1187 | for (0..2 * eta * batch_count) |_| { |
| 1188 | mask <<= eta; |
| 1189 | mask |= 1; |
| 1190 | } |
| 1191 | } |
| 1192 | |
| 1193 | var ret: Poly = undefined; |
| 1194 | for (0..comptime N / batch_count) |i| { |
| 1195 | // Read coefficients into t. In the case of η=3, |
| 1196 | // we have t = a₁ + 2a₂ + 4a₃ + 8b₁ + 16b₂ + … |
| 1197 | var t: T = 0; |
| 1198 | inline for (0..batch_bytes) |j| { |
| 1199 | t |= @as(T, buf[batch_bytes * i + j]) << (8 * j); |
| 1200 | } |
| 1201 | |
| 1202 | // Accumelate `a's and `b's together by masking them out, shifting |
| 1203 | // and adding. For η=3, we have d = a₁ + a₂ + a₃ + 8(b₁ + b₂ + b₃) + … |
| 1204 | var d: T = 0; |
| 1205 | inline for (0..eta) |j| { |
| 1206 | d += (t >> j) & mask; |
| 1207 | } |
| 1208 | |
| 1209 | // Extract each a and b separately and set coefficient in polynomial. |
| 1210 | inline for (0..batch_count) |j| { |
| 1211 | const mask2 = comptime (1 << eta) - 1; |
| 1212 | const a = @intCast(i16, (d >> (comptime (2 * j * eta))) & mask2); |
| 1213 | const b = @intCast(i16, (d >> (comptime ((2 * j + 1) * eta))) & mask2); |
| 1214 | ret.cs[batch_count * i + j] = a - b; |
| 1215 | } |
| 1216 | } |
| 1217 | |
| 1218 | return ret; |
| 1219 | } |
| 1220 | |
| 1221 | // Sample p uniformly from the given seed and x and y coordinates. |
| 1222 | fn uniform(seed: [32]u8, x: u8, y: u8) Poly { |
| 1223 | var h = sha3.Shake128.init(.{}); |
| 1224 | const suffix: [2]u8 = .{ x, y }; |
| 1225 | h.update(&seed); |
| 1226 | h.update(&suffix); |
| 1227 | |
| 1228 | const buf_len = sha3.Shake128.block_length; // rate SHAKE-128 |
| 1229 | var buf: [buf_len]u8 = undefined; |
| 1230 | |
| 1231 | var ret: Poly = undefined; |
| 1232 | var i: usize = 0; // index into ret.cs |
| 1233 | outer: while (true) { |
| 1234 | h.squeeze(&buf); |
| 1235 | |
| 1236 | var j: usize = 0; // index into buf |
| 1237 | while (j < buf_len) : (j += 3) { |
| 1238 | const b0 = @as(u16, buf[j]); |
| 1239 | const b1 = @as(u16, buf[j + 1]); |
| 1240 | const b2 = @as(u16, buf[j + 2]); |
| 1241 | |
| 1242 | const ts: [2]u16 = .{ |
| 1243 | b0 | ((b1 & 0xf) << 8), |
| 1244 | (b1 >> 4) | (b2 << 4), |
| 1245 | }; |
| 1246 | |
| 1247 | inline for (ts) |t| { |
| 1248 | if (t < Q) { |
| 1249 | ret.cs[i] = @intCast(i16, t); |
| 1250 | i += 1; |
| 1251 | |
| 1252 | if (i == N) { |
| 1253 | break :outer; |
| 1254 | } |
| 1255 | } |
| 1256 | } |
| 1257 | } |
| 1258 | } |
| 1259 | |
| 1260 | return ret; |
| 1261 | } |
| 1262 | |
| 1263 | // Packs p. |
| 1264 | // |
| 1265 | // Assumes p is normalized (and not just Barrett reduced). |
| 1266 | fn toBytes(p: Poly) [bytes_length]u8 { |
| 1267 | var ret: [bytes_length]u8 = undefined; |
| 1268 | for (0..comptime N / 2) |i| { |
| 1269 | const t0 = @intCast(u16, p.cs[2 * i]); |
| 1270 | const t1 = @intCast(u16, p.cs[2 * i + 1]); |
| 1271 | ret[3 * i] = @truncate(u8, t0); |
| 1272 | ret[3 * i + 1] = @truncate(u8, (t0 >> 8) | (t1 << 4)); |
| 1273 | ret[3 * i + 2] = @truncate(u8, t1 >> 4); |
| 1274 | } |
| 1275 | return ret; |
| 1276 | } |
| 1277 | |
| 1278 | // Unpacks a Poly from buf. |
| 1279 | // |
| 1280 | // p will not be normalized; instead 0 ≤ p[i] < 4096. |
| 1281 | fn fromBytes(buf: *const [bytes_length]u8) Poly { |
| 1282 | var ret: Poly = undefined; |
| 1283 | for (0..comptime N / 2) |i| { |
| 1284 | const b0 = @as(i16, buf[3 * i]); |
| 1285 | const b1 = @as(i16, buf[3 * i + 1]); |
| 1286 | const b2 = @as(i16, buf[3 * i + 2]); |
| 1287 | ret.cs[2 * i] = b0 | ((b1 & 0xf) << 8); |
| 1288 | ret.cs[2 * i + 1] = (b1 >> 4) | b2 << 4; |
| 1289 | } |
| 1290 | return ret; |
| 1291 | } |
| 1292 | }; |
| 1293 | |
| 1294 | // A vector of K polynomials. |
| 1295 | fn Vec(comptime K: u8) type { |
| 1296 | return struct { |
| 1297 | ps: [K]Poly, |
| 1298 | |
| 1299 | const Self = @This(); |
| 1300 | const bytes_length = K * Poly.bytes_length; |
| 1301 | |
| 1302 | fn compressedSize(comptime d: u8) usize { |
| 1303 | return Poly.compressedSize(d) * K; |
| 1304 | } |
| 1305 | |
| 1306 | fn ntt(a: Self) Self { |
| 1307 | var ret: Self = undefined; |
| 1308 | for (0..K) |i| { |
| 1309 | ret.ps[i] = a.ps[i].ntt(); |
| 1310 | } |
| 1311 | return ret; |
| 1312 | } |
| 1313 | |
| 1314 | fn invNTT(a: Self) Self { |
| 1315 | var ret: Self = undefined; |
| 1316 | for (0..K) |i| { |
| 1317 | ret.ps[i] = a.ps[i].invNTT(); |
| 1318 | } |
| 1319 | return ret; |
| 1320 | } |
| 1321 | |
| 1322 | fn normalize(a: Self) Self { |
| 1323 | var ret: Self = undefined; |
| 1324 | for (0..K) |i| { |
| 1325 | ret.ps[i] = a.ps[i].normalize(); |
| 1326 | } |
| 1327 | return ret; |
| 1328 | } |
| 1329 | |
| 1330 | fn barrettReduce(a: Self) Self { |
| 1331 | var ret: Self = undefined; |
| 1332 | for (0..K) |i| { |
| 1333 | ret.ps[i] = a.ps[i].barrettReduce(); |
| 1334 | } |
| 1335 | return ret; |
| 1336 | } |
| 1337 | |
| 1338 | fn add(a: Self, b: Self) Self { |
| 1339 | var ret: Self = undefined; |
| 1340 | for (0..K) |i| { |
| 1341 | ret.ps[i] = a.ps[i].add(b.ps[i]); |
| 1342 | } |
| 1343 | return ret; |
| 1344 | } |
| 1345 | |
| 1346 | fn sub(a: Self, b: Self) Self { |
| 1347 | var ret: Self = undefined; |
| 1348 | for (0..K) |i| { |
| 1349 | ret.ps[i] = a.ps[i].sub(b.ps[i]); |
| 1350 | } |
| 1351 | return ret; |
| 1352 | } |
| 1353 | |
| 1354 | // Samples v[i] from centered binomial distribution with the given η, |
| 1355 | // seed and nonce+i. |
| 1356 | fn noise(comptime eta: u8, nonce: u8, seed: *const [32]u8) Self { |
| 1357 | var ret: Self = undefined; |
| 1358 | for (0..K) |i| { |
| 1359 | ret.ps[i] = Poly.noise(eta, nonce + @intCast(u8, i), seed); |
| 1360 | } |
| 1361 | return ret; |
| 1362 | } |
| 1363 | |
| 1364 | // Sets p to the inner product of a and b using "pointwise" multiplication. |
| 1365 | // |
| 1366 | // See MulHat() and NTT() for a description of the multiplication. |
| 1367 | // Assumes a and b are in Montgomery form. p will be in Montgomery form, |
| 1368 | // and its coefficients will be bounded in absolute value by 2kq. |
| 1369 | // If a and b are not in Montgomery form, then the action is the same |
| 1370 | // as "pointwise" multiplication followed by multiplying by R⁻¹, the inverse |
| 1371 | // of the Montgomery factor. |
| 1372 | fn dotHat(a: Self, b: Self) Poly { |
| 1373 | var ret: Poly = Poly.zero; |
| 1374 | for (0..K) |i| { |
| 1375 | ret = ret.add(a.ps[i].mulHat(b.ps[i])); |
| 1376 | } |
| 1377 | return ret; |
| 1378 | } |
| 1379 | |
| 1380 | fn compress(v: Self, comptime d: u8) [compressedSize(d)]u8 { |
| 1381 | const cs = comptime Poly.compressedSize(d); |
| 1382 | var ret: [compressedSize(d)]u8 = undefined; |
| 1383 | inline for (0..K) |i| { |
| 1384 | mem.copy(u8, ret[i * cs .. (i + 1) * cs], &v.ps[i].compress(d)); |
| 1385 | } |
| 1386 | return ret; |
| 1387 | } |
| 1388 | |
| 1389 | fn decompress(comptime d: u8, buf: *const [compressedSize(d)]u8) Self { |
| 1390 | const cs = comptime Poly.compressedSize(d); |
| 1391 | var ret: Self = undefined; |
| 1392 | inline for (0..K) |i| { |
| 1393 | ret.ps[i] = Poly.decompress(d, buf[i * cs .. (i + 1) * cs]); |
| 1394 | } |
| 1395 | return ret; |
| 1396 | } |
| 1397 | |
| 1398 | /// Serializes the key into a byte array. |
| 1399 | fn toBytes(v: Self) [bytes_length]u8 { |
| 1400 | var ret: [bytes_length]u8 = undefined; |
| 1401 | inline for (0..K) |i| { |
| 1402 | mem.copy( |
| 1403 | u8, |
| 1404 | ret[i * Poly.bytes_length .. (i + 1) * Poly.bytes_length], |
| 1405 | &v.ps[i].toBytes(), |
| 1406 | ); |
| 1407 | } |
| 1408 | return ret; |
| 1409 | } |
| 1410 | |
| 1411 | /// Deserializes the key from a byte array. |
| 1412 | fn fromBytes(buf: *const [bytes_length]u8) Self { |
| 1413 | var ret: Self = undefined; |
| 1414 | inline for (0..K) |i| { |
| 1415 | ret.ps[i] = Poly.fromBytes( |
| 1416 | buf[i * Poly.bytes_length .. (i + 1) * Poly.bytes_length], |
| 1417 | ); |
| 1418 | } |
| 1419 | return ret; |
| 1420 | } |
| 1421 | }; |
| 1422 | } |
| 1423 | |
| 1424 | // A matrix of K vectors |
| 1425 | fn Mat(comptime K: u8) type { |
| 1426 | return struct { |
| 1427 | const Self = @This(); |
| 1428 | vs: [K]Vec(K), |
| 1429 | |
| 1430 | fn uniform(seed: [32]u8, comptime transposed: bool) Self { |
| 1431 | var ret: Self = undefined; |
| 1432 | var i: u8 = 0; |
| 1433 | while (i < K) : (i += 1) { |
| 1434 | var j: u8 = 0; |
| 1435 | while (j < K) : (j += 1) { |
| 1436 | ret.vs[i].ps[j] = Poly.uniform( |
| 1437 | seed, |
| 1438 | if (transposed) i else j, |
| 1439 | if (transposed) j else i, |
| 1440 | ); |
| 1441 | } |
| 1442 | } |
| 1443 | return ret; |
| 1444 | } |
| 1445 | |
| 1446 | // Returns transpose of A |
| 1447 | fn transpose(m: Self) Self { |
| 1448 | var ret: Self = undefined; |
| 1449 | for (0..K) |i| { |
| 1450 | for (0..K) |j| { |
| 1451 | ret.vs[i].ps[j] = m.vs[j].ps[i]; |
| 1452 | } |
| 1453 | } |
| 1454 | return ret; |
| 1455 | } |
| 1456 | }; |
| 1457 | } |
| 1458 | |
| 1459 | // Returns `true` if a ≠ b. |
| 1460 | fn ctneq(comptime len: usize, a: [len]u8, b: [len]u8) u1 { |
| 1461 | return 1 - @boolToInt(crypto.utils.timingSafeEql([len]u8, a, b)); |
| 1462 | } |
| 1463 | |
| 1464 | // Copy src into dst given b = 1. |
| 1465 | fn cmov(comptime len: usize, dst: *[len]u8, src: [len]u8, b: u1) void { |
| 1466 | const mask = @as(u8, 0) -% b; |
| 1467 | for (0..len) |i| { |
| 1468 | dst[i] ^= mask & (dst[i] ^ src[i]); |
| 1469 | } |
| 1470 | } |
| 1471 | |
| 1472 | test "MulHat" { |
| 1473 | var rnd = RndGen.init(0); |
| 1474 | |
| 1475 | for (0..100) |_| { |
| 1476 | const a = Poly.randAbsLeqQ(&rnd); |
| 1477 | const b = Poly.randAbsLeqQ(&rnd); |
| 1478 | |
| 1479 | const p2 = a.ntt().mulHat(b.ntt()).barrettReduce().invNTT().normalize(); |
| 1480 | var p: Poly = undefined; |
| 1481 | |
| 1482 | mem.set(i16, &p.cs, 0); |
| 1483 | |
| 1484 | for (0..N) |i| { |
| 1485 | for (0..N) |j| { |
| 1486 | var v = montReduce(@as(i32, a.cs[i]) * @as(i32, b.cs[j])); |
| 1487 | var k = i + j; |
| 1488 | if (k >= N) { |
| 1489 | // Recall Xᴺ = -1. |
| 1490 | k -= N; |
| 1491 | v = -v; |
| 1492 | } |
| 1493 | p.cs[k] = feBarrettReduce(v + p.cs[k]); |
| 1494 | } |
| 1495 | } |
| 1496 | |
| 1497 | p = p.toMont().normalize(); |
| 1498 | |
| 1499 | try testing.expectEqual(p, p2); |
| 1500 | } |
| 1501 | } |
| 1502 | |
| 1503 | test "NTT" { |
| 1504 | var rnd = RndGen.init(0); |
| 1505 | |
| 1506 | for (0..1000) |_| { |
| 1507 | var p = Poly.randAbsLeqQ(&rnd); |
| 1508 | const q = p.toMont().normalize(); |
| 1509 | p = p.ntt(); |
| 1510 | |
| 1511 | for (0..N) |i| { |
| 1512 | try testing.expect(p.cs[i] <= 7 * Q and -7 * Q <= p.cs[i]); |
| 1513 | } |
| 1514 | |
| 1515 | p = p.normalize().invNTT(); |
| 1516 | for (0..N) |i| { |
| 1517 | try testing.expect(p.cs[i] <= Q and -Q <= p.cs[i]); |
| 1518 | } |
| 1519 | |
| 1520 | p = p.normalize(); |
| 1521 | |
| 1522 | try testing.expectEqual(p, q); |
| 1523 | } |
| 1524 | } |
| 1525 | |
| 1526 | test "Compression" { |
| 1527 | var rnd = RndGen.init(0); |
| 1528 | inline for (.{ 1, 4, 5, 10, 11 }) |d| { |
| 1529 | for (0..1000) |_| { |
| 1530 | const p = Poly.randNormalized(&rnd); |
| 1531 | const pp = p.compress(d); |
| 1532 | const pq = Poly.decompress(d, &pp).compress(d); |
| 1533 | try testing.expectEqual(pp, pq); |
| 1534 | } |
| 1535 | } |
| 1536 | } |
| 1537 | |
| 1538 | test "noise" { |
| 1539 | var seed: [32]u8 = undefined; |
| 1540 | for (&seed, 0..) |*s, i| { |
| 1541 | s.* = @intCast(u8, i); |
| 1542 | } |
| 1543 | try testing.expectEqual(Poly.noise(3, 37, &seed).cs, .{ |
| 1544 | 0, 0, 1, -1, 0, 2, 0, -1, -1, 3, 0, 1, -2, -2, 0, 1, -2, |
| 1545 | 1, 0, -2, 3, 0, 0, 0, 1, 3, 1, 1, 2, 1, -1, -1, -1, 0, |
| 1546 | 1, 0, 1, 0, 2, 0, 1, -2, 0, -1, -1, -2, 1, -1, -1, 2, -1, |
| 1547 | 1, 1, 2, -3, -1, -1, 0, 0, 0, 0, 1, -1, -2, -2, 0, -2, 0, |
| 1548 | 0, 0, 1, 0, -1, -1, 1, -2, 2, 0, 0, 2, -2, 0, 1, 0, 1, |
| 1549 | 1, 1, 0, 1, -2, -1, -2, -1, 1, 0, 0, 0, 0, 0, 1, 0, -1, |
| 1550 | -1, 0, -1, 1, 0, 1, 0, -1, -1, 0, -2, 2, 0, -2, 1, -1, 0, |
| 1551 | 1, -1, -1, 2, 1, 0, 0, -2, -1, 2, 0, 0, 0, -1, -1, 3, 1, |
| 1552 | 0, 1, 0, 1, 0, 2, 1, 0, 0, 1, 0, 1, 0, 0, -1, -1, -1, |
| 1553 | 0, 1, 3, 1, 0, 1, 0, 1, -1, -1, -1, -1, 0, 0, -2, -1, -1, |
| 1554 | 2, 0, 1, 0, 1, 0, 2, -2, 0, 1, 1, -3, -1, -2, -1, 0, 1, |
| 1555 | 0, 1, -2, 2, 2, 1, 1, 0, -1, 0, -1, -1, 1, 0, -1, 2, 1, |
| 1556 | -1, 1, 2, -2, 1, 2, 0, 1, 2, 1, 0, 0, 2, 1, 2, 1, 0, |
| 1557 | 2, 1, 0, 0, -1, -1, 1, -1, 0, 1, -1, 2, 2, 0, 0, -1, 1, |
| 1558 | 1, 1, 1, 0, 0, -2, 0, -1, 1, 2, 0, 0, 1, 1, -1, 1, 0, |
| 1559 | 1, |
| 1560 | }); |
| 1561 | try testing.expectEqual(Poly.noise(2, 37, &seed).cs, .{ |
| 1562 | 1, 0, 1, -1, -1, -2, -1, -1, 2, 0, -1, 0, 0, -1, |
| 1563 | 1, 1, -1, 1, 0, 2, -2, 0, 1, 2, 0, 0, -1, 1, |
| 1564 | 0, -1, 1, -1, 1, 2, 1, 1, 0, -1, 1, -1, -2, -1, |
| 1565 | 1, -1, -1, -1, 2, -1, -1, 0, 0, 1, 1, -1, 1, 1, |
| 1566 | 1, 1, -1, -2, 0, 1, 0, 0, 2, 1, -1, 2, 0, 0, |
| 1567 | 1, 1, 0, -1, 0, 0, -1, -1, 2, 0, 1, -1, 2, -1, |
| 1568 | -1, -1, -1, 0, -2, 0, 2, 1, 0, 0, 0, -1, 0, 0, |
| 1569 | 0, -1, -1, 0, -1, -1, 0, -1, 0, 0, -2, 1, 1, 0, |
| 1570 | 1, 0, 1, 0, 1, 1, -1, 2, 0, 1, -1, 1, 2, 0, |
| 1571 | 0, 0, 0, -1, -1, -1, 0, 1, 0, -1, 2, 0, 0, 1, |
| 1572 | 1, 1, 0, 1, -1, 1, 2, 1, 0, 2, -1, 1, -1, -2, |
| 1573 | -1, -2, -1, 1, 0, -2, -2, -1, 1, 0, 0, 0, 0, 1, |
| 1574 | 0, 0, 0, 2, 2, 0, 1, 0, -1, -1, 0, 2, 0, 0, |
| 1575 | -2, 1, 0, 2, 1, -1, -2, 0, 0, -1, 1, 1, 0, 0, |
| 1576 | 2, 0, 1, 1, -2, 1, -2, 1, 1, 0, 2, 0, -1, 0, |
| 1577 | -1, 0, 1, 2, 0, 1, 0, -2, 1, -2, -2, 1, -1, 0, |
| 1578 | -1, 1, 1, 0, 0, 0, 1, 0, -1, 1, 1, 0, 0, 0, |
| 1579 | 0, 1, 0, 1, -1, 0, 1, -1, -1, 2, 0, 0, 1, -1, |
| 1580 | 0, 1, -1, 0, |
| 1581 | }); |
| 1582 | } |
| 1583 | |
| 1584 | test "uniform sampling" { |
| 1585 | var seed: [32]u8 = undefined; |
| 1586 | for (&seed, 0..) |*s, i| { |
| 1587 | s.* = @intCast(u8, i); |
| 1588 | } |
| 1589 | try testing.expectEqual(Poly.uniform(seed, 1, 0).cs, .{ |
| 1590 | 797, 993, 161, 6, 2608, 2385, 2096, 2661, 1676, 247, 2440, |
| 1591 | 342, 634, 194, 1570, 2848, 986, 684, 3148, 3208, 2018, 351, |
| 1592 | 2288, 612, 1394, 170, 1521, 3119, 58, 596, 2093, 1549, 409, |
| 1593 | 2156, 1934, 1730, 1324, 388, 446, 418, 1719, 2202, 1812, 98, |
| 1594 | 1019, 2369, 214, 2699, 28, 1523, 2824, 273, 402, 2899, 246, |
| 1595 | 210, 1288, 863, 2708, 177, 3076, 349, 44, 949, 854, 1371, |
| 1596 | 957, 292, 2502, 1617, 1501, 254, 7, 1761, 2581, 2206, 2655, |
| 1597 | 1211, 629, 1274, 2358, 816, 2766, 2115, 2985, 1006, 2433, 856, |
| 1598 | 2596, 3192, 1, 1378, 2345, 707, 1891, 1669, 536, 1221, 710, |
| 1599 | 2511, 120, 1176, 322, 1897, 2309, 595, 2950, 1171, 801, 1848, |
| 1600 | 695, 2912, 1396, 1931, 1775, 2904, 893, 2507, 1810, 2873, 253, |
| 1601 | 1529, 1047, 2615, 1687, 831, 1414, 965, 3169, 1887, 753, 3246, |
| 1602 | 1937, 115, 2953, 586, 545, 1621, 1667, 3187, 1654, 1988, 1857, |
| 1603 | 512, 1239, 1219, 898, 3106, 391, 1331, 2228, 3169, 586, 2412, |
| 1604 | 845, 768, 156, 662, 478, 1693, 2632, 573, 2434, 1671, 173, |
| 1605 | 969, 364, 1663, 2701, 2169, 813, 1000, 1471, 720, 2431, 2530, |
| 1606 | 3161, 733, 1691, 527, 2634, 335, 26, 2377, 1707, 767, 3020, |
| 1607 | 950, 502, 426, 1138, 3208, 2607, 2389, 44, 1358, 1392, 2334, |
| 1608 | 875, 2097, 173, 1697, 2578, 942, 1817, 974, 1165, 2853, 1958, |
| 1609 | 2973, 3282, 271, 1236, 1677, 2230, 673, 1554, 96, 242, 1729, |
| 1610 | 2518, 1884, 2272, 71, 1382, 924, 1807, 1610, 456, 1148, 2479, |
| 1611 | 2152, 238, 2208, 2329, 713, 1175, 1196, 757, 1078, 3190, 3169, |
| 1612 | 708, 3117, 154, 1751, 3225, 1364, 154, 23, 2842, 1105, 1419, |
| 1613 | 79, 5, 2013, |
| 1614 | }); |
| 1615 | } |
| 1616 | |
| 1617 | test "Polynomial packing" { |
| 1618 | var rnd = RndGen.init(0); |
| 1619 | |
| 1620 | for (0..1000) |_| { |
| 1621 | const p = Poly.randNormalized(&rnd); |
| 1622 | try testing.expectEqual(Poly.fromBytes(&p.toBytes()), p); |
| 1623 | } |
| 1624 | } |
| 1625 | |
| 1626 | test "Test inner PKE" { |
| 1627 | var seed: [32]u8 = undefined; |
| 1628 | var pt: [32]u8 = undefined; |
| 1629 | for (&seed, &pt, 0..) |*s, *p, i| { |
| 1630 | s.* = @intCast(u8, i); |
| 1631 | p.* = @intCast(u8, i + 32); |
| 1632 | } |
| 1633 | inline for (modes) |mode| { |
| 1634 | for (0..100) |i| { |
| 1635 | var pk: mode.InnerPk = undefined; |
| 1636 | var sk: mode.InnerSk = undefined; |
| 1637 | seed[0] = @intCast(u8, i); |
| 1638 | mode.innerKeyFromSeed(seed, &pk, &sk); |
| 1639 | for (0..10) |j| { |
| 1640 | seed[1] = @intCast(u8, j); |
| 1641 | try testing.expectEqual(sk.decrypt(&pk.encrypt(&pt, &seed)), pt); |
| 1642 | } |
| 1643 | } |
| 1644 | } |
| 1645 | } |
| 1646 | |
| 1647 | test "Test happy flow" { |
| 1648 | var seed: [64]u8 = undefined; |
| 1649 | for (&seed, 0..) |*s, i| { |
| 1650 | s.* = @intCast(u8, i); |
| 1651 | } |
| 1652 | inline for (modes) |mode| { |
| 1653 | for (0..100) |i| { |
| 1654 | seed[0] = @intCast(u8, i); |
| 1655 | const kp = try mode.KeyPair.create(seed); |
| 1656 | const sk = try mode.SecretKey.fromBytes(&kp.secret_key.toBytes()); |
| 1657 | try testing.expectEqual(sk, kp.secret_key); |
| 1658 | const pk = try mode.PublicKey.fromBytes(&kp.public_key.toBytes()); |
| 1659 | try testing.expectEqual(pk, kp.public_key); |
| 1660 | for (0..10) |j| { |
| 1661 | seed[1] = @intCast(u8, j); |
| 1662 | const e = pk.encaps(seed[0..32].*); |
| 1663 | try testing.expectEqual(e.shared_secret, try sk.decaps(&e.ciphertext)); |
| 1664 | } |
| 1665 | } |
| 1666 | } |
| 1667 | } |
| 1668 | |
| 1669 | // Code to test NIST Known Answer Tests (KAT), see PQCgenKAT.c. |
| 1670 | |
| 1671 | const sha2 = crypto.hash.sha2; |
| 1672 | |
| 1673 | test "NIST KAT test" { |
| 1674 | inline for (.{ |
| 1675 | .{ Kyber512, "e9c2bd37133fcb40772f81559f14b1f58dccd1c816701be9ba6214d43baf4547" }, |
| 1676 | .{ Kyber1024, "89248f2f33f7f4f7051729111f3049c409a933ec904aedadf035f30fa5646cd5" }, |
| 1677 | .{ Kyber768, "a1e122cad3c24bc51622e4c242d8b8acbcd3f618fee4220400605ca8f9ea02c2" }, |
| 1678 | }) |modeHash| { |
| 1679 | const mode = modeHash[0]; |
| 1680 | var seed: [48]u8 = undefined; |
| 1681 | for (&seed, 0..) |*s, i| { |
| 1682 | s.* = @intCast(u8, i); |
| 1683 | } |
| 1684 | var f = sha2.Sha256.init(.{}); |
| 1685 | const fw = f.writer(); |
| 1686 | var g = NistDRBG.init(seed); |
| 1687 | try std.fmt.format(fw, "# {s}\n\n", .{mode.name}); |
| 1688 | for (0..100) |i| { |
| 1689 | g.fill(&seed); |
| 1690 | try std.fmt.format(fw, "count = {}\n", .{i}); |
| 1691 | try std.fmt.format(fw, "seed = {s}\n", .{std.fmt.fmtSliceHexUpper(&seed)}); |
| 1692 | var g2 = NistDRBG.init(seed); |
| 1693 | |
| 1694 | // This is not equivalent to g2.fill(kseed[:]). As the reference |
| 1695 | // implementation calls randombytes twice generating the keypair, |
| 1696 | // we have to do that as well. |
| 1697 | var kseed: [64]u8 = undefined; |
| 1698 | var eseed: [32]u8 = undefined; |
| 1699 | g2.fill(kseed[0..32]); |
| 1700 | g2.fill(kseed[32..64]); |
| 1701 | g2.fill(&eseed); |
| 1702 | const kp = try mode.KeyPair.create(kseed); |
| 1703 | const e = kp.public_key.encaps(eseed); |
| 1704 | const ss2 = try kp.secret_key.decaps(&e.ciphertext); |
| 1705 | try testing.expectEqual(ss2, e.shared_secret); |
| 1706 | try std.fmt.format(fw, "pk = {s}\n", .{std.fmt.fmtSliceHexUpper(&kp.public_key.toBytes())}); |
| 1707 | try std.fmt.format(fw, "sk = {s}\n", .{std.fmt.fmtSliceHexUpper(&kp.secret_key.toBytes())}); |
| 1708 | try std.fmt.format(fw, "ct = {s}\n", .{std.fmt.fmtSliceHexUpper(&e.ciphertext)}); |
| 1709 | try std.fmt.format(fw, "ss = {s}\n\n", .{std.fmt.fmtSliceHexUpper(&e.shared_secret)}); |
| 1710 | } |
| 1711 | |
| 1712 | var out: [32]u8 = undefined; |
| 1713 | f.final(&out); |
| 1714 | var outHex: [64]u8 = undefined; |
| 1715 | _ = try std.fmt.bufPrint(&outHex, "{s}", .{std.fmt.fmtSliceHexLower(&out)}); |
| 1716 | try testing.expectEqual(outHex, modeHash[1].*); |
| 1717 | } |
| 1718 | } |
| 1719 | |
| 1720 | const NistDRBG = struct { |
| 1721 | key: [32]u8, |
| 1722 | v: [16]u8, |
| 1723 | |
| 1724 | fn incV(g: *NistDRBG) void { |
| 1725 | var j: usize = 15; |
| 1726 | while (j >= 0) : (j -= 1) { |
| 1727 | if (g.v[j] == 255) { |
| 1728 | g.v[j] = 0; |
| 1729 | } else { |
| 1730 | g.v[j] += 1; |
| 1731 | break; |
| 1732 | } |
| 1733 | } |
| 1734 | } |
| 1735 | |
| 1736 | // AES256_CTR_DRBG_Update(pd, &g.key, &g.v). |
| 1737 | fn update(g: *NistDRBG, pd: ?[48]u8) void { |
| 1738 | var buf: [48]u8 = undefined; |
| 1739 | const ctx = crypto.core.aes.Aes256.initEnc(g.key); |
| 1740 | var i: usize = 0; |
| 1741 | while (i < 3) : (i += 1) { |
| 1742 | g.incV(); |
| 1743 | var block: [16]u8 = undefined; |
| 1744 | ctx.encrypt(&block, &g.v); |
| 1745 | mem.copy(u8, buf[i * 16 .. (i + 1) * 16], &block); |
| 1746 | } |
| 1747 | if (pd) |p| { |
| 1748 | for (&buf, p) |*b, x| { |
| 1749 | b.* ^= x; |
| 1750 | } |
| 1751 | } |
| 1752 | mem.copy(u8, &g.key, buf[0..32]); |
| 1753 | mem.copy(u8, &g.v, buf[32..48]); |
| 1754 | } |
| 1755 | |
| 1756 | // randombytes. |
| 1757 | fn fill(g: *NistDRBG, out: []u8) void { |
| 1758 | var block: [16]u8 = undefined; |
| 1759 | var dst = out; |
| 1760 | |
| 1761 | const ctx = crypto.core.aes.Aes256.initEnc(g.key); |
| 1762 | while (dst.len > 0) { |
| 1763 | g.incV(); |
| 1764 | ctx.encrypt(&block, &g.v); |
| 1765 | if (dst.len < 16) { |
| 1766 | mem.copy(u8, dst, block[0..dst.len]); |
| 1767 | break; |
| 1768 | } |
| 1769 | mem.copy(u8, dst, &block); |
| 1770 | dst = dst[16..dst.len]; |
| 1771 | } |
| 1772 | g.update(null); |
| 1773 | } |
| 1774 | |
| 1775 | fn init(seed: [48]u8) NistDRBG { |
| 1776 | var ret: NistDRBG = .{ .key = .{0} ** 32, .v = .{0} ** 16 }; |
| 1777 | ret.update(seed); |
| 1778 | return ret; |
| 1779 | } |
| 1780 | }; |
Bas Westerbaan <bas@westerbaan.name> 2023-03-15 04:50:45+01:00
GitHub <noreply@github.com> 2023-03-15 03:50:45+00:00
Signed by PGP key 4AEE18F83AFDEB23
lib