Review comments

Co-Authored-By: Fabian Albert <[email protected]>
Co-Authored-By: Jack Lloyd <[email protected]>

src/lib/pubkey/dilithium/dilithium_common/dilithium_algos.h

@@ -17,8 +17,10 @@
 namespace Botan {
 
 DilithiumPolyMatNTT dilithium_expand_A(StrongSpan<const DilithiumSeedRho> rho, const DilithiumConstants& mode);
+
 std::pair<DilithiumPolyVec, DilithiumPolyVec> dilithium_expand_s(StrongSpan<const DilithiumSeedRhoPrime> rhoprime,
                                                                  const DilithiumConstants& mode);
+
 DilithiumPolyVec dilithium_expand_mask(StrongSpan<const DilithiumSeedRhoPrime> rhoprime,
                                        uint16_t nonce,
                                        const DilithiumConstants& mode);

src/lib/pubkey/kyber/kyber_common/kyber.cpp

@@ -210,15 +210,11 @@ Kyber_PrivateKey::Kyber_PrivateKey(RandomNumberGenerator& rng, KyberMode m) {
    auto [rho, sigma] = mode.symmetric_primitives().G(d);
    KyberPolynomialSampler ps(sigma, mode);
 
-   // TODO: Remove the need for the montgomery transformation
-   //
-   // -> When calculating A*s below, A is not in montgomery form, but s is. The
-   //    operation uses fqmul internally, which performs a montgomery reduction.
-   auto A = montgomery(kyber_sample_matrix(rho, false /* not transposed */, mode));
+   auto A = kyber_sample_matrix(rho, false /* not transposed */, mode);
    auto s = ntt(ps.sample_polynomial_vector_cbd_eta1());
    const auto e = ntt(ps.sample_polynomial_vector_cbd_eta1());
 
-   auto t = A * s;
+   auto t = montgomery(A * s);
    t += e;
    t.reduce();
 

src/lib/pubkey/kyber/kyber_common/kyber_algos.h

@@ -24,16 +24,23 @@
 namespace Botan {
 
 void kyber_encode_polynomial_vector(std::span<uint8_t> out, const KyberPolyVecNTT& p);
+
 KyberPolyVecNTT kyber_decode_polynomial_vector(std::span<const uint8_t> a, const KyberConstants& mode);
+
 KyberPoly kyber_polynomial_from_message(StrongSpan<const KyberMessage> msg);
+
 KyberMessage kyber_polynomial_to_message(const KyberPoly& p);
+
 void kyber_compress_ciphertext(StrongSpan<KyberCompressedCiphertext> out,
                                const KyberPolyVec& u,
                                const KyberPoly& v,
                                const KyberConstants& m_mode);
+
 std::pair<KyberPolyVec, KyberPoly> kyber_decompress_ciphertext(StrongSpan<const KyberCompressedCiphertext> ct,
                                                                const KyberConstants& mode);
+
 KyberPolyMat kyber_sample_matrix(StrongSpan<const KyberSeedRho> seed, bool transposed, const KyberConstants& mode);
+
 void kyber_sample_polynomial_from_cbd(KyberPoly& poly,
                                       KyberConstants::KyberEta eta,
                                       const KyberSamplingRandomness& randomness);

src/lib/pubkey/kyber/kyber_common/kyber_helpers.h

@@ -32,30 +32,29 @@ template <size_t d>
    requires(d > 0 && d < 12)
 constexpr std::make_unsigned_t<KyberConstants::T> kyber_compress(KyberConstants::T x) {
    BOTAN_DEBUG_ASSERT(x >= 0 && x < KyberConstants::Q);
-   using unsigned_T = std::make_unsigned_t<KyberConstants::T>;
+   const uint32_t n = (static_cast<uint32_t>(x) << d) + KyberConstants::Q / 2;
 
    // This is a mitigation for a potential side channel called "KyberSlash".
+   //
    // It implements the division by Q using a multiplication and a shift. Most
    // compilers would generate similar code for such a division by a constant.
    // Though, in some cases, compilers might use a variable-time int division,
    // resulting in a potential side channel.
    //
+   // The constants below work for all values that appear in Kyber with the
+   // greatest being 3328 * 2^11 + Q // 2 = 6,817,408 < 2**23 = 8,388,608.
+   //
    //   See "Hacker's Delight" (Second Edition) by Henry S. Warren, Jr.
    //   Chapter 10-9 "Unsigned Division by Divisors >= 1"
-   auto divide_by_q = [](uint32_t n) -> unsigned_T {
-      static_assert(KyberConstants::Q == 3329);
-      BOTAN_DEBUG_ASSERT(n < (1 << 23));
-
-      // These constants work for all values that appear in Kyber with the
-      // greatest being 3328 * 2^11 + Q // 2 = 6,817,408 < 2**23 = 8,388,608.
-      constexpr uint64_t m = 2580335;
-      constexpr size_t p = 33;
-      return static_cast<unsigned_T>((n * m) >> p);
-   };
+   BOTAN_DEBUG_ASSERT(n < (1 << 23));
+   static_assert(KyberConstants::Q == 3329);
+   using unsigned_T = std::make_unsigned_t<KyberConstants::T>;
 
+   constexpr uint64_t m = 2580335;
+   constexpr size_t p = 33;
    constexpr unsigned_T mask = (1 << d) - 1;
-   return divide_by_q((static_cast<uint32_t>(x) << d) + KyberConstants::Q / 2) & mask;
-}
+   return static_cast<unsigned_T>((n * m) >> p) & mask;
+};
 
 /**
  * NIST FIPS 203 IPD, Formula 4.6 (Decompress)

src/lib/pubkey/kyber/kyber_common/kyber_polynomial.h

@@ -43,6 +43,9 @@ class KyberPolyTraits final : public CRYSTALS::Trait_Base<KyberConstants, KyberP
    public:
       /**
        * NIST FIPS 203 IPD, Algorithm 8 (NTT)
+       *
+       * Produces the result of the NTT transformation without any montgomery
+       * factors in the coefficients.
        */
       constexpr static void ntt(std::span<T, N> p) {
          for(size_t len = N / 2, k = 0; len >= 2; len /= 2) {
@@ -61,6 +64,14 @@ class KyberPolyTraits final : public CRYSTALS::Trait_Base<KyberConstants, KyberP
 
       /**
        * NIST FIPS 203 IPD, Algorithm 9 (NTT^-1)
+       *
+       * The output is effectively multiplied by the montgomery parameter 2^16
+       * mod q so that the input factors 2^(-16) mod q are eliminated. Note
+       * that factors 2^(-16) mod q are introduced by multiplication and
+       * reduction of values not in montgomery domain.
+       *
+       * Produces the result of the inverse NTT transformation with a montgomery
+       * factor of (2^16 mod q) added (!). See above.
        */
       static constexpr void inverse_ntt(std::span<T, N> p) {
          for(size_t len = 2, k = 127; len <= N / 2; len *= 2) {

src/lib/pubkey/pqcrystals/pqcrystals.h

@@ -587,16 +587,6 @@ PolynomialVector<Trait, D> montgomery(PolynomialVector<Trait, D> polyvec) {
    return polyvec;
 }
 
-template <crystals_trait Trait>
-PolynomialMatrix<Trait> montgomery(PolynomialMatrix<Trait> m) {
-   for(auto& polyvec : m) {
-      for(auto& poly : polyvec) {
-         detail::montgomery(poly);
-      }
-   }
-   return m;
-}
-
 template <crystals_trait Trait>
 PolynomialVector<Trait, Domain::Normal> operator+(const PolynomialVector<Trait, Domain::Normal>& a,
                                                   const PolynomialVector<Trait, Domain::Normal>& b) {

src/lib/pubkey/pqcrystals/pqcrystals_encoding.h

@@ -12,15 +12,17 @@
 #define BOTAN_PQ_CRYSTALS_ENCODING_H_
 
 #include <limits>
+#include <numeric>
 #include <span>
 
-#include <botan/xof.h>
-
 #include <botan/internal/loadstor.h>
 #include <botan/internal/pqcrystals.h>
 #include <botan/internal/pqcrystals_helpers.h>
 #include <botan/internal/stl_util.h>
 
+#if defined(BOTAN_HAS_XOF)
+   #include <botan/xof.h>
+#endif
 namespace Botan::CRYSTALS {
 
 namespace detail {
@@ -29,9 +31,11 @@ constexpr auto as_byte_source(BufferSlicer& slicer) {
    return [&](std::span<uint8_t> out) { slicer.copy_into(out); };
 }
 
+#if defined(BOTAN_HAS_XOF)
 constexpr auto as_byte_source(Botan::XOF& xof) {
    return [&](std::span<uint8_t> out) { xof.output(out); };
 }
+#endif
 
 }  // namespace detail
 
@@ -71,7 +75,7 @@ struct BitPackingTrait final {
       constexpr static size_t bits_per_pack = [] {
          // Ensure that the bit-packing is byte-aligned and scale it
          // to utilize the collector's bit-width as much as possible.
-         size_t smallest_aligned_pack = lcm(bits_per_coeff, size_t(8));
+         size_t smallest_aligned_pack = std::lcm(bits_per_coeff, size_t(8));
          return (smallest_aligned_pack < bits_in_collector)
                    ? (bits_in_collector / smallest_aligned_pack) * smallest_aligned_pack
                    : smallest_aligned_pack;

src/lib/pubkey/pqcrystals/pqcrystals_helpers.h

@@ -39,15 +39,15 @@ using next_longer_int_t =
 
 template <std::integral T>
    requires(size_t(sizeof(T)) <= 4)
-constexpr T montgomery_R(T q) {
+consteval T montgomery_R(T q) {
    using T_unsigned = std::make_unsigned_t<T>;
    using T2 = next_longer_uint_t<T_unsigned>;
    return (T2(1) << (sizeof(T) * 8)) % q;
 }
 
 template <std::integral T>
    requires(size_t(sizeof(T)) <= 4)
-constexpr T montgomery_R2(T q) {
+consteval T montgomery_R2(T q) {
    using T2 = next_longer_int_t<T>;
    return (static_cast<T2>(montgomery_R(q)) * static_cast<T2>(montgomery_R(q))) % q;
 }
@@ -64,7 +64,7 @@ struct eea_result {
  * and b and the Bézout coefficients, u and v.
  */
 template <std::integral T>
-constexpr eea_result<T> extended_euclidean_algorithm(T a, T b) {
+consteval eea_result<T> extended_euclidean_algorithm(T a, T b) {
    if(a > b) {
       std::swap(a, b);
    }
@@ -89,15 +89,10 @@ constexpr eea_result<T> extended_euclidean_algorithm(T a, T b) {
  */
 template <std::integral T, std::integral T2 = next_longer_int_t<T>>
    requires(sizeof(T) <= 4)
-constexpr T modular_inverse(T q, T2 m = T2(1) << sizeof(T) * 8) {
+consteval T modular_inverse(T q, T2 m = T2(1) << sizeof(T) * 8) {
    return static_cast<T>(extended_euclidean_algorithm<T2>(q, m).u);
 }
 
-template <std::unsigned_integral T>
-constexpr T lcm(T a, T b) {
-   return a / extended_euclidean_algorithm(a, b).gcd * b;
-}
-
 constexpr auto bitlen(size_t x) {
    return ceil_log2(x + 1);
 };

src/lib/rng/rng.h

@@ -14,7 +14,6 @@
 #include <botan/mutex.h>
 #include <botan/secmem.h>
 
-#include <array>
 #include <chrono>
 #include <concepts>
 #include <span>