伯努利数定义：
$$\dfrac{t}{e^t - 1} = \sum_{n = 0}^{\infty} \dfrac{B_n}{n!} t^n$$
递推式：
$$\begin{align}
& \sum_{k = 0}^{n}C_{n + 1}^kB_k = 0 \\
\Longrightarrow & B_n = -\dfrac1{n+1}(C_{n+1}^0B_0 + C_{n+1}^1B_1+…+C_{n+1}^{n-1}B_{n-1})
\end{align}$$
ACdreamer有两篇博客写的很不错
应用：
求自然数的幂等和
$$
\sum_{i=1}^ni^k = \dfrac1{k+1}\sum_{i=1}^{k+1}C_{k+1}^iB_{k+1-i}(n+1)^i
$$
在k不是很大的时候，上式中除了$B_i$之外剩下的都可以轻松预处理
对于$B_i$,可以用之前的递推式子$O(k^2)$的预处理出来
然后就可以线性时间内求得自然数的幂等和
当k比较大的时候，我们重新审视一下定义：
$$\begin{align}
\sum_{n = 0}^{\infty} \dfrac{B_n}{n!} t^n =& \dfrac{t}{e^t - 1} (用泰勒公式展开e^t) \\
=& \dfrac{t}{\sum_{i=1}^{\infty}\dfrac{t^i}{i!}} \\
=& \dfrac{1}{\sum_{i=0}^{\infty}\dfrac{t^i}{(i+1)!}}
\end{align}$$
观察最后的形式，其实就是一个多项式的逆元。
那么我们把这个母函数放在模$t^{n+1}$意义下求出来，每一项就都出来了。
代码
/*
* Filename: 51nod-1258.cpp
* Created: Tuesday, October 17, 2017 03:32:24 PM
* Author: crazyX
* More:
*
*/
#include <bits/stdc++.h>
#define mp make_pair
#define pb push_back
#define fi first
#define se second
#define SZ(x) ((int) (x).size())
#define all(x) (x).begin(), (x).end()
#define sqr(x) ((x) * (x))
#define clr(a,b) (memset(a,b,sizeof(a)))
#define y0 y3487465
#define y1 y8687969
#define fastio std::ios::sync_with_stdio(false)
using namespace std;
typedef long long ll;
typedef double DB;
const int maxn = 50000 + 7;
const int maxLen = 17, maxm = 1 << maxLen | 1;
const ll maxv = 1e18 + 6; // 1e14, 1e15
const DB pi = acos(-1.0); // double is enough
ll mod = 1e9 + 7, nlim, sp, msk;
#define _ %mod
#define __ %=mod
inline ll read(){
char c=getchar(); ll x=0,f=1;
while(c<'0'||c>'9') {if(c=='-')f=-1;c=getchar();}
while(c>='0'&&c<='9') {x=x*10+c-'0';c=getchar();}
return x*f;
}
int qpow(ll x, ll p) {
ll ret = 1;
while (p) {
if (p & 1) (ret *= x) __;
(x *= x) __;
p >>= 1;
}
return ret _;
}
namespace FFT{
struct cp {
DB r, i;
cp() {}
cp(DB r, DB i) : r(r), i(i) {}
cp operator + (cp const &t) const { return cp(r + t.r, i + t.i); }
cp operator - (cp const &t) const { return cp(r - t.r, i - t.i); }
cp operator * (cp const &t) const { return cp(r * t.r - i * t.i, r * t.i + i * t.r); }
cp conj() const { return cp(r, -i); }
} w[maxm], wInv[maxm];
void init() {
for(int i = 0, ilim = 1 << maxLen; i < ilim; ++i) {
int j = i, k = ilim >> 1; // 2 pi / ilim
for( ; !(j & 1) && !(k & 1); j >>= 1, k >>= 1);
w[i] = cp(cos(pi / k * j), sin(pi / k * j));
wInv[i] = w[i].conj();
}
nlim = std::min(maxv / (mod - 1) / (mod - 1), maxn - 1LL);
for(sp = 1; 1 << (sp << 1) < mod; ++sp);
msk = (1 << sp) - 1;
}
void FFT(int n, cp a[], int flag) {
static int bitLen = 0, bitRev[maxm] = {};
if(n != (1 << bitLen)) {
for(bitLen = 0; 1 << bitLen < n; ++bitLen);
for(int i = 1; i < n; ++i)
bitRev[i] = (bitRev[i >> 1] >> 1) | ((i & 1) << (bitLen - 1));
}
for(int i = 0; i < n; ++i)
if(i < bitRev[i])
std::swap(a[i], a[bitRev[i]]);
for(int i = 1, d = 1; d < n; ++i, d <<= 1)
for(int j = 0; j < n; j += d << 1)
for(int k = 0; k < d; ++k) {
cp &AL = a[j + k], &AH = a[j + k + d];
cp TP = w[k << (maxLen - i)] * AH;
AH = AL - TP, AL = AL + TP;
}
if(flag != -1)
return;
std::reverse(a + 1, a + n);
for(int i = 0; i < n; ++i) {
a[i].r /= n;
a[i].i /= n;
}
}
void polyMul(int a[], int aLen, int b[], int bLen, int c[]) { // c not in {a, b}
static cp A[maxm], B[maxm], C[maxm], D[maxm];
int len, cLen = aLen + bLen - 1; // optional: parameter
for(len = 1; len < aLen + bLen - 1; len <<= 1);
if(std::min(aLen, bLen) <= nlim) {
for(int i = 0; i < len; ++i)
A[i] = cp(i < aLen ? a[i] : 0, i < bLen ? b[i] : 0);
FFT(len, A, 1);
cp tr(0, -0.25);
for(int i = 0, j; i < len; ++i)
j = (len - i) & (len - 1), B[i] = (A[i] * A[i] - (A[j] * A[j]).conj()) * tr;
FFT(len, B, -1);
for(int i = 0; i < cLen; ++i) c[i] = (ll)(B[i].r + 0.5) % mod;
return;
} // if min(aLen, bLen) * mod <= maxv
for(int i = 0; i < len; ++i) {
A[i] = i < aLen ? cp(a[i] & msk, a[i] >> sp) : cp(0, 0);
B[i] = i < bLen ? cp(b[i] & msk, b[i] >> sp) : cp(0, 0);
}
FFT(len, A, 1), FFT(len, B, 1);
cp trL(0.5, 0), trH(0, -0.5), tr(0, 1);
for(int i = 0, j; i < len; ++i) {
j = (len - i) & (len - 1);
cp AL = (A[i] + A[j].conj()) * trL;
cp AH = (A[i] - A[j].conj()) * trH;
cp BL = (B[i] + B[j].conj()) * trL;
cp BH = (B[i] - B[j].conj()) * trH;
C[i] = AL * (BL + BH * tr);
D[i] = AH * (BL + BH * tr);
}
FFT(len, C, -1), FFT(len, D, -1);
for(int i = 0; i < cLen; ++i) {
int v11 = (ll)(C[i].r + 0.5) % mod, v12 = (ll)(C[i].i + 0.5) % mod;
int v21 = (ll)(D[i].r + 0.5) % mod, v22 = (ll)(D[i].i + 0.5) % mod;
c[i] = (((((ll)v22 << sp) + v12 + v21) << sp) + v11) % mod;
}
}
int c[maxm], tmp[maxm];
// y should clear to 0
void polyInv(int x[], int y[], int deg) {
if (deg == 1) {
y[0] = qpow(x[0], mod - 2);
return;
}
polyInv(x, y, (deg + 1) >> 1);
copy(x, x + deg, tmp);
int p = ((deg + 1) >> 1) + deg - 1;
polyMul(y, (deg + 1) >> 1, tmp, deg, c);
for (int i = 0; i < p; i += 1) c[i] = (- c[i] + mod) _;
(c[0] += 2) __;
polyMul(y, (deg + 1) >> 1, c, deg, tmp);
copy(tmp, tmp + deg, y);
}
};
int A[maxn], B[maxn];
ll inv[maxn], fac[maxn], facInv[maxn], mi[maxn];
inline void init() {
FFT::init();
A[0] = fac[0] = fac[1] = inv[0] = inv[1] = facInv[0] = facInv[1] = 1;
for (int i = 2; i < maxn; i += 1) {
fac[i] = fac[i - 1] * i _;
facInv[i] = qpow(fac[i], mod - 2);
A[i - 1] = qpow(fac[i], mod - 2);
inv[i] = qpow(i, mod - 2);
}
FFT::polyInv(A, B, maxn - 1);
for (int i = 0; i < maxn - 1; i += 1) B[i] = fac[i] * B[i] _;
}
inline ll C(int n, int m) {return fac[n] * facInv[m] _ * facInv[n - m] _;}
inline int cal(ll n, int k) {
(++n) __;
ll ret = 0;
mi[0] = 1;
for (int i = 1; i <= k + 1; i += 1) mi[i] = mi[i - 1] * n _;
for (int i = 1; i <= k + 1; i += 1) {
if (B[k + 1 - i]) (ret += C(k + 1, i) * B[k + 1 - i] _ * mi[i] _) __;
}
return ret * inv[k + 1] _;
}
int main()
{
#ifdef AC
freopen("data.in", "r", stdin);
freopen("data.out", "w", stdout);
#endif
int T, k; ll n;
init();
T = read();
while (T--) {
n = read(); k = read();
printf("%d\n", cal(n, k));
}
return 0;
}
![hdu 5730 Shell Necklace - [分治FFT|多项式求逆]](/medias/featureimages/6.jpg)