OJ-Problems-Source/.ACM-Templates/组合数学.cpp

//组合数预处理 / 杨辉三角 O(n^2)
//C[i][i] = C[i][0] = 1
//C[i][j] = C[i - 1][j] + C[i - 1][j - 1], 0 < j < i
const int maxc = 105;
ll C[maxc][maxc];
void calC() {
  for (int i = 0; i < maxc; i++) { C[i][i] = C[i][0] = 1; }
  for (int i = 2; i < maxc; i++) {
    for (int j = 1; j < i; j++) { C[i][j] = (C[i - 1][j] + C[i - 1][j - 1]) % M; }
  }
}
//求组合数不取模 O(n)
ll Com(ll n, ll m) {
  if (m > n) { return 0; }
  if (m > n - m) { m = n - m; }
  ll ret = 1;
  for (ll i = 0, j = 1; i < m; i++) {
    for (ret *= n - i; j <= m && ret % j == 0; j++) { ret /= j; }
  }
  return ret;
}
//求组合数取模
//p <= 1e6 预处理阶乘逆元 O(min(n, p)) + O(1)
ll fac[M] = {1, 1}, invfac[M] = {1, 1};
void initFac(ll p) {
  for (int i = 2; i < p; i++) {
    fac[i] = fac[i - 1] * i % p; invfac[i] = (-invfac[p % i] * (p / i) % p + p) % p;
  }
  for (int i = 2; i < p; i++) { invfac[i] = invfac[i] * invfac[i - 1] % p; }
}
ll Com(ll n, ll m, ll p) {
  return n < m ? 0 : fac[n] * invfac[n - m] % p * invfac[m] % p;
}
//p <= 1e9 在线求逆元 O(nlogp)
ll Com(ll n, ll m, ll p) {
  if (m > n) { return 0; }
  if (m > n - m) { m = n - m; }
  ll ret = 1;
  for (ll i = 1; i <= m; i++) {
    ll a = (n + i - m) % p, b = i % p;
    ret = ret * a % p * powMod(b, p - 2, p) % p;
  }
  return ret;
}
//Lucas定理
ll Lucas(ll n, ll m, ll p) {
  if (n < m) { return 0; }
  if (m == 0 && n == 0) { return 1; }
  return Com(n % p, m % p, p) * Lucas(n / p, m / p, p) % p;
}
//第一类Stirling数 s(p, k)
//将p个物体排成k个非空循环排列的方法数
//s(p, k)的递推公式：s(p, k) = (p - 1) * s(p - 1, k) + s(p - 1, k - 1), 1 <= k <= p - 1
//边界条件：s(p, 0) = 0, p >= 1, s(p, p) = 1, p >= 0
const int maxs = 105;
ll S[maxs][maxs];
void calStir1() {
  S[0][0] = S[1][1] = 1;
  for (int i = 2; i < maxs; i++) {
    for (int j = 1; j <= i; j++) { S[i][j] = ((i - 1) * S[i - 1][j] + S[i - 1][j - 1]) % M; }
  }
}
//第二类Stirling数 S(p, k)
//将p个物体划分成k个非空的不可辨别的(可以理解为盒子没有编号)集合的方法数
//k! * S(p, k)是把p个人分进k间有差别(如被标有房号)的房间(无空房)的方法数
//S(p, k)的递推公式是：S(p, k) = k * S(p - 1, k) + S(p - 1, k - 1), 1 <= k <= p - 1
//边界条件：S(p, 0) = 0, p >= 1, S(p, p) = 1, p >= 0
const int maxs = 105;
ll S[maxs][maxs];
void calStir2() {
  S[0][0] = S[1][1] = 1;
  for (int i = 2; i < maxs; i++) {
    for (int j = 1; j <= i; j++) { S[i][j] = (j * S[i - 1][j] + S[i - 1][j - 1]) % M; }
  }
}
//Bell数
//B(n)表示基数为n的集合的划分方法的数目
//B(0) = 1, B(n + 1) = sum(C(n, k) * B(k)), 0 <= k <= n
//每个贝尔数都是第二类Stirling数的和, 即B(n) = sum(S(n, k)), 1 <= k <= n
//Bell三角形
//a[0][0] = 1
//对于n >= 1, a[n][0] = a[n - 1][n - 1]
//对于m, n >= 1, a[n][m] = a[n][m - 1] + a[n - 1][m - 1]
//每行首项是贝尔数，每行之和是第二类Stirling数
//两个重要的同余性质:
//B(p + n) = B(n) + B(n + 1) (mod p)
//B(p^m + n) = m * B(n) + B(n + 1) (mod p)
//p是不大于100的素数, 这样, 我们可以通过上面的性质来计算Bell数模小于100的素数值
//Bell数模素数p的周期为: N(p) = ((p^p) - 1) / (p - 1)
const int maxb = 105;
ll T[maxb], B[maxb];
void calBell() {
  B[0] = B[1] = T[0] = 1;
  for (int i = 2; i < maxb; i++) {
    T[i - 1] = B[i - 1];
    for (int j = i - 2; j >= 0; j--) { T[j] = (T[j] + T[j + 1]) % M; }
    B[i] = T[0];
  }
}
//计算Bell(n)对999999598 = 2 × 13 × 5281 × 7283取模 O(P^2logP)
const int N = 7284, P = 999999598;
ll n; int p[5] = {2, 13, 5281, 7283}, B[2][N], T[N];
void init() {
  T[0] = T[1] = B[0][0] = 1; B[0][1] = 2;
  for (int i = 2, crt = 1; i < N; i++, crt ^= 1) {
    T[i] = B[crt][0] = B[crt ^ 1][i - 1];
    for (int j = 1; j <= i; j++) { B[crt][j] = (B[crt ^ 1][j - 1] + B[crt][j - 1]) % P; }
  }
}
int b[N], c[N], d[70];
int cal(ll n, int mod) {
  int len = 0;
  for (int i = 0; i <= mod; i++) { b[i] = T[i] % mod; }
  while (n) { d[len++] = n % mod; n /= mod; }
  for (int i = 1; i < len; i++) {
    for (int j = 1; j <= d[i]; j++) {
      for (int k = 0; k < mod; k++) { c[k] = (b[k] * i + b[k + 1]) % mod; }
      c[mod] = (c[0] + c[1]) % mod;
      for (int k = 0; k <= mod; k++) { b[k] = c[k]; }
    }
  }
  return c[d[0]];
}
ll bell(ll n) {
  if (n < N) { return T[n]; }
  ll t = 0;
  for (int i = 0; p[i]; i++) {
    t = (t + (P / p[i]) * powMod(P / p[i], p[i] - 2, p[i]) % P * cal(n, p[i]) % P) % P;
  }
  return t;
}
//卡特兰数 Catalan Number
//Cat(1) = 1, Cat(n) = (4n - 2) * Cat(n - 1) / (n + 1) = C(2n, n) / (n + 1) = C(2n, n) - C(2n, n - 1)
//从(0, 0)点走到(n, m)点且
//不经过对角线的方法数(x > y): C(n + m - 1, m) - C(n + m - 1, m - 1)
//不穿过对角线的方法数(x >= y): C(n + m, m) - C(n + m, m - 1)
//预处理卡特兰数
int a[105][105], b[105]; //大数, 长度
void Catalan() {
  int i, j, len, carry, temp;
  a[1][0] = b[1] = len = 1;
  for (i = 2; i <= 100; i++) {
    for (j = 0; j < len; j++) { //乘法
      a[i][j] = a[i - 1][j] * (4 * (i - 1) + 2);
    }
    carry = 0;
    for (j = 0; j < len; j++) { //处理相乘结果
      temp = a[i][j] + carry;
      a[i][j] = temp % 10;
      carry = temp / 10;
    }
    while (carry) { //进位处理
      a[i][len++] = carry % 10;
      carry /= 10;
    }
    carry = 0;
    for (j = len - 1; j >= 0; j--) { //除法
      temp = carry * 10 + a[i][j];
      a[i][j] = temp / (i + 1);
      carry = temp % (i + 1);
    }
    while (!a[i][len - 1]) { len--; } //高位零处理
    b[i] = len;
  }
}
//输出大数
void printCatalan(int n) {
  for (int i = b[n] - 1; i >= 0; i--) { printf("%d", a[n][i]); }
}
//各种情况下小球放盒子的方案数
//k个球 m个盒子 是否允许有空盒子 方案数
//不同 不同 是 m^k
//不同 不同 否 m!*Stirling2(k, m)
//不同 相同 是 ∑(m, i=1)Stirling2(k, i)
//不同 相同 否 Stirling2(k, m)
//相同 不同 是 C(m + k - 1, k)
//相同 不同 否 C(k - 1, m - 1)
//相同 相同 是 1/(1-x)(1-x^2)...(1-x^m)的x^k项的系数
//相同 相同 否 x^m/(1-x)(1-x^2)...(1-x^m)的x^k项的系数
//错排公式
//D(1) = 0, D(2) = 1, D(n) = (n - 1)(D(n - 2) + D(n - 1))
//扩展 Cayley 公式
//对于n个点, m个连通块的图, 假设每个连通块有a[i]个点, 那么用s - 1条边把它连通的方案数为n^(s-2)a[1]a[2]...a[m]