|
| 1 | +#ifndef CP_ALGO_MATH_AND_CONVOLUTION_HPP |
| 2 | +#define CP_ALGO_MATH_AND_CONVOLUTION_HPP |
| 3 | +#include "../number_theory/modint.hpp" |
| 4 | +#include "../util/bit.hpp" |
| 5 | +#include "../util/checkpoint.hpp" |
| 6 | + |
| 7 | +namespace cp_algo::math { |
| 8 | + enum transform_dir { forward, inverse }; |
| 9 | + // Recursive Zeta / Möbius (AND) transform for size N (power of two) |
| 10 | + template<auto N, transform_dir direction> |
| 11 | + void and_transform(auto &&a) { |
| 12 | + if constexpr (N == 1) { |
| 13 | + return; |
| 14 | + } else { |
| 15 | + constexpr auto half = N / 2; |
| 16 | + and_transform<half, direction>(&a[0]); |
| 17 | + and_transform<half, direction>(&a[half]); |
| 18 | + for (uint32_t i = 0; i < half; i++) { |
| 19 | + if constexpr (direction == forward) { |
| 20 | + a[i] += a[i + half]; |
| 21 | + } else { |
| 22 | + a[i] -= a[i + half]; |
| 23 | + } |
| 24 | + } |
| 25 | + } |
| 26 | + } |
| 27 | + // Wrapper that deduces N at compile time via with_bit_floor |
| 28 | + template<transform_dir direction> |
| 29 | + inline void and_transform(auto &&a, auto n) { |
| 30 | + with_bit_floor(n, [&]<auto NN>() { |
| 31 | + assert(NN == n); |
| 32 | + and_transform<NN, direction>(a); |
| 33 | + }); |
| 34 | + } |
| 35 | + template<transform_dir direction = forward> |
| 36 | + inline void and_transform(auto &&a) { |
| 37 | + and_transform<direction>(a, std::size(a)); |
| 38 | + } |
| 39 | + // In-place AND convolution on sequences of equal length (power of two) |
| 40 | + void and_convolution_inplace(auto &a, auto &b) { |
| 41 | + auto N = static_cast<uint32_t>(std::size(a)); |
| 42 | + and_transform(a); |
| 43 | + and_transform(b); |
| 44 | + checkpoint("transform"); |
| 45 | + for (uint32_t i = 0; i < N; i++) { |
| 46 | + a[i] *= b[i]; |
| 47 | + } |
| 48 | + checkpoint("dot"); |
| 49 | + and_transform<inverse>(a); |
| 50 | + checkpoint("itransform"); |
| 51 | + } |
| 52 | + |
| 53 | + // Returns AND convolution of a and b; pads to next power of two |
| 54 | + auto and_convolution(auto a, auto b) { |
| 55 | + auto n = std::bit_ceil(std::max(std::size(a), std::size(b))); |
| 56 | + a.resize(n); |
| 57 | + b.resize(n); |
| 58 | + and_convolution_inplace(a, b); |
| 59 | + return a; |
| 60 | + } |
| 61 | +} |
| 62 | +#endif // CP_ALGO_MATH_AND_CONVOLUTION_HPP |
0 commit comments