String Algorithms Bundle
(string/all.hpp)
- View this file on GitHub
- Last update: 2026-08-11 13:50:43+09:00
- Include:
#include "string/all.hpp"
Overview
string/all.hpp includes the reusable string algorithms in this repository.
Use individual headers when compile time matters, or this bundle during a
contest when convenience matters more.
Included Headers
| Header | Contents |
|---|---|
string/aho_corasick.hpp |
Multi-pattern matching with failure links and occurrence counting. |
string/deque_eertree.hpp |
Double-ended palindromic tree with dynamic distinct and longest-end queries. |
string/eertree.hpp |
Online palindromic tree with suffix and series links. |
string/levenshtein_distance.hpp |
Unit-cost edit distance in linear auxiliary memory. |
string/kmp.hpp |
Prefix function and linear-time KMP occurrence search. |
string/longest_common_extension.hpp |
Static longest-common-extension queries and substring comparisons. |
string/longest_common_subsequence.hpp |
Finds one longest subsequence common to two sequences. |
string/longest_common_substring.hpp |
Finds one longest substring common to two sequences. |
string/lyndon_factorization.hpp |
Duval’s linear-time Lyndon factorization. |
string/map_trie.hpp |
Ordered-map multiset trie for large or generic alphabets. |
string/z_algorithm.hpp |
Linear-time Z array. |
string/manacher.hpp |
Odd/even palindrome radii and substring checks. |
string/minimum_rotation.hpp |
Earliest lexicographically minimum cyclic shift in linear time. |
string/palindrome_lexicographical_order.hpp |
Rank and select distinct palindromic substrings in lexicographic order. |
string/prefix_substring_lcs.hpp |
Offline LCS-length queries between prefixes and substrings. |
string/suffix_automaton.hpp |
Online suffix automaton for substring queries and occurrence classes. |
string/suffix_array.hpp |
Suffix array and LCP array. |
string/suffix_tree.hpp |
Ukkonen suffix tree with substring lookup and occurrence counts. |
string/trie.hpp |
Contiguous-alphabet multiset trie with prefix queries. |
string/wildcard_pattern_matching.hpp |
Exact wildcard matching at every text alignment. |
string/rolling_hash.hpp |
Static substring hashing, LCP, and comparison. |
string/runs.hpp |
Enumerates maximal periodic substrings and their minimum periods. |
string/string_hash.hpp |
Double whole-string hashing and constant-time hash concatenation. |
Depends on
- Convolution
(math/fps/convolution.hpp)
- math/fps/internal/ntt998_faster.hpp
- ModInt
(math/modint.hpp)
- Aho-Corasick
(string/aho_corasick.hpp)
- Deque Eertree
(string/deque_eertree.hpp)
- Eertree
(string/eertree.hpp)
- Knuth-Morris-Pratt
(string/kmp.hpp)
- Levenshtein Distance
(string/levenshtein_distance.hpp)
- Longest Common Extension
(string/longest_common_extension.hpp)
- Longest Common Subsequence
(string/longest_common_subsequence.hpp)
- Longest Common Substring
(string/longest_common_substring.hpp)
- Lyndon Factorization
(string/lyndon_factorization.hpp)
- Manacher Algorithm
(string/manacher.hpp)
- Map Trie
(string/map_trie.hpp)
- Minimum Rotation
(string/minimum_rotation.hpp)
- Palindrome Lexicographical Order
(string/palindrome_lexicographical_order.hpp)
- Prefix-Substring LCS
(string/prefix_substring_lcs.hpp)
- Static Rolling Hash
(string/rolling_hash.hpp)
- Runs
(string/runs.hpp)
- String Hash
(string/string_hash.hpp)
- Suffix Array and LCP Array
(string/suffix_array.hpp)
- Suffix Automaton
(string/suffix_automaton.hpp)
- Suffix Tree
(string/suffix_tree.hpp)
- Trie
(string/trie.hpp)
- Wildcard Pattern Matching
(string/wildcard_pattern_matching.hpp)
- Z Algorithm
(string/z_algorithm.hpp)
Verified with
Code
#ifndef M1UNE_STRING_ALL_HPP
#define M1UNE_STRING_ALL_HPP 1
#include "aho_corasick.hpp"
#include "deque_eertree.hpp"
#include "eertree.hpp"
#include "kmp.hpp"
#include "levenshtein_distance.hpp"
#include "longest_common_extension.hpp"
#include "longest_common_subsequence.hpp"
#include "longest_common_substring.hpp"
#include "lyndon_factorization.hpp"
#include "manacher.hpp"
#include "map_trie.hpp"
#include "minimum_rotation.hpp"
#include "palindrome_lexicographical_order.hpp"
#include "prefix_substring_lcs.hpp"
#include "rolling_hash.hpp"
#include "runs.hpp"
#include "string_hash.hpp"
#include "suffix_automaton.hpp"
#include "suffix_array.hpp"
#include "suffix_tree.hpp"
#include "trie.hpp"
#include "wildcard_pattern_matching.hpp"
#include "z_algorithm.hpp"
#endif // M1UNE_STRING_ALL_HPP#line 1 "string/all.hpp"
#line 1 "string/aho_corasick.hpp"
#include <array>
#include <cassert>
#include <cstddef>
#include <limits>
#include <queue>
#include <vector>
namespace m1une {
namespace string {
// Aho-Corasick automaton for a contiguous character alphabet.
template <int AlphabetSize = 26, int FirstCharacter = 'a'>
struct AhoCorasick {
static_assert(0 < AlphabetSize);
using node_id = int;
static constexpr node_id null_node = -1;
struct Node {
// Completed automaton transitions. Valid after build().
std::array<node_id, AlphabetSize> next;
node_id failure;
node_id output_link;
node_id parent;
int parent_symbol;
int depth;
std::vector<node_id> children;
std::vector<node_id> failure_children;
std::vector<int> pattern_ids;
Node(
node_id parent_value = null_node,
int parent_symbol_value = -1,
int depth_value = 0
) : failure(0),
output_link(null_node),
parent(parent_value),
parent_symbol(parent_symbol_value),
depth(depth_value) {
next.fill(null_node);
}
};
private:
std::vector<Node> _nodes;
std::vector<int> _pattern_length;
std::vector<node_id> _pattern_node;
std::vector<node_id> _bfs_order;
bool _built;
template <class Symbol>
static int symbol_index(const Symbol& symbol) {
int index = int(symbol) - FirstCharacter;
assert(0 <= index && index < AlphabetSize);
return index;
}
node_id new_node(node_id parent, int parent_symbol) {
assert(_nodes.size() < std::size_t(std::numeric_limits<int>::max()));
assert(_nodes[parent].depth < std::numeric_limits<int>::max());
_nodes.emplace_back(parent, parent_symbol, _nodes[parent].depth + 1);
return int(_nodes.size()) - 1;
}
public:
AhoCorasick() : _nodes(1), _built(false) {}
node_id root() const {
return 0;
}
bool built() const {
return _built;
}
int pattern_count() const {
return int(_pattern_length.size());
}
int pattern_length(int pattern_id) const {
assert(0 <= pattern_id && pattern_id < pattern_count());
return _pattern_length[pattern_id];
}
node_id pattern_node(int pattern_id) const {
assert(0 <= pattern_id && pattern_id < pattern_count());
return _pattern_node[pattern_id];
}
std::size_t node_count() const {
return _nodes.size();
}
const std::vector<Node>& nodes() const {
return _nodes;
}
const Node& node(node_id id) const {
assert(0 <= id && std::size_t(id) < _nodes.size());
return _nodes[id];
}
// Returns nodes in failure-link BFS order, beginning with the root.
const std::vector<node_id>& bfs_order() const {
assert(_built);
return _bfs_order;
}
void reserve(std::size_t node_capacity) {
assert(!_built);
_nodes.reserve(node_capacity);
}
void clear() {
_nodes.clear();
_nodes.emplace_back();
_pattern_length.clear();
_pattern_node.clear();
_bfs_order.clear();
_built = false;
}
// Inserts a pattern and returns its insertion-order ID.
template <class Sequence>
int insert(const Sequence& pattern) {
assert(!_built);
int pattern_id = pattern_count();
int length = 0;
node_id state = root();
for (const auto& symbol : pattern) {
assert(length < std::numeric_limits<int>::max());
int index = symbol_index(symbol);
if (_nodes[state].next[index] == null_node) {
node_id child = new_node(state, index);
_nodes[state].next[index] = child;
_nodes[state].children.push_back(child);
}
state = _nodes[state].next[index];
length++;
}
_nodes[state].pattern_ids.push_back(pattern_id);
_pattern_length.push_back(length);
_pattern_node.push_back(state);
return pattern_id;
}
// Builds failure links and completes every automaton transition.
void build() {
assert(!_built);
std::queue<node_id> queue;
_bfs_order.clear();
_bfs_order.reserve(_nodes.size());
_bfs_order.push_back(root());
for (int symbol = 0; symbol < AlphabetSize; ++symbol) {
node_id child = _nodes[root()].next[symbol];
if (child == null_node) {
_nodes[root()].next[symbol] = root();
} else {
_nodes[root()].next[symbol] = child;
_nodes[child].failure = root();
_nodes[child].output_link =
_nodes[root()].pattern_ids.empty() ? null_node : root();
_nodes[root()].failure_children.push_back(child);
queue.push(child);
}
}
while (!queue.empty()) {
node_id state = queue.front();
queue.pop();
_bfs_order.push_back(state);
for (int symbol = 0; symbol < AlphabetSize; ++symbol) {
node_id child = _nodes[state].next[symbol];
if (child == null_node) {
_nodes[state].next[symbol] =
_nodes[_nodes[state].failure].next[symbol];
continue;
}
_nodes[state].next[symbol] = child;
node_id failure =
_nodes[_nodes[state].failure].next[symbol];
_nodes[child].failure = failure;
_nodes[child].output_link =
_nodes[failure].pattern_ids.empty()
? _nodes[failure].output_link
: failure;
_nodes[failure].failure_children.push_back(child);
queue.push(child);
}
}
_built = true;
}
template <class Symbol>
node_id transition(node_id state, const Symbol& symbol) const {
assert(_built);
assert(0 <= state && std::size_t(state) < _nodes.size());
return _nodes[state].next[symbol_index(symbol)];
}
// Calls callback(pattern_id) for every pattern ending at `state`.
template <class Callback>
void for_each_output(node_id state, Callback callback) const {
assert(_built);
assert(0 <= state && std::size_t(state) < _nodes.size());
while (state != null_node) {
for (int pattern_id : _nodes[state].pattern_ids) {
callback(pattern_id);
}
state = _nodes[state].output_link;
}
}
// Calls callback(end, pattern_id) for every occurrence. `end` is the
// exclusive end position. Empty patterns occur at every text boundary.
template <class Sequence, class Callback>
void match(const Sequence& text, Callback callback) const {
assert(_built);
node_id state = root();
for_each_output(state, [&callback](int pattern_id) {
callback(0, pattern_id);
});
int end = 0;
for (const auto& symbol : text) {
state = transition(state, symbol);
end++;
for_each_output(state, [&callback, end](int pattern_id) {
callback(end, pattern_id);
});
}
}
// Counts occurrences of every inserted pattern in linear time.
template <class Sequence>
std::vector<long long> count_occurrences(const Sequence& text) const {
assert(_built);
std::vector<long long> visits(_nodes.size(), 0);
node_id state = root();
visits[root()]++;
for (const auto& symbol : text) {
state = transition(state, symbol);
visits[state]++;
}
for (std::size_t index = _bfs_order.size(); index-- > 1;) {
node_id current = _bfs_order[index];
visits[_nodes[current].failure] += visits[current];
}
std::vector<long long> result(pattern_count(), 0);
for (node_id current : _bfs_order) {
for (int pattern_id : _nodes[current].pattern_ids) {
result[pattern_id] = visits[current];
}
}
return result;
}
};
} // namespace string
} // namespace m1une
#line 1 "string/deque_eertree.hpp"
#line 7 "string/deque_eertree.hpp"
#include <deque>
#line 10 "string/deque_eertree.hpp"
namespace m1une {
namespace string {
template <int AlphabetSize = 26, int FirstCharacter = 'a'>
struct DequeEertree {
static_assert(0 < AlphabetSize);
using node_id = int;
static constexpr node_id odd_root = 0;
static constexpr node_id even_root = 1;
static constexpr node_id null_node = -1;
private:
struct Node {
std::array<node_id, AlphabetSize> next;
node_id parent;
node_id suffix_link;
node_id quick_link;
int length;
int surface_count;
int suffix_link_children;
bool active;
Node(
int length_value = 0,
node_id parent_value = null_node,
node_id suffix_link_value = null_node,
node_id quick_link_value = null_node
)
: parent(parent_value),
suffix_link(suffix_link_value),
quick_link(quick_link_value),
length(length_value),
surface_count(0),
suffix_link_children(0),
active(true) {
next.fill(null_node);
}
};
struct Position {
int symbol;
node_id prefix_surface;
node_id suffix_surface;
};
std::vector<Node> _nodes;
std::deque<Position> _text;
int _distinct_palindromes;
template <class Symbol>
static int symbol_index(const Symbol& value) {
int symbol = int(value) - FirstCharacter;
assert(0 <= symbol && symbol < AlphabetSize);
return symbol;
}
node_id new_node(node_id parent, node_id suffix_link, int length, int symbol) {
assert(_nodes.size() < std::size_t(std::numeric_limits<int>::max()));
node_id id = int(_nodes.size());
_nodes.emplace_back(length, parent, suffix_link, odd_root);
_nodes[parent].next[symbol] = id;
_nodes[suffix_link].suffix_link_children++;
_distinct_palindromes++;
return id;
}
void remove_node(node_id id, int symbol) {
Node& removed = _nodes[id];
assert(removed.active);
assert(removed.surface_count == 0);
assert(removed.suffix_link_children == 0);
assert(_nodes[removed.parent].next[symbol] == id);
_nodes[removed.parent].next[symbol] = null_node;
_nodes[removed.suffix_link].suffix_link_children--;
removed.active = false;
_distinct_palindromes--;
}
node_id back_appendable(int symbol, node_id node) const {
int n = int(_text.size());
while (true) {
int length = _nodes[node].length;
if (length == -1 || (length < n && _text[n - length - 1].symbol == symbol)) {
return node;
}
node_id suffix = _nodes[node].suffix_link;
int suffix_length = _nodes[suffix].length;
if (suffix_length == -1 || _text[n - suffix_length - 1].symbol == symbol) {
return suffix;
}
node = _nodes[node].quick_link;
}
}
node_id front_appendable(int symbol, node_id node) const {
int n = int(_text.size());
while (true) {
int length = _nodes[node].length;
if (length == -1 || (length < n && _text[length].symbol == symbol)) {
return node;
}
node_id suffix = _nodes[node].suffix_link;
int suffix_length = _nodes[suffix].length;
if (suffix_length == -1 || _text[suffix_length].symbol == symbol) {
return suffix;
}
node = _nodes[node].quick_link;
}
}
node_id prefix_node() const {
return _text.empty() ? even_root : _text.front().prefix_surface;
}
node_id suffix_node() const {
return _text.empty() ? even_root : _text.back().suffix_surface;
}
void initialize_roots() {
_nodes.clear();
_nodes.emplace_back(-1, odd_root, odd_root, odd_root);
_nodes.emplace_back(0, odd_root, odd_root, odd_root);
_distinct_palindromes = 0;
}
public:
DequeEertree() {
initialize_roots();
}
template <class Sequence>
explicit DequeEertree(const Sequence& sequence) {
initialize_roots();
build(sequence);
}
int size() const {
return _distinct_palindromes;
}
int text_length() const {
return int(_text.size());
}
bool empty() const {
return _text.empty();
}
int distinct_palindrome_count() const {
return _distinct_palindromes;
}
int longest_prefix_length() const {
return _nodes[prefix_node()].length;
}
int longest_suffix_length() const {
return _nodes[suffix_node()].length;
}
void reserve(std::size_t operation_capacity) {
_nodes.reserve(operation_capacity + 2);
}
void clear() {
_text.clear();
initialize_roots();
}
template <class Symbol>
void push_back(const Symbol& value) {
int symbol = symbol_index(value);
node_id parent = _text.empty() ? odd_root : back_appendable(symbol, suffix_node());
node_id palindrome = _nodes[parent].next[symbol];
node_id suffix = even_root;
if (palindrome == null_node) {
if (parent != odd_root) {
node_id suffix_parent = back_appendable(symbol, _nodes[parent].suffix_link);
suffix = _nodes[suffix_parent].next[symbol];
assert(suffix != null_node);
}
} else {
suffix = _nodes[palindrome].suffix_link;
}
_text.push_back(Position{symbol, even_root, even_root});
int n = int(_text.size());
if (palindrome == null_node) {
palindrome = new_node(parent, suffix, _nodes[parent].length + 2, symbol);
Node& created = _nodes[palindrome];
if (
_nodes[suffix].suffix_link != odd_root &&
_text[n - _nodes[suffix].length - 1].symbol ==
_text[n - _nodes[_nodes[suffix].suffix_link].length - 1].symbol
) {
created.quick_link = _nodes[suffix].quick_link;
} else {
created.quick_link = _nodes[suffix].suffix_link;
}
}
int left = n - _nodes[palindrome].length;
_text.back().suffix_surface = palindrome;
_text[left].prefix_surface = palindrome;
if (
_nodes[suffix].length >= 1 &&
_text[left + _nodes[suffix].length - 1].suffix_surface == suffix
) {
_text[left + _nodes[suffix].length - 1].suffix_surface = even_root;
}
_nodes[palindrome].surface_count++;
}
template <class Symbol>
void push_front(const Symbol& value) {
int symbol = symbol_index(value);
node_id parent = _text.empty() ? odd_root : front_appendable(symbol, prefix_node());
node_id palindrome = _nodes[parent].next[symbol];
node_id suffix = even_root;
if (palindrome == null_node) {
if (parent != odd_root) {
node_id suffix_parent = front_appendable(symbol, _nodes[parent].suffix_link);
suffix = _nodes[suffix_parent].next[symbol];
assert(suffix != null_node);
}
} else {
suffix = _nodes[palindrome].suffix_link;
}
_text.push_front(Position{symbol, even_root, even_root});
if (palindrome == null_node) {
palindrome = new_node(parent, suffix, _nodes[parent].length + 2, symbol);
Node& created = _nodes[palindrome];
if (
_nodes[suffix].suffix_link != odd_root &&
_text[_nodes[suffix].length].symbol ==
_text[_nodes[_nodes[suffix].suffix_link].length].symbol
) {
created.quick_link = _nodes[suffix].quick_link;
} else {
created.quick_link = _nodes[suffix].suffix_link;
}
}
_text.front().prefix_surface = palindrome;
_text[_nodes[palindrome].length - 1].suffix_surface = palindrome;
if (
_nodes[suffix].length >= 1 &&
_text[_nodes[palindrome].length - _nodes[suffix].length].prefix_surface == suffix
) {
_text[_nodes[palindrome].length - _nodes[suffix].length].prefix_surface = even_root;
}
_nodes[palindrome].surface_count++;
}
void pop_back() {
assert(!_text.empty());
node_id palindrome = suffix_node();
node_id suffix = _nodes[palindrome].suffix_link;
int left = text_length() - _nodes[palindrome].length;
int suffix_end = left + _nodes[suffix].length - 1;
if (
_nodes[palindrome].length >= 2 &&
_nodes[_text[suffix_end].suffix_surface].length < _nodes[suffix].length
) {
_text[suffix_end].suffix_surface = suffix;
_text[left].prefix_surface = suffix;
} else {
_text[left].prefix_surface = even_root;
}
_nodes[palindrome].surface_count--;
int symbol = _text.back().symbol;
if (
_nodes[palindrome].surface_count == 0 &&
_nodes[palindrome].suffix_link_children == 0
) {
remove_node(palindrome, symbol);
}
_text.pop_back();
}
void pop_front() {
assert(!_text.empty());
node_id palindrome = prefix_node();
node_id suffix = _nodes[palindrome].suffix_link;
int suffix_start = _nodes[palindrome].length - _nodes[suffix].length;
if (
_nodes[palindrome].length >= 2 &&
_nodes[_text[suffix_start].prefix_surface].length < _nodes[suffix].length
) {
_text[suffix_start].prefix_surface = suffix;
_text[_nodes[palindrome].length - 1].suffix_surface = suffix;
} else {
_text[_nodes[palindrome].length - 1].suffix_surface = even_root;
}
_nodes[palindrome].surface_count--;
int symbol = _text.front().symbol;
if (
_nodes[palindrome].surface_count == 0 &&
_nodes[palindrome].suffix_link_children == 0
) {
remove_node(palindrome, symbol);
}
_text.pop_front();
}
template <class Sequence>
void build(const Sequence& sequence) {
for (const auto& symbol : sequence) push_back(symbol);
}
};
template <int AlphabetSize = 26, int FirstCharacter = 'a'>
using DoubleEndedEertree = DequeEertree<AlphabetSize, FirstCharacter>;
template <int AlphabetSize = 26, int FirstCharacter = 'a'>
using DequePalindromicTree = DequeEertree<AlphabetSize, FirstCharacter>;
} // namespace string
} // namespace m1une
#line 1 "string/eertree.hpp"
#line 8 "string/eertree.hpp"
#include <utility>
#line 10 "string/eertree.hpp"
namespace m1une {
namespace string {
template <int AlphabetSize = 26, int FirstCharacter = 'a'>
struct Eertree {
static_assert(0 < AlphabetSize);
using node_id = int;
static constexpr node_id even_root = 0;
static constexpr node_id odd_root = 1;
static constexpr node_id null_node = -1;
struct Node {
std::array<node_id, AlphabetSize> next;
node_id suffix_link;
node_id series_link;
int length;
int diff;
int suffix_count;
int first_end;
long long suffix_occurrences;
Node(int length_value = 0, node_id suffix_link_value = even_root, node_id series_link_value = even_root)
: suffix_link(suffix_link_value),
series_link(series_link_value),
length(length_value),
diff(0),
suffix_count(0),
first_end(0),
suffix_occurrences(0) {
next.fill(null_node);
}
};
private:
std::vector<Node> _nodes;
std::vector<int> _text;
std::vector<node_id> _longest_suffix;
node_id _last;
template <class Symbol>
static int symbol_index(const Symbol& symbol) {
int index = int(symbol) - FirstCharacter;
assert(0 <= index && index < AlphabetSize);
return index;
}
node_id find_extendable(node_id node, int position, int symbol) const {
while (true) {
int length = _nodes[node].length;
int left = position - length - 1;
if (0 <= left && _text[left] == symbol) return node;
node = _nodes[node].suffix_link;
}
}
node_id new_node(int length) {
assert(_nodes.size() < std::size_t(std::numeric_limits<int>::max()));
_nodes.emplace_back(length);
return int(_nodes.size()) - 1;
}
public:
Eertree() {
clear();
}
template <class Sequence>
explicit Eertree(const Sequence& sequence) {
clear();
build(sequence);
}
int size() const {
return int(_nodes.size()) - 2;
}
bool empty() const {
return size() == 0;
}
int node_count() const {
return int(_nodes.size());
}
int text_length() const {
return int(_text.size());
}
node_id last() const {
return _last;
}
int longest_suffix_length() const {
return _nodes[_last].length;
}
const Node& node(node_id id) const {
assert(0 <= id && id < node_count());
return _nodes[id];
}
const std::vector<Node>& nodes() const {
return _nodes;
}
node_id longest_suffix_node(int prefix_length) const {
assert(1 <= prefix_length && prefix_length <= text_length());
return _longest_suffix[prefix_length - 1];
}
const std::vector<node_id>& longest_suffix_nodes() const {
return _longest_suffix;
}
template <class Callback>
void for_each_suffix(node_id id, Callback callback) const {
assert(0 <= id && id < node_count());
while (id >= 2) {
callback(id);
id = _nodes[id].suffix_link;
}
}
template <class Callback>
void for_each_suffix(Callback callback) const {
for_each_suffix(_last, callback);
}
void reserve(std::size_t text_capacity) {
_text.reserve(text_capacity);
_longest_suffix.reserve(text_capacity);
_nodes.reserve(text_capacity + 2);
}
void clear() {
_nodes.clear();
_nodes.emplace_back(0, odd_root, even_root);
_nodes.emplace_back(-1, odd_root, odd_root);
_text.clear();
_longest_suffix.clear();
_last = even_root;
}
template <class Symbol>
node_id add(const Symbol& value) {
int symbol = symbol_index(value);
int position = int(_text.size());
_text.push_back(symbol);
node_id current = find_extendable(_last, position, symbol);
node_id next = _nodes[current].next[symbol];
if (next == null_node) {
int length = _nodes[current].length + 2;
next = new_node(length);
_nodes[current].next[symbol] = next;
node_id suffix_link = even_root;
if (length != 1) {
node_id candidate = find_extendable(_nodes[current].suffix_link, position, symbol);
suffix_link = _nodes[candidate].next[symbol];
assert(suffix_link != null_node);
}
Node& created = _nodes[next];
created.suffix_link = suffix_link;
created.diff = created.length - _nodes[suffix_link].length;
created.series_link =
created.diff == _nodes[suffix_link].diff ? _nodes[suffix_link].series_link : suffix_link;
created.suffix_count = _nodes[suffix_link].suffix_count + 1;
created.first_end = position + 1;
}
_last = next;
_nodes[_last].suffix_occurrences++;
_longest_suffix.push_back(_last);
return _last;
}
template <class Sequence>
void build(const Sequence& sequence) {
for (const auto& symbol : sequence) add(symbol);
}
std::vector<long long> occurrence_counts() const {
std::vector<long long> result(_nodes.size(), 0);
for (node_id id = 0; id < node_count(); id++) {
result[id] = _nodes[id].suffix_occurrences;
}
for (node_id id = node_count() - 1; id >= 2; id--) {
result[_nodes[id].suffix_link] += result[id];
}
return result;
}
std::pair<int, int> first_occurrence(node_id id) const {
assert(2 <= id && id < node_count());
int end = _nodes[id].first_end;
return {end - _nodes[id].length, end};
}
};
template <int AlphabetSize = 26, int FirstCharacter = 'a'>
using PalindromicTree = Eertree<AlphabetSize, FirstCharacter>;
} // namespace string
} // namespace m1une
#line 1 "string/kmp.hpp"
#line 5 "string/kmp.hpp"
namespace m1une {
namespace string {
// Returns the KMP prefix function.
template <class Sequence>
std::vector<int> prefix_function(const Sequence& sequence) {
int n = int(sequence.size());
std::vector<int> prefix(n);
for (int i = 1; i < n; i++) {
int j = prefix[i - 1];
while (j > 0 && sequence[i] != sequence[j]) {
j = prefix[j - 1];
}
if (sequence[i] == sequence[j]) j++;
prefix[i] = j;
}
return prefix;
}
// Returns every starting position where pattern occurs in text.
// An empty pattern occurs at every position from 0 through text.size().
template <class Text, class Pattern>
std::vector<int> kmp_search(const Text& text, const Pattern& pattern) {
int n = int(text.size());
int m = int(pattern.size());
if (m == 0) {
std::vector<int> occurrences(n + 1);
for (int i = 0; i <= n; i++) occurrences[i] = i;
return occurrences;
}
std::vector<int> prefix = prefix_function(pattern);
std::vector<int> occurrences;
int matched = 0;
for (int i = 0; i < n; i++) {
while (matched > 0 && text[i] != pattern[matched]) {
matched = prefix[matched - 1];
}
if (text[i] == pattern[matched]) matched++;
if (matched == m) {
occurrences.push_back(i - m + 1);
matched = prefix[matched - 1];
}
}
return occurrences;
}
} // namespace string
} // namespace m1une
#line 1 "string/levenshtein_distance.hpp"
#include <algorithm>
#line 7 "string/levenshtein_distance.hpp"
namespace m1une {
namespace string {
namespace levenshtein_distance_detail {
template <class RowSequence, class ColumnSequence>
int solve(const RowSequence& rows, const ColumnSequence& columns) {
int row_count = int(rows.size());
int column_count = int(columns.size());
std::vector<int> distance(column_count + 1);
for (int column = 0; column <= column_count; column++) distance[column] = column;
for (int row = 1; row <= row_count; row++) {
int diagonal = distance[0];
distance[0] = row;
for (int column = 1; column <= column_count; column++) {
int above = distance[column];
int substitution = diagonal + (rows[row - 1] == columns[column - 1] ? 0 : 1);
distance[column] =
std::min({above + 1, distance[column - 1] + 1, substitution});
diagonal = above;
}
}
return distance[column_count];
}
template <class RowSequence, class ColumnSequence>
int solve_bounded(const RowSequence& rows, const ColumnSequence& columns,
int max_distance) {
int row_count = int(rows.size());
int column_count = int(columns.size());
assert(column_count <= row_count);
if (row_count - column_count > max_distance) return max_distance + 1;
if (max_distance >= row_count) return solve(rows, columns);
int infinity = max_distance + 1;
int previous_left = 0;
int previous_right = std::min(column_count, max_distance);
std::vector<int> previous(previous_right + 1);
for (int column = 0; column <= previous_right; column++) previous[column] = column;
std::vector<int> current;
for (int row = 1; row <= row_count; row++) {
int current_left = std::max(0, row - max_distance);
int current_right = int(std::min<long long>(column_count,
static_cast<long long>(row) + max_distance));
current.assign(current_right - current_left + 1, infinity);
for (int column = current_left; column <= current_right; column++) {
int best = infinity;
if (previous_left <= column && column <= previous_right) {
best = std::min(best, previous[column - previous_left] + 1);
}
if (current_left < column) {
best = std::min(best, current[column - current_left - 1] + 1);
}
if (0 < column && previous_left <= column - 1 && column - 1 <= previous_right) {
int substitution = previous[column - 1 - previous_left] +
(rows[row - 1] == columns[column - 1] ? 0 : 1);
best = std::min(best, substitution);
}
current[column - current_left] = std::min(best, infinity);
}
previous.swap(current);
previous_left = current_left;
previous_right = current_right;
}
return previous[column_count - previous_left];
}
} // namespace levenshtein_distance_detail
// Returns the minimum number of insertions, deletions, and substitutions
// needed to transform first into second.
template <class Sequence1, class Sequence2>
int levenshtein_distance(const Sequence1& first, const Sequence2& second) {
if (first.size() < second.size()) {
return levenshtein_distance_detail::solve(second, first);
}
return levenshtein_distance_detail::solve(first, second);
}
// Returns the exact distance when it is at most max_distance, and
// max_distance + 1 otherwise.
template <class Sequence1, class Sequence2>
int levenshtein_distance(const Sequence1& first, const Sequence2& second,
int max_distance) {
assert(0 <= max_distance);
if (first.size() < second.size()) {
return levenshtein_distance_detail::solve_bounded(second, first, max_distance);
}
return levenshtein_distance_detail::solve_bounded(first, second, max_distance);
}
} // namespace string
} // namespace m1une
#line 1 "string/longest_common_extension.hpp"
#line 6 "string/longest_common_extension.hpp"
#include <string>
#line 9 "string/longest_common_extension.hpp"
#line 1 "string/suffix_array.hpp"
#line 6 "string/suffix_array.hpp"
#include <numeric>
#line 8 "string/suffix_array.hpp"
#include <type_traits>
#line 10 "string/suffix_array.hpp"
namespace m1une {
namespace string {
namespace detail {
template <class Sequence>
std::vector<int> suffix_array_impl(const Sequence& sequence) {
int n = int(sequence.size());
if (n == 0) return {};
using Value = std::remove_cv_t<std::remove_reference_t<decltype(sequence[0])>>;
std::vector<Value> sorted(sequence.begin(), sequence.end());
std::sort(sorted.begin(), sorted.end());
sorted.erase(std::unique(sorted.begin(), sorted.end()), sorted.end());
int length = n + 1;
std::vector<int> order(length);
std::vector<int> rank(length);
std::vector<int> key(length);
key[n] = 0;
for (int i = 0; i < n; i++) {
key[i] = int(std::lower_bound(sorted.begin(), sorted.end(), sequence[i]) - sorted.begin()) + 1;
}
int alphabet = int(sorted.size()) + 1;
std::vector<int> count(std::max(length, alphabet), 0);
for (int value : key) count[value]++;
for (int i = 1; i < alphabet; i++) count[i] += count[i - 1];
for (int i = length - 1; i >= 0; i--) order[--count[key[i]]] = i;
int classes = 1;
rank[order[0]] = 0;
for (int i = 1; i < length; i++) {
if (key[order[i - 1]] != key[order[i]]) classes++;
rank[order[i]] = classes - 1;
}
std::vector<int> shifted(length);
std::vector<int> next_rank(length);
for (long long half = 1; half < length; half <<= 1) {
for (int i = 0; i < length; i++) {
long long position = order[i] - half;
if (position < 0) position += length;
shifted[i] = int(position);
}
count.assign(classes, 0);
for (int position : shifted) count[rank[position]]++;
for (int i = 1; i < classes; i++) count[i] += count[i - 1];
for (int i = length - 1; i >= 0; i--) {
int position = shifted[i];
order[--count[rank[position]]] = position;
}
int next_classes = 1;
next_rank[order[0]] = 0;
for (int i = 1; i < length; i++) {
int current = order[i];
int previous = order[i - 1];
int current_second = int((current + half) % length);
int previous_second = int((previous + half) % length);
if (
rank[current] != rank[previous] ||
rank[current_second] != rank[previous_second]
) {
next_classes++;
}
next_rank[current] = next_classes - 1;
}
rank.swap(next_rank);
classes = next_classes;
if (classes == length) break;
}
std::vector<int> suffixes(n);
for (int i = 0; i < n; i++) suffixes[i] = order[i + 1];
return suffixes;
}
} // namespace detail
template <class Sequence>
std::vector<int> suffix_array(const Sequence& sequence) {
return detail::suffix_array_impl(sequence);
}
inline std::vector<int> suffix_array(const std::string& text) {
std::vector<unsigned char> values;
values.reserve(text.size());
for (unsigned char character : text) values.push_back(character);
return detail::suffix_array_impl(values);
}
template <class Sequence>
std::vector<int> lcp_array(const Sequence& sequence, const std::vector<int>& suffixes) {
int n = int(sequence.size());
assert(int(suffixes.size()) == n);
if (n == 0) return {};
std::vector<int> rank(n);
for (int i = 0; i < n; i++) {
assert(0 <= suffixes[i] && suffixes[i] < n);
rank[suffixes[i]] = i;
}
std::vector<int> lcp(n - 1);
int common = 0;
for (int i = 0; i < n; i++) {
int position = rank[i];
if (position == n - 1) {
common = 0;
continue;
}
int j = suffixes[position + 1];
while (
i + common < n &&
j + common < n &&
sequence[i + common] == sequence[j + common]
) {
common++;
}
lcp[position] = common;
if (common > 0) common--;
}
return lcp;
}
} // namespace string
} // namespace m1une
#line 11 "string/longest_common_extension.hpp"
namespace m1une {
namespace string {
template <class Sequence = std::string>
struct LongestCommonExtension {
private:
Sequence _sequence;
std::vector<int> _suffix_array;
std::vector<int> _rank;
std::vector<int> _lcp;
std::vector<int> _log;
std::vector<std::vector<int>> _table;
int range_min(int left, int right) const {
assert(0 <= left && left < right && right <= int(_lcp.size()));
int k = _log[right - left];
return std::min(_table[k][left], _table[k][right - (1 << k)]);
}
void build() {
int n = int(_sequence.size());
_suffix_array = m1une::string::suffix_array(_sequence);
_rank.assign(n, 0);
for (int i = 0; i < n; i++) {
_rank[_suffix_array[i]] = i;
}
_lcp = m1une::string::lcp_array(_sequence, _suffix_array);
int m = int(_lcp.size());
_log.assign(m + 1, 0);
for (int i = 2; i <= m; i++) {
_log[i] = _log[i >> 1] + 1;
}
_table.clear();
if (m == 0) return;
_table.assign(_log[m] + 1, std::vector<int>());
_table[0] = _lcp;
for (int k = 1; k < int(_table.size()); k++) {
int width = 1 << k;
int half = width >> 1;
_table[k].resize(m - width + 1);
for (int i = 0; i + width <= m; i++) {
_table[k][i] = std::min(_table[k - 1][i], _table[k - 1][i + half]);
}
}
}
public:
LongestCommonExtension() = default;
explicit LongestCommonExtension(const Sequence& sequence) : _sequence(sequence) {
build();
}
explicit LongestCommonExtension(Sequence&& sequence) : _sequence(std::move(sequence)) {
build();
}
int size() const {
return int(_sequence.size());
}
bool empty() const {
return _sequence.empty();
}
const Sequence& sequence() const {
return _sequence;
}
const std::vector<int>& suffix_array() const {
return _suffix_array;
}
const std::vector<int>& rank() const {
return _rank;
}
const std::vector<int>& lcp_array() const {
return _lcp;
}
int longest_common_extension(int i, int j) const {
int n = size();
assert(0 <= i && i <= n);
assert(0 <= j && j <= n);
if (i == j) return n - i;
if (i == n || j == n) return 0;
int left = _rank[i];
int right = _rank[j];
if (left > right) std::swap(left, right);
return range_min(left, right);
}
int longest_common_extension(int i, int j, int limit) const {
assert(0 <= limit);
return std::min(longest_common_extension(i, j), limit);
}
int lcp(int i, int j) const {
return longest_common_extension(i, j);
}
int operator()(int i, int j) const {
return longest_common_extension(i, j);
}
int compare_suffix(int i, int j) const {
int n = size();
assert(0 <= i && i <= n);
assert(0 <= j && j <= n);
if (i == j) return 0;
int common = longest_common_extension(i, j);
if (i + common == n && j + common == n) return 0;
if (i + common == n) return -1;
if (j + common == n) return 1;
return _sequence[i + common] < _sequence[j + common] ? -1 : 1;
}
int compare(int l1, int r1, int l2, int r2) const {
int n = size();
assert(0 <= l1 && l1 <= r1 && r1 <= n);
assert(0 <= l2 && l2 <= r2 && r2 <= n);
int len1 = r1 - l1;
int len2 = r2 - l2;
int common = longest_common_extension(l1, l2, std::min(len1, len2));
if (common == len1 && common == len2) return 0;
if (common == len1) return -1;
if (common == len2) return 1;
return _sequence[l1 + common] < _sequence[l2 + common] ? -1 : 1;
}
};
} // namespace string
} // namespace m1une
#line 1 "string/longest_common_subsequence.hpp"
#line 9 "string/longest_common_subsequence.hpp"
namespace m1une {
namespace string {
struct LongestCommonSubsequence {
std::vector<std::pair<int, int>> matches;
int length() const {
return int(matches.size());
}
bool empty() const {
return matches.empty();
}
std::vector<int> first_indices() const {
std::vector<int> result;
result.reserve(matches.size());
for (auto [i, j] : matches) {
(void)j;
result.push_back(i);
}
return result;
}
std::vector<int> second_indices() const {
std::vector<int> result;
result.reserve(matches.size());
for (auto [i, j] : matches) {
(void)i;
result.push_back(j);
}
return result;
}
template <class Sequence>
std::vector<std::remove_cv_t<std::remove_reference_t<decltype(std::declval<const Sequence&>()[0])>>>
values_from_first(const Sequence& first) const {
using Value = std::remove_cv_t<std::remove_reference_t<decltype(std::declval<const Sequence&>()[0])>>;
std::vector<Value> result;
result.reserve(matches.size());
for (auto [i, j] : matches) {
(void)j;
result.push_back(first[i]);
}
return result;
}
template <class Sequence>
std::vector<std::remove_cv_t<std::remove_reference_t<decltype(std::declval<const Sequence&>()[0])>>>
values_from_second(const Sequence& second) const {
using Value = std::remove_cv_t<std::remove_reference_t<decltype(std::declval<const Sequence&>()[0])>>;
std::vector<Value> result;
result.reserve(matches.size());
for (auto [i, j] : matches) {
(void)i;
result.push_back(second[j]);
}
return result;
}
};
template <class FirstSequence, class SecondSequence>
int longest_common_subsequence_length(const FirstSequence& first, const SecondSequence& second) {
int n = int(first.size());
int m = int(second.size());
if (m <= n) {
std::vector<int> dp(m + 1, 0);
for (int i = 0; i < n; i++) {
int diagonal = 0;
for (int j = 0; j < m; j++) {
int up = dp[j + 1];
if (first[i] == second[j]) {
dp[j + 1] = diagonal + 1;
} else {
dp[j + 1] = std::max(dp[j + 1], dp[j]);
}
diagonal = up;
}
}
return dp[m];
} else {
std::vector<int> dp(n + 1, 0);
for (int j = 0; j < m; j++) {
int diagonal = 0;
for (int i = 0; i < n; i++) {
int up = dp[i + 1];
if (first[i] == second[j]) {
dp[i + 1] = diagonal + 1;
} else {
dp[i + 1] = std::max(dp[i + 1], dp[i]);
}
diagonal = up;
}
}
return dp[n];
}
}
template <class FirstSequence, class SecondSequence>
LongestCommonSubsequence longest_common_subsequence(
const FirstSequence& first,
const SecondSequence& second
) {
int n = int(first.size());
int m = int(second.size());
std::vector<std::vector<int>> dp(n + 1, std::vector<int>(m + 1, 0));
for (int i = 0; i < n; i++) {
for (int j = 0; j < m; j++) {
if (first[i] == second[j]) {
dp[i + 1][j + 1] = dp[i][j] + 1;
} else {
dp[i + 1][j + 1] = std::max(dp[i][j + 1], dp[i + 1][j]);
}
}
}
LongestCommonSubsequence result;
result.matches.reserve(dp[n][m]);
int i = n;
int j = m;
while (i > 0 && j > 0) {
if (first[i - 1] == second[j - 1]) {
result.matches.emplace_back(i - 1, j - 1);
i--;
j--;
} else if (dp[i - 1][j] >= dp[i][j - 1]) {
i--;
} else {
j--;
}
}
std::reverse(result.matches.begin(), result.matches.end());
return result;
}
} // namespace string
} // namespace m1une
#line 1 "string/longest_common_substring.hpp"
#line 9 "string/longest_common_substring.hpp"
#line 11 "string/longest_common_substring.hpp"
namespace m1une {
namespace string {
struct LongestCommonSubstring {
int first_left = 0;
int first_right = 0;
int second_left = 0;
int second_right = 0;
int length() const {
assert(first_right - first_left == second_right - second_left);
return first_right - first_left;
}
bool empty() const {
return length() == 0;
}
std::pair<int, int> first_interval() const {
return {first_left, first_right};
}
std::pair<int, int> second_interval() const {
return {second_left, second_right};
}
};
namespace detail {
template <class Sequence>
std::vector<int> compressed_join_with_separator(const Sequence& first, const Sequence& second) {
using Value = std::remove_cv_t<std::remove_reference_t<decltype(first[0])>>;
std::vector<Value> values;
values.reserve(first.size() + second.size());
for (const auto& value : first) values.push_back(value);
for (const auto& value : second) values.push_back(value);
std::sort(values.begin(), values.end());
values.erase(std::unique(values.begin(), values.end()), values.end());
std::vector<int> joined;
joined.reserve(first.size() + second.size() + 1);
for (const auto& value : first) {
joined.push_back(int(std::lower_bound(values.begin(), values.end(), value) - values.begin()) + 2);
}
joined.push_back(1);
for (const auto& value : second) {
joined.push_back(int(std::lower_bound(values.begin(), values.end(), value) - values.begin()) + 2);
}
return joined;
}
} // namespace detail
template <class Sequence>
LongestCommonSubstring longest_common_substring(const Sequence& first, const Sequence& second) {
int n = int(first.size());
int m = int(second.size());
std::vector<int> joined = detail::compressed_join_with_separator(first, second);
std::vector<int> suffixes = suffix_array(joined);
std::vector<int> lcp = lcp_array(joined, suffixes);
LongestCommonSubstring result;
for (int i = 0; i + 1 < int(suffixes.size()); i++) {
int a = suffixes[i];
int b = suffixes[i + 1];
if (a == n || b == n) continue;
bool a_first = a < n;
bool b_first = b < n;
if (a_first == b_first) continue;
int first_left = a_first ? a : b;
int second_left = a_first ? b - n - 1 : a - n - 1;
int length = lcp[i];
length = std::min(length, n - first_left);
length = std::min(length, m - second_left);
if (length > result.length()) {
result.first_left = first_left;
result.first_right = first_left + length;
result.second_left = second_left;
result.second_right = second_left + length;
}
}
return result;
}
} // namespace string
} // namespace m1une
#line 1 "string/lyndon_factorization.hpp"
#line 6 "string/lyndon_factorization.hpp"
#line 1 "string/minimum_rotation.hpp"
namespace m1une {
namespace string {
// Returns the smallest starting index of a lexicographically minimum cyclic shift.
template <class Sequence>
int minimum_cyclic_shift(const Sequence& sequence) {
const int size = int(sequence.size());
if (size == 0) return 0;
auto less = [&](int left, int right) {
return sequence[left < size ? left : left - size] <
sequence[right < size ? right : right - size];
};
int answer = 0;
int start = 0;
while (start < size) {
answer = start;
int scan = start + 1;
int matched = start;
while (scan < 2 * size && !less(scan, matched)) {
if (less(matched, scan)) {
matched = start;
} else {
matched++;
}
scan++;
}
const int period = scan - matched;
while (start <= matched) start += period;
}
return answer;
}
} // namespace string
} // namespace m1une
#line 8 "string/lyndon_factorization.hpp"
namespace m1une {
namespace string {
// Returns boundaries 0 = a[0] < a[1] < ... < a[k] = sequence.size()
// of the Lyndon factorization.
template <class Sequence>
std::vector<int> lyndon_factor_boundaries(const Sequence& sequence) {
int n = int(sequence.size());
std::vector<int> boundaries;
boundaries.push_back(0);
int i = 0;
while (i < n) {
int j = i + 1;
int k = i;
while (j < n && !(sequence[j] < sequence[k])) {
if (sequence[k] < sequence[j]) {
k = i;
} else {
k++;
}
j++;
}
int length = j - k;
while (i <= k) {
i += length;
boundaries.push_back(i);
}
}
return boundaries;
}
// Returns half-open intervals [left, right) of the Lyndon factorization.
template <class Sequence>
std::vector<std::pair<int, int>> lyndon_factorization(const Sequence& sequence) {
std::vector<int> boundaries = lyndon_factor_boundaries(sequence);
std::vector<std::pair<int, int>> factors;
factors.reserve(boundaries.size() - 1);
for (int i = 0; i + 1 < int(boundaries.size()); i++) {
factors.emplace_back(boundaries[i], boundaries[i + 1]);
}
return factors;
}
} // namespace string
} // namespace m1une
#line 1 "string/manacher.hpp"
#line 7 "string/manacher.hpp"
namespace m1une {
namespace string {
struct ManacherResult {
// odd[i] is the radius including center i.
// The palindrome is [i - odd[i] + 1, i + odd[i]).
std::vector<int> odd;
// even[i] is the radius centered between i - 1 and i.
// The palindrome is [i - even[i], i + even[i]).
std::vector<int> even;
int size() const {
return int(odd.size());
}
bool empty() const {
return odd.empty();
}
bool is_palindrome(int left, int right) const {
int n = size();
assert(0 <= left && left <= right && right <= n);
int length = right - left;
if (length == 0) return true;
if (length & 1) {
int center = (left + right) / 2;
return length / 2 + 1 <= odd[center];
}
int center = (left + right) / 2;
return length / 2 <= even[center];
}
int longest_length() const {
int result = 0;
for (int radius : odd) result = std::max(result, 2 * radius - 1);
for (int radius : even) result = std::max(result, 2 * radius);
return result;
}
};
template <class Sequence>
ManacherResult manacher(const Sequence& sequence) {
int n = int(sequence.size());
ManacherResult result;
result.odd.assign(n, 0);
result.even.assign(n, 0);
int left = 0;
int right = -1;
for (int i = 0; i < n; i++) {
int radius = i > right ? 1 : std::min(result.odd[left + right - i], right - i + 1);
while (
0 <= i - radius &&
i + radius < n &&
sequence[i - radius] == sequence[i + radius]
) {
radius++;
}
result.odd[i] = radius;
if (right < i + radius - 1) {
left = i - radius + 1;
right = i + radius - 1;
}
}
left = 0;
right = -1;
for (int i = 0; i < n; i++) {
int radius = i > right ? 0 : std::min(result.even[left + right - i + 1], right - i + 1);
while (
0 <= i - radius - 1 &&
i + radius < n &&
sequence[i - radius - 1] == sequence[i + radius]
) {
radius++;
}
result.even[i] = radius;
if (right < i + radius - 1) {
left = i - radius;
right = i + radius - 1;
}
}
return result;
}
} // namespace string
} // namespace m1une
#line 1 "string/map_trie.hpp"
#line 6 "string/map_trie.hpp"
#include <functional>
#line 8 "string/map_trie.hpp"
#include <map>
#line 10 "string/map_trie.hpp"
namespace m1une {
namespace string {
// A multiset trie whose outgoing edges are stored in ordered maps.
template <class Symbol, class Compare = std::less<Symbol>>
struct MapTrie {
using node_id = int;
static constexpr node_id null_node = -1;
struct Node {
std::map<Symbol, node_id, Compare> child;
int subtree_count = 0;
int terminal_count = 0;
};
private:
std::vector<Node> _nodes;
int _distinct_size;
node_id new_node() {
assert(_nodes.size() < std::size_t(std::numeric_limits<int>::max()));
_nodes.emplace_back();
return int(_nodes.size()) - 1;
}
template <class Sequence>
node_id find_node(const Sequence& sequence) const {
node_id node = 0;
for (const auto& symbol : sequence) {
auto iterator = _nodes[node].child.find(symbol);
if (iterator == _nodes[node].child.end()) return null_node;
node = iterator->second;
if (_nodes[node].subtree_count == 0) return null_node;
}
return node;
}
public:
MapTrie() : _nodes(1), _distinct_size(0) {}
int size() const {
return _nodes[0].subtree_count;
}
int distinct_size() const {
return _distinct_size;
}
bool empty() const {
return size() == 0;
}
node_id root() const {
return 0;
}
const Node& node(node_id id) const {
assert(0 <= id && std::size_t(id) < _nodes.size());
return _nodes[id];
}
template <class Sequence>
node_id find(const Sequence& sequence) const {
return find_node(sequence);
}
std::size_t node_count() const {
return _nodes.size();
}
void reserve(std::size_t node_capacity) {
_nodes.reserve(node_capacity);
}
void clear() {
_nodes.clear();
_nodes.emplace_back();
_distinct_size = 0;
}
template <class Sequence>
node_id insert(const Sequence& sequence, int multiplicity = 1) {
assert(0 < multiplicity);
node_id node = 0;
_nodes[node].subtree_count += multiplicity;
for (const auto& symbol : sequence) {
auto iterator = _nodes[node].child.find(symbol);
node_id child;
if (iterator == _nodes[node].child.end()) {
child = new_node();
_nodes[node].child.emplace(symbol, child);
} else {
child = iterator->second;
}
node = child;
_nodes[node].subtree_count += multiplicity;
}
if (_nodes[node].terminal_count == 0) _distinct_size++;
_nodes[node].terminal_count += multiplicity;
return node;
}
template <class Sequence>
int count(const Sequence& sequence) const {
node_id node = find_node(sequence);
return node == null_node ? 0 : _nodes[node].terminal_count;
}
template <class Sequence>
bool contains(const Sequence& sequence) const {
return count(sequence) != 0;
}
// Returns the number of stored sequences beginning with prefix.
template <class Sequence>
int prefix_count(const Sequence& prefix) const {
node_id node = find_node(prefix);
return node == null_node ? 0 : _nodes[node].subtree_count;
}
template <class Sequence>
bool starts_with(const Sequence& prefix) const {
return prefix_count(prefix) != 0;
}
template <class Sequence>
bool erase_one(const Sequence& sequence) {
node_id terminal = find_node(sequence);
if (terminal == null_node || _nodes[terminal].terminal_count == 0) {
return false;
}
node_id node = 0;
_nodes[node].subtree_count--;
for (const auto& symbol : sequence) {
node = _nodes[node].child.find(symbol)->second;
_nodes[node].subtree_count--;
}
_nodes[node].terminal_count--;
if (_nodes[node].terminal_count == 0) _distinct_size--;
return true;
}
template <class Sequence>
bool erase(const Sequence& sequence) {
return erase_one(sequence);
}
template <class Sequence>
int erase_all(const Sequence& sequence) {
int multiplicity = count(sequence);
if (multiplicity == 0) return 0;
node_id node = 0;
_nodes[node].subtree_count -= multiplicity;
for (const auto& symbol : sequence) {
node = _nodes[node].child.find(symbol)->second;
_nodes[node].subtree_count -= multiplicity;
}
_nodes[node].terminal_count = 0;
_distinct_size--;
return multiplicity;
}
// Calls callback(length, multiplicity) for every stored prefix.
// The empty prefix is reported with length 0 when it is stored.
template <class Sequence, class Callback>
void for_each_prefix(const Sequence& sequence, Callback callback) const {
node_id node = 0;
if (_nodes[node].terminal_count != 0) {
callback(0, _nodes[node].terminal_count);
}
int length = 0;
for (const auto& symbol : sequence) {
auto iterator = _nodes[node].child.find(symbol);
if (iterator == _nodes[node].child.end()) return;
node = iterator->second;
if (_nodes[node].subtree_count == 0) return;
length++;
if (_nodes[node].terminal_count != 0) {
callback(length, _nodes[node].terminal_count);
}
}
}
// Returns the length of the longest stored sequence that is a prefix.
// Returns -1 when no stored prefix exists.
template <class Sequence>
int longest_prefix(const Sequence& sequence) const {
int result = _nodes[0].terminal_count == 0 ? -1 : 0;
for_each_prefix(sequence, [&result](int length, int) {
result = length;
});
return result;
}
};
} // namespace string
} // namespace m1une
#line 1 "string/palindrome_lexicographical_order.hpp"
#line 9 "string/palindrome_lexicographical_order.hpp"
#line 12 "string/palindrome_lexicographical_order.hpp"
namespace m1une {
namespace string {
// Indexes the distinct nonempty palindromic substrings of one sequence.
template <
class Sequence = std::string,
int AlphabetSize = 26,
int FirstCharacter = 'a'
>
struct PalindromeLexicographicalOrder {
static_assert(0 < AlphabetSize);
using eertree_type = Eertree<AlphabetSize, FirstCharacter>;
using node_id = typename eertree_type::node_id;
private:
Sequence _sequence;
eertree_type _eertree;
std::vector<node_id> _nodes_in_order;
std::vector<int> _order_of_node;
template <class Symbol>
static int symbol_index(const Symbol& symbol) {
int index = int(symbol) - FirstCharacter;
assert(0 <= index && index < AlphabetSize);
return index;
}
void build_order() {
const int node_count = _eertree.node_count();
std::vector<std::vector<node_id>> suffix_children(node_count);
for (node_id id = 0; id < node_count; id++) {
if (id == eertree_type::odd_root) continue;
suffix_children[_eertree.node(id).suffix_link].push_back(id);
}
std::vector<int> enter(node_count);
std::vector<int> leave(node_count);
std::vector<std::pair<node_id, bool>> stack;
stack.reserve(2 * node_count);
stack.emplace_back(eertree_type::odd_root, false);
int timer = 0;
while (!stack.empty()) {
auto [id, exiting] = stack.back();
stack.pop_back();
if (exiting) {
leave[id] = timer;
continue;
}
enter[id] = timer++;
stack.emplace_back(id, true);
const auto& children = suffix_children[id];
for (int i = int(children.size()) - 1; i >= 0; i--) {
stack.emplace_back(children[i], false);
}
}
std::vector<int> suffixes = suffix_array(_sequence);
std::vector<int> suffix_rank(_sequence.size());
for (int rank = 0; rank < int(suffixes.size()); rank++) {
suffix_rank[suffixes[rank]] = rank;
}
_nodes_in_order.resize(_eertree.size());
for (int i = 0; i < _eertree.size(); i++) {
_nodes_in_order[i] = i + 2;
}
auto is_ancestor = [&](node_id ancestor, node_id descendant) {
return
enter[ancestor] <= enter[descendant] &&
leave[descendant] <= leave[ancestor];
};
std::sort(
_nodes_in_order.begin(),
_nodes_in_order.end(),
[&](node_id first, node_id second) {
if (first == second) return false;
// A palindromic prefix is also a palindromic suffix, so prefix
// cases are exactly the ancestor cases in the suffix-link tree.
if (is_ancestor(first, second)) return true;
if (is_ancestor(second, first)) return false;
// Otherwise the first mismatch occurs inside both substrings,
// and the ranks of representative suffixes give their order.
int first_start = _eertree.first_occurrence(first).first;
int second_start = _eertree.first_occurrence(second).first;
return suffix_rank[first_start] < suffix_rank[second_start];
}
);
_order_of_node.assign(node_count, -1);
for (int order = 0; order < size(); order++) {
_order_of_node[_nodes_in_order[order]] = order;
}
}
public:
PalindromeLexicographicalOrder() : _order_of_node(2, -1) {}
explicit PalindromeLexicographicalOrder(const Sequence& sequence)
: _sequence(sequence), _eertree(_sequence) {
build_order();
}
explicit PalindromeLexicographicalOrder(Sequence&& sequence)
: _sequence(std::move(sequence)), _eertree(_sequence) {
build_order();
}
int size() const {
return int(_nodes_in_order.size());
}
bool empty() const {
return _nodes_in_order.empty();
}
int text_length() const {
return int(_sequence.size());
}
const Sequence& sequence() const {
return _sequence;
}
const eertree_type& eertree() const {
return _eertree;
}
const std::vector<node_id>& nodes_in_order() const {
return _nodes_in_order;
}
int order_of_node(node_id id) const {
assert(2 <= id && id < _eertree.node_count());
return _order_of_node[id];
}
node_id node_by_order(int order) const {
assert(0 <= order && order < size());
return _nodes_in_order[order];
}
template <class Palindrome>
node_id find(const Palindrome& palindrome) const {
const int length = int(palindrome.size());
if (length == 0) return eertree_type::null_node;
for (int i = 0; i < length / 2; i++) {
if (palindrome[i] != palindrome[length - 1 - i]) {
return eertree_type::null_node;
}
}
node_id id =
length & 1 ? eertree_type::odd_root : eertree_type::even_root;
for (int i = (length - 1) / 2; i >= 0; i--) {
int symbol = symbol_index(palindrome[i]);
id = _eertree.node(id).next[symbol];
if (id == eertree_type::null_node) return id;
}
return id;
}
template <class Palindrome>
bool contains(const Palindrome& palindrome) const {
return find(palindrome) != eertree_type::null_node;
}
template <class Palindrome>
int order_of_palindrome(const Palindrome& palindrome) const {
node_id id = find(palindrome);
return id == eertree_type::null_node ? -1 : order_of_node(id);
}
std::pair<int, int> representative_occurrence(int order) const {
return _eertree.first_occurrence(node_by_order(order));
}
Sequence palindrome(int order) const {
auto [left, right] = representative_occurrence(order);
return Sequence(_sequence.begin() + left, _sequence.begin() + right);
}
Sequence kth(int order) const {
return palindrome(order);
}
};
} // namespace string
} // namespace m1une
#line 1 "string/prefix_substring_lcs.hpp"
#line 9 "string/prefix_substring_lcs.hpp"
namespace m1une {
namespace string {
// Answers LCS-length queries between a prefix of the first sequence and a
// substring of the second sequence. Queries are evaluated as one offline batch.
template <class FirstSequence, class SecondSequence>
class PrefixSubstringLcs {
private:
struct Query {
int first_prefix;
int second_left;
int second_right;
};
FirstSequence _first;
SecondSequence _second;
std::vector<Query> _queries;
public:
PrefixSubstringLcs(FirstSequence first, SecondSequence second)
: _first(std::move(first)), _second(std::move(second)) {}
int first_size() const {
return int(_first.size());
}
int second_size() const {
return int(_second.size());
}
int query_count() const {
return int(_queries.size());
}
bool empty() const {
return _queries.empty();
}
void reserve(int query_capacity) {
assert(0 <= query_capacity);
_queries.reserve(query_capacity);
}
void clear() {
_queries.clear();
}
// Adds LCS(first[0..first_prefix), second[second_left..second_right)) and
// returns its insertion-order ID.
int add_query(int first_prefix, int second_left, int second_right) {
assert(0 <= first_prefix && first_prefix <= first_size());
assert(0 <= second_left && second_left <= second_right);
assert(second_right <= second_size());
const int id = query_count();
_queries.push_back(Query{first_prefix, second_left, second_right});
return id;
}
std::vector<int> calculate() const {
const int first_length = first_size();
const int second_length = second_size();
const int count = query_count();
std::vector<int> answers(count, 0);
if (count == 0 || first_length == 0 || second_length == 0) {
return answers;
}
std::vector<std::vector<int>> queries_by_prefix(first_length + 1);
for (int id = 0; id < count; id++) {
const Query& query = _queries[id];
if (query.first_prefix > 0 &&
query.second_left < query.second_right) {
queries_by_prefix[query.first_prefix].push_back(id);
}
}
// seaweed[j] is the bottom endpoint of the seaweed entering at j for
// the current prefix of the first sequence. -1 denotes the left edge.
std::vector<int> seaweed(second_length);
for (int j = 0; j < second_length; j++) seaweed[j] = j;
std::vector<int> heads(second_length + 1, -1);
std::vector<int> next(count, -1);
std::vector<int> fenwick(second_length + 1, 0);
for (int i = 0; i < first_length; i++) {
int displaced = -1;
for (int j = 0; j < second_length; j++) {
if (_first[i] == _second[j] || seaweed[j] < displaced) {
std::swap(seaweed[j], displaced);
}
}
const std::vector<int>& prefix_queries = queries_by_prefix[i + 1];
if (prefix_queries.empty()) continue;
std::fill(heads.begin(), heads.end(), -1);
for (int id : prefix_queries) {
const int right = _queries[id].second_right;
next[id] = heads[right];
heads[right] = id;
}
std::fill(fenwick.begin(), fenwick.end(), 0);
int inserted = 0;
for (int right = 1; right <= second_length; right++) {
const int endpoint = seaweed[right - 1];
if (endpoint >= 0) {
inserted++;
for (int position = endpoint + 1;
position <= second_length;
position += position & -position) {
fenwick[position]++;
}
}
for (int id = heads[right]; id != -1; id = next[id]) {
const int left = _queries[id].second_left;
int below_left = 0;
for (int position = left; position > 0;
position -= position & -position) {
below_left += fenwick[position];
}
const int crossing = inserted - below_left;
answers[id] = (right - left) - crossing;
}
}
}
return answers;
}
};
template <class FirstSequence, class SecondSequence>
PrefixSubstringLcs(FirstSequence&&, SecondSequence&&)
-> PrefixSubstringLcs<
std::decay_t<FirstSequence>,
std::decay_t<SecondSequence>
>;
} // namespace string
} // namespace m1une
#line 1 "string/rolling_hash.hpp"
#line 8 "string/rolling_hash.hpp"
namespace m1une {
namespace string {
// Standard Rolling Hash for static strings.
// Precomputes hashes to answer substring queries in O(1).
// Provides advanced operations like LCP, lexicographical comparison, and string repetition in O(log N).
template <long long Base = 10007, long long Mod = (1LL << 61) - 1>
struct RollingHash {
std::string s;
std::vector<long long> hash;
std::vector<long long> power;
RollingHash() = default;
// Constructs the rolling hash table for the given string.
explicit RollingHash(const std::string& str) : s(str) {
int n = s.size();
hash.assign(n + 1, 0);
power.assign(n + 1, 1);
for (int i = 0; i < n; ++i) {
// Use __int128_t to prevent overflow during multiplication
hash[i + 1] = (static_cast<__int128_t>(hash[i]) * Base + s[i]) % Mod;
power[i + 1] = (static_cast<__int128_t>(power[i]) * Base) % Mod;
}
}
// Returns the hash of the substring S[l..r) in O(1).
long long get(int l, int r) const {
long long res = hash[r] - (static_cast<__int128_t>(hash[l]) * power[r - l]) % Mod;
if (res < 0) res += Mod;
return res;
}
// Returns the hash of the concatenated substrings S[l1..r1) and S[l2..r2).
long long concat(int l1, int r1, int l2, int r2) const {
long long h1 = get(l1, r1);
long long h2 = get(l2, r2);
return combine(h1, h2, power[r2 - l2]);
}
// Calculates the Longest Common Prefix (LCP) length of S[l1..r1) and S[l2..r2) in O(log N).
int lcp(int l1, int r1, int l2, int r2) const {
int len = std::min(r1 - l1, r2 - l2);
int low = 0, high = len + 1;
while (high - low > 1) {
int mid = low + (high - low) / 2;
if (get(l1, l1 + mid) == get(l2, l2 + mid)) {
low = mid;
} else {
high = mid;
}
}
return low;
}
// Lexicographically compares S[l1..r1) and S[l2..r2) in O(log N).
// Returns -1 if S[l1..r1) < S[l2..r2), 0 if equal, and 1 if S[l1..r1) > S[l2..r2).
int compare(int l1, int r1, int l2, int r2) const {
int l = lcp(l1, r1, l2, r2);
bool end1 = (l1 + l == r1);
bool end2 = (l2 + l == r2);
if (end1 && end2) return 0;
if (end1) return -1;
if (end2) return 1;
return s[l1 + l] < s[l2 + l] ? -1 : 1;
}
// Returns the hash of the substring S[l..r) repeated 'k' times.
long long repeat(int l, int r, long long k) const {
long long h = get(l, r);
long long p = power[r - l];
return repeat_hash(h, p, k);
}
// --- Static Helpers for dynamic processing and Monoid integration ---
// Computes the hash of a single string in O(N) time and O(1) space.
static long long compute_hash(const std::string& str) {
long long h = 0;
for (char c : str) {
h = (static_cast<__int128_t>(h) * Base + c) % Mod;
}
return h;
}
// Combines two hashes. Equivalent to concatenating string 'b' to the right of string 'a'.
static constexpr long long combine(long long h1, long long h2, long long base_power2) {
return (static_cast<__int128_t>(h1) * base_power2 + h2) % Mod;
}
// Returns the hash of a string (with hash 'h' and base_power 'p') repeated 'k' times.
static constexpr long long repeat_hash(long long h, long long p, long long k) {
long long res_h = 0;
long long res_p = 1;
long long cur_h = h;
long long cur_p = p;
while (k > 0) {
if (k & 1) {
res_h = combine(res_h, cur_h, cur_p);
res_p = (static_cast<__int128_t>(res_p) * cur_p) % Mod;
}
cur_h = combine(cur_h, cur_h, cur_p);
cur_p = (static_cast<__int128_t>(cur_p) * cur_p) % Mod;
k >>= 1;
}
return res_h;
}
// Creates the state pair {hash_value, base_power} for a single character.
static constexpr std::pair<long long, long long> make_single(long long c) {
return {c % Mod, Base % Mod};
}
};
} // namespace string
} // namespace m1une
#line 1 "string/runs.hpp"
#line 5 "string/runs.hpp"
#include <set>
#line 8 "string/runs.hpp"
namespace m1une {
namespace string {
struct Run {
int period;
int left;
int right;
bool operator==(const Run&) const = default;
};
namespace internal {
template <class Sequence>
class RunEnumerator {
private:
const Sequence& _sequence;
int _size;
std::vector<std::vector<std::pair<int, int>>> _candidates;
template <class Access>
static std::vector<int> z_algorithm(int length, Access access) {
std::vector<int> z(length + 1, 0);
if (length == 0) return z;
z[0] = length;
int left = 0;
int right = 0;
for (int i = 1; i < length; i++) {
if (i < right) z[i] = std::min(right - i, z[i - left]);
while (
i + z[i] < length &&
access(z[i]) == access(i + z[i])
) {
z[i]++;
}
if (right < i + z[i]) {
left = i;
right = i + z[i];
}
}
return z;
}
decltype(auto) element(int index, bool reversed) const {
int original_index = reversed ? _size - 1 - index : index;
return _sequence[original_index];
}
void add_candidate(int period, int left, int right, bool reversed) {
if (reversed) {
left = _size - left;
right = _size - right;
std::swap(left, right);
}
_candidates[period].emplace_back(left, right);
}
void collect(int range_left, int range_right, int phase, bool reversed) {
if (range_right - range_left <= 1) return;
int middle = (range_left + range_right + phase) / 2;
collect(range_left, middle, phase, reversed);
collect(middle, range_right, phase, reversed);
int left_length = middle - range_left;
int right_length = range_right - middle;
std::vector<int> left_z = z_algorithm(left_length, [&](int index) -> decltype(auto) {
return element(middle - 1 - index, reversed);
});
int combined_length = right_length + range_right - range_left;
std::vector<int> right_z = z_algorithm(combined_length, [&](int index) -> decltype(auto) {
if (index < right_length) return element(middle + index, reversed);
return element(range_left + index - right_length, reversed);
});
for (int start = middle - 1; start >= range_left; start--) {
int period = middle - start;
int extend_left = std::min(start - range_left, left_z[period]);
int extend_right = std::min(
range_right - middle,
right_z[range_right - range_left - period]
);
int left = start - extend_left;
int right = middle + extend_right;
if (right - left >= 2 * period) {
add_candidate(period, left, right, reversed);
}
}
}
public:
explicit RunEnumerator(const Sequence& sequence)
: _sequence(sequence),
_size(int(sequence.size())),
_candidates(_size / 2 + 1) {}
std::vector<Run> enumerate() {
collect(0, _size, 0, true);
collect(0, _size, 1, false);
std::set<std::pair<int, int>> used_intervals;
std::vector<Run> result;
for (int period = 1; period <= _size / 2; period++) {
std::vector<std::pair<int, int>>& candidates = _candidates[period];
std::sort(
candidates.begin(),
candidates.end(),
[](const auto& first, const auto& second) {
if (first.first != second.first) {
return first.first < second.first;
}
return first.second > second.second;
}
);
int farthest_right = -1;
for (const auto& interval : candidates) {
if (interval.second <= farthest_right) continue;
farthest_right = interval.second;
if (!used_intervals.insert(interval).second) continue;
result.push_back(Run{period, interval.first, interval.second});
}
}
return result;
}
};
} // namespace internal
// Returns all runs as (minimum period, maximal half-open interval),
// sorted lexicographically by (period, left, right).
template <class Sequence>
std::vector<Run> enumerate_runs(const Sequence& sequence) {
return internal::RunEnumerator<Sequence>(sequence).enumerate();
}
} // namespace string
} // namespace m1une
#line 1 "string/string_hash.hpp"
#line 5 "string/string_hash.hpp"
#include <cstdint>
#line 7 "string/string_hash.hpp"
#include <string_view>
namespace m1une {
namespace string {
struct StringHash {
std::uint32_t first;
std::uint32_t second;
std::uint32_t first_power;
std::uint32_t second_power;
std::size_t length;
friend constexpr bool operator==(const StringHash& left, const StringHash& right) {
return left.length == right.length && left.first == right.first && left.second == right.second;
}
};
namespace string_hash_detail {
inline constexpr std::uint64_t first_mod = 1'000'000'007;
inline constexpr std::uint64_t second_mod = 1'000'000'009;
inline constexpr std::uint64_t base = 911'382'323;
} // namespace string_hash_detail
// Computes a double polynomial hash. Bytes are interpreted as unsigned.
constexpr StringHash hash_string(std::string_view value) {
using namespace string_hash_detail;
std::uint64_t first = 0;
std::uint64_t second = 0;
std::uint64_t first_power = 1;
std::uint64_t second_power = 1;
for (char character : value) {
std::uint64_t symbol = static_cast<unsigned char>(character) + std::uint64_t(1);
first = (first * base + symbol) % first_mod;
second = (second * base + symbol) % second_mod;
first_power = first_power * base % first_mod;
second_power = second_power * base % second_mod;
}
return StringHash{
static_cast<std::uint32_t>(first),
static_cast<std::uint32_t>(second),
static_cast<std::uint32_t>(first_power),
static_cast<std::uint32_t>(second_power),
value.size(),
};
}
constexpr StringHash hash_string(const std::string& value) {
return hash_string(std::string_view(value));
}
constexpr StringHash hash_string(const char* value) {
return hash_string(std::string_view(value));
}
// Returns the hash of the concatenation represented by `left` and `right`.
constexpr StringHash concat_string_hash(const StringHash& left, const StringHash& right) {
using namespace string_hash_detail;
return StringHash{
static_cast<std::uint32_t>((std::uint64_t(left.first) * right.first_power + right.first) % first_mod),
static_cast<std::uint32_t>((std::uint64_t(left.second) * right.second_power + right.second) % second_mod),
static_cast<std::uint32_t>(std::uint64_t(left.first_power) * right.first_power % first_mod),
static_cast<std::uint32_t>(std::uint64_t(left.second_power) * right.second_power % second_mod),
left.length + right.length,
};
}
// Hash adapter for std::unordered_map and std::unordered_set.
struct StringHasher {
using is_transparent = void;
constexpr std::size_t operator()(std::string_view value) const {
return operator()(hash_string(value));
}
constexpr std::size_t operator()(const std::string& value) const {
return operator()(std::string_view(value));
}
constexpr std::size_t operator()(const char* value) const {
return operator()(std::string_view(value));
}
constexpr std::size_t operator()(const StringHash& value) const {
std::uint64_t combined = (std::uint64_t(value.first) << 32) | value.second;
combined ^= std::uint64_t(value.length) + 0x9e3779b97f4a7c15ULL;
combined ^= combined >> 30;
combined *= 0xbf58476d1ce4e5b9ULL;
combined ^= combined >> 27;
combined *= 0x94d049bb133111ebULL;
combined ^= combined >> 31;
return static_cast<std::size_t>(combined);
}
};
} // namespace string
} // namespace m1une
#line 1 "string/suffix_automaton.hpp"
#line 11 "string/suffix_automaton.hpp"
namespace m1une {
namespace string {
template <int AlphabetSize = 26, int FirstCharacter = 'a'>
struct SuffixAutomaton {
static_assert(0 < AlphabetSize);
using state_id = int;
static constexpr state_id root_state = 0;
static constexpr state_id null_state = -1;
struct State {
std::array<state_id, AlphabetSize> next;
state_id suffix_link;
int length;
int first_end;
int direct_occurrences;
bool clone;
State(int length_value = 0)
: suffix_link(null_state),
length(length_value),
first_end(0),
direct_occurrences(0),
clone(false) {
next.fill(null_state);
}
};
private:
std::vector<State> _states;
state_id _last;
int _text_length;
template <class Symbol>
static int symbol_index(const Symbol& symbol) {
int index = int(symbol) - FirstCharacter;
assert(0 <= index && index < AlphabetSize);
return index;
}
state_id new_state(int length) {
assert(_states.size() < std::size_t(std::numeric_limits<int>::max()));
_states.emplace_back(length);
return int(_states.size()) - 1;
}
public:
SuffixAutomaton() {
clear();
}
template <class Sequence>
explicit SuffixAutomaton(const Sequence& sequence) {
clear();
build(sequence);
}
int state_count() const {
return int(_states.size());
}
int size() const {
return state_count();
}
bool empty() const {
return _text_length == 0;
}
int text_length() const {
return _text_length;
}
state_id root() const {
return root_state;
}
state_id last() const {
return _last;
}
const State& state(state_id id) const {
assert(0 <= id && id < state_count());
return _states[id];
}
const std::vector<State>& states() const {
return _states;
}
int minimum_length(state_id id) const {
assert(0 <= id && id < state_count());
return id == root_state ? 0 : _states[_states[id].suffix_link].length + 1;
}
template <class Symbol>
state_id transition(state_id id, const Symbol& symbol) const {
assert(0 <= id && id < state_count());
return _states[id].next[symbol_index(symbol)];
}
void reserve(std::size_t text_capacity) {
_states.reserve(2 * text_capacity);
}
void clear() {
_states.clear();
_states.emplace_back();
_last = root_state;
_text_length = 0;
}
template <class Symbol>
state_id add(const Symbol& value) {
int symbol = symbol_index(value);
assert(_text_length < std::numeric_limits<int>::max());
_text_length++;
state_id current = new_state(_states[_last].length + 1);
_states[current].first_end = _text_length;
_states[current].direct_occurrences = 1;
state_id p = _last;
while (p != null_state && _states[p].next[symbol] == null_state) {
_states[p].next[symbol] = current;
p = _states[p].suffix_link;
}
if (p == null_state) {
_states[current].suffix_link = root_state;
} else {
state_id q = _states[p].next[symbol];
if (_states[p].length + 1 == _states[q].length) {
_states[current].suffix_link = q;
} else {
state_id clone = new_state(_states[p].length + 1);
_states[clone] = _states[q];
_states[clone].length = _states[p].length + 1;
_states[clone].direct_occurrences = 0;
_states[clone].clone = true;
while (p != null_state && _states[p].next[symbol] == q) {
_states[p].next[symbol] = clone;
p = _states[p].suffix_link;
}
_states[q].suffix_link = clone;
_states[current].suffix_link = clone;
}
}
_last = current;
return current;
}
template <class Sequence>
void build(const Sequence& sequence) {
for (const auto& symbol : sequence) add(symbol);
}
template <class Sequence>
state_id find(const Sequence& sequence) const {
state_id current = root_state;
for (const auto& symbol : sequence) {
current = transition(current, symbol);
if (current == null_state) return null_state;
}
return current;
}
template <class Sequence>
bool contains(const Sequence& sequence) const {
return find(sequence) != null_state;
}
std::vector<state_id> length_order() const {
std::vector<int> count(_text_length + 1, 0);
for (const State& current : _states) count[current.length]++;
for (int length = 1; length <= _text_length; length++) count[length] += count[length - 1];
std::vector<state_id> order(state_count());
for (state_id id = state_count() - 1; id >= 0; id--) {
order[--count[_states[id].length]] = id;
}
return order;
}
std::vector<long long> occurrence_counts() const {
std::vector<long long> result(state_count(), 0);
for (state_id id = 0; id < state_count(); id++) {
result[id] = _states[id].direct_occurrences;
}
std::vector<state_id> order = length_order();
for (int i = int(order.size()) - 1; i > 0; i--) {
state_id id = order[i];
result[_states[id].suffix_link] += result[id];
}
return result;
}
std::vector<bool> terminal_states() const {
std::vector<bool> result(state_count(), false);
for (state_id id = _last; id != null_state; id = _states[id].suffix_link) {
result[id] = true;
}
return result;
}
long long distinct_substring_count() const {
long long result = 0;
for (state_id id = 1; id < state_count(); id++) {
result += _states[id].length - _states[_states[id].suffix_link].length;
}
return result;
}
std::pair<int, int> longest_representative(state_id id) const {
assert(0 <= id && id < state_count());
int end = _states[id].first_end;
return {end - _states[id].length, end};
}
template <class Sequence>
std::pair<int, int> representative_occurrence(const Sequence& sequence) const {
state_id id = root_state;
int length = 0;
for (const auto& symbol : sequence) {
id = transition(id, symbol);
if (id == null_state) return {-1, -1};
length++;
}
int end = _states[id].first_end;
return {end - length, end};
}
template <class Sequence>
std::pair<int, int> longest_common_substring(const Sequence& sequence) const {
state_id current = root_state;
int current_length = 0;
int best_length = 0;
int best_end = 0;
int end = 0;
for (const auto& value : sequence) {
int symbol = symbol_index(value);
while (current != root_state && _states[current].next[symbol] == null_state) {
current = _states[current].suffix_link;
current_length = std::min(current_length, _states[current].length);
}
state_id next = _states[current].next[symbol];
if (next == null_state) {
current = root_state;
current_length = 0;
} else {
current = next;
current_length++;
}
end++;
if (best_length < current_length) {
best_length = current_length;
best_end = end;
}
}
return {best_end - best_length, best_end};
}
};
} // namespace string
} // namespace m1une
#line 1 "string/suffix_tree.hpp"
#line 11 "string/suffix_tree.hpp"
namespace m1une {
namespace string {
template <int AlphabetSize = 26, int FirstCharacter = 'a'>
struct SuffixTree {
static_assert(0 < AlphabetSize);
using node_id = int;
static constexpr node_id root_node = 0;
static constexpr node_id null_node = -1;
static constexpr int terminal_symbol = AlphabetSize;
struct Node {
std::array<node_id, AlphabetSize + 1> next;
node_id suffix_link;
node_id parent;
int left;
int right;
int suffix_start;
int representative_suffix;
int leaf_count;
int incoming_symbol;
node_id first_child;
node_id next_sibling;
int child_count;
Node(int left_value = 0, int right_value = 0, node_id parent_value = null_node)
: suffix_link(null_node),
parent(parent_value),
left(left_value),
right(right_value),
suffix_start(-1),
representative_suffix(-1),
leaf_count(0),
incoming_symbol(-1),
first_child(null_node),
next_sibling(null_node),
child_count(0) {
next.fill(null_node);
}
};
struct Locus {
node_id node;
int offset;
explicit operator bool() const {
return node != null_node;
}
friend bool operator==(const Locus&, const Locus&) = default;
};
private:
struct ActivePoint {
node_id node;
int offset;
};
std::vector<Node> _nodes;
std::vector<int> _text;
ActivePoint _active;
int _text_length;
template <class Symbol>
static int symbol_index(const Symbol& symbol) {
int index = int(symbol) - FirstCharacter;
assert(0 <= index && index < AlphabetSize);
return index;
}
int edge_length_unchecked(node_id id) const {
return _nodes[id].right - _nodes[id].left;
}
node_id new_node(int left, int right, node_id parent) {
assert(_nodes.size() < std::size_t(std::numeric_limits<int>::max()));
_nodes.emplace_back(left, right, parent);
return int(_nodes.size()) - 1;
}
ActivePoint go(ActivePoint point, int left, int right) const {
while (left < right) {
if (point.offset == edge_length_unchecked(point.node)) {
point = {_nodes[point.node].next[_text[left]], 0};
if (point.node == null_node) return point;
} else {
if (_text[_nodes[point.node].left + point.offset] != _text[left]) {
return {null_node, 0};
}
int remaining = edge_length_unchecked(point.node) - point.offset;
if (right - left < remaining) {
point.offset += right - left;
return point;
}
left += remaining;
point.offset = edge_length_unchecked(point.node);
}
}
return point;
}
node_id split(ActivePoint point) {
if (point.offset == edge_length_unchecked(point.node)) return point.node;
if (point.offset == 0) return _nodes[point.node].parent;
node_id child = point.node;
node_id parent = _nodes[child].parent;
int left = _nodes[child].left;
node_id middle = new_node(left, left + point.offset, parent);
_nodes[parent].next[_text[left]] = middle;
_nodes[middle].next[_text[left + point.offset]] = child;
_nodes[child].parent = middle;
_nodes[child].left += point.offset;
return middle;
}
node_id get_suffix_link(node_id id) {
if (_nodes[id].suffix_link != null_node) return _nodes[id].suffix_link;
node_id parent = _nodes[id].parent;
if (parent == null_node) return root_node;
node_id parent_link = get_suffix_link(parent);
ActivePoint point = {
parent_link,
edge_length_unchecked(parent_link)
};
int left = _nodes[id].left + (parent == root_node);
point = go(point, left, _nodes[id].right);
assert(point.node != null_node);
return _nodes[id].suffix_link = split(point);
}
void extend(int position) {
while (true) {
ActivePoint next = go(_active, position, position + 1);
if (next.node != null_node) {
_active = next;
return;
}
node_id middle = split(_active);
node_id leaf = new_node(position, int(_text.size()), middle);
_nodes[middle].next[_text[position]] = leaf;
_active.node = get_suffix_link(middle);
_active.offset = edge_length_unchecked(_active.node);
if (middle == root_node) return;
}
}
void finish_metadata() {
std::vector<node_id> order;
order.reserve(_nodes.size());
order.push_back(root_node);
std::vector<int> depth(_nodes.size(), 0);
for (std::size_t i = 0; i < order.size(); i++) {
node_id id = order[i];
node_id previous_child = null_node;
for (int symbol = 0; symbol <= terminal_symbol; symbol++) {
node_id child = _nodes[id].next[symbol];
if (child == null_node) continue;
_nodes[child].incoming_symbol = symbol;
if (previous_child == null_node) {
_nodes[id].first_child = child;
} else {
_nodes[previous_child].next_sibling = child;
}
previous_child = child;
_nodes[id].child_count++;
depth[child] = depth[id] + edge_length_unchecked(child);
order.push_back(child);
}
}
for (int i = int(order.size()) - 1; i >= 0; i--) {
node_id id = order[i];
bool leaf = true;
for (node_id child : _nodes[id].next) {
if (child == null_node) continue;
leaf = false;
_nodes[id].leaf_count += _nodes[child].leaf_count;
if (_nodes[id].representative_suffix == -1) {
_nodes[id].representative_suffix = _nodes[child].representative_suffix;
}
}
if (leaf) {
_nodes[id].suffix_start = int(_text.size()) - depth[id];
_nodes[id].representative_suffix = _nodes[id].suffix_start;
_nodes[id].leaf_count = 1;
}
}
}
void initialize() {
_nodes.clear();
_nodes.reserve(2 * _text.size() + 1);
_nodes.emplace_back();
_nodes[root_node].suffix_link = root_node;
_active = {root_node, 0};
for (int position = 0; position < int(_text.size()); position++) extend(position);
finish_metadata();
}
public:
SuffixTree() {
clear();
}
template <class Sequence>
explicit SuffixTree(const Sequence& sequence) {
build(sequence);
}
int size() const {
return node_count();
}
bool empty() const {
return _text_length == 0;
}
int node_count() const {
return int(_nodes.size());
}
int text_length() const {
return _text_length;
}
node_id root() const {
return root_node;
}
const Node& node(node_id id) const {
assert(0 <= id && id < node_count());
return _nodes[id];
}
const std::vector<Node>& nodes() const {
return _nodes;
}
int edge_length(node_id id) const {
assert(0 <= id && id < node_count());
return edge_length_unchecked(id);
}
bool is_leaf(node_id id) const {
assert(0 <= id && id < node_count());
return _nodes[id].suffix_start != -1;
}
template <class Symbol>
node_id child(node_id id, const Symbol& symbol) const {
assert(0 <= id && id < node_count());
return _nodes[id].next[symbol_index(symbol)];
}
node_id child_by_index(node_id id, int symbol) const {
assert(0 <= id && id < node_count());
assert(0 <= symbol && symbol <= terminal_symbol);
return _nodes[id].next[symbol];
}
template <class Callback>
void for_each_child(node_id id, Callback callback) const {
assert(0 <= id && id < node_count());
for (
node_id child_id = _nodes[id].first_child;
child_id != null_node;
child_id = _nodes[child_id].next_sibling
) {
callback(_nodes[child_id].incoming_symbol, child_id);
}
}
void clear() {
_text.clear();
_text.push_back(terminal_symbol);
_text_length = 0;
initialize();
}
template <class Sequence>
void build(const Sequence& sequence) {
_text.clear();
for (const auto& symbol : sequence) _text.push_back(symbol_index(symbol));
assert(_text.size() < std::size_t(std::numeric_limits<int>::max()));
_text_length = int(_text.size());
_text.push_back(terminal_symbol);
initialize();
}
template <class Sequence>
Locus find(const Sequence& sequence) const {
ActivePoint point = {root_node, 0};
for (const auto& value : sequence) {
int symbol = symbol_index(value);
if (point.offset == edge_length_unchecked(point.node)) {
point = {_nodes[point.node].next[symbol], 0};
if (point.node == null_node) return {null_node, 0};
}
if (_text[_nodes[point.node].left + point.offset] != symbol) {
return {null_node, 0};
}
point.offset++;
}
return {point.node, point.offset};
}
template <class Sequence>
bool contains(const Sequence& sequence) const {
return bool(find(sequence));
}
template <class Sequence>
int count_occurrences(const Sequence& sequence) const {
Locus locus = find(sequence);
return locus ? _nodes[locus.node].leaf_count : 0;
}
template <class Sequence>
std::pair<int, int> representative_occurrence(const Sequence& sequence) const {
Locus locus = {root_node, 0};
int length = 0;
for (const auto& value : sequence) {
int symbol = symbol_index(value);
if (locus.offset == edge_length_unchecked(locus.node)) {
locus = {_nodes[locus.node].next[symbol], 0};
if (locus.node == null_node) return {-1, -1};
}
if (_text[_nodes[locus.node].left + locus.offset] != symbol) return {-1, -1};
locus.offset++;
length++;
}
int left = _nodes[locus.node].representative_suffix;
return {left, left + length};
}
long long distinct_substring_count() const {
long long result = 0;
for (node_id id = 1; id < node_count(); id++) {
result += std::max(0, std::min(_nodes[id].right, _text_length) - _nodes[id].left);
}
return result;
}
};
} // namespace string
} // namespace m1une
#line 1 "string/trie.hpp"
#line 9 "string/trie.hpp"
namespace m1une {
namespace string {
// A multiset trie for a contiguous character alphabet.
template <int AlphabetSize = 26, int FirstCharacter = 'a'>
struct Trie {
static_assert(0 < AlphabetSize);
using node_id = int;
static constexpr node_id null_node = -1;
struct Node {
std::array<node_id, AlphabetSize> child;
int subtree_count;
int terminal_count;
Node() : subtree_count(0), terminal_count(0) {
child.fill(null_node);
}
};
private:
std::vector<Node> _nodes;
int _distinct_size;
template <class Symbol>
static int symbol_index(const Symbol& symbol) {
int index = int(symbol) - FirstCharacter;
assert(0 <= index && index < AlphabetSize);
return index;
}
node_id new_node() {
assert(_nodes.size() < std::size_t(std::numeric_limits<int>::max()));
_nodes.emplace_back();
return int(_nodes.size()) - 1;
}
template <class Sequence>
node_id find_node(const Sequence& sequence) const {
node_id node = 0;
for (const auto& symbol : sequence) {
node = _nodes[node].child[symbol_index(symbol)];
if (node == null_node || _nodes[node].subtree_count == 0) {
return null_node;
}
}
return node;
}
public:
Trie() : _nodes(1), _distinct_size(0) {}
int size() const {
return _nodes[0].subtree_count;
}
int distinct_size() const {
return _distinct_size;
}
bool empty() const {
return size() == 0;
}
node_id root() const {
return 0;
}
const Node& node(node_id id) const {
assert(0 <= id && std::size_t(id) < _nodes.size());
return _nodes[id];
}
template <class Sequence>
node_id find(const Sequence& sequence) const {
return find_node(sequence);
}
std::size_t node_count() const {
return _nodes.size();
}
void reserve(std::size_t node_capacity) {
_nodes.reserve(node_capacity);
}
void clear() {
_nodes.clear();
_nodes.emplace_back();
_distinct_size = 0;
}
template <class Sequence>
node_id insert(const Sequence& sequence, int multiplicity = 1) {
assert(0 < multiplicity);
node_id node = 0;
_nodes[node].subtree_count += multiplicity;
for (const auto& symbol : sequence) {
int index = symbol_index(symbol);
node_id child = _nodes[node].child[index];
if (child == null_node) {
child = new_node();
_nodes[node].child[index] = child;
}
node = child;
_nodes[node].subtree_count += multiplicity;
}
if (_nodes[node].terminal_count == 0) _distinct_size++;
_nodes[node].terminal_count += multiplicity;
return node;
}
template <class Sequence>
int count(const Sequence& sequence) const {
node_id node = find_node(sequence);
return node == null_node ? 0 : _nodes[node].terminal_count;
}
template <class Sequence>
bool contains(const Sequence& sequence) const {
return count(sequence) != 0;
}
// Returns the number of stored strings beginning with prefix.
template <class Sequence>
int prefix_count(const Sequence& prefix) const {
node_id node = find_node(prefix);
return node == null_node ? 0 : _nodes[node].subtree_count;
}
template <class Sequence>
bool starts_with(const Sequence& prefix) const {
return prefix_count(prefix) != 0;
}
template <class Sequence>
bool erase_one(const Sequence& sequence) {
node_id terminal = find_node(sequence);
if (terminal == null_node || _nodes[terminal].terminal_count == 0) {
return false;
}
int node = 0;
_nodes[node].subtree_count--;
for (const auto& symbol : sequence) {
node = _nodes[node].child[symbol_index(symbol)];
_nodes[node].subtree_count--;
}
_nodes[node].terminal_count--;
if (_nodes[node].terminal_count == 0) _distinct_size--;
return true;
}
template <class Sequence>
bool erase(const Sequence& sequence) {
return erase_one(sequence);
}
template <class Sequence>
int erase_all(const Sequence& sequence) {
int multiplicity = count(sequence);
if (multiplicity == 0) return 0;
int node = 0;
_nodes[node].subtree_count -= multiplicity;
for (const auto& symbol : sequence) {
node = _nodes[node].child[symbol_index(symbol)];
_nodes[node].subtree_count -= multiplicity;
}
_nodes[node].terminal_count = 0;
_distinct_size--;
return multiplicity;
}
// Calls callback(length, multiplicity) for every stored prefix.
// The empty prefix is reported with length 0 when it is stored.
template <class Sequence, class Callback>
void for_each_prefix(const Sequence& sequence, Callback callback) const {
int node = 0;
if (_nodes[node].terminal_count != 0) {
callback(0, _nodes[node].terminal_count);
}
int length = 0;
for (const auto& symbol : sequence) {
node = _nodes[node].child[symbol_index(symbol)];
if (node == null_node || _nodes[node].subtree_count == 0) return;
length++;
if (_nodes[node].terminal_count != 0) {
callback(length, _nodes[node].terminal_count);
}
}
}
// Returns the length of the longest stored string that is a prefix.
// Returns -1 when no stored prefix exists.
template <class Sequence>
int longest_prefix(const Sequence& sequence) const {
int result = _nodes[0].terminal_count == 0 ? -1 : 0;
for_each_prefix(sequence, [&result](int length, int) {
result = length;
});
return result;
}
};
} // namespace string
} // namespace m1une
#line 1 "string/wildcard_pattern_matching.hpp"
#line 7 "string/wildcard_pattern_matching.hpp"
#line 1 "math/fps/convolution.hpp"
#line 8 "math/fps/convolution.hpp"
#include <cstring>
#include <new>
#line 13 "math/fps/convolution.hpp"
#if defined(__GNUC__) && !defined(__clang__) && \
(defined(__x86_64__) || defined(__i386__)) && \
!defined(M1UNE_FPS_DISABLE_X86_SIMD)
#include <immintrin.h>
#define M1UNE_FPS_HAS_X86_SIMD 1
#pragma GCC push_options
#pragma GCC target("avx2,bmi")
#endif
#line 1 "math/fps/internal/ntt998_faster.hpp"
#ifdef M1UNE_FPS_HAS_X86_SIMD
#line 9 "math/fps/internal/ntt998_faster.hpp"
#include <immintrin.h>
namespace m1une {
namespace fps {
namespace internal {
namespace fast998_v2 {
// Fixed-modulus AVX2 transform with an in-register degree-8 residue product.
using u32=unsigned;
using u64=unsigned long long;
using idt=std::size_t;
using I256=__m256i;
inline void store256(void*p,I256 x){
_mm256_store_si256((I256*)p,x);
}
inline I256 load256(const void*p){
return _mm256_load_si256((const I256*)p);
}
constexpr u32 shrk(u32 x,u32 M){
return std::min(x,x-M);
}
constexpr u32 dilt(u32 x,u32 M){
return std::min(x,x+M);
}
constexpr u32 reduce(u64 x,u32 niv,u32 M){
return (x+u64(u32(x)*niv)*M)>>32;
}
constexpr u32 mul(u32 x,u32 y,u32 niv,u32 M){
return reduce(u64(x)*y,niv,M);
}
constexpr u32 mul_s(u32 x,u32 y,u32 niv,u32 M){
return shrk(reduce(u64(x)*y,niv,M),M);
}
constexpr u32 qpw(u32 a,u32 b,u32 niv,u32 M,u32 r){
for(;b;b>>=1,a=mul(a,a,niv,M)){
if(b&1){
r=mul(r,a,niv,M);
}
}
return r;
}
constexpr u32 qpw_s(u32 a,u32 b,u32 niv,u32 M,u32 r){
return shrk(qpw(a,b,niv,M,r),M);
}
inline I256 shrk32(I256 x,I256 M){
return _mm256_min_epu32(x,_mm256_sub_epi32(x,M));
}
inline I256 dilt32(I256 x,I256 M){
return _mm256_min_epu32(x,_mm256_add_epi32(x,M));
}
inline I256 Ladd32(I256 x,I256 y,I256){
return _mm256_add_epi32(x,y);
}
inline I256 Lsub32(I256 x,I256 y,I256 M){
return _mm256_add_epi32(_mm256_sub_epi32(x,y),M);
}
inline I256 add32(I256 x,I256 y,I256 M){
return shrk32(_mm256_add_epi32(x,y),M);
}
inline I256 sub32(I256 x,I256 y,I256 M){
return dilt32(_mm256_sub_epi32(x,y),M);
}
template<int msk>inline I256 neg32_m(I256 x,I256 M){
return _mm256_blend_epi32(x,_mm256_sub_epi32(M,x),msk);
}
inline I256 reduce(I256 a,I256 b,I256 niv,I256 M){
I256 c=_mm256_mul_epu32(a,niv),d=_mm256_mul_epu32(b,niv);
c=_mm256_mul_epu32(c,M),d=_mm256_mul_epu32(d,M);
return _mm256_blend_epi32(_mm256_srli_epi64(_mm256_add_epi64(a,c),32),_mm256_add_epi64(b,d),0xaa);
}
inline I256 mul(I256 a,I256 b,I256 niv,I256 M){
return reduce(_mm256_mul_epu32(a,b),_mm256_mul_epu32(_mm256_srli_epi64(a,32),_mm256_srli_epi64(b,32)),niv,M);
}
inline I256 mul_s(I256 a,I256 b,I256 niv,I256 M){
return shrk32(mul(a,b,niv,M),M);
}
inline I256 mul_bsm(I256 a,I256 b,I256 niv,I256 M){
return reduce(_mm256_mul_epu32(a,b),_mm256_mul_epu32(_mm256_srli_epi64(a,32),b),niv,M);
}
inline I256 mul_bsmfxd(I256 a,I256 b,I256 bniv,I256 M){
I256 cc=_mm256_mul_epu32(a,bniv),dd=_mm256_mul_epu32(_mm256_srli_epi64(a,32),bniv);
I256 c=_mm256_mul_epu32(a,b),d=_mm256_mul_epu32(_mm256_srli_epi64(a,32),b);
cc=_mm256_mul_epu32(cc,M),dd=_mm256_mul_epu32(dd,M);
return _mm256_blend_epi32(_mm256_srli_epi64(_mm256_add_epi64(c,cc),32),_mm256_add_epi64(d,dd),0xaa);
}
inline I256 mul_bfxd(I256 a,I256 b,I256 bniv,I256 M){
I256 cc=_mm256_mul_epu32(a,bniv),dd=_mm256_mul_epu32(_mm256_srli_epi64(a,32),_mm256_srli_epi64(bniv,32));
I256 c=_mm256_mul_epu32(a,b),d=_mm256_mul_epu32(_mm256_srli_epi64(a,32),_mm256_srli_epi64(b,32));
cc=_mm256_mul_epu32(cc,M),dd=_mm256_mul_epu32(dd,M);
return _mm256_blend_epi32(_mm256_srli_epi64(_mm256_add_epi64(c,cc),32),_mm256_add_epi64(d,dd),0xaa);
}
inline I256 mul_upd_rt(I256 a,I256 bu,I256 M){
I256 cc=_mm256_mul_epu32(a,bu),c=_mm256_mul_epu32(a,_mm256_srli_epi64(bu,32));
cc=_mm256_mul_epu32(cc,M);
return shrk32(_mm256_srli_epi64(_mm256_add_epi64(c,cc),32),M);
}
constexpr auto _mxlg=26,_lg_itth=6;
constexpr auto _itth=idt(1)<<_lg_itth;
static_assert(_lg_itth%2==0);
struct FNTT32_info{
u32 mod,mod2,niv,one,r2,r3,img,imgniv,RT1[_mxlg];
alignas(32) std::array<u32,8> rt3[_mxlg-2],rt3i[_mxlg-2],bwbr,bwb,bwbi,rt4[_mxlg-3],rt4niv[_mxlg-3],rt4i[_mxlg-3],rt4iniv[_mxlg-3],pr2,pr4,pr2niv,pr4niv,pr2i,pr2iniv,pr4i,pr4iniv;
constexpr FNTT32_info(const u32 m):mod(m),mod2(m*2),niv([&]{u32 n=2+m;for(int i=0;i<4;++i){n*=2+m*n;}return n;}()),one((-m)%m),r2((-u64(m))%m),r3(mul_s(r2,r2,niv,m)),img{},imgniv{},RT1{},rt3{},rt3i{},bwbr{},bwb{},bwbi{},rt4{},rt4niv{},rt4i{},rt4iniv{},pr2{},pr4{},pr2niv{},pr4niv{},pr2i{},pr2iniv{},pr4i{},pr4iniv{}{
const int k=__builtin_ctz(m-1);
u32 _g=mul(3,r2,niv,mod);
for(;;++_g){
if(qpw_s(_g,mod>>1,niv,mod,one)!=one){
break;
}
}
_g=qpw(_g,mod>>k,niv,mod,one);
u32 rt1[_mxlg-1],rt1i[_mxlg-1];
rt1[k-2]=_g,rt1i[k-2]=qpw(_g,mod-2,niv,mod,one);
for(int i=k-2;i>0;--i){
rt1[i-1]=mul(rt1[i],rt1[i],niv,mod);
rt1i[i-1]=mul(rt1i[i],rt1i[i],niv,mod);
}
RT1[k-1]=qpw_s(_g,3,niv,mod,one);
for(int i=k-1;i>0;--i){
RT1[i-1]=mul_s(RT1[i],RT1[i],niv,mod);
}
img=rt1[0],imgniv=img*niv;
bwbr={one,0,one,0,one};
bwb={rt1[1],0,rt1[0],0,mod-mul_s(rt1[0],rt1[1],niv,mod)};
bwbi={rt1i[1],0,rt1i[0],0,mul_s(rt1i[0],rt1i[1],niv,mod)};
u32 pr=one,pri=one;
for(int i=0;i<k-2;++i){
const u32 r=mul_s(pr,rt1[i+1],niv,mod),ri=mul_s(pri,rt1i[i+1],niv,mod);
const u32 r2=mul_s(r,r,niv,mod),r2i=mul_s(ri,ri,niv,mod);
const u32 r3=mul_s(r,r2,niv,mod),r3i=mul_s(ri,r2i,niv,mod);
rt3[i]={r*niv,r,r2*niv,r2,r3*niv,r3};
rt3i[i]={ri*niv,ri,r2i*niv,r2i,r3i*niv,r3i};
pr=mul(pr,rt1i[i+1],niv,mod),pri=mul(pri,rt1[i+1],niv,mod);
}
pr=one,pri=one;
for(int i=0;i<k-3;++i){
const u32 r=mul_s(pr,rt1[i+2],niv,mod),ri=mul_s(pri,rt1i[i+2],niv,mod);
rt4[i][0]=rt4i[i][0]=one;
for(int j=1;j<8;++j){
rt4[i][j]=mul_s(rt4[i][j-1],r,niv,mod);
rt4i[i][j]=mul_s(rt4i[i][j-1],ri,niv,mod);
}
for(int j=0;j<8;++j){
rt4niv[i][j]=rt4[i][j]*niv;
rt4iniv[i][j]=rt4i[i][j]*niv;
}
pr=mul(pr,rt1i[i+2],niv,mod),pri=mul(pri,rt1[i+2],niv,mod);
}
pr2={one,one,one,img,one,one,one,img};
pr4={one,one,one,one,one,rt1[1],img,mul_s(img,rt1[1],niv,mod)};
const u32 nr2=mod-r2,imgr2=mul_s(img,r2,niv,mod);
pr2i={nr2,nr2,nr2,imgr2,nr2,nr2,nr2,imgr2};
pr4i={one,one,one,one,one,rt1i[1],rt1i[0],mul_s(rt1i[0],rt1i[1],niv,mod)};
for(int j=0;j<8;++j){
pr2niv[j]=pr2[j]*niv,pr4niv[j]=pr4[j]*niv;
pr2iniv[j]=pr2i[j]*niv,pr4iniv[j]=pr4i[j]*niv;
}
}
};
inline void vector_dif(I256*const f,const idt n,const FNTT32_info*info){
alignas(32) std::array<u32,8> st_1[_mxlg>>1];
const I256 Mod=_mm256_set1_epi32(info->mod),Mod2=_mm256_set1_epi32(info->mod2),Niv=_mm256_set1_epi32(info->niv);
const I256 Img=_mm256_set1_epi32(info->img),ImgNiv=_mm256_set1_epi32(info->imgniv),id=_mm256_setr_epi32(0,2,0,4,0,2,0,4);
const int lgn=__builtin_ctzll(n);
std::fill(st_1,st_1+(lgn>>1),info->bwb);
const idt nn=n>>(lgn&1),m=std::min(n,_itth),mm=std::min(nn,_itth);
// I256 rr=_mm256_set1_epi32(info->one);
if(nn!=n){
for(idt i=0;i<nn;++i){
auto const p0=f+i,p1=f+nn+i;
const auto f0=load256(p0),f1=load256(p1);
const auto g0=add32(f0,f1,Mod2),g1=Lsub32(f0,f1,Mod2);
store256(p0,g0),store256(p1,g1);
}
}
for(idt L=nn>>2;L>0;L>>=2){
for(idt i=0;i<L;++i){
auto const p0=f+i,p1=p0+L,p2=p1+L,p3=p2+L;
const auto f1=load256(p1),f3=load256(p3),f2=load256(p2),f0=load256(p0);
const auto g3=mul_bsmfxd(Lsub32(f1,f3,Mod2),Img,ImgNiv,Mod),g1=add32(f1,f3,Mod2);
const auto g0=add32(f0,f2,Mod2),g2=sub32(f0,f2,Mod2);
const auto h0=add32(g0,g1,Mod2),h1=Lsub32(g0,g1,Mod2);
const auto h2=Ladd32(g2,g3,Mod2),h3=Lsub32(g2,g3,Mod2);
store256(p0,h0),store256(p1,h1),store256(p2,h2),store256(p3,h3);
}
}
for(idt j=0;j<n;j+=m){
int t=((j==0)?std::min(_lg_itth,lgn):__builtin_ctzll(j))&-2,p=(t-2)>>1;
for(idt L=(idt(1)<<t)>>2;L>=_itth;L>>=2,t-=2,--p){
auto rt=load256(st_1+p);
const auto r1=_mm256_permutevar8x32_epi32(rt,id);
const auto r1Niv=_mm256_permutevar8x32_epi32(_mm256_mul_epu32(rt,Niv),id);
rt=mul_upd_rt(rt,load256(info->rt3+__builtin_ctzll(~j>>t)),Mod);
const auto r2=_mm256_shuffle_epi32(r1,_MM_PERM_BBBB),nr3=_mm256_shuffle_epi32(r1,_MM_PERM_DDDD);
const auto r2Niv=_mm256_shuffle_epi32(r1Niv,_MM_PERM_BBBB),nr3Niv=_mm256_shuffle_epi32(r1Niv,_MM_PERM_DDDD);
store256(st_1+p,rt);
for(idt i=0;i<L;++i){
auto const p0=f+i+j,p1=p0+L,p2=p1+L,p3=p2+L;
const auto f1=load256(p1),f3=load256(p3),f2=load256(p2),f0=load256(p0);
const auto g1=mul_bsmfxd(f1,r1,r1Niv,Mod),ng3=mul_bsmfxd(f3,nr3,nr3Niv,Mod);
const auto g2=mul_bsmfxd(f2,r2,r2Niv,Mod),g0=shrk32(f0,Mod2);
const auto h3=mul_bsmfxd(Ladd32(g1,ng3,Mod2),Img,ImgNiv,Mod),h1=sub32(g1,ng3,Mod2);
const auto h0=add32(g0,g2,Mod2),h2=sub32(g0,g2,Mod2);
const auto u0=Ladd32(h0,h1,Mod2),u1=Lsub32(h0,h1,Mod2);
const auto u2=Ladd32(h2,h3,Mod2),u3=Lsub32(h2,h3,Mod2);
store256(p0,u0),store256(p1,u1),store256(p2,u2),store256(p3,u3);
}
}
I256*const g=f+j;
for(idt l=mm,L=mm>>2;L;l=L,L>>=2,t-=2,--p){
auto rt=load256(st_1+p);
for(idt i=(j==0?l:0),k=(j+i)>>t;i<m;i+=l,++k){
const auto r1=_mm256_permutevar8x32_epi32(rt,id);
const auto r2=_mm256_shuffle_epi32(r1,_MM_PERM_BBBB);
const auto nr3=_mm256_shuffle_epi32(r1,_MM_PERM_DDDD);
for(idt j=0;j<L;++j){
auto const p0=g+i+j,p1=p0+L,p2=p1+L,p3=p2+L;
const auto f1=load256(p1),f3=load256(p3),f2=load256(p2),f0=load256(p0);
const auto g1=mul_bsm(f1,r1,Niv,Mod),ng3=mul_bsm(f3,nr3,Niv,Mod);
const auto g2=mul_bsm(f2,r2,Niv,Mod),g0=shrk32(f0,Mod2);
const auto h3=mul_bsmfxd(Ladd32(g1,ng3,Mod2),Img,ImgNiv,Mod),h1=sub32(g1,ng3,Mod2);
const auto h0=add32(g0,g2,Mod2),h2=sub32(g0,g2,Mod2);
const auto u0=Ladd32(h0,h1,Mod2),u1=Lsub32(h0,h1,Mod2);
const auto u2=Ladd32(h2,h3,Mod2),u3=Lsub32(h2,h3,Mod2);
store256(p0,u0),store256(p1,u1),store256(p2,u2),store256(p3,u3);
}
rt=mul_upd_rt(rt,load256(info->rt3+__builtin_ctzll(~k)),Mod);
}
store256(st_1+p,rt);
}
// const auto pr2=load256(&info->pr2),pr4=load256(&info->pr4);
// const auto pr2Niv=load256(&info->pr2niv),pr4Niv=load256(&info->pr4niv);
// for(idt i=j;i<j+m;++i){
// auto fi=load256(f+i);
// fi=mul(fi,rr,Niv,Mod);
// rr=shrk32(mul_bfxd(rr,load256(info->rt4+__builtin_ctzll(~i)),load256(info->rt4niv+__builtin_ctzll(~i)),Mod),Mod);
// fi=mul_bfxd(Ladd32(neg32_m<0xf0>(fi,Mod2),_mm256_permute2x128_si256(fi,fi,1),Mod2),pr4,pr4Niv,Mod);
// fi=mul_bfxd(Ladd32(neg32_m<0xcc>(fi,Mod2),_mm256_shuffle_epi32(fi,0x4e),Mod2),pr2,pr2Niv,Mod);
// fi=sub32(_mm256_shuffle_epi32(fi,0xb1),neg32_m<0x55>(fi,Mod2),Mod2);
// store256(f+i,fi);
// }
}
}
template<bool shrk=false>inline void vector_dit(I256*const f,idt n,const FNTT32_info*const info){
alignas(32) std::array<u32,8> st_1[_mxlg>>1];
const I256 Mod=_mm256_set1_epi32(info->mod),Mod2=_mm256_set1_epi32(info->mod2),Niv=_mm256_set1_epi32(info->niv);
const I256 Img=_mm256_set1_epi32(info->img),ImgNiv=_mm256_set1_epi32(info->imgniv),id=_mm256_setr_epi32(0,2,0,4,0,2,0,4);
const int lgn=__builtin_ctzll(n);
std::fill(st_1,st_1+(_lg_itth>>1),info->bwbr);
std::fill(st_1+(_lg_itth>>1),st_1+(_mxlg>>1),info->bwbi);
const idt nn=n>>(lgn&1),mm=std::min(nn,_itth);
// I256 rr=_mm256_set1_epi32((info->mod-1)>>(lgn+3));
for(idt j=0;j<n;j+=mm){
// const auto pr2=load256(&info->pr2i),pr4=load256(&info->pr4i);
// const auto pr2Niv=load256(&info->pr2iniv),pr4Niv=load256(&info->pr4iniv);
// for(idt i=j;i<j+mm;++i){
// auto fi=load256(f+i);
// const auto rt=rr;
// rr=shrk32(mul_bfxd(rr,load256(info->rt4i+__builtin_ctzll(~i)),load256(info->rt4iniv+__builtin_ctzll(~i)),Mod),Mod);
// fi=mul_bfxd(Ladd32(neg32_m<0xaa>(fi,Mod2),_mm256_shuffle_epi32(fi,0xb1),Mod2),pr2,pr2Niv,Mod);
// fi=mul_bfxd(Ladd32(neg32_m<0xcc>(fi,Mod2),_mm256_shuffle_epi32(fi,0x4e),Mod2),pr4,pr4Niv,Mod);
// fi=mul(Ladd32(neg32_m<0xf0>(fi,Mod2),_mm256_permute2x128_si256(fi,fi,1),Mod2),rt,Niv,Mod);
// store256(f+i,fi);
// }
I256*const g=f+j;
int t=2,p=0;
for(idt l=4,L=1;l<=mm;L=l,l<<=2,t+=2,++p){
auto rt=load256(st_1+p);
for(idt i=0,k=j>>t;i<mm;i+=l,++k){
const auto r1=_mm256_permutevar8x32_epi32(rt,id);
const auto r2=_mm256_shuffle_epi32(r1,_MM_PERM_BBBB);
const auto r3=_mm256_shuffle_epi32(r1,_MM_PERM_DDDD);
for(idt j=0;j<L;++j){
auto const p0=g+i+j,p1=p0+L,p2=p1+L,p3=p2+L;
const auto f0=load256(p0),f1=load256(p1),f2=load256(p2),f3=load256(p3);
const auto g0=add32(f0,f1,Mod2),g1=sub32(f0,f1,Mod2);
const auto g2=add32(f2,f3,Mod2),g3=mul_bsmfxd(Lsub32(f3,f2,Mod2),Img,ImgNiv,Mod);
const auto h0=Ladd32(g0,g2,Mod2),h1=Ladd32(g1,g3,Mod2);
const auto h2=Lsub32(g0,g2,Mod2),h3=Lsub32(g1,g3,Mod2);
const auto u0=shrk32(h0,Mod2),u1=mul_bsm(h1,r1,Niv,Mod);
const auto u2=mul_bsm(h2,r2,Niv,Mod),u3=mul_bsm(h3,r3,Niv,Mod);
store256(p0,u0),store256(p1,u1),store256(p2,u2),store256(p3,u3);
}
rt=mul_upd_rt(rt,load256(info->rt3i+__builtin_ctzll(~k)),Mod);
}
store256(st_1+p,rt);
}
int tt=std::min(__builtin_ctzll(~(j>>_lg_itth))+_lg_itth,lgn);
for(idt L=_itth,l=L<<2;t<=tt;L=l,l<<=2,t+=2,++p){
if((j+_itth)==l){
if(shrk && l==n){
for(idt i=0;i<L;++i){
auto const p0=f+i,p1=p0+L,p2=p1+L,p3=p2+L;
const auto f2=load256(p2),f3=load256(p3),f0=load256(p0),f1=load256(p1);
const auto g3=mul_bsmfxd(Lsub32(f3,f2,Mod2),Img,ImgNiv,Mod),g2=add32(f2,f3,Mod2);
const auto g0=add32(f0,f1,Mod2),g1=sub32(f0,f1,Mod2);
const auto h0=add32(g0,g2,Mod2),h1=add32(g1,g3,Mod2);
const auto h2=sub32(g0,g2,Mod2),h3=sub32(g1,g3,Mod2);
const auto u0=shrk32(h0,Mod),u1=shrk32(h1,Mod);
const auto u2=shrk32(h2,Mod),u3=shrk32(h3,Mod);
store256(p0,u0),store256(p1,u1),store256(p2,u2),store256(p3,u3);
}
}
else{
for(idt i=0;i<L;++i){
auto const p0=f+i,p1=p0+L,p2=p1+L,p3=p2+L;
const auto f2=load256(p2),f3=load256(p3),f0=load256(p0),f1=load256(p1);
const auto g3=mul_bsmfxd(Lsub32(f3,f2,Mod2),Img,ImgNiv,Mod),g2=add32(f2,f3,Mod2);
const auto g0=add32(f0,f1,Mod2),g1=sub32(f0,f1,Mod2);
const auto h0=add32(g0,g2,Mod2),h1=add32(g1,g3,Mod2);
const auto h2=sub32(g0,g2,Mod2),h3=sub32(g1,g3,Mod2);
store256(p0,h0),store256(p1,h1),store256(p2,h2),store256(p3,h3);
}
}
}
else{
auto rt=load256(st_1+p);
const auto r1=_mm256_permutevar8x32_epi32(rt,id);
const auto r1Niv=_mm256_permutevar8x32_epi32(_mm256_mul_epu32(rt,Niv),id);
rt=mul_upd_rt(rt,load256(info->rt3i+__builtin_ctzll(~j>>t)),Mod);
const auto r2=_mm256_shuffle_epi32(r1,_MM_PERM_BBBB),r3=_mm256_shuffle_epi32(r1,_MM_PERM_DDDD);
const auto r2Niv=_mm256_shuffle_epi32(r1Niv,_MM_PERM_BBBB),r3Niv=_mm256_shuffle_epi32(r1Niv,_MM_PERM_DDDD);
store256(st_1+p,rt);
for(idt i=0;i<L;++i){
auto const p0=f+j+_itth-l+i,p1=p0+L,p2=p1+L,p3=p2+L;
const auto f0=load256(p0),f1=load256(p1),f2=load256(p2),f3=load256(p3);
const auto g0=add32(f0,f1,Mod2),g1=sub32(f0,f1,Mod2);
const auto g2=add32(f2,f3,Mod2),g3=mul_bsmfxd(Lsub32(f3,f2,Mod2),Img,ImgNiv,Mod);
const auto h0=Ladd32(g0,g2,Mod2),h1=Ladd32(g1,g3,Mod2);
const auto h2=Lsub32(g0,g2,Mod2),h3=Lsub32(g1,g3,Mod2);
const auto u0=shrk32(h0,Mod2),u1=mul_bsmfxd(h1,r1,r1Niv,Mod);
const auto u2=mul_bsmfxd(h2,r2,r2Niv,Mod),u3=mul_bsmfxd(h3,r3,r3Niv,Mod);
store256(p0,u0),store256(p1,u1),store256(p2,u2),store256(p3,u3);
}
}
}
}
if(shrk && nn==n && n<=_itth){
for(idt i=0;i<n;++i){
const auto f0=load256(f+i);
store256(f+i,shrk32(f0,Mod));
}
}
if(nn!=n){
for(idt i=0;i<nn;++i){
auto const p0=f+i,p1=f+nn+i;
const auto f0=load256(p0),f1=load256(p1);
const auto g0=add32(f0,f1,Mod2),g1=sub32(f0,f1,Mod2);
if constexpr(shrk){
const auto h0=shrk32(g0,Mod),h1=shrk32(g1,Mod);
store256(p0,h0),store256(p1,h1);
}
else{
store256(p0,g0),store256(p1,g1);
}
}
}
}
// Returns fx * f[0,8) * g[0,8) (mod x^8 - ww).
[[gnu::always_inline]] inline I256 convolve8(const I256*f,const I256*g,I256 ww,I256 fx,I256 Niv,I256 Mod,I256 Mod2){
const auto raa=load256(f),rbb=load256(g);
const auto taa=shrk32(raa,Mod2),bb=shrk32(mul_bsm(rbb,fx,Niv,Mod),Mod);
const auto aw=shrk32(mul_bsm(taa,ww,Niv,Mod),Mod);
const auto aa=shrk32(taa,Mod);
const auto awa=_mm256_permute2x128_si256(aa,aw,3);
const auto b0=_mm256_permute4x64_epi64(bb,0x00),b1=_mm256_shuffle_epi32(b0,_MM_PERM_CDAB);
const auto a0=aa,a1=_mm256_srli_epi64(a0,32);
const auto aw7=_mm256_alignr_epi8(aa,awa,12);
auto res00=_mm256_mul_epu32(a0,b0);
auto res01=_mm256_mul_epu32(a1,b0);
auto res10=_mm256_mul_epu32(aw7,b1);
auto res11=_mm256_mul_epu32(a0,b1);
const auto b2=_mm256_permute4x64_epi64(bb,0x55),b3=_mm256_shuffle_epi32(b2,_MM_PERM_CDAB);
const auto aw6=_mm256_alignr_epi8(aa,awa,8);
const auto aw5=_mm256_alignr_epi8(aa,awa,4);
res00=_mm256_add_epi64(res00,_mm256_mul_epu32(aw6,b2));
res01=_mm256_add_epi64(res01,_mm256_mul_epu32(aw7,b2));
res10=_mm256_add_epi64(res10,_mm256_mul_epu32(aw5,b3));
res11=_mm256_add_epi64(res11,_mm256_mul_epu32(aw6,b3));
const auto b4=_mm256_permute4x64_epi64(bb,0xaa),b5=_mm256_shuffle_epi32(b4,_MM_PERM_CDAB);
const auto aw3=_mm256_alignr_epi8(awa,aw,12);
res00=_mm256_add_epi64(res00,_mm256_mul_epu32(awa,b4));
res01=_mm256_add_epi64(res01,_mm256_mul_epu32(aw5,b4));
res10=_mm256_add_epi64(res10,_mm256_mul_epu32(aw3,b5));
res11=_mm256_add_epi64(res11,_mm256_mul_epu32(awa,b5));
const auto b6=_mm256_permute4x64_epi64(bb,0xff),b7=_mm256_shuffle_epi32(b6,_MM_PERM_CDAB);
const auto aw2=_mm256_alignr_epi8(awa,aw,8);
const auto aw1=_mm256_alignr_epi8(awa,aw,4);
res00=_mm256_add_epi64(res00,_mm256_mul_epu32(aw2,b6));
res01=_mm256_add_epi64(res01,_mm256_mul_epu32(aw3,b6));
res10=_mm256_add_epi64(res10,_mm256_mul_epu32(aw1,b7));
res11=_mm256_add_epi64(res11,_mm256_mul_epu32(aw2,b7));
res00=_mm256_add_epi64(res00,res10);
res01=_mm256_add_epi64(res01,res11);
return shrk32(reduce(res00,res01,Niv,Mod),Mod2);
}
inline void vector_convolution_direct(I256*f,const I256*g,idt lm,const FNTT32_info*const info){
u32 RR=info->one;
const auto mod=info->mod,niv=info->niv;
const auto Fx=_mm256_set1_epi32(mul_s((mod-((mod-1)>>(__builtin_ctzll(lm)))),info->r3,niv,mod));
const auto Niv=_mm256_set1_epi32(niv),Mod=_mm256_set1_epi32(mod),Mod2=_mm256_set1_epi32(info->mod2);
for(idt i=0;i<lm;++i){
store256(f+i,convolve8(f+i,g+i,_mm256_set1_epi32(RR),Fx,Niv,Mod,Mod2));
RR=mul(RR,info->RT1[__builtin_ctzll(~i)],niv,mod);
}
}
inline void vector_convolution_accumulate(I256*const result,const I256*const f,
const I256*const g,idt lm,
const FNTT32_info*const info){
u32 RR=info->one;
const auto mod=info->mod,niv=info->niv;
const auto Fx=_mm256_set1_epi32(mul_s((mod-((mod-1)>>(__builtin_ctzll(lm)))),info->r3,niv,mod));
const auto Niv=_mm256_set1_epi32(niv),Mod=_mm256_set1_epi32(mod),Mod2=_mm256_set1_epi32(info->mod2);
for(idt i=0;i<lm;++i){
const auto product=convolve8(f+i,g+i,_mm256_set1_epi32(RR),Fx,Niv,Mod,Mod2);
store256(result+i,add32(load256(result+i),product,Mod2));
RR=mul(RR,info->RT1[__builtin_ctzll(~i)],niv,mod);
}
}
} // namespace fast998_v2
} // namespace internal
} // namespace fps
} // namespace m1une
#endif // M1UNE_FPS_HAS_X86_SIMD
#line 24 "math/fps/convolution.hpp"
#ifdef M1UNE_FPS_HAS_X86_SIMD
#pragma GCC pop_options
#endif
#line 1 "math/modint.hpp"
#line 6 "math/modint.hpp"
#include <iostream>
#line 9 "math/modint.hpp"
namespace m1une {
namespace math {
template <uint32_t Modulus>
struct ModInt {
static_assert(0 < Modulus, "Modulus must be positive");
private:
uint32_t _v;
public:
static constexpr uint32_t mod() {
return Modulus;
}
static constexpr ModInt raw(uint32_t v) noexcept {
ModInt x;
x._v = v;
return x;
}
constexpr ModInt() noexcept : _v(0) {}
template <class Integer, std::enable_if_t<std::is_integral_v<Integer>, int> = 0>
constexpr ModInt(Integer v) noexcept {
if constexpr (std::is_signed_v<Integer>) {
int64_t x = static_cast<int64_t>(v) % static_cast<int64_t>(Modulus);
if (x < 0) x += Modulus;
_v = static_cast<uint32_t>(x);
} else {
_v = static_cast<uint32_t>(static_cast<uint64_t>(v) % Modulus);
}
}
constexpr uint32_t val() const noexcept {
return _v;
}
constexpr ModInt& operator++() noexcept {
_v++;
if (_v == Modulus) _v = 0;
return *this;
}
constexpr ModInt& operator--() noexcept {
if (_v == 0) _v = Modulus;
_v--;
return *this;
}
constexpr ModInt operator++(int) noexcept {
ModInt res = *this;
++*this;
return res;
}
constexpr ModInt operator--(int) noexcept {
ModInt res = *this;
--*this;
return res;
}
constexpr ModInt& operator+=(const ModInt& rhs) noexcept {
_v += rhs._v;
if (_v >= Modulus) _v -= Modulus;
return *this;
}
constexpr ModInt& operator-=(const ModInt& rhs) noexcept {
_v -= rhs._v;
if (_v >= Modulus) _v += Modulus;
return *this;
}
constexpr ModInt& operator*=(const ModInt& rhs) noexcept {
uint64_t z = _v;
z *= rhs._v;
_v = static_cast<uint32_t>(z % Modulus);
return *this;
}
constexpr ModInt& operator/=(const ModInt& rhs) noexcept {
return *this *= rhs.inv();
}
constexpr ModInt operator+(const ModInt& rhs) const noexcept {
return ModInt(*this) += rhs;
}
constexpr ModInt operator-(const ModInt& rhs) const noexcept {
return ModInt(*this) -= rhs;
}
constexpr ModInt operator*(const ModInt& rhs) const noexcept {
return ModInt(*this) *= rhs;
}
constexpr ModInt operator/(const ModInt& rhs) const noexcept {
return ModInt(*this) /= rhs;
}
constexpr bool operator==(const ModInt& rhs) const noexcept {
return _v == rhs._v;
}
constexpr bool operator!=(const ModInt& rhs) const noexcept {
return _v != rhs._v;
}
constexpr ModInt pow(long long n) const noexcept {
ModInt res = raw(1 % Modulus);
ModInt x = n < 0 ? inv() : *this;
uint64_t exponent = n < 0 ? uint64_t(-(n + 1)) + 1 : uint64_t(n);
while (exponent > 0) {
if (exponent & 1) res *= x;
x *= x;
exponent >>= 1;
}
return res;
}
constexpr ModInt inv() const noexcept {
int64_t a = _v, b = Modulus, u = 1, v = 0;
while (b) {
int64_t t = a / b;
a -= t * b;
std::swap(a, b);
u -= t * v;
std::swap(u, v);
}
assert(a == 1);
u %= Modulus;
if (u < 0) u += Modulus;
return raw(static_cast<uint32_t>(u));
}
friend std::ostream& operator<<(std::ostream& os, const ModInt& rhs) {
return os << rhs._v;
}
friend std::istream& operator>>(std::istream& is, ModInt& rhs) {
long long v;
is >> v;
rhs = ModInt(v);
return is;
}
};
using modint998244353 = ModInt<998244353>;
using modint1000000007 = ModInt<1000000007>;
template <int Id = 0>
struct DynamicModInt {
private:
uint32_t _v;
inline static uint32_t _mod = 1;
public:
static uint32_t mod() noexcept {
return _mod;
}
static void set_mod(uint32_t modulus) noexcept {
assert(modulus > 0);
assert(modulus <= uint32_t(1) << 31);
_mod = modulus;
}
static DynamicModInt raw(uint32_t v) noexcept {
assert(v < _mod);
DynamicModInt x;
x._v = v;
return x;
}
DynamicModInt() noexcept : _v(0) {}
template <class Integer, std::enable_if_t<std::is_integral_v<Integer>, int> = 0>
DynamicModInt(Integer v) noexcept {
if constexpr (std::is_signed_v<Integer>) {
int64_t x = static_cast<int64_t>(v) % static_cast<int64_t>(_mod);
if (x < 0) x += _mod;
_v = static_cast<uint32_t>(x);
} else {
_v = static_cast<uint32_t>(static_cast<uint64_t>(v) % _mod);
}
}
uint32_t val() const noexcept {
return _v;
}
DynamicModInt& operator++() noexcept {
_v++;
if (_v == _mod) _v = 0;
return *this;
}
DynamicModInt& operator--() noexcept {
if (_v == 0) _v = _mod;
_v--;
return *this;
}
DynamicModInt operator++(int) noexcept {
DynamicModInt result = *this;
++*this;
return result;
}
DynamicModInt operator--(int) noexcept {
DynamicModInt result = *this;
--*this;
return result;
}
DynamicModInt& operator+=(const DynamicModInt& rhs) noexcept {
_v += rhs._v;
if (_v >= _mod) _v -= _mod;
return *this;
}
DynamicModInt& operator-=(const DynamicModInt& rhs) noexcept {
_v -= rhs._v;
if (_v >= _mod) _v += _mod;
return *this;
}
DynamicModInt& operator*=(const DynamicModInt& rhs) noexcept {
_v = static_cast<uint32_t>(uint64_t(_v) * rhs._v % _mod);
return *this;
}
DynamicModInt& operator/=(const DynamicModInt& rhs) noexcept {
return *this *= rhs.inv();
}
DynamicModInt operator+(const DynamicModInt& rhs) const noexcept {
return DynamicModInt(*this) += rhs;
}
DynamicModInt operator-(const DynamicModInt& rhs) const noexcept {
return DynamicModInt(*this) -= rhs;
}
DynamicModInt operator*(const DynamicModInt& rhs) const noexcept {
return DynamicModInt(*this) *= rhs;
}
DynamicModInt operator/(const DynamicModInt& rhs) const noexcept {
return DynamicModInt(*this) /= rhs;
}
bool operator==(const DynamicModInt& rhs) const noexcept {
return _v == rhs._v;
}
bool operator!=(const DynamicModInt& rhs) const noexcept {
return _v != rhs._v;
}
DynamicModInt pow(long long exponent) const noexcept {
DynamicModInt result = raw(1 % _mod);
DynamicModInt base = exponent < 0 ? inv() : *this;
uint64_t magnitude =
exponent < 0 ? uint64_t(-(exponent + 1)) + 1 : uint64_t(exponent);
while (magnitude > 0) {
if (magnitude & 1) result *= base;
base *= base;
magnitude >>= 1;
}
return result;
}
DynamicModInt inv() const noexcept {
int64_t a = _v, b = _mod, u = 1, v = 0;
while (b) {
int64_t quotient = a / b;
a -= quotient * b;
std::swap(a, b);
u -= quotient * v;
std::swap(u, v);
}
assert(a == 1);
u %= _mod;
if (u < 0) u += _mod;
return raw(static_cast<uint32_t>(u));
}
friend std::ostream& operator<<(std::ostream& os, const DynamicModInt& rhs) {
return os << rhs._v;
}
friend std::istream& operator>>(std::istream& is, DynamicModInt& rhs) {
long long value;
is >> value;
rhs = DynamicModInt(value);
return is;
}
};
} // namespace math
} // namespace m1une
#line 29 "math/fps/convolution.hpp"
namespace m1une {
namespace fps {
namespace internal {
template <class Mint, class = void>
struct has_static_modulus : std::false_type {};
template <class Mint>
struct has_static_modulus<
Mint, std::void_t<decltype(std::integral_constant<uint32_t, Mint::mod()>{})>>
: std::true_type {};
constexpr uint32_t primitive_root_constexpr(uint32_t mod) {
if (mod == 2) return 1;
if (mod == 167772161) return 3;
if (mod == 469762049) return 3;
if (mod == 754974721) return 11;
if (mod == 998244353) return 3;
if (mod == 1224736769) return 3;
uint32_t divisors[32] = {};
int count = 0;
uint32_t x = mod - 1;
for (uint32_t p = 2; uint64_t(p) * p <= x; p++) {
if (x % p != 0) continue;
divisors[count++] = p;
while (x % p == 0) x /= p;
}
if (x > 1) divisors[count++] = x;
for (uint32_t g = 2;; g++) {
bool ok = true;
for (int i = 0; i < count; i++) {
uint64_t value = 1;
uint64_t base = g;
uint32_t exponent = (mod - 1) / divisors[i];
while (exponent > 0) {
if (exponent & 1) value = value * base % mod;
base = base * base % mod;
exponent >>= 1;
}
if (value == 1) {
ok = false;
break;
}
}
if (ok) return g;
}
}
constexpr int two_adic_order(uint32_t x) {
int result = 0;
while ((x & 1) == 0) {
x >>= 1;
result++;
}
return result;
}
template <class Mint>
struct NttRoots {
static constexpr int max_base = two_adic_order(Mint::mod() - 1);
std::array<Mint, max_base + 1> root;
std::array<Mint, max_base + 1> inverse_root;
std::array<Mint, max_base> rate;
std::array<Mint, max_base> inverse_rate;
std::array<Mint, max_base> rate_radix4;
std::array<Mint, max_base> inverse_rate_radix4;
NttRoots() {
constexpr uint32_t primitive_root = primitive_root_constexpr(Mint::mod());
for (int level = 1; level <= max_base; level++) {
root[level] = Mint(primitive_root).pow((Mint::mod() - 1) >> level);
inverse_root[level] = root[level].inv();
}
Mint product = 1;
Mint inverse_product = 1;
for (int i = 0; i + 1 < max_base; i++) {
rate[i] = root[i + 2] * product;
inverse_rate[i] = inverse_root[i + 2] * inverse_product;
product *= inverse_root[i + 2];
inverse_product *= root[i + 2];
}
product = 1;
inverse_product = 1;
for (int i = 0; i + 2 < max_base; i++) {
rate_radix4[i] = root[i + 3] * product;
inverse_rate_radix4[i] = inverse_root[i + 3] * inverse_product;
product *= inverse_root[i + 3];
inverse_product *= root[i + 3];
}
}
};
template <class Mint>
const NttRoots<Mint>& ntt_roots() {
static const NttRoots<Mint> roots;
return roots;
}
template <class Mint>
void ntt(std::vector<Mint>& a, bool inverse, bool normalize = true) {
const int n = int(a.size());
assert(n > 0 && (n & (n - 1)) == 0);
assert((Mint::mod() - 1) % uint32_t(n) == 0);
const auto& roots = ntt_roots<Mint>();
const int height = two_adic_order(uint32_t(n));
if (!inverse) {
int phase = 0;
while (phase < height) {
if (height - phase == 1) {
const int width = 1 << (height - phase - 1);
Mint twiddle = 1;
for (int block = 0; block < (1 << phase); block++) {
const int offset = block << (height - phase);
for (int i = 0; i < width; i++) {
const Mint left = a[offset + i];
const Mint right = a[offset + i + width] * twiddle;
a[offset + i] = left + right;
a[offset + i + width] = left - right;
}
if (block + 1 != (1 << phase))
twiddle *= roots.rate[__builtin_ctz(~uint32_t(block))];
}
phase++;
continue;
}
const int width = 1 << (height - phase - 2);
Mint twiddle = 1;
const Mint imaginary = roots.root[2];
for (int block = 0; block < (1 << phase); block++) {
const Mint twiddle2 = twiddle * twiddle;
const Mint twiddle3 = twiddle2 * twiddle;
const int offset = block << (height - phase);
for (int i = 0; i < width; i++) {
const uint64_t mod2 = uint64_t(Mint::mod()) * Mint::mod();
const uint64_t a0 = a[offset + i].val();
const uint64_t a1 = uint64_t(a[offset + i + width].val()) * twiddle.val();
const uint64_t a2 =
uint64_t(a[offset + i + 2 * width].val()) * twiddle2.val();
const uint64_t a3 =
uint64_t(a[offset + i + 3 * width].val()) * twiddle3.val();
const uint64_t a1na3i =
uint64_t(Mint(a1 + mod2 - a3).val()) * imaginary.val();
const uint64_t negative_a2 = mod2 - a2;
a[offset + i] = Mint(a0 + a2 + a1 + a3);
a[offset + i + width] = Mint(a0 + a2 + 2 * mod2 - a1 - a3);
a[offset + i + 2 * width] = Mint(a0 + negative_a2 + a1na3i);
a[offset + i + 3 * width] = Mint(a0 + negative_a2 + mod2 - a1na3i);
}
if (block + 1 != (1 << phase))
twiddle *= roots.rate_radix4[__builtin_ctz(~uint32_t(block))];
}
phase += 2;
}
} else {
int phase = height;
while (phase > 0) {
if (phase == 1) {
const int width = 1 << (height - phase);
Mint twiddle = 1;
for (int block = 0; block < (1 << (phase - 1)); block++) {
const int offset = block << (height - phase + 1);
for (int i = 0; i < width; i++) {
const Mint left = a[offset + i];
const Mint right = a[offset + i + width];
a[offset + i] = left + right;
a[offset + i + width] = (left - right) * twiddle;
}
if (block + 1 != (1 << (phase - 1)))
twiddle *= roots.inverse_rate[__builtin_ctz(~uint32_t(block))];
}
phase--;
continue;
}
const int width = 1 << (height - phase);
Mint twiddle = 1;
const Mint inverse_imaginary = roots.inverse_root[2];
for (int block = 0; block < (1 << (phase - 2)); block++) {
const Mint twiddle2 = twiddle * twiddle;
const Mint twiddle3 = twiddle2 * twiddle;
const int offset = block << (height - phase + 2);
for (int i = 0; i < width; i++) {
const uint64_t a0 = a[offset + i].val();
const uint64_t a1 = a[offset + i + width].val();
const uint64_t a2 = a[offset + i + 2 * width].val();
const uint64_t a3 = a[offset + i + 3 * width].val();
const uint64_t a2na3i =
uint64_t(Mint((Mint::mod() + a2 - a3) * inverse_imaginary.val()).val());
a[offset + i] = Mint(a0 + a1 + a2 + a3);
a[offset + i + width] =
Mint((a0 + Mint::mod() - a1 + a2na3i) * twiddle.val());
a[offset + i + 2 * width] = Mint(
(a0 + a1 + 2ULL * Mint::mod() - a2 - a3) * twiddle2.val());
a[offset + i + 3 * width] = Mint(
(a0 + Mint::mod() - a1 + Mint::mod() - a2na3i) * twiddle3.val());
}
if (block + 1 != (1 << (phase - 2)))
twiddle *= roots.inverse_rate_radix4[__builtin_ctz(~uint32_t(block))];
}
phase -= 2;
}
if (normalize) {
const Mint inverse_n = Mint(n).inv();
for (Mint& value : a) value *= inverse_n;
}
}
}
#ifdef M1UNE_FPS_HAS_X86_SIMD
#pragma GCC push_options
#pragma GCC target("avx2,bmi")
template <class Mint>
__attribute__((target("avx2,bmi"), hot))
std::vector<Mint> convolution_998244353_simd(const std::vector<Mint>& a,
const std::vector<Mint>& b) {
const int result_size = int(a.size() + b.size() - 1);
int n = 1;
while (n < result_size) n <<= 1;
const bool squaring = &a == &b;
auto* transformed_a = static_cast<uint32_t*>(
::operator new[](sizeof(uint32_t) * n, std::align_val_t(32)));
auto* transformed_b = squaring
? transformed_a
: static_cast<uint32_t*>(::operator new[](
sizeof(uint32_t) * n, std::align_val_t(32)));
if constexpr (std::is_same_v<Mint, math::ModInt<998244353>>) {
static_assert(sizeof(Mint) == sizeof(uint32_t) && std::is_trivially_copyable_v<Mint>);
std::memcpy(transformed_a, a.data(), sizeof(uint32_t) * a.size());
if (!squaring)
std::memcpy(transformed_b, b.data(), sizeof(uint32_t) * b.size());
} else {
for (int i = 0; i < int(a.size()); i++) transformed_a[i] = a[i].val();
if (!squaring)
for (int i = 0; i < int(b.size()); i++) transformed_b[i] = b[i].val();
}
std::memset(transformed_a + a.size(), 0, sizeof(uint32_t) * (n - a.size()));
if (!squaring)
std::memset(transformed_b + b.size(), 0, sizeof(uint32_t) * (n - b.size()));
static constexpr fast998_v2::FNTT32_info transform(998244353);
const std::size_t vector_size = std::size_t(n) >> 3;
fast998_v2::vector_dif(reinterpret_cast<__m256i*>(transformed_a), vector_size, &transform);
if (!squaring)
fast998_v2::vector_dif(reinterpret_cast<__m256i*>(transformed_b), vector_size,
&transform);
fast998_v2::vector_convolution_direct(
reinterpret_cast<__m256i*>(transformed_a),
reinterpret_cast<const __m256i*>(transformed_b), vector_size, &transform);
fast998_v2::vector_dit<true>(reinterpret_cast<__m256i*>(transformed_a), vector_size,
&transform);
std::vector<Mint> result(result_size);
for (int j = 0; j < result_size; j++) result[j] = Mint::raw(transformed_a[j]);
::operator delete[](transformed_a, std::align_val_t(32));
if (!squaring) ::operator delete[](transformed_b, std::align_val_t(32));
return result;
}
#pragma GCC pop_options
#endif
} // namespace internal
template <class Mint>
std::vector<Mint> convolution_naive(const std::vector<Mint>& a, const std::vector<Mint>& b) {
if (a.empty() || b.empty()) return {};
std::vector<Mint> result(a.size() + b.size() - 1);
if (a.size() < b.size()) {
for (int i = 0; i < int(a.size()); i++) {
for (int j = 0; j < int(b.size()); j++) result[i + j] += a[i] * b[j];
}
} else {
for (int j = 0; j < int(b.size()); j++) {
for (int i = 0; i < int(a.size()); i++) result[i + j] += a[i] * b[j];
}
}
return result;
}
template <class Mint>
std::vector<Mint> convolution_ntt(const std::vector<Mint>& a, const std::vector<Mint>& b) {
const int result_size = int(a.size() + b.size() - 1);
int n = 1;
while (n < result_size) n <<= 1;
assert((Mint::mod() - 1) % uint32_t(n) == 0);
#ifdef M1UNE_FPS_HAS_X86_SIMD
if constexpr (Mint::mod() == 998244353) {
if (n >= 64 && __builtin_cpu_supports("avx2"))
return internal::convolution_998244353_simd(a, b);
}
#endif
// Allocate the padded buffers directly. Constructing from the inputs and
// then resizing used to allocate and copy both large operands twice.
const bool squaring = &a == &b;
std::vector<Mint> fa(n);
std::copy(a.begin(), a.end(), fa.begin());
internal::ntt(fa, false);
const Mint inverse_n = Mint(n).inv();
if (squaring) {
for (int i = 0; i < n; i++) fa[i] *= fa[i] * inverse_n;
} else {
std::vector<Mint> fb(n);
std::copy(b.begin(), b.end(), fb.begin());
internal::ntt(fb, false);
for (int i = 0; i < n; i++) fa[i] *= fb[i] * inverse_n;
}
internal::ntt(fa, true, false);
fa.resize(result_size);
return fa;
}
namespace internal {
template <class Mint>
std::vector<Mint> convolution_998244353_blocked_scalar(const std::vector<Mint>& a,
const std::vector<Mint>& b,
int transform_size) {
assert(Mint::mod() == 998244353);
assert(transform_size >= 2 && (transform_size & (transform_size - 1)) == 0);
assert((Mint::mod() - 1) % uint32_t(transform_size) == 0);
const int block_size = transform_size / 2;
const int a_blocks = int((a.size() + block_size - 1) / block_size);
const int b_blocks = int((b.size() + block_size - 1) / block_size);
auto transform_blocks = [&](const std::vector<Mint>& values, int block_count) {
std::vector<std::vector<Mint>> blocks;
blocks.reserve(block_count);
for (int block = 0; block < block_count; block++) {
const int begin = block * block_size;
const int count = std::min(block_size, int(values.size()) - begin);
std::vector<Mint> transformed(transform_size);
std::copy_n(values.begin() + begin, count, transformed.begin());
ntt(transformed, false);
blocks.emplace_back(std::move(transformed));
}
return blocks;
};
std::vector<std::vector<Mint>> transformed_a = transform_blocks(a, a_blocks);
std::vector<std::vector<Mint>> transformed_b = transform_blocks(b, b_blocks);
const int result_size = int(a.size() + b.size() - 1);
std::vector<Mint> result(result_size);
std::vector<Mint> transformed_result(transform_size);
for (int diagonal = 0; diagonal < a_blocks + b_blocks - 1; diagonal++) {
std::fill(transformed_result.begin(), transformed_result.end(), Mint(0));
const int first_a = std::max(0, diagonal - (b_blocks - 1));
const int last_a = std::min(a_blocks - 1, diagonal);
for (int a_block = first_a; a_block <= last_a; a_block++) {
const int b_block = diagonal - a_block;
for (int i = 0; i < transform_size; i++)
transformed_result[i] +=
transformed_a[a_block][i] * transformed_b[b_block][i];
}
ntt(transformed_result, true);
const int output_offset = diagonal * block_size;
const int output_count = std::min(transform_size, result_size - output_offset);
for (int i = 0; i < output_count; i++)
result[output_offset + i] += transformed_result[i];
}
return result;
}
#ifdef M1UNE_FPS_HAS_X86_SIMD
class AlignedUint32Buffer {
private:
uint32_t* data_;
public:
explicit AlignedUint32Buffer(std::size_t size)
: data_(static_cast<uint32_t*>(
::operator new[](sizeof(uint32_t) * size, std::align_val_t(32)))) {}
AlignedUint32Buffer(const AlignedUint32Buffer&) = delete;
AlignedUint32Buffer& operator=(const AlignedUint32Buffer&) = delete;
AlignedUint32Buffer(AlignedUint32Buffer&& other) noexcept : data_(other.data_) {
other.data_ = nullptr;
}
AlignedUint32Buffer& operator=(AlignedUint32Buffer&& other) noexcept {
if (this == &other) return *this;
::operator delete[](data_, std::align_val_t(32));
data_ = other.data_;
other.data_ = nullptr;
return *this;
}
~AlignedUint32Buffer() {
::operator delete[](data_, std::align_val_t(32));
}
uint32_t* data() {
return data_;
}
const uint32_t* data() const {
return data_;
}
};
template <class Mint>
__attribute__((target("avx2,bmi"), hot))
std::vector<Mint> convolution_998244353_blocked_simd(const std::vector<Mint>& a,
const std::vector<Mint>& b,
int transform_size) {
assert(Mint::mod() == 998244353);
assert(transform_size >= 64 && (transform_size & (transform_size - 1)) == 0);
assert((Mint::mod() - 1) % uint32_t(transform_size) == 0);
const int block_size = transform_size / 2;
const int a_blocks = int((a.size() + block_size - 1) / block_size);
const int b_blocks = int((b.size() + block_size - 1) / block_size);
static constexpr fast998_v2::FNTT32_info transform(998244353);
const std::size_t vector_size = std::size_t(transform_size) / 8;
auto transform_blocks = [&](const std::vector<Mint>& values, int block_count) {
std::vector<AlignedUint32Buffer> blocks;
blocks.reserve(block_count);
for (int block = 0; block < block_count; block++) {
const int begin = block * block_size;
const int count = std::min(block_size, int(values.size()) - begin);
AlignedUint32Buffer transformed(transform_size);
if constexpr (std::is_same_v<Mint, math::ModInt<998244353>>) {
static_assert(sizeof(Mint) == sizeof(uint32_t) &&
std::is_trivially_copyable_v<Mint>);
std::memcpy(transformed.data(), values.data() + begin,
sizeof(uint32_t) * count);
} else {
for (int i = 0; i < count; i++)
transformed.data()[i] = values[begin + i].val();
}
std::memset(transformed.data() + count, 0,
sizeof(uint32_t) * (transform_size - count));
fast998_v2::vector_dif(reinterpret_cast<__m256i*>(transformed.data()),
vector_size, &transform);
blocks.emplace_back(std::move(transformed));
}
return blocks;
};
std::vector<AlignedUint32Buffer> transformed_a = transform_blocks(a, a_blocks);
std::vector<AlignedUint32Buffer> transformed_b = transform_blocks(b, b_blocks);
const int result_size = int(a.size() + b.size() - 1);
std::vector<Mint> result(result_size);
AlignedUint32Buffer transformed_result(transform_size);
for (int diagonal = 0; diagonal < a_blocks + b_blocks - 1; diagonal++) {
std::memset(transformed_result.data(), 0, sizeof(uint32_t) * transform_size);
const int first_a = std::max(0, diagonal - (b_blocks - 1));
const int last_a = std::min(a_blocks - 1, diagonal);
for (int a_block = first_a; a_block <= last_a; a_block++) {
const int b_block = diagonal - a_block;
fast998_v2::vector_convolution_accumulate(
reinterpret_cast<__m256i*>(transformed_result.data()),
reinterpret_cast<const __m256i*>(transformed_a[a_block].data()),
reinterpret_cast<const __m256i*>(transformed_b[b_block].data()),
vector_size, &transform);
}
fast998_v2::vector_dit<true>(
reinterpret_cast<__m256i*>(transformed_result.data()), vector_size,
&transform);
const int output_offset = diagonal * block_size;
const int output_count = std::min(transform_size, result_size - output_offset);
for (int i = 0; i < output_count; i++) {
uint32_t value = result[output_offset + i].val() + transformed_result.data()[i];
if (value >= Mint::mod()) value -= Mint::mod();
result[output_offset + i] = Mint::raw(value);
}
}
return result;
}
#endif
template <class Mint>
std::vector<Mint> convolution_998244353_blocked(const std::vector<Mint>& a,
const std::vector<Mint>& b,
int transform_size = 1 << 23) {
#ifdef M1UNE_FPS_HAS_X86_SIMD
if (transform_size >= 64 && __builtin_cpu_supports("avx2"))
return convolution_998244353_blocked_simd(a, b, transform_size);
#endif
return convolution_998244353_blocked_scalar(a, b, transform_size);
}
} // namespace internal
template <class Mint>
std::vector<Mint> convolution(const std::vector<Mint>& a, const std::vector<Mint>& b) {
if (a.empty() || b.empty()) return {};
if (std::min(a.size(), b.size()) <= 32) return convolution_naive(a, b);
const int result_size = int(a.size() + b.size() - 1);
int n = 1;
while (n < result_size) n <<= 1;
if constexpr (internal::has_static_modulus<Mint>::value) {
if constexpr (Mint::mod() == 998244353) {
if (n > (1 << 23))
return internal::convolution_998244353_blocked(a, b);
}
if ((Mint::mod() - 1) % uint32_t(n) == 0) return convolution_ntt(a, b);
}
using Mint1 = math::ModInt<167772161>;
using Mint2 = math::ModInt<469762049>;
using Mint3 = math::ModInt<754974721>;
assert(n <= (1 << 24));
[[maybe_unused]] const unsigned __int128 coefficient_bound =
static_cast<unsigned __int128>(std::min(a.size(), b.size())) * (Mint::mod() - 1) *
(Mint::mod() - 1);
[[maybe_unused]] const unsigned __int128 crt_modulus =
static_cast<unsigned __int128>(Mint1::mod()) * Mint2::mod() * Mint3::mod();
assert(coefficient_bound < crt_modulus);
auto converted_convolution = [&]<class OtherMint>() {
std::vector<OtherMint> converted_a(a.size());
std::vector<OtherMint> converted_b(b.size());
for (int i = 0; i < int(a.size()); i++) converted_a[i] = OtherMint(a[i].val());
for (int i = 0; i < int(b.size()); i++) converted_b[i] = OtherMint(b[i].val());
return convolution_ntt(converted_a, converted_b);
};
std::vector<Mint1> c1 = converted_convolution.template operator()<Mint1>();
std::vector<Mint2> c2 = converted_convolution.template operator()<Mint2>();
std::vector<Mint3> c3 = converted_convolution.template operator()<Mint3>();
static const uint64_t inverse_mod1_mod2 = Mint2(Mint1::mod()).inv().val();
static const uint64_t mod1_mod3 = Mint1::mod() % Mint3::mod();
static const uint64_t mod1_mod2_mod3 =
mod1_mod3 * (Mint2::mod() % Mint3::mod()) % Mint3::mod();
static const uint64_t inverse_mod1_mod2_mod3 = Mint3(uint32_t(mod1_mod2_mod3)).inv().val();
const uint64_t target_mod = Mint::mod();
const uint64_t mod1_target = Mint1::mod() % target_mod;
const uint64_t mod1_mod2_target = mod1_target * (Mint2::mod() % target_mod) % target_mod;
std::vector<Mint> result(result_size);
for (int i = 0; i < result_size; i++) {
const uint64_t r1 = c1[i].val();
const uint64_t r2 = c2[i].val();
const uint64_t r3 = c3[i].val();
const uint64_t first =
(r2 + Mint2::mod() - r1 % Mint2::mod()) % Mint2::mod() * inverse_mod1_mod2 %
Mint2::mod();
const uint64_t combined_mod3 =
(r1 % Mint3::mod() + mod1_mod3 * (first % Mint3::mod())) % Mint3::mod();
const uint64_t second =
(r3 + Mint3::mod() - combined_mod3) % Mint3::mod() * inverse_mod1_mod2_mod3 %
Mint3::mod();
uint64_t value = r1 % target_mod;
value = (value + mod1_target * (first % target_mod)) % target_mod;
value = (value + mod1_mod2_target * (second % target_mod)) % target_mod;
result[i] = Mint::raw(uint32_t(value));
}
return result;
}
} // namespace fps
} // namespace m1une
#ifdef M1UNE_FPS_HAS_X86_SIMD
#undef M1UNE_FPS_HAS_X86_SIMD
#endif
#line 9 "string/wildcard_pattern_matching.hpp"
namespace m1une {
namespace string {
namespace internal {
template <class Mint>
std::vector<Mint> wildcard_mismatch_scores(
const std::string& text,
const std::string& pattern,
char wildcard
) {
int n = int(text.size());
int m = int(pattern.size());
std::array<std::vector<Mint>, 3> text_powers;
std::array<std::vector<Mint>, 3> pattern_powers;
for (auto& powers : text_powers) powers.resize(n);
for (auto& powers : pattern_powers) powers.resize(m);
auto encode = [wildcard](char character) {
if (character == wildcard) return 0;
return int(static_cast<unsigned char>(character)) + 1;
};
for (int i = 0; i < n; i++) {
Mint value = encode(text[i]);
text_powers[0][i] = value;
text_powers[1][i] = value * value;
text_powers[2][i] = value * value * value;
}
for (int i = 0; i < m; i++) {
Mint value = encode(pattern[m - 1 - i]);
pattern_powers[0][i] = value;
pattern_powers[1][i] = value * value;
pattern_powers[2][i] = value * value * value;
}
std::vector<Mint> scores(n - m + 1);
for (int power = 0; power < 3; power++) {
std::vector<Mint> product = fps::convolution(
text_powers[power],
pattern_powers[2 - power]
);
Mint multiplier = power == 1 ? Mint(-2) : Mint(1);
for (int start = 0; start <= n - m; start++) {
scores[start] += multiplier * product[start + m - 1];
}
}
return scores;
}
} // namespace internal
// result[i] is true exactly when pattern matches text[i, i + pattern.size()).
// The wildcard character matches every character on either side.
inline std::vector<bool> wildcard_pattern_matching(
const std::string& text,
const std::string& pattern,
char wildcard = '*'
) {
int n = int(text.size());
int m = int(pattern.size());
if (m == 0) return std::vector<bool>(n + 1, true);
if (n < m) return {};
using Mint1 = math::ModInt<469762049>;
using Mint2 = math::ModInt<754974721>;
std::vector<Mint1> first_scores =
internal::wildcard_mismatch_scores<Mint1>(text, pattern, wildcard);
std::vector<Mint2> second_scores =
internal::wildcard_mismatch_scores<Mint2>(text, pattern, wildcard);
std::vector<bool> result(n - m + 1, false);
for (int start = 0; start <= n - m; start++) {
result[start] = first_scores[start].val() == 0 &&
second_scores[start].val() == 0;
}
return result;
}
} // namespace string
} // namespace m1une
#line 1 "string/z_algorithm.hpp"
#line 6 "string/z_algorithm.hpp"
namespace m1une {
namespace string {
// Returns z[i] = LCP(sequence, sequence[i..]).
template <class Sequence>
std::vector<int> z_algorithm(const Sequence& sequence) {
int n = int(sequence.size());
if (n == 0) return {};
std::vector<int> z(n);
z[0] = n;
int left = 0;
int right = 0;
for (int i = 1; i < n; i++) {
if (i < right) z[i] = std::min(right - i, z[i - left]);
while (i + z[i] < n && sequence[z[i]] == sequence[i + z[i]]) {
z[i]++;
}
if (right < i + z[i]) {
left = i;
right = i + z[i];
}
}
return z;
}
} // namespace string
} // namespace m1une
#line 27 "string/all.hpp"