伯努利数（详解 + 例题：P3711 仓鼠的数学题）

2024-03-20 16:21:28

伯努利数

定义 S k ( n ) = ∑ i = 0 n − 1 i k S_k(n) = \sum\limits_{i = 0} ^{n - 1} i ^ k Sk(n)=i=0∑n−1ik。

从二项式出发

( 0 + 1 ) k + 1 = ∑ i = 0 k C k + 1 i 0 i + 0 k + 1 ⋮ ( n − 1 + 1 ) k + 1 = ∑ i = 0 k C k + 1 i ( n − 1 ) i + ( n − 1 ) k + 1 把次方 k + 1 的移项，再整体相加，得 n k + 1 = ∑ i = 0 k C k + 1 i S i ( n ) n k + 1 = ∑ i = 0 k − 1 C k + 1 i S i ( n ) + ( k + 1 ) S k ( n ) S k ( n ) = 1 k + 1 ( n k + 1 − ∑ i = 0 k − 1 C k + 1 i S i ( n ) ) (0 + 1) ^ {k + 1} = \sum_{i = 0} ^{k} C_{k + 1} ^{i} 0 ^ {i} + 0 ^{k + 1}\\ \vdots\\ (n - 1 + 1) ^{k + 1} = \sum_{i = 0} ^{k} C_{k + 1} ^{i} (n - 1) ^{i} + (n - 1) ^{k + 1}\\ 把次方k + 1的移项，再整体相加，得\\ n ^{k + 1} = \sum_{i = 0} ^{k} C_{k + 1} ^{i} S_i(n)\\ n ^{k + 1} = \sum_{i = 0} ^{k - 1} C_{k + 1} ^{i} S_i(n) + (k + 1) S_k(n)\\ S_k(n) = \frac{1}{k + 1}\left(n ^{k + 1} - \sum_{i = 0} ^{k - 1} C_{k + 1} ^{i} S_i(n) \right)\\ (0+1)k+1=i=0∑kCk+1i0i+0k+1⋮(n−1+1)k+1=i=0∑kCk+1i(n−1)i+(n−1)k+1把次方k+1的移项，再整体相加，得nk+1=i=0∑kCk+1iSi(n)nk+1=i=0∑k−1Ck+1iSi(n)+(k+1)Sk(n)Sk(n)=k+11(nk+1−i=0∑k−1Ck+1iSi(n))
于是我们有了一个 O ( n 2 ) O(n ^ 2) O(n2)递推求 S i ( n ) S_i(n) Si(n)得算法，

对二项式再进一步考虑：

( 0 + m ) k = ∑ i = 0 k C k i 0 i m k − i ⋮ ( n − 1 + m ) k = ∑ i = 0 k C k i ( n − 1 ) i m k − i 左右两边同时相加 S k ( n + m ) − S k ( m ) = ∑ i = 0 k C k i S i ( n ) m k − i 左右两边对 n 求导 S k ′ ( n + m ) = ∑ i = 0 k C k i S i ′ ( n ) m k − i n = 0 ，代入 S k ′ ( m ) = ∑ i = 0 k C k i S i ′ ( 0 ) m k − i 两边对 m 从 0 − > n 求积分 S k ( n ) − S k ( 0 ) = ∑ i = 0 k C k i S i ′ ( 0 ) n k − i + 1 k − i + 1 − ∑ i = 0 k C k i S i ′ ( 0 ) 0 k − i + 1 k − i + 1 S k ( 0 ) = 0 , 对于最右侧 i ≠ k + 1 , 所以 0 k − i + 1 = 0 S k ( n ) = ∑ i = 0 k C k i S i ′ ( 0 ) n k − i + 1 k − i + 1 = 1 k + 1 ∑ i = 0 k C k + 1 i S i ′ ( 0 ) n k + i − 1 (0 + m) ^{k} = \sum_{i = 0} ^{k} C_{k} ^{i} 0 ^{i} m ^{k - i}\\ \vdots\\ (n - 1 + m) ^{k} = \sum_{i = 0} ^{k} C_{k} ^{i}(n - 1) ^{i} m ^{k - i}\\ 左右两边同时相加\\ S_k(n + m) - S_k(m) = \sum_{i = 0} ^{k} C_{k} ^{i} S_{i} (n) m ^{k - i}\\ 左右两边对n求导\\ S'_k(n + m) = \sum_{i = 0} ^{k} C_{k} ^{i} S'_i(n) m^{k - i}\\ \\ n = 0，代入\\ S'_k(m) = \sum_{i = 0} ^{k} C_{k} ^{i} S'_i(0) m^{k - i}\\ 两边对m从0->n求积分\\ S_k(n) - S_k(0) = \sum_{i = 0} ^{k} C_{k} ^{i} S'_i(0) \frac{n ^{k - i + 1}}{k - i + 1} - \sum_{i = 0} ^{k} C_{k} ^{i} S'_i(0) \frac{0 ^{k - i + 1}}{k - i + 1}\\ S_k(0) = 0, 对于最右侧i \neq k + 1,所以0 ^{k - i + 1} = 0 \\ S_k(n) = \sum_{i = 0} ^{k} C_{k} ^{i} S'_i(0) \frac{n ^{k - i + 1}}{k - i + 1} = \frac{1}{k + 1} \sum_{i = 0} ^{k} C_{k + 1} ^{i} S'_i(0) n ^{k + i - 1}\\ (0+m)k=i=0∑kCki0imk−i⋮(n−1+m)k=i=0∑kCki(n−1)imk−i左右两边同时相加Sk(n+m)−Sk(m)=i=0∑kCkiSi(n)mk−i左右两边对n求导Sk′(n+m)=i=0∑kCkiSi′(n)mk−in=0，代入Sk′(m)=i=0∑kCkiSi′(0)mk−i两边对m从0−>n求积分Sk(n)−Sk(0)=i=0∑kCkiSi′(0)k−i+1nk−i+1−i=0∑kCkiSi′(0)k−i+10k−i+1Sk(0)=0,对于最右侧i=k+1,所以0k−i+1=0Sk(n)=i=0∑kCkiSi′(0)k−i+1nk−i+1=k+11i=0∑kCk+1iSi′(0)nk+i−1
其实，这里的 S i ( 0 ) S_i(0) Si(0)就是伯努利数，为了方便书写， B i = S i ′ ( 0 ) B_i = S'_i(0) Bi=Si′(0)。
我们先将 n = 1 代入上式 S k ( 1 ) = 1 k + 1 ∑ i = 0 k C k + 1 i B i 同样的我们把右边的最高项给提出来 S k ( 1 ) = 1 k + 1 ∑ i = 0 k − 1 C k + 1 i B i + B k 则有 B k − S k ( 1 ) = − 1 k + 1 ∑ i = 0 k − 1 C k + 1 i B i 特殊地，考虑 k = 0 , 有 S 0 ( 1 ) = 1 , 所以 B 0 = 1 k ≥ 1 , S k ( 1 ) = 0 , 所以 B k = − 1 k + 1 ∑ i = 0 k − 1 C k + 1 i B i 我们先将n = 1代入上式\\ S_k(1) = \frac{1}{k + 1} \sum_{i = 0} ^{k} C_{k + 1} ^{i} B_i\\ 同样的我们把右边的最高项给提出来\\ S_k(1) = \frac{1}{k + 1} \sum_{i = 0} ^{k - 1} C_{k + 1} ^{i} B_i + B_k\\ 则有B_k - S_k(1) = - \frac{1}{k + 1} \sum_{i = 0} ^{k - 1} C_{k + 1} ^{i} B_i\\ 特殊地，考虑k = 0,有S_0(1) = 1, 所以B_0 = 1\\ k \geq 1, S_k(1) = 0, 所以B_k = - \frac{1}{k + 1} \sum_{i = 0} ^{k - 1} C_{k + 1} ^{i} B_i\\ 我们先将n=1代入上式Sk(1)=k+11i=0∑kCk+1iBi同样的我们把右边的最高项给提出来Sk(1)=k+11i=0∑k−1Ck+1iBi+Bk则有Bk−Sk(1)=−k+11i=0∑k−1Ck+1iBi特殊地，考虑k=0,有S0(1)=1,所以B0=1k≥1,Sk(1)=0,所以Bk=−k+11i=0∑k−1Ck+1iBi
由此我们可以 O ( k 2 ) O(k ^ 2) O(k2)递推出伯努利数出来了。

F a u l h a b e r Faulhaber Faulhaber公式

S k ( x ) = 1 k + 1 ∑ i = 0 k C k + 1 i B i x k − i + 1 S_k(x) = \frac{1}{k + 1} \sum\limits_{i = 0} ^{k}C_{k + 1} ^{i}B_i x ^{k - i + 1} Sk(x)=k+11i=0∑kCk+1iBixk−i+1

一个自然数幂次求和的公式，其实就是我们推导过程中的一步。

#include <bits/stdc++.h>

using namespace std;

const int N = 3e3 +10, mod = 1e9 + 7;

int C[N][N], B[N], inv[N];

void init() {
  for (int i = 0; i < N; i++) {
    C[i][0] = C[i][i] = 1;
    for (int j = 1; j < i; j++) {
      C[i][j] = (C[i - 1][j] + C[i - 1][j - 1]) % mod;
    }
  }
  inv[1] = 1;
  for (int i = 2; i < N; i++) {
    inv[i] = 1ll * (mod - mod / i) * inv[mod % i] % mod;
  }
  B[0] = 1;
  for (int i = 1; i < N - 1; i++) {
    int cur = 0;
    for (int j = 0; j < i; j++) {
      cur = (cur + 1ll * C[i + 1][j] * B[j] % mod) % mod;
    }
    cur = 1ll * cur * inv[i + 1] % mod;
    B[i] = (mod - cur) % mod;
  }
}

int main() {
  // freopen("in.txt", "r", stdin);
  // freopen("out.txt", "w", stdout);
  // ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
  init();
  int T, k;
  scanf("%d", &T);
  while (T--) {
    long long n, ans = 0, cur = 1;
    scanf("%lld %d", &n, &k);
    n++;
    n %= mod;
    for (int i = k; i >= 0; i--) {
      cur = cur * n % mod;
      ans = (ans + 1ll * C[k + 1][i] * B[i] % mod * cur % mod) % mod;
    }
    printf("%lld\n", ans * inv[k + 1] % mod);
  }
  return 0;
}

考虑生成函数

∑ i = 0 n C n + 1 i B i = [ n = 0 ] 对两边加上 B n + 1 ∑ i = 0 n + 1 C n + 1 i B i = [ n = 0 ] + B n + 1 ∑ i = 0 n C n i B i = [ n = 1 ] + B n ∑ i = 0 n 1 i ! ( n − i ) ! B i = [ n = 1 ] + B n n ! 设指数型生成函数 B ( x ) = ∑ n ≥ 0 B n x n 有 B ( x ) e x = x + B ( x ) 从而有 B ( x ) = x e x − 1 B ( x ) = x ∑ n ≥ 0 1 n ! x n − 1 x ∑ n ≥ 1 1 n ! x n 1 ∑ n ≥ 0 1 ( n + 1 ) ! x n \sum_{i = 0} ^{n} C_{n + 1} ^{i} B_i = [n = 0]\\ 对两边加上B_{n + 1}\\ \sum_{i = 0} ^{n + 1} C_{n + 1} ^{i} B_i = [n = 0] + B_{n + 1}\\ \sum_{i = 0} ^{n} C_{n} ^{i} B_i = [n = 1] + B_{n}\\ \sum_{i = 0} ^{n} \frac{1}{i!(n - i)!} B_{i} = [n = 1] + \frac{B_{n}}{n!}\\ 设指数型生成函数B(x) = \sum_{n \geq 0} B_n x ^ n\\ 有B(x) e ^x = x + B(x)\\ 从而有B(x) = \frac{x}{e ^ x - 1}\\ B(x) = \frac{x}{\sum\limits_{n \geq 0} \frac{1}{n!} x ^ n - 1}\\ \frac{x}{\sum\limits_{n \geq 1} \frac{1}{n!} x ^ n}\\ \frac{1}{\sum\limits_{n \geq 0} \frac{1}{(n + 1)!}x ^ n}\\ i=0∑nCn+1iBi=[n=0]对两边加上Bn+1i=0∑n+1Cn+1iBi=[n=0]+Bn+1i=0∑nCniBi=[n=1]+Bni=0∑ni!(n−i)!1Bi=[n=1]+n!Bn设指数型生成函数B(x)=n≥0∑Bnxn有B(x)ex=x+B(x)从而有B(x)=ex−1xB(x)=n≥0∑n!1xn−1xn≥1∑n!1xnxn≥0∑(n+1)!1xn1
只需多项式求逆即可解决。

#include <bits/stdc++.h>

using namespace std;

const int N = 1e6 + 10, mod = 998244353;

int r[N], b[N], t[N];

int quick_pow(int a, int n) {
  int ans = 1;
  while (n) {
    if (n & 1) {
      ans = 1ll * a * ans % mod;
    }
    a = 1ll * a * a % mod;
    n >>= 1;
  }
  return ans;
}

void get_r(int lim) {
  for (int i = 0; i < lim; i++) {
    r[i] = (i & 1) * (lim >> 1) + (r[i >> 1] >> 1);
  }
}

void NTT(int *f, int lim, int rev) {
  for (int i = 0; i < lim; i++) {
    if (i < r[i]) {
      swap(f[i], f[r[i]]);
    }
  }
  for (int mid = 1; mid < lim; mid <<= 1) {
    int wn = quick_pow(3, (mod - 1) / (mid << 1));
    for (int len = mid << 1, cur = 0; cur < lim; cur += len) {
      int w = 1;
      for (int k = 0; k < mid; k++, w = 1ll * w * wn % mod) {
        int x = f[cur + k], y = 1ll * w * f[cur + mid + k] % mod;
        f[cur + k] = (x + y) % mod, f[cur + mid + k] = (x - y + mod) % mod;
      }
    }
  }
  if (rev == -1) {
    int inv = quick_pow(lim, mod - 2);
    reverse(f + 1, f + lim);
    for (int i = 0; i < lim; i++) {
      f[i] = 1ll * f[i] * inv % mod;
    }
  }
}

void polyinv(int *f, int *g, int n) {
  if (n == 1) {
    g[0] = quick_pow(f[0], mod - 2);
    return ;
  }
  polyinv(f, g, n + 1 >> 1);
  for (int i = 0; i < n; i++) {
    t[i] = f[i];
  }
  int lim = 1;
  while (lim < 2 * n) {
    lim <<= 1;
  }
  get_r(lim);
  NTT(t, lim, 1);
  NTT(g, lim, 1);
  for (int i = 0; i < lim; i++) {
    int cur = (2 - 1ll * g[i] * t[i] % mod + mod) % mod;
    g[i] = 1ll * g[i] * cur % mod;
    t[i] = 0;
  }
  NTT(g, lim, -1);
  for (int i = n; i < lim; i++) {
    g[i] = 0;
  }
}

int fac[N], ifac[N], B[N], n;

void init() {
  fac[0] = 1;
  for (int i = 1; i < N; i++) {
    fac[i] = 1ll * fac[i - 1] * i % mod;
  }
  ifac[N - 1] = quick_pow(fac[N - 1], mod - 2);
  for (int i = N - 2; i >= 0; i--) {
    ifac[i] = 1ll * ifac[i + 1] * (i + 1) % mod;
  }
}

int main() {
  // freopen("in.txt", "r", stdin);
  // freopen("out.txt", "w", stdout);
  // ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
  init();
  n = 100000;
  for (int i = 0; i <= n + 1; i++) {
    B[i] = ifac[i + 1];
  }
  polyinv(B, b, n + 2);
  for (int i = 0; i <= n + 1; i++) {
    B[i] = 1ll * b[i] * fac[i] % mod;
  }
  for (int i = 0; i <= n + 1; i++) {
    printf("%d%c", B[i], i == n + 1 ? '\n' : ' ');
  }
  return 0;
}

P3711 仓鼠的数学题

F ( x ) = ∑ k = 0 n S k ( x ) a k 原本定义的 S k ( x ) = ∑ i = 0 x i k 根据伯努利数的定义 S k ′ ( x ) = ∑ i = 0 x − 1 i k 则我们求 F ( x ) = ∑ k = 0 n S k ′ ( x ) a k , 答案即为 F ( x + 1 ) 考虑先求 F ( x ) = ∑ k = 0 n a k 1 k + 1 ∑ i = 0 k C k + 1 i B i x k − i + 1 ∑ k = 0 n a k k + 1 ∑ i = 0 k ( k + 1 ) ! i ! ( k + 1 − i ) ! B i x k − i + 1 ∑ k = 0 n a k k ! ∑ i = 0 k B i i ! x k − i + 1 ( k − i + 1 ) ! ∑ i = 1 n + 1 x i i ! ∑ k = i − 1 n ( a k k ! ) ( B k + 1 − i ( k + 1 − i ) ! ) 设 f ( n ) = a k k ! , g ( n ) = B n n ! , 另 g ( n ) = g ( n − i ) 即翻转一下 ∑ i = 1 n + 1 x i i ! ∑ k = i − 1 n f ( k ) g ( n + i − k ) 容易发现后面是一个卷积，所以可以优化 F(x) = \sum_{k = 0} ^{n} S_k(x) a_k\\ 原本定义的S_k(x) = \sum_{i = 0} ^{x} i ^ k\\ 根据伯努利数的定义S'_k(x) = \sum_{i = 0} ^{x - 1} i ^ k\\ 则我们求F(x) = \sum_{k = 0} ^{n} S'_k(x) a_k, 答案即为F(x + 1)\\ 考虑先求F(x) = \sum_{k = 0} ^{n} a_k \frac{1}{k + 1} \sum\limits_{i = 0} ^{k}C_{k + 1} ^{i}B_i x ^{k - i + 1}\\ \sum_{k = 0} ^{n} \frac{a_k}{k + 1} \sum_{i = 0} ^{k} \frac{(k + 1)!}{i!(k + 1 - i)!} B_i x ^{k - i + 1}\\ \sum_{k = 0} ^{n} a_k k! \sum_{i = 0} ^{k} \frac{B_i}{i!}\frac{x ^{k - i + 1}}{(k - i + 1)!}\\ \sum_{i = 1} ^{n + 1} \frac{x ^{i}}{i!} \sum_{k = i - 1} ^{n} \left(a_k k!\right) \left( \frac{B_{k + 1 - i}}{(k + 1 - i)!} \right)\\ 设f(n) = a_k k!, g(n) = \frac{B_{n}}{n!}, 另g(n) = g(n - i)即翻转一下\\ \sum_{i = 1} ^{n + 1} \frac{x ^ i}{i!} \sum_{k = i - 1} ^{n} f(k) g(n + i - k)\\ 容易发现后面是一个卷积，所以可以优化\\ F(x)=k=0∑nSk(x)ak原本定义的Sk(x)=i=0∑xik根据伯努利数的定义Sk′(x)=i=0∑x−1ik则我们求F(x)=k=0∑nSk′(x)ak,答案即为F(x+1)考虑先求F(x)=k=0∑nakk+11i=0∑kCk+1iBixk−i+1k=0∑nk+1aki=0∑ki!(k+1−i)!(k+1)!Bixk−i+1k=0∑nakk!i=0∑ki!Bi(k−i+1)!xk−i+1i=1∑n+1i!xik=i−1∑n(akk!)((k+1−i)!Bk+1−i)设f(n)=akk!,g(n)=n!Bn,另g(n)=g(n−i)即翻转一下i=1∑n+1i!xik=i−1∑nf(k)g(n+i−k)容易发现后面是一个卷积，所以可以优化
下面的 n = n + 1 n = n + 1 n=n+1，由上易知 F ( x ) F(x) F(x)是一个 n + 1 n + 1 n+1次多项式。
设 F ( x ) = ∑ i = 0 n a i x i F ( x + 1 ) = ∑ i = 0 n a i ( x + 1 ) i ∑ i = 0 n a i ∑ j = 0 i C j i x j ∑ i = 0 n a i ∑ j = 0 i i ! j ! ( i − j ) ! x j ∑ i = 0 n x i i ! ∑ j = i n a j j ! 1 ( j − i ) ! 另 f ( n ) = a n n ! , g ( n ) = 1 n ! , 再翻转 g ( n ) ∑ i = 0 n x i i ! ∑ j = i n f ( j ) g ( n + i − j ) 后面也同样是一个卷积，所以可以优化设F(x) = \sum_{i = 0} ^{n} a_i x ^ i\\ F(x + 1) = \sum_{i = 0} ^{n} a_i (x + 1) ^ i\\ \sum_{i = 0} ^{n} a_i \sum_{j = 0} ^{i} C_{j} ^{i} x ^ j\\ \sum_{i = 0} ^{n} a_i \sum_{j = 0} ^{i} \frac{i!}{j!(i - j)!} x ^ j\\ \sum_{i = 0} ^{n} \frac{x ^ i}{i!} \sum_{j = i} ^{n} a_j j! \frac{1}{(j - i)!}\\ 另f(n) = a_n n!, g(n) = \frac{1}{n!},再翻转g(n)\\ \sum_{i = 0} ^{n} \frac{x ^{i}}{i!} \sum_{j = i} ^{n} f(j) g(n + i - j)\\ 后面也同样是一个卷积，所以可以优化\\ 设F(x)=i=0∑naixiF(x+1)=i=0∑nai(x+1)ii=0∑naij=0∑iCjixji=0∑naij=0∑ij!(i−j)!i!xji=0∑ni!xij=i∑najj!(j−i)!1另f(n)=ann!,g(n)=n!1,再翻转g(n)i=0∑ni!xij=i∑nf(j)g(n+i−j)后面也同样是一个卷积，所以可以优化

#include <bits/stdc++.h>

using namespace std;

const int N = 1e6 + 10, mod = 998244353;

int r[N], t[N];

int quick_pow(int a, int n) {
  int ans = 1;
  while (n) {
    if (n & 1) {
      ans = 1ll * a * ans % mod;
    }
    a = 1ll * a * a % mod;
    n >>= 1;
  }
  return ans;
}

void get_r(int lim) {
  for (int i = 0; i < lim; i++) {
    r[i] = (i & 1) * (lim >> 1) + (r[i >> 1] >> 1);
  }
}

void NTT(int *f, int lim, int rev) {
  for (int i = 0; i < lim; i++) {
    if (i < r[i]) {
      swap(f[i], f[r[i]]);
    }
  }
  for (int mid = 1; mid < lim; mid <<= 1) {
    int wn = quick_pow(3, (mod - 1) / (mid << 1));
    for (int len = mid << 1, cur = 0; cur < lim; cur += len) {
      int w = 1;
      for (int k = 0; k < mid; k++, w = 1ll * w * wn % mod) {
        int x = f[cur + k], y = 1ll * w * f[cur + mid + k] % mod;
        f[cur + k] = (x + y) % mod, f[cur + mid + k] = (x - y + mod) % mod;
      }
    }
  }
  if (rev == -1) {
    int inv = quick_pow(lim, mod - 2);
    reverse(f + 1, f + lim);
    for (int i = 0; i < lim; i++) {
      f[i] = 1ll * f[i] * inv % mod;
    }
  }
}

void polyinv(int *f, int *g, int n) {
  if (n == 1) {
    g[0] = quick_pow(f[0], mod - 2);
    return ;
  }
  polyinv(f, g, n + 1 >> 1);
  for (int i = 0; i < n; i++) {
    t[i] = f[i];
  }
  int lim = 1;
  while (lim < 2 * n) {
    lim <<= 1;
  }
  get_r(lim);
  NTT(t, lim, 1);
  NTT(g, lim, 1);
  for (int i = 0; i < lim; i++) {
    int cur = (2 - 1ll * g[i] * t[i] % mod + mod) % mod;
    g[i] = 1ll * g[i] * cur % mod;
    t[i] = 0;
  }
  NTT(g, lim, -1);
  for (int i = n; i < lim; i++) {
    g[i] = 0;
  }
}

int fac[N], ifac[N], a[N], f[N], g[N], B[N], n;

void init() {
  fac[0] = 1;
  for (int i = 1; i < N; i++) {
    fac[i] = 1ll * fac[i - 1] * i % mod;
  }
  ifac[N - 1] = quick_pow(fac[N - 1], mod - 2);
  for (int i = N - 2; i >= 0; i--) {
    ifac[i] = 1ll * ifac[i + 1] * (i + 1) % mod;
  }
}

int main() {
  // freopen("in.txt", "r", stdin);
  // freopen("out.txt", "w", stdout);
  // ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
  init();
  scanf("%d", &n);
  for (int i = 0; i <= n; i++) {
    scanf("%d", &a[i]);
  }
  for (int i = 0; i <= n + 1; i++) {
    B[i] = ifac[i + 1];
  }
  polyinv(B, g, n + 2);
  for (int i = 0; i <= n; i++) {
    f[i] = 1ll * a[i] * fac[i] % mod;
  }
  reverse(g, g + n + 2);
  int lim = 1;
  while (lim < 2 * n + 10) {
    lim <<= 1;
  }
  get_r(lim);
  NTT(f, lim, 1), NTT(g, lim, 1);
  for (int i = 0; i < lim; i++) {
    f[i] = 1ll * f[i] * g[i] % mod;
  }
  NTT(f, lim, -1);
  for (int i = 1; i <= n + 1; i++) {
    a[i] = 1ll * ifac[i] * f[n + i] % mod;
  }
  a[0] = 0;
  n++;
  for (int i = 0; i < lim; i++) {
    f[i] = g[i] = 0;
  }
  for (int i = 0; i <= n; i++) {
    f[i] = 1ll * a[i] * fac[i] % mod;
    g[i] = ifac[n - i];
  }
  NTT(f, lim, 1), NTT(g, lim, 1);
  for (int i = 0; i < lim; i++) {
    f[i] = 1ll * f[i] * g[i] % mod;
  }
  NTT(f, lim, -1);
  for (int i = 0; i <= n; i++) {
    printf("%lld%c", 1ll * ifac[i] * f[n + i] % mod, i == n ? '\n' : ' ');
  }
  return 0;
}

码农公寓

伯努利数

从二项式出发

考虑生成函数

P3711 仓鼠的数学题

相关文章