add

Author: KumaCS

convolution/intmod.hpp

@@ -5,11 +5,11 @@
 
 namespace ConvolutionIntMod {
 using ll = long long;
-const ll Mod1 = 754974721;
-const ll Mod2 = 167772161;
-const ll Mod3 = 469762049;
-const ll M1invM2 = 95869806;
-const ll M12invM3 = 187290749;
+static constexpr ll Mod1 = 754974721;
+static constexpr ll Mod2 = 167772161;
+static constexpr ll Mod3 = 469762049;
+static constexpr ll M1invM2 = 95869806;
+static constexpr ll M12invM3 = 187290749;
 using M1 = ModInt<Mod1>;
 using M2 = ModInt<Mod2>;
 using M3 = ModInt<Mod3>;

convolution/mod2_64.hpp

@@ -0,0 +1,94 @@
+#pragma once
+
+#include "modint/modint.hpp"
+#include "fft/ntt.hpp"
+#include "math/util.hpp"
+
+namespace ConvolutionMod2_64 {
+using ull = unsigned long long;
+static constexpr ull M1 = 645922817;
+static constexpr ull M2 = 754974721;
+static constexpr ull M3 = 880803841;
+static constexpr ull M4 = 897581057;
+static constexpr ull M5 = 998244353;
+static constexpr ull M12M4 = M1 * M2 % M4;
+static constexpr ull M12M5 = M1 * M2 % M5;
+static constexpr ull M123M5 = M12M5 * M3 % M5;
+static constexpr ull M12 = M1 * M2;
+static constexpr ull M123 = M12 * M3;
+static constexpr ull M1234 = M123 * M4;
+static constexpr ull I2 = Math::inv_mod(M1, M2);
+static constexpr ull I3 = Math::inv_mod(M1 * M2 % M3, M3);
+static constexpr ull I4 = Math::inv_mod(M1 * M2 % M4 * M3 % M4, M4);
+static constexpr ull I5 = Math::inv_mod(M1 * M2 % M5 * M3 % M5 * M4 % M5, M5);
+
+using mint1 = ModInt<M1>;
+using mint2 = ModInt<M2>;
+using mint3 = ModInt<M3>;
+using mint4 = ModInt<M4>;
+using mint5 = ModInt<M5>;
+
+NTT<mint1> ntt1;
+NTT<mint2> ntt2;
+NTT<mint3> ntt3;
+NTT<mint4> ntt4;
+NTT<mint5> ntt5;
+
+template <class mint>
+vector<mint> inner_mult(const vector<ull>& a, const vector<ull>& b, NTT<mint>& ntt) {
+  constexpr unsigned int mod = mint::get_mod();
+  vector<mint> a1(a.size()), b1(b.size());
+  for (int i = 0; i < a.size(); i++) a1[i] = a[i] % mod;
+  for (int i = 0; i < b.size(); i++) b1[i] = b[i] % mod;
+  mint c = ntt.multiply(a1, b1)[0];
+  return ntt.multiply(a1, b1);
+}
+template <class mint>
+vector<mint> inner_middle_prod(const vector<ull>& a, const vector<ull>& b, NTT<mint>& ntt) {
+  constexpr unsigned int mod = mint::get_mod();
+  vector<mint> a1(a.size()), b1(b.size());
+  for (int i = 0; i < a.size(); i++) a1[i] = a[i] % mod;
+  for (int i = 0; i < b.size(); i++) b1[i] = b[i] % mod;
+  return ntt.middle_product(a1, b1);
+}
+vector<ull> multiply(const vector<ull>& a, const vector<ull>& b) {
+  if (a.empty() || b.empty()) return {};
+  auto c1 = inner_mult(a, b, ntt1);
+  auto c2 = inner_mult(a, b, ntt2);
+  auto c3 = inner_mult(a, b, ntt3);
+  auto c4 = inner_mult(a, b, ntt4);
+  auto c5 = inner_mult(a, b, ntt5);
+  vector<ull> c(a.size() + b.size() - 1, 0);
+  for (int i = 0; i < c.size(); i++) {
+    ull y1 = c1[i].val();
+    ull y2 = (c2[i].val() + M2 - y1) * I2 % M2;
+    ull y3 = (c3[i].val() + M3 - (y1 + y2 * M1) % M3) * I3 % M3;
+    ull y4 = (c4[i].val() + M4 - (y1 + y2 * M1 + y3 * M12M4) % M4) * I4 % M4;
+    ull y5 = (c5[i].val() + M5 - (y1 + y2 * M1 + y3 * M12M5 + y4 * M123M5) % M5) * I5 % M5;
+    c[i] = y1 + y2 * M1 + y3 * M12 + y4 * M123 + y5 * M1234;
+  }
+  return c;
+}
+vector<ull> middle_product(const vector<ull>& a, const vector<ull>& b) {
+  if (b.empty() || a.size() > b.size()) return {};
+  auto c1 = inner_middle_prod(a, b, ntt1);
+  auto c2 = inner_middle_prod(a, b, ntt2);
+  auto c3 = inner_middle_prod(a, b, ntt3);
+  auto c4 = inner_middle_prod(a, b, ntt4);
+  auto c5 = inner_middle_prod(a, b, ntt5);
+  vector<ull> c(c1.size(), 0);
+  for (int i = 0; i < c.size(); i++) {
+    ull y1 = c1[i].val();
+    ull y2 = (c2[i].val() + M2 - y1) * I2 % M2;
+    ull y3 = (c3[i].val() + M3 * 2 - (y1 + y2 * M1 % M3)) * I3 % M3;
+    ull y4 = (c4[i].val() + M4 * 3 - (y1 + y2 * M1 + y3 * M12M4) % M4) * I4 % M4;
+    ull y5 = (c5[i].val() + M5 * 4 - (y1 + y2 * M1 + y3 * M12M5 + y4 * M123M5) % M5) * I5 % M5;
+    c[i] = y1 + y2 * M1 + y3 * M12 + y4 * M123 + y5 * M1234;
+  }
+  return c;
+}
+};  // namespace ConvolutionMod2_64
+
+/**
+ * @brief 畳み込み mod 2^64
+ */

data-structure/integer-set.hpp

@@ -3,7 +3,7 @@
 template <class T = int, class U = unsigned long long>
 class IntegerSet {
  private:
-  const static T B = 6, W = 64, MASK = W - 1;
+  static constexpr T B = 6, W = 64, MASK = W - 1;
   T size;
   vector<T> start;
   vector<U> data;

docs/fps/famous-sequences.md

@@ -8,6 +8,22 @@
 $$\prod_{k=1}^{\infty}(1-x^k)=\sum_{k=-\infty}^{\infty}(-1)^kx^{k(3k-1)/2}$$
 であることを用いれば $O(N\log N)$ 時間で $p_0,p_1,\dots,p_N$ が列挙できる．
 
+## ベル数
+
+$n$ 元集合を空でない部分集合に分割する方法の数 $B_n$ をベル数という．
+指数型母関数が下のように計算できる．
+$$\sum_{n}\frac{B_n}{n!}x^n
+=\prod_{i=1}^{\infty}\sum_{j=0}^{\infty}\frac{1}{(i!)^jj!}x^{ij}
+=\prod_{i=1}^{\infty}\exp\left(\frac{x^i}{i!}\right)
+=\exp\left(\sum_{i=1}^{\infty}\frac{x^i}{i!}\right)
+=\exp(e^x-1)$$
+
+## モンモール数
+
+長さ $n$ の撹乱順列，すなわち $(1,2,\dots,n)$ の順列 $(p_1,p_2,\dots,p_n)$ で $p_i\neq i$ を満たすものの個数．
+
+$a_n$ とおくと $a_n=(n-1)(a_{n-1}+a_{n-2})$ であるから $O(n)$ 時間で列挙できる．
+
 ## 第一種スターリング数
 
 $s(n,k)$ を以下で定める．

fps/famous-sequences.hpp

@@ -12,6 +12,22 @@ FormalPowerSeries<mint> PartitionFunction(int n) {
   return g.inv(n + 1);
 }
 template <class mint>
+FormalPowerSeries<mint> BellNumber(int n) {
+  using fact = Factorial<mint>;
+  FormalPowerSeries<mint> f(n + 1);
+  for (int i = 1; i < f.size(); i++) f[i] = fact::fact_inv(i);
+  f = f.exp();
+  for (int i = 0; i < f.size(); i++) f[i] *= fact::fact(i);
+  return f;
+}
+template <class mint>
+vector<mint> MontmortNumber(int n) {
+  vector<mint> f(n + 1);
+  f[0] = 1, f[1] = 0;
+  for (int i = 2; i < f.size(); i++) f[i] = (i - 1) * (f[i - 1] + f[i - 2]);
+  return f;
+}
+template <class mint>
 FormalPowerSeries<mint> FirstKindStirlingNumbers(int n) {
   FormalPowerSeries<mint> f{1};
   for (int l = 30; l >= 0; l--) {

fps/formal-power-series.hpp

@@ -94,8 +94,10 @@ struct FormalPowerSeries : vector<mint> {
   }
   FPS operator>>=(int sz) {
     assert(sz >= 0);
-    if ((int)this->size() <= sz) return {};
-    this->erase(this->begin(), this->begin() + sz);
+    if ((int)this->size() <= sz)
+      this->clear();
+    else
+      this->erase(this->begin(), this->begin() + sz);
     return *this;
   }
   FPS operator>>(int sz) const {

fps/product-of-polynomials.hpp

@@ -7,7 +7,7 @@ FormalPowerSeries<mint> ProductOfPolynomials(vector<FormalPowerSeries<mint>> fs)
   for (auto& f : fs) f.shrink();
   for (auto& f : fs)
     if (f.empty()) return {};
-  const int B = 1 << 6;
+  static constexpr int B = 1 << 6;
   for (int i = 0, j = -1; i < fs.size(); i++) {
     if (fs[i].size() > B) continue;
     if (j == -1 || fs[i].size() + fs[j].size() - 1 > B) {

fps/relaxed.hpp

@@ -5,7 +5,7 @@
 
 template <class mint>
 class RelaxedMultiply {
-  const int B = 6;
+  static constexpr int B = 6;
   using fps = FormalPowerSeries<mint>;
   int n;
   fps f, g, h;
@@ -61,7 +61,7 @@ class RelaxedMultiply {
 
 template <class mint>
 class SemiRelaxedMultiply {
-  const int B = 6;
+  static constexpr int B = 6;
   using fps = FormalPowerSeries<mint>;
   int n, m0;
   fps f, g, h;

fps/sum-of-rationals.hpp

@@ -4,12 +4,12 @@
 
 template <class mint>
 FPSRational<mint> SumOfRationals(vector<FPSRational<mint>> rs) {
-  if (ps.empty()) return {};
+  if (rs.empty()) return {};
   for (auto& r : rs) {
     r.num.shrink(), r.den.shrink();
     if (r.den.size() < r.num.size()) r.num.resize(r.den.size());
   }
-  const int B = 1 << 5;
+  static constexpr int B = 1 << 5;
   for (int i = 0, j = -1; i < rs.size(); i++) {
     if (rs[i].den.size() > B) continue;
     if (j == -1 || rs[i].den.size() + rs[j].den.size() - 1 > B) {
@@ -23,7 +23,7 @@ FPSRational<mint> SumOfRationals(vector<FPSRational<mint>> rs) {
   if (rs.size() == 1) return rs[0];
   for (auto& r : rs) {
     int sz = B;
-    while (sz < r.size()) sz <<= 1;
+    while (sz < r.num.size() || sz < r.den.size()) sz <<= 1;
     r.num.resize(sz);
     r.num.ntt();
     r.den.resize(sz);

math/garner-online.hpp

@@ -9,7 +9,7 @@ struct GarnerOnline {
   void push(T r, T m) {
     T x = 0, p = 1;
     for (int j = 0; j < ms.size(); j++) {
-      x = (x + y[j] * p) % ms[i];
+      x = (x + y[j] * p) % m;
       p = p * ms[j] % m;
     }
     ms.push_back(m);
@@ -23,12 +23,12 @@ struct GarnerOnline {
   }
   T get() {
     T res = 0;
-    for (int i = (int)ms.size() - 1; i >= 0; i--) res = res * m[i] + y[i];
+    for (int i = (int)ms.size() - 1; i >= 0; i--) res = res * ms[i] + y[i];
     return res;
   }
   T get(T mod) {
     T res = 0;
-    for (int i = (int)ms.size() - 1; i >= 0; i--) res = (res * m[i] + y[i]) % mod;
+    for (int i = (int)ms.size() - 1; i >= 0; i--) res = (res * ms[i] + y[i]) % mod;
     return res;
   }
 };