Small prime ntt

libs/math/include/nil/crypto3/math/detail/small_prime_fft_arithmetic.hpp

@@ -0,0 +1,182 @@
+//---------------------------------------------------------------------------//
+// Copyright (c) 2026
+//
+// MIT License
+//
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included in all
+// copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+// SOFTWARE.
+//---------------------------------------------------------------------------//
+
+#ifndef CRYPTO3_MATH_DETAIL_SMALL_PRIME_FFT_ARITHMETIC_HPP
+#define CRYPTO3_MATH_DETAIL_SMALL_PRIME_FFT_ARITHMETIC_HPP
+
+#include <cassert>
+#include <cfenv>
+#include <cfloat>
+#include <cmath>
+#include <cstdint>
+#include <limits>
+#include <stdexcept>
+
+#include <nil/crypto3/multiprecision/cpp_int_modular.hpp>
+
+namespace nil::crypto3::math::detail {
+
+    /**
+     * Scalar residue arithmetic for an odd 50-bit modulus q. Values are integer-valued binary64
+     * numbers, not approximate coefficients. Primality is a caller concern (the prime catalogue
+     * is separately certified); these arithmetic operations also work for odd composite moduli.
+     *
+     * Select the kernel once per operation with fma_supported(). The FMA methods require this
+     * predicate to remain true until the operation finishes, including FE_TONEAREST on the calling
+     * thread. The integer methods are the exact fallback for other rounding modes/compiler flags.
+     * Keeping this dispatch outside each scalar operation avoids a rounding-mode query per butterfly.
+     * No method changes the floating-point environment.
+     */
+    class small_prime_fft_arithmetic {
+        using wide_type = boost::multiprecision::uint128_modular_t;
+
+        static_assert(std::numeric_limits<double>::radix == 2 && std::numeric_limits<double>::digits == 53 &&
+                          std::numeric_limits<double>::max_exponent == 1024,
+                      "small-prime residues require binary64 doubles");
+
+    public:
+        explicit small_prime_fft_arithmetic(std::uint64_t modulus) : modulus_(modulus) {
+            if (modulus < (std::uint64_t(1) << 49) || modulus >= (std::uint64_t(1) << 50) || modulus % 2 == 0) {
+                throw std::invalid_argument("small_prime_fft_arithmetic: expected an odd 50-bit modulus");
+            }
+
+            // Since 2^49 < q < 2^50, the nearest binary64 reciprocal is S/2^102, where
+            // S = round(2^102/q) fits exactly in 53 bits. An odd q cannot produce a halfway tie.
+            // Construct S using integers so that setup is independent of the active rounding mode.
+            const wide_type scale = wide_type(1) << 102;
+            const wide_type significand = (scale + modulus / 2) / modulus;
+            reciprocal_ = std::ldexp(static_cast<double>(significand.convert_to<std::uint64_t>()), -102);
+
+            // For h = RN(a*b), l = fma(a,b,-h), and k = nearbyint(h*reciprocal), let
+            // epsilon = reciprocal - 1/q and |a*b| < 2^e. Then h+l is exact and
+            // |a*b/q - k| <= |a*b*epsilon| + 2^(e-54)*reciprocal + rounding_error.
+            // For |a*b| < 2*q^2 the rounding_error is at most 5/8. With
+            // D = |q*S - 2^102|, the remaining numerator over 2^102 is
+            // T = 2*q*D + 2^(e-54)*S. Requiring T < 3*2^99 makes the total < 1.
+            // For |a*b| < 4*q^2, e increases by one and the rounding_error is 3/4:
+            // the condition becomes 2*T < 3*2^100, giving an output strictly inside (-1.5*q,1.5*q).
+            // All certificate intermediates fit within 128 bits.
+            const wide_type scaled_reciprocal = significand * modulus;
+            const wide_type error = scaled_reciprocal > scale ? scaled_reciprocal - scale : scale - scaled_reciprocal;
+            const wide_type product_limit = 2 * wide_type(modulus) * modulus;
+            const unsigned exponent = boost::multiprecision::msb(product_limit) + 1;
+            const wide_type numerator = 2 * wide_type(modulus) * error + (significand << (exponent - 54));
+            reciprocal_satisfies_bounds_ = numerator < (wide_type(3) << 99);
+        }
+
+        std::uint64_t modulus() const {
+            return modulus_;
+        }
+
+        double reciprocal() const {
+            return reciprocal_;
+        }
+
+        /** Check once before calling FMA kernels; use the integer kernels when this returns false. */
+        bool fma_supported() const {
+#if defined(__FAST_MATH__) || defined(__ASSOCIATIVE_MATH__) || defined(__RECIPROCAL_MATH__) || defined(_M_FP_FAST) || \
+    (defined(__FINITE_MATH_ONLY__) && __FINITE_MATH_ONLY__ > 0) || FLT_EVAL_METHOD != 0
+            return false;
+#else
+            return std::numeric_limits<double>::is_iec559 && reciprocal_satisfies_bounds_ &&
+                   std::fegetround() == FE_TONEAREST;
+#endif
+        }
+
+        /**
+         * Return an integer representative of value modulo q strictly inside (-q,q).
+         * @pre value is an integer with |value| < 2^53, and fma_supported() is true.
+         */
+        double reduce_fma(double value) const {
+#if defined(__clang__)
+            // Clang does not expose every unsafe floating-point option through predefined macros.
+#pragma clang fp reassociate(off)
+#pragma clang fp contract(off)
+#endif
+            assert(valid_integer(value));
+            assert(fma_supported());
+            // |value/q| < 16, so quotient rounding error is < 1. The FMA remainder is an
+            // exactly representable integer of magnitude < q, even for a large input representative.
+            const double quotient = std::nearbyint(value * reciprocal_);
+            return std::fma(-quotient, static_cast<double>(modulus_), value);
+        }
+
+        /**
+         * Return an integer representative of left*right modulo q. The output is strictly inside
+         * (-q,q) when |left*right| < 2*q^2, and inside (-1.5*q,1.5*q) when |left*right| < 4*q^2.
+         * @pre both operands are integers of magnitude < 2^53, |left*right| < 4*q^2, and
+         *      fma_supported() is true. These are integer product bounds, not rounded double tests.
+         */
+        double multiply_fma(double left, double right) const {
+#if defined(__clang__)
+            // Preserve the rounded high product even when the caller enables reassociation/contraction.
+#pragma clang fp reassociate(off)
+#pragma clang fp contract(off)
+#endif
+            assert(valid_integer(left) && valid_integer(right));
+            assert(product_magnitude(left, right) < 4 * wide_type(modulus_) * modulus_);
+            assert(fma_supported());
+            const double high = left * right;
+            const double low = std::fma(left, right, -high);
+            const double quotient = std::nearbyint(high * reciprocal_);
+            // Both the compensated subtraction and the final remainder are small enough to be
+            // exact integers in binary64. Preserve the split product: reassociation is not allowed.
+            return std::fma(-quotient, static_cast<double>(modulus_), high) + low;
+        }
+
+        /** Exact fallback for reduce_fma, with the same input/output bounds and any rounding mode. */
+        double reduce_integer(double value) const {
+            assert(valid_integer(value));
+            return static_cast<double>(static_cast<std::int64_t>(value) % static_cast<std::int64_t>(modulus_));
+        }
+
+        /**
+         * Exact fallback: both operands must be integers of magnitude < 2^53. Their full product
+         * fits in 106 bits, so the existing 128-bit integer type suffices. The result is inside (-q,q).
+         */
+        double multiply_integer(double left, double right) const {
+            assert(valid_integer(left) && valid_integer(right));
+            const wide_type residue = product_magnitude(left, right) % modulus_;
+            const double result = static_cast<double>(residue.convert_to<std::uint64_t>());
+            return (left < 0) != (right < 0) ? -result : result;
+        }
+
+    private:
+        static bool valid_integer(double value) {
+            return std::isfinite(value) && std::abs(value) < 0x1p53 && std::trunc(value) == value;
+        }
+
+        // Call only after checking the integer/range preconditions, so both conversions are exact.
+        static wide_type product_magnitude(double left, double right) {
+            return wide_type(static_cast<std::uint64_t>(std::abs(left))) * static_cast<std::uint64_t>(std::abs(right));
+        }
+
+        std::uint64_t modulus_;
+        double reciprocal_;
+        bool reciprocal_satisfies_bounds_;
+    };
+
+}    // namespace nil::crypto3::math::detail
+
+#endif    // CRYPTO3_MATH_DETAIL_SMALL_PRIME_FFT_ARITHMETIC_HPP

libs/math/include/nil/crypto3/math/detail/small_prime_fft_params.hpp

@@ -0,0 +1,188 @@
+//---------------------------------------------------------------------------//
+// Copyright (c) 2026
+//
+// MIT License
+//
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included in all
+// copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+// SOFTWARE.
+//---------------------------------------------------------------------------//
+
+#ifndef CRYPTO3_MATH_DETAIL_SMALL_PRIME_FFT_PARAMS_HPP
+#define CRYPTO3_MATH_DETAIL_SMALL_PRIME_FFT_PARAMS_HPP
+
+#include <algorithm>
+#include <array>
+#include <bit>
+#include <cstddef>
+#include <cstdint>
+#include <limits>
+#include <optional>
+#include <utility>
+
+#include <boost/multiprecision/cpp_int.hpp>
+
+#include <nil/crypto3/algebra/fields/alt_bn128/base_field.hpp>
+#include <nil/crypto3/algebra/fields/field_order.hpp>
+
+namespace nil::crypto3::math::detail {
+
+    struct small_prime_fft_prime {
+        std::uint64_t modulus;
+        std::uint64_t root;    // Exact order 2^two_adicity; shorter roots are obtained by squaring.
+        unsigned two_adicity;
+    };
+
+    // These distinct primes q = k*2^s + 1 (k odd) balance residue arithmetic, FFT length and CRT capacity:
+    //  - 2^49 < q < 2^50 matches the FMA kernel's certified range and leaves three headroom bits in
+    //    binary64's 53-bit significand for exact integer-valued lazy sums with magnitude below 8*q.
+    //    Smaller primes contribute fewer CRT bits, generally requiring more primes and transforms
+    //    to reconstruct the same full-width BN254 Fp12 product.
+    //  - s ranges from 41 to 44, providing roots for every power-of-two transform length up to 2^s.
+    //    Entries are ordered by decreasing s, then decreasing q, favoring root capacity first and
+    //    more CRT bits among primes with the same root capacity.
+    //  - Twelve entries provide spare CRT capacity for larger convolutions; the 32768-point
+    //    full-width BN254 Fp12 workload needs eleven. Selection takes only the shortest prefix whose
+    //    product M exceeds twice the signed coefficient bound, so unused entries incur no transforms.
+    // Other prime/root choices meeting these conditions can also work. Tests certify this catalogue's
+    // primality, exact root orders and FMA reciprocal error bounds; arithmetic checks remain separate
+    // from these integer reconstruction parameters, which use no floating-point arithmetic.
+    inline constexpr std::array<small_prime_fft_prime, 12> small_prime_fft_primes = {{
+        {1108307720798209ULL, 194751219211145ULL, 44},
+        {659706976665601ULL, 422061946989193ULL, 43},
+        {1086317488242689ULL, 813873581740013ULL, 42},
+        {910395627798529ULL, 212465101045679ULL, 42},
+        {699289395265537ULL, 643876812462936ULL, 42},
+        {1022545813831681ULL, 136965991847555ULL, 41},
+        {1013749720809473ULL, 549289857321204ULL, 41},
+        {868614185943041ULL, 322646709537195ULL, 41},
+        {864216139431937ULL, 46940759627233ULL, 41},
+        {837827860365313ULL, 128509036413366ULL, 41},
+        {815837627809793ULL, 584069832576029ULL, 41},
+        {600333348765697ULL, 295748699958414ULL, 41},
+    }};
+
+    struct small_prime_fft_crt_params {
+        std::size_t prime_count;
+        boost::multiprecision::cpp_int modulus;
+    };
+
+    /**
+     * Select the shortest nonempty catalogue prefix with M > 2*coefficient_bound, where every
+     * prime supports transform_size. Then each integer in [-coefficient_bound, coefficient_bound]
+     * has a unique centered CRT representative. The bound must cover the integer expression being
+     * reconstructed, before reduction modulo the target field characteristic.
+     *
+     * Return nullopt for a negative bound, a non-power-of-two size (including zero), insufficient
+     * root order, or insufficient CRT capacity. Arbitrary-precision setup arithmetic avoids overflow.
+     */
+    inline std::optional<small_prime_fft_crt_params>
+        select_small_prime_fft_crt_params(std::size_t transform_size,
+                                          const boost::multiprecision::cpp_int &coefficient_bound) {
+        if (!std::has_single_bit(transform_size) || coefficient_bound < 0) {
+            return std::nullopt;
+        }
+        const unsigned depth = std::countr_zero(transform_size);
+        const boost::multiprecision::cpp_int twice_bound = 2 * coefficient_bound;
+        boost::multiprecision::cpp_int modulus = 1;
+        for (std::size_t i = 0; i < small_prime_fft_primes.size(); ++i) {
+            if (depth > small_prime_fft_primes[i].two_adicity) {
+                return std::nullopt;
+            }
+            modulus *= small_prime_fft_primes[i].modulus;
+            if (modulus > twice_bound) {
+                return small_prime_fft_crt_params {i + 1, std::move(modulus)};
+            }
+        }
+        return std::nullopt;
+    }
+
+    /**
+     * Bound one coordinate of a sum of term_count BN254 Fp12 products, with input coordinates
+     * lifted canonically into [0, p-1]. In the tower u^2=-1, v^3=9+u, w^2=v, use the equivalent
+     * basis w^j*u^k (0<=j<6, 0<=k<2). At output w^j there are j+1 unwrapped pairs of w powers,
+     * with absolute coordinate weight 2 each, and 5-j wrapped pairs, with weight 2*(9+1) each.
+     * Thus each output coordinate has weight 2*(j+1)+20*(5-j) <= 102. The triangle inequality
+     * gives C = 102*term_count*(p-1)^2, including the signs introduced by u^2=-1.
+     *
+     * This is a bound on the final integer tower product, not on lazy arithmetic intermediates.
+     * Cyclic input folding, digit splitting and other input lifts require their own bounds.
+     */
+    inline boost::multiprecision::cpp_int bn254_fp12_convolution_bound(std::size_t term_count) {
+        const boost::multiprecision::cpp_int largest_coordinate =
+            algebra::fields::field_characteristic<algebra::fields::alt_bn128_base_field<254>>() - 1;
+        return 102 * boost::multiprecision::cpp_int(term_count) * largest_coordinate * largest_coordinate;
+    }
+
+    struct small_prime_fft_params {
+        std::size_t left_size;
+        std::size_t right_size;
+        std::size_t result_size;
+        std::size_t transform_size;
+        std::size_t spectrum_bytes;    // One spectrum: 12 coordinates * prime_count * transform_size 64-bit cells.
+        boost::multiprecision::cpp_int coefficient_bound;
+        small_prime_fft_crt_params crt;
+    };
+
+    /**
+     * Plan a full BN254 Fp12 product, or its first coefficient_count coefficients. Only input prefixes
+     * that can affect that output are used, but the transform still contains their COMPLETE linear
+     * product: shortening it to the output count would permit cyclic aliasing into low coefficients.
+     * Square uses the same plan with equal operand lengths. No transform buffers are allocated here.
+     *
+     * Return nullopt when either input length or coefficient_count is zero (the caller handles a zero
+     * result without a transform), when size arithmetic would overflow, when the catalogue cannot
+     * supply the roots/CRT capacity, or when one spectrum exceeds max_spectrum_bytes. This byte limit
+     * covers only one spectrum, not all future operation scratch, twiddle tables or prepared operands.
+     */
+    inline std::optional<small_prime_fft_params>
+        make_bn254_small_prime_fft_params(std::size_t left_size, std::size_t right_size,
+                                          std::size_t coefficient_count = std::numeric_limits<std::size_t>::max(),
+                                          std::size_t max_spectrum_bytes = std::numeric_limits<std::size_t>::max()) {
+        left_size = std::min(left_size, coefficient_count);
+        right_size = std::min(right_size, coefficient_count);
+        if (left_size == 0 || right_size == 0 ||
+            left_size > std::numeric_limits<std::size_t>::max() - (right_size - 1)) {
+            return std::nullopt;
+        }
+        const std::size_t product_size = left_size + (right_size - 1);
+        // std::bit_ceil requires the result to be representable in size_t.
+        if (product_size > (std::size_t(1) << (std::numeric_limits<std::size_t>::digits - 1))) {
+            return std::nullopt;
+        }
+        const std::size_t transform_size = std::bit_ceil(product_size);
+        boost::multiprecision::cpp_int bound = bn254_fp12_convolution_bound(std::min(left_size, right_size));
+        auto crt = select_small_prime_fft_crt_params(transform_size, bound);
+        if (!crt) {
+            return std::nullopt;
+        }
+        const std::size_t bytes_per_index = 12 * crt->prime_count * sizeof(std::uint64_t);
+        if (transform_size > max_spectrum_bytes / bytes_per_index) {
+            return std::nullopt;
+        }
+        return small_prime_fft_params {left_size,
+                                       right_size,
+                                       std::min(product_size, coefficient_count),
+                                       transform_size,
+                                       transform_size * bytes_per_index,
+                                       std::move(bound),
+                                       std::move(*crt)};
+    }
+
+}    // namespace nil::crypto3::math::detail
+
+#endif    // CRYPTO3_MATH_DETAIL_SMALL_PRIME_FFT_PARAMS_HPP