From 25b4848155772908d8e3979d7b43806617689156 Mon Sep 17 00:00:00 2001 From: Riccardo Abbate Date: Thu, 22 Oct 2026 12:20:14 -0400 Subject: [PATCH] Accelerate BN254 CRT export with AVX-512 IFMA --- .../detail/element/fp12_fast/limb_ops_x86.hpp | 20 ++ .../math/detail/bn254_small_prime_fft.hpp | 199 +++++++++++++++++ libs/math/test/small_prime_fft_crt.cpp | 211 +++++++++++++++++- 3 files changed, 429 insertions(+), 1 deletion(-) diff --git a/libs/algebra/include/nil/crypto3/algebra/fields/detail/element/fp12_fast/limb_ops_x86.hpp b/libs/algebra/include/nil/crypto3/algebra/fields/detail/element/fp12_fast/limb_ops_x86.hpp index 1ccbec8a3..560565e43 100644 --- a/libs/algebra/include/nil/crypto3/algebra/fields/detail/element/fp12_fast/limb_ops_x86.hpp +++ b/libs/algebra/include/nil/crypto3/algebra/fields/detail/element/fp12_fast/limb_ops_x86.hpp @@ -1427,6 +1427,16 @@ namespace nil::crypto3::algebra::fields::detail::fp12_fast { static value_type shift52(value_type a) { return {_mm256_srli_epi64(a.data, 52)}; } + template + static value_type shift_right(value_type a) { + static_assert(Bits < 64); + return {_mm256_srli_epi64(a.data, Bits)}; + } + template + static value_type shift_left(value_type a) { + static_assert(Bits < 64); + return {_mm256_slli_epi64(a.data, Bits)}; + } static __mmask8 less(value_type a, value_type b) { return _mm256_cmp_epu64_mask(a.data, b.data, _MM_CMPINT_LT); } @@ -1468,6 +1478,16 @@ namespace nil::crypto3::algebra::fields::detail::fp12_fast { static value_type shift52(value_type a) { return {_mm512_srli_epi64(a.data, 52)}; } + template + static value_type shift_right(value_type a) { + static_assert(Bits < 64); + return {_mm512_srli_epi64(a.data, Bits)}; + } + template + static value_type shift_left(value_type a) { + static_assert(Bits < 64); + return {_mm512_slli_epi64(a.data, Bits)}; + } static __mmask8 less(value_type a, value_type b) { return _mm512_cmp_epu64_mask(a.data, b.data, _MM_CMPINT_LT); } diff --git a/libs/math/include/nil/crypto3/math/detail/bn254_small_prime_fft.hpp b/libs/math/include/nil/crypto3/math/detail/bn254_small_prime_fft.hpp index d17eb5397..33a919c63 100644 --- a/libs/math/include/nil/crypto3/math/detail/bn254_small_prime_fft.hpp +++ b/libs/math/include/nil/crypto3/math/detail/bn254_small_prime_fft.hpp @@ -54,8 +54,108 @@ static_assert(CRYPTO3_SMALL_PRIME_FFT_RADIX == 4 || CRYPTO3_SMALL_PRIME_FFT_RADI CRYPTO3_SMALL_PRIME_FFT_RADIX == 16, "small-prime FFT stage radix must be 4, 8 or 16"); +#if defined(CRYPTO3_HAS_AVX512IFMA) && defined(CRYPTO3_SMALL_PRIME_FFT_AVX512_AVAILABLE) && \ + !defined(CRYPTO3_SMALL_PRIME_CRT_DISABLE_IFMA) +#define CRYPTO3_SMALL_PRIME_CRT_IFMA_SELECTED +#endif + namespace nil::crypto3::math::detail { + // Independent coefficient lanes, using the existing IFMA operation policy. A CRT digit is + // below 2^50 and each five-digit export weight is below p. With at most twelve primes, + // every unnormalized column is below 24*2^52 and the complete sum is below 2^54*p. + template + struct small_prime_fft_crt_accumulator { + using radix = algebra::fields::detail::fp12_fast::montgomery_radix52; + using vector = typename Operations::value_type; + std::array columns; + + small_prime_fft_crt_accumulator() { + columns.fill(Operations::broadcast(0)); + } + + void add(const typename radix::value_type &weight, vector digits) { + for (std::size_t limb = 0; limb < weight.size(); ++limb) { + const auto w = Operations::broadcast(weight[limb]); + columns[limb] = Operations::multiply_add_low(columns[limb], digits, w); + columns[limb + 1] = Operations::multiply_add_high(columns[limb + 1], digits, w); + } + } + + // The certificate reciprocal has up to 61 bits, so it cannot be an IFMA operand + // without splitting. Compute floor(d*r/2^50) exactly from four 52-bit split products. + static vector certificate_term(vector digits, std::uint64_t reciprocal) { + const auto zero = Operations::broadcast(0); + const auto low = Operations::broadcast(reciprocal & radix::mask); + const auto high = Operations::broadcast(reciprocal >> 52); + const auto a = Operations::multiply_add_low(zero, digits, low); + const auto b = Operations::multiply_add_high(zero, digits, low); + const auto c = Operations::multiply_add_low(zero, digits, high); + const auto d = Operations::multiply_add_high(zero, digits, high); + return Operations::add( + Operations::add(Operations::template shift_right<50>(a), Operations::template shift_left<2>(b)), + Operations::add(Operations::template shift_left<2>(c), Operations::template shift_left<54>(d))); + } + + struct certificate { + vector quotients; + vector accepted; + }; + + static certificate certify(vector lower, std::uint64_t scaled_bound, std::size_t prime_count) { + constexpr std::uint64_t scale = std::uint64_t(1) << 60; + const auto zero = Operations::broadcast(0); + const auto one = Operations::broadcast(1); + const auto quotient = + Operations::template shift_right<60>(Operations::add(lower, Operations::broadcast(scale / 2))); + const auto center = Operations::template shift_left<60>(quotient); + const auto bound = Operations::broadcast(scaled_bound); + const auto upper = Operations::add(lower, Operations::broadcast(2 * prime_count)); + const auto lower_ok = + Operations::select(Operations::less(Operations::add(lower, bound), center), zero, one); + const auto upper_ok = + Operations::select(Operations::less(Operations::add(center, bound), upper), zero, one); + return {quotient, Operations::bit_and(lower_ok, upper_ok)}; + } + + // Add the selected centering correction and divide by S=2^260. The caller prepares + // weights and corrections with R*S (R=2^256), so the result is already R-Montgomery. + std::array reduce(const std::array &correction) const { + const auto zero = Operations::broadcast(0); + const auto mask = Operations::broadcast(radix::mask); + const auto inverse = Operations::broadcast(radix::p_dash()); + std::array t; + t.fill(zero); + std::copy(columns.begin(), columns.end(), t.begin()); + for (std::size_t i = 0; i < correction.size(); ++i) + t[i] = Operations::add(t[i], correction[i]); + for (std::size_t i = 0; i < 9; ++i) { + t[i + 1] = Operations::add(t[i + 1], Operations::shift52(t[i])); + t[i] = Operations::bit_and(t[i], mask); + } + // Each cancellation adds at most five split products to a column. All columns + // stay below 2^57; the initial sum plus correction is < 2^54*p < S*p. + for (std::size_t i = 0; i < 5; ++i) { + const auto m = Operations::multiply_add_low(zero, t[i], inverse); + for (std::size_t j = 0; j < 5; ++j) { + const auto p = Operations::broadcast(radix::modulus()[j]); + t[i + j] = Operations::multiply_add_low(t[i + j], m, p); + t[i + j + 1] = Operations::multiply_add_high(t[i + j + 1], m, p); + } + t[i + 1] = Operations::add(t[i + 1], Operations::shift52(t[i])); + } + std::array reduced; + for (std::size_t i = 5; i < 9; ++i) { + reduced[i - 5] = Operations::bit_and(t[i], mask); + t[i + 1] = Operations::add(t[i + 1], Operations::shift52(t[i])); + } + reduced[4] = t[9]; + // REDC(T)<2*p::subtract_modulus(reduced); + return reduced; + } + }; + /** * Convert canonical BN254 Fp12 coordinates to a prefix of the small-prime catalogue. * Coordinate order is Crypto3's [1,u,v,uv,v^2,uv^2,w,uw,vw,uvw,v^2w,uv^2w]. @@ -239,6 +339,9 @@ namespace nil::crypto3::math::detail { static_assert(50 * small_prime_fft_primes.size() <= 600, "CRT export folding requires M < 2^600"); using fast_params = bn254_small_prime_fft_converter::field_type::extension_policy::fast_params; using direct_weight_type = std::array; +#if defined(CRYPTO3_SMALL_PRIME_CRT_IFMA_SELECTED) + using radix52 = algebra::fields::detail::fp12_fast::montgomery_radix52; +#endif static constexpr bool direct_export_supported = base_field_type::integral_type::backend_type::limb_bits == 64; static constexpr unsigned certificate_bits = 60; static constexpr std::uint64_t certificate_scale = std::uint64_t(1) << certificate_bits; @@ -252,6 +355,9 @@ namespace nil::crypto3::math::detail { double inverse; std::uint64_t certificate_reciprocal; direct_weight_type export_weight; +#if defined(CRYPTO3_SMALL_PRIME_CRT_IFMA_SELECTED) + radix52::value_type export_weight52; +#endif prime_context(std::uint64_t modulus, const integer_type &cofactor, const direct_weight_type &weight) : arithmetic(modulus), cofactor(cofactor), export_weight(weight) { @@ -306,6 +412,13 @@ namespace nil::crypto3::math::detail { direct_weight_type limbs {}; boost::multiprecision::export_bits(weight, limbs.begin(), 64, false); primes_.emplace_back(modulus, multiples_[0] / modulus, limbs); +#if defined(CRYPTO3_SMALL_PRIME_CRT_IFMA_SELECTED) + // S/R=16. Export after radix-52 REDC must retain the field's R=2^256 scaling. + const boost::multiprecision::cpp_int vector_weight = (weight * 16) % characteristic; + limbs.fill(0); + boost::multiprecision::export_bits(vector_weight, limbs.begin(), 64, false); + primes_.back().export_weight52 = radix52::pack(limbs.data()); +#endif } // The centering quotient is in [0,prime_count]. Store canonical -k*M*R^2 mod p. for (std::size_t k = 0; k <= params.prime_count; ++k) { @@ -314,6 +427,14 @@ namespace nil::crypto3::math::detail { correction = characteristic - correction; } boost::multiprecision::export_bits(correction, export_corrections_[k].begin(), 64, false); +#if defined(CRYPTO3_SMALL_PRIME_CRT_IFMA_SELECTED) + const boost::multiprecision::cpp_int vector_correction = (correction * 16) % characteristic; + direct_weight_type limbs {}; + boost::multiprecision::export_bits(vector_correction, limbs.begin(), 64, false); + const auto packed = radix52::pack(limbs.data()); + for (std::size_t digit = 0; digit < packed.size(); ++digit) + export_corrections52_[digit][k] = packed[digit]; +#endif } } @@ -382,6 +503,10 @@ namespace nil::crypto3::math::detail { return prime.arithmetic.fma_supported(); }); if (use_fma) { +#if defined(CRYPTO3_SMALL_PRIME_CRT_IFMA_SELECTED) + reconstruct_ifma(output, residues, stride); + return; +#endif reconstruct_impl(output, residues, stride); } else { reconstruct_impl(output, residues, stride); @@ -574,6 +699,77 @@ namespace nil::crypto3::math::detail { } } +#if defined(CRYPTO3_SMALL_PRIME_CRT_IFMA_SELECTED) + template + void reconstruct_ifma(std::span output, std::span residues, + std::size_t stride) const { + using operations = algebra::fields::detail::fp12_fast::ifma52_operations<8>; + for (std::size_t coordinate = 0; coordinate < field_type::arity; ++coordinate) { + for (std::size_t offset = 0; offset < output.size();) { + const std::size_t count = std::min(operations::lanes, output.size() - offset); + const auto active = static_cast<__mmask8>((1u << count) - 1); + using accumulator_type = small_prime_fft_crt_accumulator; + accumulator_type accumulator; + auto lower = operations::broadcast(0); + std::array, small_prime_fft_primes.size()> + saved_digits; + for (std::size_t prime = 0; prime < primes_.size(); ++prime) { + const auto &context = primes_[prime]; + const small_prime_fft_arithmetic::wide_fma_kernel packed(context.arithmetic); + const double *source = + residues.data() + (prime * field_type::arity + coordinate) * stride + offset; + auto digits = packed.reduce_fma(_mm512_maskz_loadu_pd(active, source)); + if constexpr (!Weighted) { + digits = packed.multiply_fma(digits, _mm512_set1_pd(context.inverse)); + } + const auto negative = _mm512_cmp_pd_mask(digits, _mm512_setzero_pd(), _CMP_LT_OQ); + digits = _mm512_mask_add_pd(digits, negative, digits, + _mm512_set1_pd(static_cast(context.arithmetic.modulus()))); + const operations::value_type canonical {_mm512_cvttpd_epu64(digits)}; + operations::store(saved_digits[prime].data(), canonical); + accumulator.add(context.export_weight52, canonical); + lower = operations::add( + lower, accumulator_type::certificate_term(canonical, context.certificate_reciprocal)); + } + const auto certificate = accumulator_type::certify(lower, scaled_bound_, primes_.size()); + const auto accepted = operations::less(operations::broadcast(0), certificate.accepted) & active; + std::array, 5> reduced; + if (accepted != 0) { + std::array corrections; + for (std::size_t digit = 0; digit < corrections.size(); ++digit) { + // Every quotient is in [0,prime_count], including uncertified lanes. + corrections[digit] = {_mm512_i64gather_epi64(certificate.quotients.data, + export_corrections52_[digit].data(), 8)}; + } + const auto values = accumulator.reduce(corrections); + for (std::size_t digit = 0; digit < values.size(); ++digit) + operations::store(reduced[digit].data(), values[digit]); + } + for (std::size_t lane = 0; lane < count; ++lane) { + auto &destination = output[offset + lane].coordinate(coordinate); + if ((accepted & (1u << lane)) != 0) { + // Numeric repacking only: the R*S weights already account for S/R. + base_field_type::value_type value; + auto *words = value.data.backend().base_data().limbs(); + words[0] = reduced[0][lane] | (reduced[1][lane] << 52); + words[1] = (reduced[1][lane] >> 12) | (reduced[2][lane] << 40); + words[2] = (reduced[2][lane] >> 24) | (reduced[3][lane] << 28); + words[3] = (reduced[3][lane] >> 36) | (reduced[4][lane] << 16); + destination = value; + } else { + integer_type exact_sum = 0; + for (std::size_t prime = 0; prime < primes_.size(); ++prime) { + add_cofactor_product(exact_sum, primes_[prime].cofactor, saved_digits[prime][lane]); + } + destination = export_exact(exact_sum); + } + } + offset += count; + } + } + } +#endif + std::vector primes_; std::array multiples_; integer_type half_modulus_; @@ -582,6 +778,9 @@ namespace nil::crypto3::math::detail { export_type export_radix_; std::uint64_t scaled_bound_; std::array export_corrections_ {}; +#if defined(CRYPTO3_SMALL_PRIME_CRT_IFMA_SELECTED) + std::array, 5> export_corrections52_ {}; +#endif }; /** diff --git a/libs/math/test/small_prime_fft_crt.cpp b/libs/math/test/small_prime_fft_crt.cpp index 6fb091d26..25d46a055 100644 --- a/libs/math/test/small_prime_fft_crt.cpp +++ b/libs/math/test/small_prime_fft_crt.cpp @@ -49,6 +49,92 @@ namespace { using value_type = reconstructor_type::value_type; using coordinates_type = std::array; + template + struct scalar_ifma52_operations { + using limb = nil::crypto3::algebra::fields::detail::fp12_fast::limb; + using wide = nil::crypto3::algebra::fields::detail::fp12_fast::wide_limb; + static constexpr std::size_t lanes = Lanes; + static constexpr limb mask = (limb(1) << 52) - 1; + using value_type = std::array; + using mask_type = std::array; + static value_type broadcast(limb x) { + value_type r; + r.fill(x); + return r; + } + static value_type load(const limb *x) { + value_type r; + std::copy_n(x, lanes, r.begin()); + return r; + } + static void store(limb *x, const value_type &a) { + std::copy(a.begin(), a.end(), x); + } + static value_type add(value_type a, const value_type &b) { + for (std::size_t i = 0; i < lanes; ++i) { + if (a[i] > std::numeric_limits::max() - b[i]) { + throw std::overflow_error("IFMA lane overflow"); + } + a[i] += b[i]; + } + return a; + } + static value_type subtract(value_type a, const value_type &b) { + for (std::size_t i = 0; i < lanes; ++i) + a[i] -= b[i]; + return a; + } + static value_type bit_and(value_type a, const value_type &b) { + for (std::size_t i = 0; i < lanes; ++i) + a[i] &= b[i]; + return a; + } + static value_type shift52(value_type a) { + for (auto &x : a) + x >>= 52; + return a; + } + template + static value_type shift_right(value_type a) { + for (auto &x : a) + x >>= Bits; + return a; + } + template + static value_type shift_left(value_type a) { + for (auto &x : a) { + if (x > (std::numeric_limits::max() >> Bits)) + throw std::overflow_error("IFMA shift overflow"); + x <<= Bits; + } + return a; + } + static mask_type less(const value_type &a, const value_type &b) { + mask_type r; + for (std::size_t i = 0; i < lanes; ++i) + r[i] = a[i] < b[i]; + return r; + } + static value_type select(const mask_type &mask, value_type yes, const value_type &no) { + for (std::size_t i = 0; i < lanes; ++i) + if (!mask[i]) + yes[i] = no[i]; + return yes; + } + static value_type multiply_add_low(value_type sum, const value_type &a, const value_type &b) { + value_type product; + for (std::size_t i = 0; i < lanes; ++i) + product[i] = limb(wide(a[i] & mask) * (b[i] & mask)) & mask; + return add(sum, product); + } + static value_type multiply_add_high(value_type sum, const value_type &a, const value_type &b) { + value_type product; + for (std::size_t i = 0; i < lanes; ++i) + product[i] = limb((wide(a[i] & mask) * (b[i] & mask)) >> 52); + return add(sum, product); + } + }; + struct restore_rounding { int saved = std::fegetround(); ~restore_rounding() { @@ -155,10 +241,133 @@ namespace { } return result; } + template + void check_crt_accumulation() { + using params = reconstructor_type::field_type::extension_policy::fast_params; + using accumulator_type = detail::small_prime_fft_crt_accumulator; + using radix = typename accumulator_type::radix; + const cpp_int p(fq::modulus.backend().to_cpp_int()); + std::mt19937_64 random(0x52c47); + for (std::size_t count = 1; count <= detail::small_prime_fft_primes.size(); ++count) { + accumulator_type accumulator; + std::array expected {}; + for (std::size_t prime = 0; prime < count; ++prime) { + const auto q = detail::small_prime_fft_primes[prime].modulus; + cpp_int weight = p - 1; + if (prime % 3 != 0) { + weight = 0; + for (unsigned i = 0; i < 4; ++i) + weight = (weight << 64) + random(); + weight %= p; + } + std::array words {}; + boost::multiprecision::export_bits(weight, words.begin(), 64, false); + const auto packed = radix::pack(words.data()); + std::array digits; + for (std::size_t lane = 0; lane < Operations::lanes; ++lane) { + digits[lane] = lane == 0 ? 0 : lane == 1 ? q - 1 : random() % q; + expected[lane] += weight * digits[lane]; + } + accumulator.add(packed, Operations::load(digits.data())); + } + std::array, 6> columns; + for (std::size_t digit = 0; digit < columns.size(); ++digit) + Operations::store(columns[digit].data(), accumulator.columns[digit]); + std::array correction; + std::array, 5> correction_digits; + for (std::size_t lane = 0; lane < Operations::lanes; ++lane) { + cpp_int actual = 0; + for (std::size_t i = columns.size(); i-- > 0;) + actual = (actual << 52) + columns[i][lane]; + BOOST_CHECK_EQUAL(actual, expected[lane]); + correction[lane] = lane % 2 ? p - 1 : cpp_int(lane); + cpp_int c = correction[lane]; + for (std::size_t i = 0; i < 5; ++i) { + correction_digits[i][lane] = (c & radix::mask).template convert_to(); + c >>= 52; + } + } + std::array vector_correction; + for (std::size_t i = 0; i < 5; ++i) + vector_correction[i] = Operations::load(correction_digits[i].data()); + const auto reduced = accumulator.reduce(vector_correction); + std::array, 5> digits; + for (std::size_t i = 0; i < 5; ++i) + Operations::store(digits[i].data(), reduced[i]); + const cpp_int s_inverse = boost::multiprecision::powm(cpp_int((cpp_int(1) << 260) % p), cpp_int(p - 2), p); + for (std::size_t lane = 0; lane < Operations::lanes; ++lane) { + cpp_int actual = 0; + for (std::size_t i = 5; i-- > 0;) { + BOOST_CHECK_LE(digits[i][lane], radix::mask); + actual = (actual << 52) + digits[i][lane]; + } + BOOST_CHECK_LT(actual, p); + BOOST_CHECK_EQUAL(actual, ((expected[lane] + correction[lane]) * s_inverse) % p); + } + } + } + template + void check_crt_certificates() { + using params = reconstructor_type::field_type::extension_policy::fast_params; + using accumulator_type = detail::small_prime_fft_crt_accumulator; + constexpr std::uint64_t scale = std::uint64_t(1) << 60; + std::mt19937_64 random(0xc471f1); + for (std::size_t count = 1; count <= detail::small_prime_fft_primes.size(); ++count) { + const auto plan = make_params(count); + for (unsigned trial = 0; trial < 32; ++trial) { + auto lower = Operations::broadcast(0); + std::array sums {}, expected_lower {}; + for (std::size_t prime = 0; prime < count; ++prime) { + const auto q = detail::small_prime_fft_primes[prime].modulus; + const auto reciprocal = ((cpp_int(1) << 110) / q).convert_to(); + std::array input; + for (std::size_t lane = 0; lane < Operations::lanes; ++lane) { + input[lane] = trial == 0 ? 0 : trial == 1 ? q - 1 : random() % q; + sums[lane] += cpp_int(input[lane]) * (plan.modulus / q); + expected_lower[lane] += (cpp_int(input[lane]) * reciprocal) >> 50; + } + lower = Operations::add( + lower, accumulator_type::certificate_term(Operations::load(input.data()), reciprocal)); + } + std::array actual_lower; + Operations::store(actual_lower.data(), lower); + for (std::size_t lane = 0; lane < Operations::lanes; ++lane) + BOOST_CHECK_EQUAL(actual_lower[lane], expected_lower[lane]); + for (std::uint64_t bound : {std::uint64_t(0), std::uint64_t(1), scale / 8, scale / 2 - 1}) { + const auto certificate = accumulator_type::certify(lower, bound, count); + std::array quotients, accepted; + Operations::store(quotients.data(), certificate.quotients); + Operations::store(accepted.data(), certificate.accepted); + for (std::size_t lane = 0; lane < Operations::lanes; ++lane) { + const cpp_int quotient = (expected_lower[lane] + scale / 2) / scale; + const cpp_int center = quotient * scale; + const bool expected = expected_lower[lane] + bound >= center && + expected_lower[lane] + 2 * count <= center + bound; + BOOST_CHECK_EQUAL(quotients[lane], quotient); + BOOST_CHECK_EQUAL(bool(accepted[lane]), expected); + if (accepted[lane]) { + const cpp_int difference = sums[lane] - quotient * plan.modulus; + BOOST_CHECK_LE(boost::multiprecision::abs(difference) * scale, + cpp_int(bound) * plan.modulus); + } + } + } + } + } + } } // namespace BOOST_AUTO_TEST_SUITE(small_prime_fft_crt_test_suite) +BOOST_AUTO_TEST_CASE(radix52_accumulation_matches_full_integer_sums) { + check_crt_accumulation>(); + check_crt_certificates>(); +#if defined(CRYPTO3_SMALL_PRIME_CRT_IFMA_SELECTED) + check_crt_accumulation>(); + check_crt_certificates>(); +#endif +} + BOOST_AUTO_TEST_CASE(combined_normalization_is_centered_and_exact_in_every_rounding_mode) { restore_rounding restore; for (int mode : {FE_DOWNWARD, FE_UPWARD, FE_TOWARDZERO, FE_TONEAREST}) { @@ -301,7 +510,7 @@ BOOST_AUTO_TEST_CASE(bound_failure_preserves_completed_prefixes_across_block_bou const reconstructor_type reconstructor(plan->crt, bound); for (int mode : {FE_TONEAREST, FE_UPWARD}) { BOOST_REQUIRE_EQUAL(std::fesetround(mode), 0); - for (std::size_t bad : {0, 1, 31, 32, 33}) { + for (std::size_t bad : {0, 1, 7, 8, 9, 15, 16, 17, 31, 32, 33}) { for (int sign : {-1, 1}) { std::vector coordinates(35); coordinates[bad][0] = sign * (bound + 1);