Combinatorial Sequences
(math/combinatorial_sequences.hpp)
- View this file on GitHub
- Last update: 2026-08-10 17:30:05+09:00
- Include:
#include "math/combinatorial_sequences.hpp"
Overview
This header generates several standard counting sequences:
- Catalan numbers;
- Bernoulli numbers;
- Bell numbers;
- Stirling numbers of the second kind;
- integer partition numbers;
- derangement numbers.
Every function returns a vector of modular values. Except for the Stirling
function, requesting maximum returns every value from index 0 through
maximum.
The Bernoulli, Bell, Stirling, and partition implementations use formal power series and fast convolution. They are suitable for hundreds of thousands of terms, rather than only the small sizes supported by quadratic dynamic programming.
Catalan Numbers
The n-th Catalan number is
The first values are:
n: 0 1 2 3 4 5
C_n: 1 1 2 5 14 42
Catalan numbers count many equivalent structures, including:
- balanced parenthesis strings containing
npairs; - binary tree shapes with
ninternal vertices; - ways to triangulate a convex polygon with
n + 2vertices; - monotone grid paths that never cross above the diagonal.
For example, the two balanced strings with two pairs are:
()()
(())
catalan_numbers<Mint>(maximum) computes the sequence in $O(\text{maximum})$
time using
If a Combinatorics<Mint> object already contains factorials through 2 * n,
its catalan(n) method returns one Catalan number in $O(1)$ time.
Bernoulli Numbers
Bernoulli numbers are rational numbers defined by the exponential generating function
\[\frac{x}{e^x-1} = \sum_{n=0}^{\infty} B_n \frac{x^n}{n!}.\]This header uses the convention
\[B_1 = -\frac{1}{2}.\]The first values are:
B_0 = 1
B_1 = -1/2
B_2 = 1/6
B_3 = 0
B_4 = -1/30
B_5 = 0
B_6 = 1/42
Bernoulli numbers do not directly count a simple family of objects. Their most common contest use is evaluating sums of powers. For example,
\[1^p + 2^p + \cdots + n^p = \frac{1}{p+1} \sum_{j=0}^{p} (-1)^j \binom{p+1}{j} B_j n^{p+1-j}.\]They also appear in finite differences, polynomial interpolation, and advanced number theory.
The returned values are modular representations of the rational numbers.
For example, B_1 is -1 * inverse(2) modulo Mint::mod().
bernoulli_numbers<Mint>(maximum) computes all terms in
$O(\text{maximum} \log \text{maximum})$ time by inverting the series
The function is defined in math/bernoulli.hpp, which also
provides Bernoulli polynomials, power sums, and polynomial discrete
integration. It remains available through this header for compatibility.
Stirling Numbers of the Second Kind
The Stirling number of the second kind $S(n,k)$ counts ways to divide a set of
n distinct elements into exactly k nonempty, unlabeled groups.
For the set containing a, b, and c:
- $S(3,1)=1$: all elements are in one group;
- $S(3,2)=3$: one element is alone and the other two are together;
- $S(3,3)=1$: every element is alone.
The row for n = 5 is:
k: 0 1 2 3 4 5
S(5,k): 0 1 15 25 10 1
stirling_numbers_second_kind<Mint>(n) returns this entire row. It uses the
convolution identity
and runs in $O(n \log n)$ time.
Bell Numbers
The Bell number counts ways to divide a set of n distinct elements into any
number of nonempty, unlabeled groups.
The first values are:
n: 0 1 2 3 4 5
Bell_n: 1 1 2 5 15 52
Bell numbers and Stirling numbers are related by
\[\operatorname{Bell}_n = \sum_{k=0}^{n} S(n,k).\]Their exponential generating function is
\[\exp(e^x-1) = \sum_{n=0}^{\infty} \operatorname{Bell}_n \frac{x^n}{n!}.\]bell_numbers<Mint>(maximum) evaluates this series in
$O(\text{maximum} \log \text{maximum})$ time.
Bernoulli numbers and Bell numbers are both traditionally written as B_n,
but they are unrelated sequences. The API uses distinct function names to
avoid that ambiguity.
Integer Partition Numbers
The partition number $p(n)$ counts ways to write n as a sum of positive
integers when order does not matter.
For example, $p(5)=7$:
5
4 + 1
3 + 2
3 + 1 + 1
2 + 2 + 1
2 + 1 + 1 + 1
1 + 1 + 1 + 1 + 1
This is different from a Bell number: integer partitions split an integer, while Bell numbers split a set of distinct elements.
The generating function is
\[\prod_{k=1}^{\infty}\frac{1}{1-x^k} = \sum_{n=0}^{\infty}p(n)x^n.\]Euler’s pentagonal number theorem gives a sparse denominator, which the
implementation inverts with formal power series. Therefore
partition_function<Mint>(maximum) and its compatibility alias
partition_numbers<Mint>(maximum) run in
$O(\text{maximum} \log \text{maximum})$ time.
Derangements
A derangement is a permutation in which no element remains in its original position.
For three elements, the two derangements are:
2 3 1
3 1 2
The first values are:
n: 0 1 2 3 4 5
D_n: 1 0 1 2 9 44
The sequence satisfies
\[D_n = (n-1)(D_{n-1}+D_{n-2}).\]derangement_numbers<Mint>(maximum) uses this recurrence and runs in
$O(\text{maximum})$ time.
API
template <class Mint>
std::vector<Mint> catalan_numbers(int maximum);
template <class Mint>
std::vector<Mint> bernoulli_numbers(int maximum);
template <class Mint>
std::vector<Mint> bell_numbers(int maximum);
template <class Mint>
std::vector<Mint> stirling_numbers_second_kind(int n);
template <class Mint>
std::vector<Mint> partition_numbers(int maximum);
template <class Mint>
std::vector<Mint> derangement_numbers(int maximum);
The template argument Mint is the element type of the returned
std::vector<Mint>. The function argument is an int and must be
nonnegative. Except for Stirling numbers, a call with maximum returns a
vector of length maximum + 1. The Stirling function returns a vector of
length n + 1.
| Function | Returned values | Complexity |
|---|---|---|
catalan_numbers<Mint>(maximum) |
C_0 through C_maximum
|
$O(n)$ |
bernoulli_numbers<Mint>(maximum) |
B_0 through B_maximum
|
$O(n \log n)$ |
bell_numbers<Mint>(maximum) |
Bell numbers 0 through maximum
|
$O(n \log n)$ |
stirling_numbers_second_kind<Mint>(n) |
S(n, 0) through S(n, n)
|
$O(n \log n)$ |
partition_numbers<Mint>(maximum) |
p(0) through p(maximum)
|
$O(n \log n)$ |
derangement_numbers<Mint>(maximum) |
D_0 through D_maximum
|
$O(n)$ |
Every function uses $O(n)$ memory.
Requirements
The fast formulas are intended for a static modular integer such as
modint998244353.
- The modulus should be prime for Catalan, Bernoulli, Bell, and Stirling numbers because the algorithms use modular inverses.
-
maximum + 1must be smaller than the modulus for Catalan and Bernoulli numbers. -
maximummust be smaller than the modulus for Bell and Stirling numbers. - The formal-power-series transform length must be supported by the convolution
implementation. Modulus
998244353is the standard fast choice.
These functions return residues modulo Mint::mod(), not arbitrary-precision
integers.
Example
#include "math/combinatorial_sequences.hpp"
#include "math/modint.hpp"
#include <iostream>
using Mint = m1une::math::modint998244353;
int main() {
auto catalan = m1une::math::catalan_numbers<Mint>(10);
auto bernoulli = m1une::math::bernoulli_numbers<Mint>(10);
auto partitions = m1une::math::partition_numbers<Mint>(10);
std::cout << catalan[5] << "\n"; // 42
std::cout << bernoulli[1] * 2 << "\n"; // -1 modulo 998244353
std::cout << partitions[5] << "\n"; // 7
}
Depends on
- Bernoulli Numbers and Power Sums
(math/bernoulli.hpp)
- Combinatorics
(math/combinatorics.hpp)
- Convolution
(math/fps/convolution.hpp)
- Formal Power Series
(math/fps/formal_power_series.hpp)
- math/fps/internal/ntt998_faster.hpp
- ModInt
(math/modint.hpp)
- Modular Square Root
(math/modular_square_root.hpp)
- Partition Function
(math/partition_function.hpp)
Required by
Verified with
- verify/math/bell_number.test.cpp
- verify/math/math_algorithms.test.cpp
- verify/math/stirling_number_of_the_second_kind.test.cpp
Code
#ifndef M1UNE_MATH_COMBINATORIAL_SEQUENCES_HPP
#define M1UNE_MATH_COMBINATORIAL_SEQUENCES_HPP 1
#include <cassert>
#include <cstdint>
#include <vector>
#include "fps/formal_power_series.hpp"
#include "bernoulli.hpp"
#include "combinatorics.hpp"
#include "partition_function.hpp"
namespace m1une {
namespace math {
template <class Mint>
std::vector<Mint> catalan_numbers(int maximum) {
assert(maximum >= 0);
assert(static_cast<uint64_t>(maximum) + 1 < Mint::mod());
std::vector<Mint> inverse(maximum + 2);
inverse[1] = 1;
for (int i = 2; i <= maximum + 1; i++) {
inverse[i] = Mint(0) - Mint(Mint::mod() / uint32_t(i)) * inverse[Mint::mod() % uint32_t(i)];
}
std::vector<Mint> result(maximum + 1);
result[0] = 1;
for (int n = 0; n < maximum; n++) {
result[n + 1] = result[n] * Mint(2) * Mint(2LL * n + 1) * inverse[n + 2];
}
return result;
}
template <class Mint>
std::vector<Mint> bell_numbers(int maximum) {
assert(maximum >= 0);
assert(static_cast<uint64_t>(maximum) < Mint::mod());
using Fps = fps::FormalPowerSeries<Mint>;
Combinatorics<Mint> combinations(maximum);
Fps exponent(maximum + 1);
for (int i = 1; i <= maximum; i++) {
exponent[i] = combinations.inverse_factorial(i);
}
Fps generating_function = exponent.exp(maximum + 1);
std::vector<Mint> result(maximum + 1);
for (int i = 0; i <= maximum; i++) {
result[i] = generating_function[i] * combinations.factorial(i);
}
return result;
}
template <class Mint>
std::vector<Mint> stirling_numbers_second_kind(int n) {
assert(n >= 0);
assert(static_cast<uint64_t>(n) < Mint::mod());
Combinatorics<Mint> combinations(n);
std::vector<Mint> powers(n + 1);
std::vector<Mint> signs(n + 1);
for (int i = 0; i <= n; i++) {
powers[i] = Mint(i).pow(n) * combinations.inverse_factorial(i);
signs[i] = combinations.inverse_factorial(i);
if (i & 1) signs[i] = Mint(0) - signs[i];
}
std::vector<Mint> result = fps::convolution(powers, signs);
result.resize(n + 1);
return result;
}
template <class Mint>
std::vector<Mint> derangement_numbers(int maximum) {
assert(maximum >= 0);
std::vector<Mint> result(maximum + 1);
result[0] = 1;
if (maximum >= 1) result[1] = 0;
for (int n = 2; n <= maximum; n++) {
result[n] = Mint(n - 1) * (result[n - 1] + result[n - 2]);
}
return result;
}
} // namespace math
} // namespace m1une
#endif // M1UNE_MATH_COMBINATORIAL_SEQUENCES_HPP#line 1 "math/combinatorial_sequences.hpp"
#include <cassert>
#include <cstdint>
#include <vector>
#line 1 "math/fps/formal_power_series.hpp"
#include <algorithm>
#line 7 "math/fps/formal_power_series.hpp"
#include <optional>
#include <utility>
#line 10 "math/fps/formal_power_series.hpp"
#line 1 "math/modular_square_root.hpp"
#line 7 "math/modular_square_root.hpp"
namespace m1une {
namespace math {
namespace internal {
inline uint64_t modular_square_root_multiply(uint64_t lhs, uint64_t rhs, uint64_t mod) {
return static_cast<uint64_t>(static_cast<unsigned __int128>(lhs) * rhs % mod);
}
inline uint64_t modular_square_root_power(uint64_t base, uint64_t exponent, uint64_t mod) {
uint64_t result = 1 % mod;
while (exponent > 0) {
if (exponent & 1) result = modular_square_root_multiply(result, base, mod);
base = modular_square_root_multiply(base, base, mod);
exponent >>= 1;
}
return result;
}
} // namespace internal
// Returns x such that x * x = value (mod prime), or nullopt when no such x exists.
// The modulus must be prime.
inline std::optional<uint64_t> modular_square_root(uint64_t value, uint64_t prime) {
assert(prime >= 2);
value %= prime;
if (value == 0 || prime == 2) return value;
if (internal::modular_square_root_power(value, (prime - 1) / 2, prime) != 1) {
return std::nullopt;
}
if (prime % 4 == 3) {
return internal::modular_square_root_power(value, prime / 4 + 1, prime);
}
uint64_t odd_part = prime - 1;
int power_of_two = 0;
while ((odd_part & 1) == 0) {
odd_part >>= 1;
power_of_two++;
}
uint64_t non_residue = 2;
while (internal::modular_square_root_power(non_residue, (prime - 1) / 2, prime) == 1) {
non_residue++;
}
uint64_t c = internal::modular_square_root_power(non_residue, odd_part, prime);
uint64_t root = internal::modular_square_root_power(value, odd_part / 2 + 1, prime);
uint64_t remainder = internal::modular_square_root_power(value, odd_part, prime);
int remaining_power = power_of_two;
while (remainder != 1) {
int exponent = 1;
uint64_t squared = internal::modular_square_root_multiply(remainder, remainder, prime);
while (squared != 1) {
squared = internal::modular_square_root_multiply(squared, squared, prime);
exponent++;
}
uint64_t correction = c;
for (int i = 0; i < remaining_power - exponent - 1; i++) {
correction = internal::modular_square_root_multiply(correction, correction, prime);
}
root = internal::modular_square_root_multiply(root, correction, prime);
c = internal::modular_square_root_multiply(correction, correction, prime);
remainder = internal::modular_square_root_multiply(remainder, c, prime);
remaining_power = exponent;
}
return root;
}
template <class Mint>
std::optional<Mint> modular_square_root(Mint value) {
auto root = modular_square_root(static_cast<uint64_t>(value.val()),
static_cast<uint64_t>(Mint::mod()));
if (!root.has_value()) return std::nullopt;
return Mint(*root);
}
} // namespace math
} // namespace m1une
#line 1 "math/fps/convolution.hpp"
#line 5 "math/fps/convolution.hpp"
#include <array>
#line 8 "math/fps/convolution.hpp"
#include <cstring>
#include <new>
#include <type_traits>
#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 13 "math/fps/formal_power_series.hpp"
namespace m1une {
namespace fps {
template <class Mint>
struct FormalPowerSeries : std::vector<Mint> {
using std::vector<Mint>::vector;
using Fps = FormalPowerSeries;
FormalPowerSeries() = default;
FormalPowerSeries(const std::vector<Mint>& values) : std::vector<Mint>(values) {}
FormalPowerSeries(std::vector<Mint>&& values) : std::vector<Mint>(std::move(values)) {}
Fps& shrink() {
while (!this->empty() && this->back() == Mint(0)) this->pop_back();
return *this;
}
Fps pre(int degree) const {
assert(degree >= 0);
Fps result(this->begin(), this->begin() + std::min<int>(degree, this->size()));
result.resize(degree);
return result;
}
Fps reversed(int size = -1) const {
Fps result = *this;
if (size >= 0) result.resize(size);
std::reverse(result.begin(), result.end());
return result;
}
Fps& operator+=(const Fps& rhs) {
if (this->size() < rhs.size()) this->resize(rhs.size());
for (int i = 0; i < int(rhs.size()); i++) (*this)[i] += rhs[i];
return *this;
}
Fps& operator-=(const Fps& rhs) {
if (this->size() < rhs.size()) this->resize(rhs.size());
for (int i = 0; i < int(rhs.size()); i++) (*this)[i] -= rhs[i];
return *this;
}
Fps& operator*=(const Fps& rhs) {
std::vector<Mint> lhs(this->begin(), this->end());
*this = convolution(lhs, rhs);
return *this;
}
Fps& operator*=(Mint rhs) {
for (Mint& value : *this) value *= rhs;
return *this;
}
Fps& operator/=(Mint rhs) {
return *this *= rhs.inv();
}
Fps& operator<<=(int shift) {
assert(shift >= 0);
this->insert(this->begin(), shift, Mint(0));
return *this;
}
Fps& operator>>=(int shift) {
assert(shift >= 0);
if (shift >= int(this->size())) {
this->clear();
} else {
this->erase(this->begin(), this->begin() + shift);
}
return *this;
}
Fps operator+() const {
return *this;
}
Fps operator-() const {
Fps result = *this;
for (Mint& value : result) value = Mint(0) - value;
return result;
}
friend Fps operator+(Fps lhs, const Fps& rhs) {
return lhs += rhs;
}
friend Fps operator-(Fps lhs, const Fps& rhs) {
return lhs -= rhs;
}
friend Fps operator*(Fps lhs, const Fps& rhs) {
return lhs *= rhs;
}
friend Fps operator*(Fps lhs, Mint rhs) {
return lhs *= rhs;
}
friend Fps operator*(Mint lhs, Fps rhs) {
return rhs *= lhs;
}
friend Fps operator/(Fps lhs, Mint rhs) {
return lhs /= rhs;
}
friend Fps operator<<(Fps lhs, int shift) {
return lhs <<= shift;
}
friend Fps operator>>(Fps lhs, int shift) {
return lhs >>= shift;
}
Fps derivative() const {
if (this->empty()) return {};
Fps result(this->size() - 1);
for (int i = 1; i < int(this->size()); i++) result[i - 1] = (*this)[i] * Mint(i);
return result;
}
Fps integral() const {
Fps result(this->size() + 1);
if (this->empty()) return result;
assert(this->size() < Mint::mod());
std::vector<Mint> inverse(this->size() + 1);
inverse[1] = 1;
for (int i = 2; i <= int(this->size()); i++) {
inverse[i] = Mint(0) - Mint(Mint::mod() / uint32_t(i)) * inverse[Mint::mod() % uint32_t(i)];
}
for (int i = 0; i < int(this->size()); i++) result[i + 1] = (*this)[i] * inverse[i + 1];
return result;
}
Mint evaluate(Mint x) const {
Mint result = 0;
for (auto it = this->rbegin(); it != this->rend(); ++it) result = result * x + *it;
return result;
}
Fps inv(int degree = -1) const {
if (degree < 0) degree = int(this->size());
assert(degree >= 0);
if (degree == 0) return {};
assert(!this->empty() && (*this)[0] != Mint(0));
Fps result(1, (*this)[0].inv());
for (int size = 1; size < degree; size <<= 1) {
const int next_size = std::min(size << 1, degree);
const int transform_size = size << 1;
if (size >= 32 && (Mint::mod() - 1) % uint32_t(transform_size) == 0) {
// Newton's g <- g(2-fg), restricted to the newly determined
// half. Keeping g in the frequency domain avoids two general
// convolutions and their 2x larger padding.
std::vector<Mint> transformed_f(transform_size);
std::copy_n(this->begin(), std::min<int>(this->size(), next_size),
transformed_f.begin());
std::vector<Mint> transformed_g(transform_size);
std::copy(result.begin(), result.end(), transformed_g.begin());
internal::ntt(transformed_f, false);
internal::ntt(transformed_g, false);
std::vector<Mint> error(transform_size);
for (int i = 0; i < transform_size; i++)
error[i] = transformed_f[i] * transformed_g[i];
internal::ntt(error, true);
std::fill(error.begin(), error.begin() + size, Mint(0));
internal::ntt(error, false);
for (int i = 0; i < transform_size; i++) error[i] *= transformed_g[i];
internal::ntt(error, true);
result.resize(next_size);
for (int i = size; i < next_size; i++) result[i] = Mint(0) - error[i];
continue;
}
Fps product = this->pre(next_size) * result;
product.resize(next_size);
for (Mint& value : product) value = Mint(0) - value;
product[0] += Mint(2);
result = (result * product).pre(next_size);
}
return result.pre(degree);
}
Fps log(int degree = -1) const {
if (degree < 0) degree = int(this->size());
assert(degree >= 0);
if (degree == 0) return {};
assert(!this->empty() && (*this)[0] == Mint(1));
return (derivative() * inv(degree)).pre(degree - 1).integral();
}
Fps exp(int degree = -1) const {
if (degree < 0) degree = int(this->size());
assert(degree >= 0);
if (degree == 0) return {};
assert(this->empty() || (*this)[0] == Mint(0));
Fps result(1, Mint(1));
for (int size = 1; size < degree; size <<= 1) {
const int next_size = std::min(size << 1, degree);
Fps correction = this->pre(next_size) - result.log(next_size);
correction[0] += Mint(1);
result = (result * correction).pre(next_size);
}
return result.pre(degree);
}
Fps pow(long long exponent, int degree = -1) const {
if (degree < 0) degree = int(this->size());
assert(exponent >= 0 && degree >= 0);
if (degree == 0) return {};
if (exponent == 0) {
Fps result(degree);
result[0] = 1;
return result;
}
int first = 0;
while (first < int(this->size()) && (*this)[first] == Mint(0)) first++;
if (first == int(this->size()) || first > (degree - 1) / exponent) return Fps(degree);
const int shift = int(first * exponent);
const Mint leading = (*this)[first];
Fps normalized = (*this >> first) / leading;
Fps result = (normalized.log(degree - shift) * Mint(exponent)).exp(degree - shift);
result *= leading.pow(exponent);
result <<= shift;
result.resize(degree);
return result;
}
std::optional<Fps> sqrt(int degree = -1) const {
if (degree < 0) degree = int(this->size());
assert(degree >= 0);
if (degree == 0) return Fps();
int first = 0;
while (first < int(this->size()) && (*this)[first] == Mint(0)) first++;
if (first == int(this->size())) return Fps(degree);
if (first >= degree) return Fps(degree);
if (first & 1) return std::nullopt;
const int shift = first / 2;
auto leading_root = m1une::math::modular_square_root((*this)[first]);
if (!leading_root.has_value()) return std::nullopt;
const int result_degree = degree - shift;
Fps normalized = (*this >> first) / (*this)[first];
Fps result = (normalized.log(result_degree) / Mint(2)).exp(result_degree);
result *= *leading_root;
result <<= shift;
result.resize(degree);
return result;
}
std::pair<Fps, Fps> divmod(const Fps& divisor) const {
Fps dividend = *this;
Fps normalized_divisor = divisor;
dividend.shrink();
normalized_divisor.shrink();
assert(!normalized_divisor.empty());
if (dividend.size() < normalized_divisor.size()) return std::make_pair(Fps(), dividend);
const int quotient_size = int(dividend.size() - normalized_divisor.size() + 1);
Fps quotient =
(dividend.reversed().pre(quotient_size) * normalized_divisor.reversed().inv(quotient_size))
.pre(quotient_size)
.reversed();
quotient.shrink();
Fps remainder = dividend - normalized_divisor * quotient;
remainder.resize(normalized_divisor.size() - 1);
remainder.shrink();
return std::make_pair(std::move(quotient), std::move(remainder));
}
Fps& operator/=(const Fps& rhs) {
*this = divmod(rhs).first;
return *this;
}
Fps& operator%=(const Fps& rhs) {
*this = divmod(rhs).second;
return *this;
}
friend Fps operator/(Fps lhs, const Fps& rhs) {
return lhs /= rhs;
}
friend Fps operator%(Fps lhs, const Fps& rhs) {
return lhs %= rhs;
}
Fps taylor_shift(Mint shift) const {
const int n = int(this->size());
if (n == 0) return {};
assert(uint32_t(n) < Mint::mod());
std::vector<Mint> factorial(n, Mint(1));
std::vector<Mint> inverse_factorial(n, Mint(1));
for (int i = 1; i < n; i++) factorial[i] = factorial[i - 1] * Mint(i);
inverse_factorial[n - 1] = factorial[n - 1].inv();
for (int i = n - 1; i > 0; i--) inverse_factorial[i - 1] = inverse_factorial[i] * Mint(i);
Fps left(n);
Fps right(n);
Mint power = 1;
for (int i = 0; i < n; i++) {
left[n - 1 - i] = (*this)[i] * factorial[i];
right[i] = power * inverse_factorial[i];
power *= shift;
}
Fps product = left * right;
Fps result(n);
for (int i = 0; i < n; i++) result[i] = product[n - 1 - i] * inverse_factorial[i];
return result;
}
};
} // namespace fps
} // namespace m1une
#line 1 "math/bernoulli.hpp"
#line 7 "math/bernoulli.hpp"
#line 1 "math/combinatorics.hpp"
#line 7 "math/combinatorics.hpp"
namespace m1une {
namespace math {
template <class Mint>
struct Combinatorics {
private:
std::vector<Mint> _factorial;
std::vector<Mint> _inverse_factorial;
public:
explicit Combinatorics(int maximum = 0) : _factorial(1, Mint(1)), _inverse_factorial(1, Mint(1)) {
ensure(maximum);
}
int maximum() const {
return int(_factorial.size()) - 1;
}
void ensure(int maximum) {
assert(maximum >= 0);
assert(static_cast<uint64_t>(maximum) < Mint::mod());
if (maximum <= this->maximum()) return;
const int old_maximum = this->maximum();
_factorial.resize(maximum + 1);
_inverse_factorial.resize(maximum + 1);
for (int i = old_maximum + 1; i <= maximum; i++) {
_factorial[i] = _factorial[i - 1] * Mint(i);
}
_inverse_factorial[maximum] = _factorial[maximum].inv();
for (int i = maximum; i > old_maximum; i--) {
_inverse_factorial[i - 1] = _inverse_factorial[i] * Mint(i);
}
}
Mint factorial(int n) const {
assert(0 <= n && n <= maximum());
return _factorial[n];
}
Mint inverse_factorial(int n) const {
assert(0 <= n && n <= maximum());
return _inverse_factorial[n];
}
Mint inverse(int n) const {
assert(1 <= n && n <= maximum());
return _factorial[n - 1] * _inverse_factorial[n];
}
Mint binom(int n, int k) const {
if (k < 0 || k > n) return Mint(0);
assert(n <= maximum());
return _factorial[n] * _inverse_factorial[k] * _inverse_factorial[n - k];
}
Mint perm(int n, int k) const {
if (k < 0 || k > n) return Mint(0);
assert(n <= maximum());
return _factorial[n] * _inverse_factorial[n - k];
}
Mint multiset(int types, int count) const {
if (types < 0 || count < 0) return Mint(0);
if (types == 0) return Mint(count == 0);
const long long total = static_cast<long long>(types) + count - 1;
assert(total <= maximum());
return binom(static_cast<int>(total), count);
}
Mint catalan(int n) const {
assert(n >= 0);
const long long doubled = 2LL * n;
assert(doubled <= maximum());
return binom(int(doubled), n) - binom(int(doubled), n + 1);
}
};
} // namespace math
} // namespace m1une
#line 10 "math/bernoulli.hpp"
namespace m1une {
namespace math {
namespace bernoulli_detail {
template <class Mint>
std::vector<Mint> numbers(
int maximum,
const Combinatorics<Mint>& combinations
) {
using Fps = fps::FormalPowerSeries<Mint>;
Fps denominator(maximum + 1);
for (int index = 0; index <= maximum; ++index) {
denominator[index] = combinations.inverse_factorial(index + 1);
}
Fps generating_function = denominator.inv(maximum + 1);
std::vector<Mint> result(maximum + 1);
for (int index = 0; index <= maximum; ++index) {
result[index] =
generating_function[index] * combinations.factorial(index);
}
return result;
}
template <class Mint>
Mint evaluate_polynomial(const std::vector<Mint>& coefficients, Mint x) {
Mint result = 0;
for (int index = int(coefficients.size()) - 1; index >= 0; --index) {
result = result * x + coefficients[index];
}
return result;
}
} // namespace bernoulli_detail
// Uses x / (exp(x) - 1), so B_1 = -1/2.
template <class Mint>
std::vector<Mint> bernoulli_numbers(int maximum) {
assert(maximum >= 0);
assert(static_cast<uint64_t>(maximum) + 1 < Mint::mod());
Combinatorics<Mint> combinations(maximum + 1);
return bernoulli_detail::numbers(maximum, combinations);
}
template <class Mint>
class Bernoulli {
public:
explicit Bernoulli(int maximum)
: combinations_(checked_maximum(maximum) + 1),
numbers_(bernoulli_detail::numbers(maximum, combinations_)) {}
int maximum() const {
return int(numbers_.size()) - 1;
}
const std::vector<Mint>& numbers() const {
return numbers_;
}
Mint number(int degree) const {
assert(0 <= degree && degree <= maximum());
return numbers_[degree];
}
// Coefficients of B_degree(x), in increasing order of powers of x.
std::vector<Mint> polynomial_coefficients(int degree) const {
assert(0 <= degree && degree <= maximum());
std::vector<Mint> result(degree + 1);
for (int power = 0; power <= degree; ++power) {
result[power] =
combinations_.binom(degree, power) *
numbers_[degree - power];
}
return result;
}
Mint polynomial(int degree, Mint x) const {
assert(0 <= degree && degree <= maximum());
std::vector<Mint> powers(degree + 1, Mint(1));
for (int power = 0; power < degree; ++power) {
powers[power + 1] = powers[power] * x;
}
Mint result = 0;
for (int index = 0; index <= degree; ++index) {
result += combinations_.binom(degree, index) * numbers_[index] *
powers[degree - index];
}
return result;
}
// Returns sum_{i=0}^{n-1} i^degree, evaluated as a polynomial in n.
Mint power_sum(Mint n, int degree) const {
assert(0 <= degree && degree <= maximum());
std::vector<Mint> powers(degree + 2, Mint(1));
for (int power = 0; power <= degree; ++power) {
powers[power + 1] = powers[power] * n;
}
Mint result = 0;
for (int index = 0; index <= degree; ++index) {
result += combinations_.binom(degree + 1, index) *
numbers_[index] * powers[degree + 1 - index];
}
return result * combinations_.inverse(degree + 1);
}
// Returns sum_{i=left}^{right-1} i^degree.
Mint power_sum(Mint left, Mint right, int degree) const {
return power_sum(right, degree) - power_sum(left, degree);
}
// Coefficients of sum_{i=0}^{n-1} i^degree as a polynomial in n.
std::vector<Mint> power_sum_polynomial(int degree) const {
assert(0 <= degree && degree <= maximum());
std::vector<Mint> result(degree + 2);
Mint inverse = combinations_.inverse(degree + 1);
for (int index = 0; index <= degree; ++index) {
result[degree + 1 - index] +=
combinations_.binom(degree + 1, index) * numbers_[index] *
inverse;
}
return result;
}
// If P is given by coefficients, returns coefficients of the unique Q
// with Q(0) = 0 and Q(n) = sum_{i=0}^{n-1} P(i).
std::vector<Mint> polynomial_prefix_sum(
const std::vector<Mint>& coefficients
) const {
if (coefficients.empty()) return std::vector<Mint>{Mint(0)};
int degree = int(coefficients.size()) - 1;
assert(degree <= maximum());
std::vector<Mint> result(degree + 2);
for (int source_degree = 0;
source_degree <= degree;
++source_degree) {
Mint inverse = combinations_.inverse(source_degree + 1);
for (int index = 0; index <= source_degree; ++index) {
result[source_degree + 1 - index] +=
coefficients[source_degree] *
combinations_.binom(source_degree + 1, index) *
numbers_[index] * inverse;
}
}
return result;
}
// Returns sum_{i=left}^{right-1} P(i).
Mint polynomial_sum(
const std::vector<Mint>& coefficients,
Mint left,
Mint right
) const {
std::vector<Mint> prefix = polynomial_prefix_sum(coefficients);
return bernoulli_detail::evaluate_polynomial(prefix, right) -
bernoulli_detail::evaluate_polynomial(prefix, left);
}
// Returns sum_{i=0}^{count-1} (start + step*i)^degree.
Mint arithmetic_progression_power_sum(
Mint start,
Mint step,
Mint count,
int degree
) const {
assert(0 <= degree && degree <= maximum());
std::vector<Mint> start_powers(degree + 1, Mint(1));
std::vector<Mint> step_powers(degree + 1, Mint(1));
for (int power = 0; power < degree; ++power) {
start_powers[power + 1] = start_powers[power] * start;
step_powers[power + 1] = step_powers[power] * step;
}
Mint result = 0;
for (int power = 0; power <= degree; ++power) {
result += combinations_.binom(degree, power) *
start_powers[degree - power] * step_powers[power] *
power_sum(count, power);
}
return result;
}
private:
static int checked_maximum(int maximum) {
assert(maximum >= 0);
assert(static_cast<uint64_t>(maximum) + 1 < Mint::mod());
return maximum;
}
Combinatorics<Mint> combinations_;
std::vector<Mint> numbers_;
};
} // namespace math
} // namespace m1une
#line 1 "math/partition_function.hpp"
#line 5 "math/partition_function.hpp"
#line 8 "math/partition_function.hpp"
namespace m1une {
namespace math {
// Returns p(0), p(1), ..., p(maximum), where p(n) is the number of integer
// partitions of n.
template <class Mint>
std::vector<Mint> partition_function(int maximum) {
assert(maximum >= 0);
using Fps = fps::FormalPowerSeries<Mint>;
Fps denominator(maximum + 1);
denominator[0] = 1;
for (long long k = 1;; k++) {
long long first = k * (3 * k - 1) / 2;
long long second = k * (3 * k + 1) / 2;
if (first > maximum) break;
Mint sign = (k & 1) ? Mint(-1) : Mint(1);
denominator[int(first)] += sign;
if (second <= maximum) denominator[int(second)] += sign;
}
return denominator.inv(maximum + 1);
}
template <class Mint>
std::vector<Mint> partition_numbers(int maximum) {
return partition_function<Mint>(maximum);
}
} // namespace math
} // namespace m1une
#line 12 "math/combinatorial_sequences.hpp"
namespace m1une {
namespace math {
template <class Mint>
std::vector<Mint> catalan_numbers(int maximum) {
assert(maximum >= 0);
assert(static_cast<uint64_t>(maximum) + 1 < Mint::mod());
std::vector<Mint> inverse(maximum + 2);
inverse[1] = 1;
for (int i = 2; i <= maximum + 1; i++) {
inverse[i] = Mint(0) - Mint(Mint::mod() / uint32_t(i)) * inverse[Mint::mod() % uint32_t(i)];
}
std::vector<Mint> result(maximum + 1);
result[0] = 1;
for (int n = 0; n < maximum; n++) {
result[n + 1] = result[n] * Mint(2) * Mint(2LL * n + 1) * inverse[n + 2];
}
return result;
}
template <class Mint>
std::vector<Mint> bell_numbers(int maximum) {
assert(maximum >= 0);
assert(static_cast<uint64_t>(maximum) < Mint::mod());
using Fps = fps::FormalPowerSeries<Mint>;
Combinatorics<Mint> combinations(maximum);
Fps exponent(maximum + 1);
for (int i = 1; i <= maximum; i++) {
exponent[i] = combinations.inverse_factorial(i);
}
Fps generating_function = exponent.exp(maximum + 1);
std::vector<Mint> result(maximum + 1);
for (int i = 0; i <= maximum; i++) {
result[i] = generating_function[i] * combinations.factorial(i);
}
return result;
}
template <class Mint>
std::vector<Mint> stirling_numbers_second_kind(int n) {
assert(n >= 0);
assert(static_cast<uint64_t>(n) < Mint::mod());
Combinatorics<Mint> combinations(n);
std::vector<Mint> powers(n + 1);
std::vector<Mint> signs(n + 1);
for (int i = 0; i <= n; i++) {
powers[i] = Mint(i).pow(n) * combinations.inverse_factorial(i);
signs[i] = combinations.inverse_factorial(i);
if (i & 1) signs[i] = Mint(0) - signs[i];
}
std::vector<Mint> result = fps::convolution(powers, signs);
result.resize(n + 1);
return result;
}
template <class Mint>
std::vector<Mint> derangement_numbers(int maximum) {
assert(maximum >= 0);
std::vector<Mint> result(maximum + 1);
result[0] = 1;
if (maximum >= 1) result[1] = 0;
for (int n = 2; n <= maximum; n++) {
result[n] = Mint(n - 1) * (result[n - 1] + result[n - 2]);
}
return result;
}
} // namespace math
} // namespace m1une