/*

 * Copyright 2024-2025 The OpenSSL Project Authors. All Rights Reserved.

 *

 * Licensed under the Apache License 2.0 (the "License").  You may not use

 * this file except in compliance with the License.  You can obtain a copy

 * in the file LICENSE in the source distribution or at

 * https://www.openssl.org/source/license.html

 */

#include <openssl/byteorder.h>

#include <openssl/rand.h>

#include <openssl/proverr.h>

#include "crypto/ml_kem.h"

#include "internal/common.h"

#include "internal/constant_time.h"

#include "internal/sha3.h"

#if defined(OPENSSL_CONSTANT_TIME_VALIDATION)

#include <valgrind/memcheck.h>

#endif

#if ML_KEM_SEED_BYTES != ML_KEM_SHARED_SECRET_BYTES + ML_KEM_RANDOM_BYTES

#error "ML-KEM keygen seed length != shared secret + random bytes length"

#endif

#if ML_KEM_SHARED_SECRET_BYTES != ML_KEM_RANDOM_BYTES

#error "Invalid unequal lengths of ML-KEM shared secret and random inputs"

#endif

#if UINT_MAX < UINT32_MAX

#error "Unsupported compiler: sizeof(unsigned int) < sizeof(uint32_t)"

#endif

/* Handy function-like bit-extraction macros */

#define bit0(b) ((b) & 1)

#define bitn(n, b) (((b) >> n) & 1)

/*

 * 12 bits are sufficient to losslessly represent values in [0, q-1].

 * INVERSE_DEGREE is (n/2)^-1 mod q; used in inverse NTT.

 */

#define DEGREE ML_KEM_DEGREE

#define INVERSE_DEGREE (ML_KEM_PRIME - 2 * 13)

#define LOG2PRIME 12

#define BARRETT_SHIFT (2 * LOG2PRIME)

#ifdef SHA3_BLOCKSIZE

#define SHAKE128_BLOCKSIZE SHA3_BLOCKSIZE(128)

#endif

/*

 * Return whether a value that can only be 0 or 1 is non-zero, in constant time

 * in practice!  The return value is a mask that is all ones if true, and all

 * zeros otherwise (twos-complement arithmetic assumed for unsigned values).

 *

 * Although this is used in constant-time selects, we omit a value barrier

 * here.  Value barriers impede auto-vectorization (likely because it forces

 * the value to transit through a general-purpose register). On AArch64, this

 * is a difference of 2x.

 *

 * We usually add value barriers to selects because Clang turns consecutive

 * selects with the same condition into a branch instead of CMOV/CSEL. This

 * condition does not occur in Kyber, so omitting it seems to be safe so far,

 * but see |cbd_2|, |cbd_3|, where reduction needs to be specialised to the

 * sign of the input, rather than adding |q| in advance, and using the generic

 * |reduce_once|.  (David Benjamin, Chromium)

 */

#if 0

#define constish_time_non_zero(b) (~constant_time_is_zero(b));

#else

#define constish_time_non_zero(b) (0u - (b))

#endif

/*

 * The scalar rejection-sampling buffer size needs to be a multiple of 12, but

 * is otherwise arbitrary, the preferred block size matches the internal buffer

 * size of SHAKE128, avoiding internal buffering and copying in SHAKE128. That

 * block size of (1600 - 256)/8 bytes, or 168, just happens to divide by 12!

 *

 * If the blocksize is unknown, or is not divisible by 12, 168 is used as a

 * fallback.

 */

#if defined(SHAKE128_BLOCKSIZE) && (SHAKE128_BLOCKSIZE) % 12 == 0

#define SCALAR_SAMPLING_BUFSIZE (SHAKE128_BLOCKSIZE)

#else

#define SCALAR_SAMPLING_BUFSIZE 168

#endif

/*

 * Structure of keys

 */

typedef struct ossl_ml_kem_scalar_st {

    /* On every function entry and exit, 0 <= c[i] < ML_KEM_PRIME. */

    uint16_t c[ML_KEM_DEGREE];

} scalar;

/* Key material allocation layout */

#define DECLARE_ML_KEM_KEYDATA(name, rank, private_sz)               \

    struct name##_alloc {                                            \

        /* Public vector |t| */                                      \

        scalar tbuf[(rank)];                                         \

        /* Pre-computed matrix |m| (FIPS 203 |A| transpose) */       \

        scalar mbuf[(rank) * (rank)] /* optional private key data */ \

            private_sz                                               \

    }

/* Declare variant-specific public and private storage */

#define DECLARE_ML_KEM_VARIANT_KEYDATA(bits)                        \

    DECLARE_ML_KEM_KEYDATA(pubkey_##bits, ML_KEM_##bits##_RANK, ;); \

    DECLARE_ML_KEM_KEYDATA(prvkey_##bits, ML_KEM_##bits##_RANK, ; scalar sbuf[ML_KEM_##bits##_RANK]; uint8_t zbuf[2 * ML_KEM_RANDOM_BYTES];)

DECLARE_ML_KEM_VARIANT_KEYDATA(512);

DECLARE_ML_KEM_VARIANT_KEYDATA(768);

DECLARE_ML_KEM_VARIANT_KEYDATA(1024);

#undef DECLARE_ML_KEM_VARIANT_KEYDATA

#undef DECLARE_ML_KEM_KEYDATA

typedef __owur int (*CBD_FUNC)(scalar *out, uint8_t in[ML_KEM_RANDOM_BYTES + 1],

    EVP_MD_CTX *mdctx, const ML_KEM_KEY *key);

static void scalar_encode(uint8_t *out, const scalar *s, int bits);

/*

 * The wire-form of a losslessly encoded vector uses 12-bits per element.

 *

 * The wire-form public key consists of the lossless encoding of the public

 * vector |t|, followed by the public seed |rho|.

 *

 * Our serialised private key concatenates serialisations of the private vector

 * |s|, the public key, the public key hash, and the failure secret |z|.

 */

#define VECTOR_BYTES(b) ((3 * DEGREE / 2) * ML_KEM_##b##_RANK)

#define PUBKEY_BYTES(b) (VECTOR_BYTES(b) + ML_KEM_RANDOM_BYTES)

#define PRVKEY_BYTES(b) (2 * PUBKEY_BYTES(b) + ML_KEM_PKHASH_BYTES)

/*

 * Encapsulation produces a vector "u" and a scalar "v", whose coordinates

 * (numbers modulo the ML-KEM prime "q") are lossily encoded using as "du" and

 * "dv" bits, respectively.  This encoding is the ciphertext input for

 * decapsulation.

 */

#define U_VECTOR_BYTES(b) ((DEGREE / 8) * ML_KEM_##b##_DU * ML_KEM_##b##_RANK)

#define V_SCALAR_BYTES(b) ((DEGREE / 8) * ML_KEM_##b##_DV)

#define CTEXT_BYTES(b) (U_VECTOR_BYTES(b) + V_SCALAR_BYTES(b))

#if defined(OPENSSL_CONSTANT_TIME_VALIDATION)

/*

 * CONSTTIME_SECRET takes a pointer and a number of bytes and marks that region

 * of memory as secret. Secret data is tracked as it flows to registers and

 * other parts of a memory. If secret data is used as a condition for a branch,

 * or as a memory index, it will trigger warnings in valgrind.

 */

#define CONSTTIME_SECRET(ptr, len) VALGRIND_MAKE_MEM_UNDEFINED(ptr, len)

/*

 * CONSTTIME_DECLASSIFY takes a pointer and a number of bytes and marks that

 * region of memory as public. Public data is not subject to constant-time

 * rules.

 */

#define CONSTTIME_DECLASSIFY(ptr, len) VALGRIND_MAKE_MEM_DEFINED(ptr, len)

#else

#define CONSTTIME_SECRET(ptr, len)

#define CONSTTIME_DECLASSIFY(ptr, len)

#endif

/*

 * Indices of slots in the vinfo tables below

 */

#define ML_KEM_512_VINFO 0

#define ML_KEM_768_VINFO 1

#define ML_KEM_1024_VINFO 2

/*

 * Per-variant fixed parameters

 */

static const ML_KEM_VINFO vinfo_map[3] = {

    { "ML-KEM-512",

        PRVKEY_BYTES(512),

        sizeof(struct prvkey_512_alloc),

        PUBKEY_BYTES(512),

        sizeof(struct pubkey_512_alloc),

        CTEXT_BYTES(512),

        VECTOR_BYTES(512),

        U_VECTOR_BYTES(512),

        EVP_PKEY_ML_KEM_512,

        ML_KEM_512_BITS,

        ML_KEM_512_RANK,

        ML_KEM_512_DU,

        ML_KEM_512_DV,

        ML_KEM_512_SECBITS },

    { "ML-KEM-768",

        PRVKEY_BYTES(768),

        sizeof(struct prvkey_768_alloc),

        PUBKEY_BYTES(768),

        sizeof(struct pubkey_768_alloc),

        CTEXT_BYTES(768),

        VECTOR_BYTES(768),

        U_VECTOR_BYTES(768),

        EVP_PKEY_ML_KEM_768,

        ML_KEM_768_BITS,

        ML_KEM_768_RANK,

        ML_KEM_768_DU,

        ML_KEM_768_DV,

        ML_KEM_768_SECBITS },

    { "ML-KEM-1024",

        PRVKEY_BYTES(1024),

        sizeof(struct prvkey_1024_alloc),

        PUBKEY_BYTES(1024),

        sizeof(struct pubkey_1024_alloc),

        CTEXT_BYTES(1024),

        VECTOR_BYTES(1024),

        U_VECTOR_BYTES(1024),

        EVP_PKEY_ML_KEM_1024,

        ML_KEM_1024_BITS,

        ML_KEM_1024_RANK,

        ML_KEM_1024_DU,

        ML_KEM_1024_DV,

        ML_KEM_1024_SECBITS }

};

/*

 * Remainders modulo `kPrime`, for sufficiently small inputs, are computed in

 * constant time via Barrett reduction, and a final call to reduce_once(),

 * which reduces inputs that are at most 2*kPrime and is also constant-time.

 */

static const int kPrime = ML_KEM_PRIME;

static const unsigned int kBarrettShift = BARRETT_SHIFT;

static const size_t kBarrettMultiplier = (1 << BARRETT_SHIFT) / ML_KEM_PRIME;

static const uint16_t kHalfPrime = (ML_KEM_PRIME - 1) / 2;

static const uint16_t kInverseDegree = INVERSE_DEGREE;

/*

 * Python helper:

 *

 * p = 3329

 * def bitreverse(i):

 *     ret = 0

 *     for n in range(7):

 *         bit = i & 1

 *         ret <<= 1

 *         ret |= bit

 *         i >>= 1

 *     return ret

 */

/*-

 * First precomputed array from Appendix A of FIPS 203, or else Python:

 * kNTTRoots = [pow(17, bitreverse(i), p) for i in range(128)]

 */

static const uint16_t kNTTRoots[128] = {

    1,

    1729,

    2580,

    3289,

    2642,

    630,

    1897,

    848,

    1062,

    1919,

    193,

    797,

    2786,

    3260,

    569,

    1746,

    296,

    2447,

    1339,

    1476,

    3046,

    56,

    2240,

    1333,

    1426,

    2094,

    535,

    2882,

    2393,

    2879,

    1974,

    821,

    289,

    331,

    3253,

    1756,

    1197,

    2304,

    2277,

    2055,

    650,

    1977,

    2513,

    632,

    2865,

    33,

    1320,

    1915,

    2319,

    1435,

    807,

    452,

    1438,

    2868,

    1534,

    2402,

    2647,

    2617,

    1481,

    648,

    2474,

    3110,

    1227,

    910,

    17,

    2761,

    583,

    2649,

    1637,

    723,

    2288,

    1100,

    1409,

    2662,

    3281,

    233,

    756,

    2156,

    3015,

    3050,

    1703,

    1651,

    2789,

    1789,

    1847,

    952,

    1461,

    2687,

    939,

    2308,

    2437,

    2388,

    733,

    2337,

    268,

    641,

    1584,

    2298,

    2037,

    3220,

    375,

    2549,

    2090,

    1645,

    1063,

    319,

    2773,

    757,

    2099,

    561,

    2466,

    2594,

    2804,

    1092,

    403,

    1026,

    1143,

    2150,

    2775,

    886,

    1722,

    1212,

    1874,

    1029,

    2110,

    2935,

    885,

    2154,

};

/*

 * InverseNTTRoots = [pow(17, -bitreverse(i), p) for i in range(128)]

 * Listed in order of use in the inverse NTT loop (index 0 is skipped):

 *

 *  0, 64, 65, ..., 127, 32, 33, ..., 63, 16, 17, ..., 31, 8, 9, ...

 */

static const uint16_t kInverseNTTRoots[128] = {

    1,

    1175,

    2444,

    394,

    1219,

    2300,

    1455,

    2117,

    1607,

    2443,

    554,

    1179,

    2186,

    2303,

    2926,

    2237,

    525,

    735,

    863,

    2768,

    1230,

    2572,

    556,

    3010,

    2266,

    1684,

    1239,

    780,

    2954,

    109,

    1292,

    1031,

    1745,

    2688,

    3061,

    992,

    2596,

    941,

    892,

    1021,

    2390,

    642,

    1868,

    2377,

    1482,

    1540,

    540,

    1678,

    1626,

    279,

    314,

    1173,

    2573,

    3096,

    48,

    667,

    1920,

    2229,

    1041,

    2606,

    1692,

    680,

    2746,

    568,

    3312,

    2419,

    2102,

    219,

    855,

    2681,

    1848,

    712,

    682,

    927,

    1795,

    461,

    1891,

    2877,

    2522,

    1894,

    1010,

    1414,

    2009,

    3296,

    464,

    2697,

    816,

    1352,

    2679,

    1274,

    1052,

    1025,

    2132,

    1573,

    76,

    2998,

    3040,

    2508,

    1355,

    450,

    936,

    447,

    2794,

    1235,

    1903,

    1996,

    1089,

    3273,

    283,

    1853,

    1990,

    882,

    3033,

    1583,

    2760,

    69,

    543,

    2532,

    3136,

    1410,

    2267,

    2481,

    1432,

    2699,

    687,

    40,

    749,

    1600,

};

/*

 * Second precomputed array from Appendix A of FIPS 203 (normalised positive),

 * or else Python:

 * ModRoots = [pow(17, 2*bitreverse(i) + 1, p) for i in range(128)]

 */

static const uint16_t kModRoots[128] = {

    17,

    3312,

    2761,

    568,

    583,

    2746,

    2649,

    680,

    1637,

    1692,

    723,

    2606,

    2288,

    1041,

    1100,

    2229,

    1409,

    1920,

    2662,

    667,

    3281,

    48,

    233,

    3096,

    756,

    2573,

    2156,

    1173,

    3015,

    314,

    3050,

    279,

    1703,

    1626,

    1651,

    1678,

    2789,

    540,

    1789,

    1540,

    1847,

    1482,

    952,

    2377,

    1461,

    1868,

    2687,

    642,

    939,

    2390,

    2308,

    1021,

    2437,

    892,

    2388,

    941,

    733,

    2596,

    2337,

    992,

    268,

    3061,

    641,

    2688,

    1584,

    1745,

    2298,

    1031,

    2037,

    1292,

    3220,

    109,

    375,

    2954,

    2549,

    780,

    2090,

    1239,

    1645,

    1684,

    1063,

    2266,

    319,

    3010,

    2773,

    556,

    757,

    2572,

    2099,

    1230,

    561,

    2768,

    2466,

    863,

    2594,

    735,

    2804,

    525,

    1092,

    2237,

    403,

    2926,

    1026,

    2303,

    1143,

    2186,

    2150,

    1179,

    2775,

    554,

    886,

    2443,

    1722,

    1607,

    1212,

    2117,

    1874,

    1455,

    1029,

    2300,

    2110,

    1219,

    2935,

    394,

    885,

    2444,

    2154,

    1175,

};

/*

 * single_keccak hashes |inlen| bytes from |in| and writes |outlen| bytes of

 * output to |out|. If the |md| specifies a fixed-output function, like

 * SHA3-256, then |outlen| must be the correct length for that function.

 */

static __owur int single_keccak(uint8_t *out, size_t outlen, const uint8_t *in, size_t inlen,

    EVP_MD_CTX *mdctx)

{

    unsigned int sz = (unsigned int)outlen;

    if (!EVP_DigestUpdate(mdctx, in, inlen))

        return 0;

    if (EVP_MD_xof(EVP_MD_CTX_get0_md(mdctx)))

        return EVP_DigestFinalXOF(mdctx, out, outlen);

    return EVP_DigestFinal_ex(mdctx, out, &sz)

        && ossl_assert((size_t)sz == outlen);

}

/*

 * FIPS 203, Section 4.1, equation (4.3): PRF. Takes 32+1 input bytes, and uses

 * SHAKE256 to produce the input to SamplePolyCBD_eta: FIPS 203, algorithm 8.

 */

static __owur int prf(uint8_t *out, size_t len, const uint8_t in[ML_KEM_RANDOM_BYTES + 1],

    EVP_MD_CTX *mdctx, const ML_KEM_KEY *key)

{

    return EVP_DigestInit_ex(mdctx, key->shake256_md, NULL)

        && single_keccak(out, len, in, ML_KEM_RANDOM_BYTES + 1, mdctx);

}

/*

 * FIPS 203, Section 4.1, equation (4.4): H.  SHA3-256 hash of a variable

 * length input, producing 32 bytes of output.

 */

static __owur int hash_h(uint8_t out[ML_KEM_PKHASH_BYTES], const uint8_t *in, size_t len,

    EVP_MD_CTX *mdctx, const ML_KEM_KEY *key)

{

    return EVP_DigestInit_ex(mdctx, key->sha3_256_md, NULL)

        && single_keccak(out, ML_KEM_PKHASH_BYTES, in, len, mdctx);

}

/* Incremental hash_h of expanded public key */

static int

hash_h_pubkey(uint8_t pkhash[ML_KEM_PKHASH_BYTES],

    EVP_MD_CTX *mdctx, ML_KEM_KEY *key)

{

    const ML_KEM_VINFO *vinfo = key->vinfo;

    const scalar *t = key->t, *end = t + vinfo->rank;

    unsigned int sz;

    if (!EVP_DigestInit_ex(mdctx, key->sha3_256_md, NULL))

        return 0;

    do {

        uint8_t buf[3 * DEGREE / 2];

        scalar_encode(buf, t++, 12);

        if (!EVP_DigestUpdate(mdctx, buf, sizeof(buf)))

            return 0;

    } while (t < end);

    if (!EVP_DigestUpdate(mdctx, key->rho, ML_KEM_RANDOM_BYTES))

        return 0;

    return EVP_DigestFinal_ex(mdctx, pkhash, &sz)

        && ossl_assert(sz == ML_KEM_PKHASH_BYTES);

}

/*

 * FIPS 203, Section 4.1, equation (4.5): G.  SHA3-512 hash of a variable

 * length input, producing 64 bytes of output, in particular the seeds

 * (d,z) for key generation.

 */

static __owur int hash_g(uint8_t out[ML_KEM_SEED_BYTES], const uint8_t *in, size_t len,

    EVP_MD_CTX *mdctx, const ML_KEM_KEY *key)

{

    return EVP_DigestInit_ex(mdctx, key->sha3_512_md, NULL)

        && single_keccak(out, ML_KEM_SEED_BYTES, in, len, mdctx);

}

/*

 * FIPS 203, Section 4.1, equation (4.4): J. SHAKE256 taking a variable length

 * input to compute a 32-byte implicit rejection shared secret, of the same

 * length as the expected shared secret.  (Computed even on success to avoid

 * side-channel leaks).

 */

static __owur int kdf(uint8_t out[ML_KEM_SHARED_SECRET_BYTES],

    const uint8_t z[ML_KEM_RANDOM_BYTES],

    const uint8_t *ctext, size_t len,

    EVP_MD_CTX *mdctx, const ML_KEM_KEY *key)

{

    return EVP_DigestInit_ex(mdctx, key->shake256_md, NULL)

        && EVP_DigestUpdate(mdctx, z, ML_KEM_RANDOM_BYTES)

        && EVP_DigestUpdate(mdctx, ctext, len)

        && EVP_DigestFinalXOF(mdctx, out, ML_KEM_SHARED_SECRET_BYTES);

}

/*

 * FIPS 203, Section 4.2.2, Algorithm 7: "SampleNTT" (steps 3-17, steps 1, 2

 * are performed by the caller). Rejection-samples a Keccak stream to get

 * uniformly distributed elements in the range [0,q). This is used for matrix

 * expansion and only operates on public inputs.

 */

static __owur int sample_scalar(scalar *out, EVP_MD_CTX *mdctx)

{

    uint16_t *curr = out->c, *endout = curr + DEGREE;

    uint8_t buf[SCALAR_SAMPLING_BUFSIZE], *in;

    uint8_t *endin = buf + sizeof(buf);

    uint16_t d;

    uint8_t b1, b2, b3;

    do {

        if (!EVP_DigestSqueeze(mdctx, in = buf, sizeof(buf)))

            return 0;

        do {

            b1 = *in++;

            b2 = *in++;

            b3 = *in++;

            if (curr >= endout)

                break;

            if ((d = ((b2 & 0x0f) << 8) + b1) < kPrime)

                *curr++ = d;

            if (curr >= endout)

                break;

            if ((d = (b3 << 4) + (b2 >> 4)) < kPrime)

                *curr++ = d;

        } while (in < endin);

    } while (curr < endout);

    return 1;

}

/*-

 * reduce_once reduces 0 <= x < 2*kPrime, mod kPrime.

 *

 * Subtract |q| if the input is larger, without exposing a side-channel,

 * avoiding the "clangover" attack.  See |constish_time_non_zero| for a

 * discussion on why the value barrier is by default omitted.

 */

static __owur uint16_t reduce_once(uint16_t x)

{

    const uint16_t subtracted = x - kPrime;

    uint16_t mask = constish_time_non_zero(subtracted >> 15);

    return (mask & x) | (~mask & subtracted);

}

/*

 * Constant-time reduce x mod kPrime using Barrett reduction. x must be less

 * than kPrime + 2 * kPrime^2.  This is sufficient to reduce a product of

 * two already reduced u_int16 values, in fact it is sufficient for each

 * to be less than 2^12, because (kPrime * (2 * kPrime + 1)) > 2^24.

 */

static __owur uint16_t reduce(uint32_t x)

{

    uint64_t product = (uint64_t)x * kBarrettMultiplier;

    uint32_t quotient = (uint32_t)(product >> kBarrettShift);

    uint32_t remainder = x - quotient * kPrime;

    return reduce_once(remainder);

}

/* Multiply a scalar by a constant. */

static void scalar_mult_const(scalar *s, uint16_t a)

{

    uint16_t *curr = s->c, *end = curr + DEGREE, tmp;

    do {

        tmp = reduce(*curr * a);

        *curr++ = tmp;

    } while (curr < end);

}

/*-

 * FIPS 203, Section 4.3, Algorithm 9: "NTT".

 * In-place number theoretic transform of a given scalar.  Note that ML-KEM's

 * kPrime 3329 does not have a 512th root of unity, so this transform leaves

 * off the last iteration of the usual FFT code, with the 128 relevant roots of

 * unity being stored in NTTRoots.  This means the output should be seen as 128

 * elements in GF(3329^2), with the coefficients of the elements being

 * consecutive entries in |s->c|.

 */

static void scalar_ntt(scalar *s)

{

    const uint16_t *roots = kNTTRoots;

    uint16_t *end = s->c + DEGREE;

    int offset = DEGREE / 2;

    do {

        uint16_t *curr = s->c, *peer;

        do {

            uint16_t *pause = curr + offset, even, odd;

            uint32_t zeta = *++roots;

            peer = pause;

            do {

                even = *curr;

                odd = reduce(*peer * zeta);

                *peer++ = reduce_once(even - odd + kPrime);

                *curr++ = reduce_once(odd + even);

            } while (curr < pause);

        } while ((curr = peer) < end);

    } while ((offset >>= 1) >= 2);

}

/*-

 * FIPS 203, Section 4.3, Algorithm 10: "NTT^(-1)".

 * In-place inverse number theoretic transform of a given scalar, with pairs of

 * entries of s->v being interpreted as elements of GF(3329^2). Just as with

 * the number theoretic transform, this leaves off the first step of the normal

 * iFFT to account for the fact that 3329 does not have a 512th root of unity,

 * using the precomputed 128 roots of unity stored in InverseNTTRoots.

 */

static void scalar_inverse_ntt(scalar *s)

{

    const uint16_t *roots = kInverseNTTRoots;

    uint16_t *end = s->c + DEGREE;

    int offset = 2;

    do {

        uint16_t *curr = s->c, *peer;

        do {

            uint16_t *pause = curr + offset, even, odd;

            uint32_t zeta = *++roots;

            peer = pause;

            do {

                even = *curr;

                odd = *peer;

                *peer++ = reduce(zeta * (even - odd + kPrime));

                *curr++ = reduce_once(odd + even);

            } while (curr < pause);

        } while ((curr = peer) < end);

    } while ((offset <<= 1) < DEGREE);

    scalar_mult_const(s, kInverseDegree);

}

/* Addition updating the LHS scalar in-place. */

static void scalar_add(scalar *lhs, const scalar *rhs)

{

    int i;

    for (i = 0; i < DEGREE; i++)

        lhs->c[i] = reduce_once(lhs->c[i] + rhs->c[i]);

}

/* Subtraction updating the LHS scalar in-place. */

static void scalar_sub(scalar *lhs, const scalar *rhs)

{

    int i;

    for (i = 0; i < DEGREE; i++)

        lhs->c[i] = reduce_once(lhs->c[i] - rhs->c[i] + kPrime);

}

/*

 * Multiplying two scalars in the number theoretically transformed state. Since

 * 3329 does not have a 512th root of unity, this means we have to interpret

 * the 2*ith and (2*i+1)th entries of the scalar as elements of

 * GF(3329)[X]/(X^2 - 17^(2*bitreverse(i)+1)).

 *

 * The value of 17^(2*bitreverse(i)+1) mod 3329 is stored in the precomputed

 * ModRoots table. Note that our Barrett transform only allows us to multiply

 * two reduced numbers together, so we need some intermediate reduction steps,

 * even if an uint64_t could hold 3 multiplied numbers.

 */

static void scalar_mult(scalar *out, const scalar *lhs,

    const scalar *rhs)

{

    uint16_t *curr = out->c, *end = curr + DEGREE;

    const uint16_t *lc = lhs->c, *rc = rhs->c;

    const uint16_t *roots = kModRoots;

    do {

        uint32_t l0 = *lc++, r0 = *rc++;

        uint32_t l1 = *lc++, r1 = *rc++;

        uint32_t zetapow = *roots++;

        *curr++ = reduce(l0 * r0 + reduce(l1 * r1) * zetapow);

        *curr++ = reduce(l0 * r1 + l1 * r0);

    } while (curr < end);

}

/* Above, but add the result to an existing scalar */

static ossl_inline void scalar_mult_add(scalar *out, const scalar *lhs,

    const scalar *rhs)

{

    uint16_t *curr = out->c, *end = curr + DEGREE;

    const uint16_t *lc = lhs->c, *rc = rhs->c;

    const uint16_t *roots = kModRoots;

    do {

        uint32_t l0 = *lc++, r0 = *rc++;

        uint32_t l1 = *lc++, r1 = *rc++;

        uint16_t *c0 = curr++;

        uint16_t *c1 = curr++;

        uint32_t zetapow = *roots++;

        *c0 = reduce(*c0 + l0 * r0 + reduce(l1 * r1) * zetapow);

        *c1 = reduce(*c1 + l0 * r1 + l1 * r0);

    } while (curr < end);

}

/*-

 * FIPS 203, Section 4.2.1, Algorithm 5: "ByteEncode_d", for 2<=d<=12.

 * Here |bits| is |d|.  For efficiency, we handle the d=1 case separately.

 */

static void scalar_encode(uint8_t *out, const scalar *s, int bits)

{

    const uint16_t *curr = s->c, *end = curr + DEGREE;

    uint64_t accum = 0, element;

    int used = 0;

    do {

        element = *curr++;

        if (used + bits < 64) {

            accum |= element << used;

            used += bits;

        } else if (used + bits > 64) {

            out = OPENSSL_store_u64_le(out, accum | (element << used));

            accum = element >> (64 - used);

            used = (used + bits) - 64;

        } else {

            out = OPENSSL_store_u64_le(out, accum | (element << used));

            accum = 0;

            used = 0;

        }

    } while (curr < end);

}

/*

 * scalar_encode_1 is |scalar_encode| specialised for |bits| == 1.

 */

static void scalar_encode_1(uint8_t out[DEGREE / 8], const scalar *s)

{

    int i, j;

    uint8_t out_byte;

    for (i = 0; i < DEGREE; i += 8) {

        out_byte = 0;

        for (j = 0; j < 8; j++)

            out_byte |= bit0(s->c[i + j]) << j;

        *out = out_byte;

        out++;

    }

}

/*-

 * FIPS 203, Section 4.2.1, Algorithm 6: "ByteDecode_d", for 2<=d<12.

 * Here |bits| is |d|.  For efficiency, we handle the d=1 and d=12 cases

 * separately.

 *

 * scalar_decode parses |DEGREE * bits| bits from |in| into |DEGREE| values in

 * |out|.

 */

static void scalar_decode(scalar *out, const uint8_t *in, int bits)

{

    uint16_t *curr = out->c, *end = curr + DEGREE;

    uint64_t accum = 0;

    int accum_bits = 0, todo = bits;

    uint16_t bitmask = (((uint16_t)1) << bits) - 1, mask = bitmask;

    uint16_t element = 0;

    do {

        if (accum_bits == 0) {

            in = OPENSSL_load_u64_le(&accum, in);

            accum_bits = 64;

        }

        if (todo == bits && accum_bits >= bits) {

            /* No partial "element", and all the required bits available */

            *curr++ = ((uint16_t)accum) & mask;

            accum >>= bits;

            accum_bits -= bits;

        } else if (accum_bits >= todo) {

            /* A partial "element", and all the required bits available */

            *curr++ = element | ((((uint16_t)accum) & mask) << (bits - todo));

            accum >>= todo;

            accum_bits -= todo;

            element = 0;

            todo = bits;

            mask = bitmask;

        } else {

            /*

             * Only some of the requisite bits accumulated, store |accum_bits|

             * of these in |element|.  The accumulated bitcount becomes 0, but

             * as soon as we have more bits we'll want to merge accum_bits

             * fewer of them into the final |element|.

             *

             * Note that with a 64-bit accumulator and |bits| always 12 or

             * less, if we're here, the previous iteration had all the

             * requisite bits, and so there are no kept bits in |element|.

             */

            element = ((uint16_t)accum) & mask;

            todo -= accum_bits;

            mask = bitmask >> accum_bits;

            accum_bits = 0;

        }

    } while (curr < end);

}

static __owur int scalar_decode_12(scalar *out, const uint8_t in[3 * DEGREE / 2])

{

    int i;

    uint16_t *c = out->c;

    for (i = 0; i < DEGREE / 2; ++i) {

        uint8_t b1 = *in++;

        uint8_t b2 = *in++;

        uint8_t b3 = *in++;

        int outOfRange1 = (*c++ = b1 | ((b2 & 0x0f) << 8)) >= kPrime;

        int outOfRange2 = (*c++ = (b2 >> 4) | (b3 << 4)) >= kPrime;

        if (outOfRange1 | outOfRange2)

            return 0;

    }

    return 1;

}

/*-

 * scalar_decode_decompress_add is a combination of decoding and decompression

 * both specialised for |bits| == 1, with the result added (and sum reduced) to

 * the output scalar.

 *

 * NOTE: this function MUST not leak an input-data-depedennt timing signal.

 * A timing leak in a related function in the reference Kyber implementation

 * made the "clangover" attack (CVE-2024-37880) possible, giving key recovery

 * for ML-KEM-512 in minutes, provided the attacker has access to precise

 * timing of a CPU performing chosen-ciphertext decap.  Admittedly this is only

 * a risk when private keys are reused (perhaps KEMTLS servers).

 */

static void

scalar_decode_decompress_add(scalar *out, const uint8_t in[DEGREE / 8])

{

    static const uint16_t half_q_plus_1 = (ML_KEM_PRIME >> 1) + 1;

    uint16_t *curr = out->c, *end = curr + DEGREE;

    uint16_t mask;

    uint8_t b;

    /*

     * Add |half_q_plus_1| if the bit is set, without exposing a side-channel,

     * avoiding the "clangover" attack.  See |constish_time_non_zero| for a

     * discussion on why the value barrier is by default omitted.

     */

#define decode_decompress_add_bit                        \

    mask = constish_time_non_zero(bit0(b));              \

    *curr = reduce_once(*curr + (mask & half_q_plus_1)); \

    curr++;                                              \

    b >>= 1