斯特林数的四种求法

几乎是抄神仙yyc的blog
P5395 第二类斯特林数·行
m n = ∑ i = 0 m { n i } ( m i ) i ! m^n=\sum_{i=0}^m\begin{Bmatrix}n \\ i\end{Bmatrix} \binom{m}{i}i! mn=i=0m{ni}(im)i!
二项式反演得到
{ n m } = 1 m ! ∑ i = 0 m ( − 1 ) m − i ( m i ) i n \begin{Bmatrix}n \\ m\end{Bmatrix}=\frac{1}{m!}\sum_{i=0}^m(-1)^{m-i}\binom{m}{i}i^n {nm}=m!1i=0m(1)mi(im)in
把这个式子化简一下可以得到
∑ i = 0 m ( − 1 ) m − i ( m − i ) ! × i n i ! \sum_{i=0}^m\frac{(-1)^{m-i}}{(m-i)!}\times\frac{i^n}{i!} i=0m(mi)!(1)mi×i!in
大力卷积即可
code:

#include<bits/stdc++.h>
#define mod 167772161
#define N 4000500
using namespace std;
int add(int x, int y) { x += y;
    if(x >= mod) x -= mod;
    return x;
}
int sub(int x, int y) { x -= y;
    if(x < 0) x += mod;
    return x;
}
int mul(int x, int y) {
    return 1ll * x * y % mod;
}
int qpow(int x, int y) {
    int ret = 1;
    for(; y; y >>= 1, x = mul(x, x)) if(y & 1) ret = mul(ret, x);
    return ret;
}
const int G = 3;
const int Ginv = qpow(G, mod - 2);
int rev[N];
void ntt(int *a, int n, int o) {
    for(int i = 1; i < n; i ++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) * (n >> 1));
    for(int i = 1; i < n; i ++) if(i > rev[i]) swap(a[i], a[rev[i]]);
    for(int len = 2; len <= n; len <<= 1) {
        int w0 = qpow(o == 1? G : Ginv, (mod - 1) / len);
        for(int j = 0; j < n; j += len) {
            int wn = 1;
            for(int k = j; k < j + (len >> 1); k ++, wn = mul(wn, w0)) {
                int X = a[k], Y = mul(a[k + (len >> 1)], wn);
                a[k] = add(X, Y); a[k + (len >> 1)] = sub(X, Y);
            }
        }
    }
    int ninv = qpow(n, mod - 2);
    if(o == -1)
        for(int i = 0; i < n; i ++) a[i] = mul(a[i], ninv);
}
int n, fac[N], ifac[N], a[N], b[N];
int main() {
    scanf("%d", &n);
    fac[0] = 1;
    for(int i = 1; i <= n; i ++) fac[i] = mul(fac[i - 1], i);
    ifac[n] = qpow(fac[n], mod - 2);
    for(int i = n - 1; i >= 0; i --) ifac[i] = mul(ifac[i + 1], i + 1);
    for(int i = 0; i <= n; i ++) a[i] = mul(qpow(mod - 1, i), ifac[i]);
    for(int i = 0; i <= n; i ++) b[i] = mul(qpow(i, n), ifac[i]);
    
    int len = 1;
    for(; len <= n + n; ) len <<= 1;
    ntt(a, len, 1), ntt(b, len, 1);
    for(int i = 0; i < len; i ++) a[i] = mul(a[i], b[i]);
    ntt(a, len, - 1);
    for(int i = 0; i <= n; i ++) printf("%d ", a[i]);
    return 0;
}

P5408 第一类斯特林数·行
发现 x n ‾ x^{\overline{n}} xn的生成函数各项系数就是要求的,所以我们考虑计算这个即可
发现可以倍增,再加上多项式平移的技巧即可
即由 f ( x ) f(x) f(x)得到 f ( x + n ) f(x+n) f(x+n)
考虑用 x + n x+n x+n替换原来的 x x x,再用二项式定理展开

f ( x + n ) = ∑ i a [ i ] ( x + n ) i = ∑ i x i ∑ j = i a [ i ] ( j i ) n j − i f(x+n)=\sum_ia[i](x+n)^i\\=\sum_ix^i\sum_{j=i}a[i]\binom{j}{i}n^{j-i} f(x+n)=ia[i](x+n)i=ixij=ia[i](ij)nji
右边那块一看就很能卷,按照套路去推即可
code:

#include<bits/stdc++.h>
#define mod 167772161
#define N 4000500
using namespace std;
int add(int x, int y) { x += y;
    if(x >= mod) x -= mod;
    return x;
}
int sub(int x, int y) { x -= y;
    if(x < 0) x += mod;
    return x;
}
int mul(int x, int y) {
    return 1ll * x * y % mod;
}
int qpow(int x, int y) {
    int ret = 1;
    for(; y; y >>= 1, x = mul(x, x)) if(y & 1) ret = mul(ret, x);
    return ret;
}
const int G = 3;
const int Ginv = qpow(G, mod - 2);
int rev[N];
void ntt(int *a, int n, int o) {
    for(int i = 1; i < n; i ++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) * (n >> 1));
    for(int i = 1; i < n; i ++) if(i > rev[i]) swap(a[i], a[rev[i]]);
    for(int len = 2; len <= n; len <<= 1) {
        int w0 = qpow(o == 1? G : Ginv, (mod - 1) / len);
        for(int j = 0; j < n; j += len) {
            int wn = 1;
            for(int k = j; k < j + (len >> 1); k ++, wn = mul(wn, w0)) {
                int X = a[k], Y = mul(a[k + (len >> 1)], wn);
                a[k] = add(X, Y); a[k + (len >> 1)] = sub(X, Y);
            }
        }
    }
    int ninv = qpow(n, mod - 2);
    if(o == -1)
        for(int i = 0; i < n; i ++) a[i] = mul(a[i], ninv);
}
int n, fac[N], ifac[N], a[N], b[N], f[N], g[N];
void init(int n) {
    fac[0] = 1;
    for(int i = 1; i <= n; i ++) fac[i] = mul(fac[i - 1], i);
    ifac[n] = qpow(fac[n], mod - 2);
    for(int i = n - 1; i >= 0; i --) ifac[i] = mul(ifac[i + 1], i + 1);
}
void Mul(int *a, int *b, int n, int m) {
    int len = 1;
    for(; len <= n + m; ) len <<= 1;
    ntt(a, len, 1), ntt(b, len, 1);
    for(int i = 0; i < len; i ++) a[i] = mul(a[i], b[i]);
    ntt(a, len, - 1);
    for(int i = n + m + 1; i <= len; i ++) a[i] = 0;
}
void solve(int *f, int n) {
    if(n == 1) {f[1] = 1; return ;}
    if(n & 1) {
        solve(f, n - 1);
        for(int i = n; i >= 1; i --) {
            f[i] = add(f[i - 1], mul(f[i], n - 1));
        }
        f[0] = mul(f[0], n - 1);
    } else {
        solve(f, n >> 1);
        int m = n / 2, mi = 1;
        for(int i = 0; i <= m; i ++, mi = mul(mi, m))
            a[i] = mul(f[i], fac[i]), b[i] = mul(mi, ifac[i]);
        reverse(a, a + 1 + m);
        Mul(a, b, m, m);
        for(int i = 0; i <= m; i ++) g[i] = mul(a[m - i], ifac[i]);
        Mul(f, g, m, m);
        int len = 1;
        for( ;len <= m + m; ) len <<= 1;
        for(int i = 0; i <= len; i ++) g[i] = a[i] = b[i] = 0;
     //   for(int i = 0; i <= n; i ++) printf("%d ", f[i]); printf("    %d\n", n);
    }
}
int main() {
    scanf("%d", &n);
    init(n);
    solve(f, n);
    for(int i = 0; i <= n; i ++) printf("%d ", f[i]);
    return 0;
}

P5396 第二类斯特林数·列

考虑组合意义:一个盒子的EGF是 ( e x − 1 ) (e^x-1) (ex1) k个盒子的无序组合就是
1 k ! ( e x − 1 ) k \frac{1}{k!}(e^x-1)^k k!1(ex1)k
这个就是第二类斯特林数的生成函数
直接多项式快速幂走起
注意最后要 ∗ i ! *i! i!
code:


#include<bits/stdc++.h>
#define int long long
#define mod 167772161
#define G 3
#define N 800005
using namespace std;
int qpow(int x, int y){
    int ret = 1;
    for(; y; y >>= 1, x = x * x % mod) if(y & 1) ret = ret * x % mod;
    return ret;
}
int rev[N], G_inv, len_inv;
void ntt(int *a, int len, int o){
    len_inv = qpow(len, mod - 2);
    for(int i = 0; i <= len; i ++) rev[i] = (rev[i >> 1] >> 1) | ((i&1) * len >> 1);
    for(int i = 0; i <= len; i ++) if(i < rev[i]) swap(a[i], a[rev[i]]);
    for(int i = 2; i <= len; i <<= 1){
        int wn = qpow((o == 1)? G:G_inv, (mod - 1) / i);
        for(int j = 0, p = i / 2; j + i - 1 <= len; j += i){
            int w0 = 1;
            for(int k = j; k < j + p; k ++, w0 = w0 * wn % mod){
                int X = a[k], Y = w0 * a[k + p] % mod;
                a[k] = (X + Y) % mod;
                a[k + p] = (X - Y + mod) % mod;
            }
        }
    }
    if(o == -1)
        for(int i = 0; i <= len; i ++) a[i] = a[i] * len_inv % mod;
}
int c[N];
void inv(int *a, int *b, int sz){
    if(sz == 0) {b[0] = qpow(a[0], mod - 2); return;}
    inv(a, b, sz / 2);
    int len = 1;
    for(; len <= sz + sz; len <<= 1);
    for(int i = 0; i <= sz; i ++) c[i] = a[i];
    for(int i = sz + 1; i <= len; i ++) c[i] = 0;
    ntt(c, len, 1), ntt(b, len, 1);
    for(int i = 0; i <= len; i ++) b[i] = (b[i] * 2 % mod - b[i] * b[i] % mod * c[i] % mod + mod) % mod;
    ntt(b, len, -1);
    for(int i = sz + 1; i <= len; i ++) b[i] = 0;
}
void qiudao(int *a, int sz) {
    for(int i = 0; i < sz; i ++) a[i] = a[i + 1] * (i + 1) % mod;
    a[sz] = 0;
}
void jifen(int *a, int sz) {
    for(int i = sz; i >= 1; i --) a[i] = a[i - 1] * qpow(i, mod - 2) % mod;
    a[0] = 0;
}
int Ad[N], An[N];
void ln(int *A, int n) {
    for(int i = 0; i <= n; i ++) Ad[i] = A[i];
    qiudao(Ad, n);
    inv(A, An, n);
    int len = 1;
    for(; len <= n + n;) len <<= 1;
    ntt(Ad, len, 1), ntt(An, len, 1);
    for(int i = 0; i <= len; i ++) Ad[i] = Ad[i] * An[i] % mod;
    ntt(Ad, len, -1);
    jifen(Ad, n);
    for(int i = 0; i <= n; i ++) A[i] = Ad[i];
    for(int i = 0; i <= len; i ++) An[i] = Ad[i] = 0;
}
int fln[N];
void exp(int *a, int *b, int n) {
    if(n == 0) {b[0] = 1; return;}
    exp(a, b, n / 2);
    for(int i = 0; i <= n; i ++) fln[i] = b[i]; ln(fln, n);
    fln[0] = 1;
    for(int i = 1; i <= n; i ++) fln[i] = (a[i] - fln[i] + mod ) % mod;
    int len = 1;
    for(; len <= n + n;) len <<= 1;
    ntt(b, len, 1), ntt(fln, len, 1);
    for(int i = 0; i <= len; i ++) b[i] = b[i] * fln[i] % mod;
    ntt(b, len, -1);
    for(int i = 0; i <= len; i ++) fln[i] = 0;
}
int a[N], b[N], n, m, fac[N], ifac[N];
void init(int n) {
    fac[0] = 1;
    for(int i = 1; i <= n; i ++) fac[i] = fac[i - 1] * i % mod;
    ifac[n] = qpow(fac[n], mod - 2);
    for(int i = n - 1; i >= 0; i --) ifac[i] = 1ll * ifac[i + 1] * (i + 1) % mod;
}
signed main(){
    init(N - 10);
    G_inv = qpow(G, mod - 2);
    scanf("%lld%lld", &n, &m);
    for(int i = 0; i < n; i ++) a[i] = ifac[i + 1];
 //   for(int i = 0; i < n; i ++) printf("%lld ", a[i]); printf("\n");
    ln(a, n);
    for(int i = 0; i <= n; i ++) a[i] = a[i] * m % mod;
    exp(a, b, n);
    for(int i = 1; i <= min(n + 1, m); i ++) printf("0 ");
    for(int i = 0; i <= n - m; i ++) printf("%lld ", b[i] * fac[i + m] % mod * ifac[m] % mod);
    return 0;
}

P5409 第一类斯特林数·列

同上,考虑组合意义:一个盒子的生成函数是 ∑ i ( i − 1 ) ! x i i ! = ∑ i 1 i x i \large \sum_i(i-1)!\frac{x^i}{i!}=\sum_i \frac{1}{i}x^i i(i1)!i!xi=ii1xi
通过luogu P4389 付公主的背包
我们知道这个就是 ln ⁡ ( 1 1 − x ) = ∑ 1 i x i \ln(\frac{1}{1-x})=\sum\frac{1}{i}x^i ln(1x1)=i1xi
于是第一类斯特林数的指数型生成函数就是
1 k ! ln ⁡ ( 1 1 − x ) k \frac{1}{k!}\ln(\frac{1}{1-x})^k k!1ln(1x1)k
同样的冲一波多项式快速幂即可
code:


#include<bits/stdc++.h>
#define int long long
#define mod 167772161
#define G 3
#define N 800005
using namespace std;
int qpow(int x, int y){
    int ret = 1;
    for(; y; y >>= 1, x = x * x % mod) if(y & 1) ret = ret * x % mod;
    return ret;
}
int rev[N], G_inv, len_inv;
void ntt(int *a, int len, int o){
    len_inv = qpow(len, mod - 2);
    for(int i = 0; i <= len; i ++) rev[i] = (rev[i >> 1] >> 1) | ((i&1) * len >> 1);
    for(int i = 0; i <= len; i ++) if(i < rev[i]) swap(a[i], a[rev[i]]);
    for(int i = 2; i <= len; i <<= 1){
        int wn = qpow((o == 1)? G:G_inv, (mod - 1) / i);
        for(int j = 0, p = i / 2; j + i - 1 <= len; j += i){
            int w0 = 1;
            for(int k = j; k < j + p; k ++, w0 = w0 * wn % mod){
                int X = a[k], Y = w0 * a[k + p] % mod;
                a[k] = (X + Y) % mod;
                a[k + p] = (X - Y + mod) % mod;
            }
        }
    }
    if(o == -1)
        for(int i = 0; i <= len; i ++) a[i] = a[i] * len_inv % mod;
}
int c[N];
void inv(int *a, int *b, int sz){
    if(sz == 0) {b[0] = qpow(a[0], mod - 2); return;}
    inv(a, b, sz / 2);
    int len = 1;
    for(; len <= sz + sz; len <<= 1);
    for(int i = 0; i <= sz; i ++) c[i] = a[i];
    for(int i = sz + 1; i <= len; i ++) c[i] = 0;
    ntt(c, len, 1), ntt(b, len, 1);
    for(int i = 0; i <= len; i ++) b[i] = (b[i] * 2 % mod - b[i] * b[i] % mod * c[i] % mod + mod) % mod;
    ntt(b, len, -1);
    for(int i = sz + 1; i <= len; i ++) b[i] = 0;
}
void qiudao(int *a, int sz) {
    for(int i = 0; i < sz; i ++) a[i] = a[i + 1] * (i + 1) % mod;
    a[sz] = 0;
}
void jifen(int *a, int sz) {
    for(int i = sz; i >= 1; i --) a[i] = a[i - 1] * qpow(i, mod - 2) % mod;
    a[0] = 0;
}
int Ad[N], An[N];
void ln(int *A, int n) {
    for(int i = 0; i <= n; i ++) Ad[i] = A[i];
    qiudao(Ad, n);
    inv(A, An, n);
    int len = 1;
    for(; len <= n + n;) len <<= 1;
    ntt(Ad, len, 1), ntt(An, len, 1);
    for(int i = 0; i <= len; i ++) Ad[i] = Ad[i] * An[i] % mod;
    ntt(Ad, len, -1);
    jifen(Ad, n);
    for(int i = 0; i <= n; i ++) A[i] = Ad[i];
    for(int i = 0; i <= len; i ++) An[i] = Ad[i] = 0;
}
int fln[N];
void exp(int *a, int *b, int n) {
    if(n == 0) {b[0] = 1; return;}
    exp(a, b, n / 2);
    for(int i = 0; i <= n; i ++) fln[i] = b[i]; ln(fln, n);
    fln[0] = 1;
    for(int i = 1; i <= n; i ++) fln[i] = (a[i] - fln[i] + mod ) % mod;
    int len = 1;
    for(; len <= n + n;) len <<= 1;
    ntt(b, len, 1), ntt(fln, len, 1);
    for(int i = 0; i <= len; i ++) b[i] = b[i] * fln[i] % mod;
    ntt(b, len, -1);
    for(int i = 0; i <= len; i ++) fln[i] = 0;
}
int a[N], b[N], n, m, fac[N], ifac[N];
void init(int n) {
    fac[0] = 1;
    for(int i = 1; i <= n; i ++) fac[i] = fac[i - 1] * i % mod;
    ifac[n] = qpow(fac[n], mod - 2);
    for(int i = n - 1; i >= 0; i --) ifac[i] = 1ll * ifac[i + 1] * (i + 1) % mod;
}
signed main(){
    init(N - 10);
    G_inv = qpow(G, mod - 2);
    scanf("%lld%lld", &n, &m);
    for(int i = 0; i < n; i ++) a[i] = qpow(i + 1, mod - 2);
    ln(a, n);
    for(int i = 0; i <= n; i ++) a[i] = a[i] * m % mod;
    exp(a, b, n);
    for(int i = 1; i <= min(n + 1, m); i ++) printf("0 ");
    for(int i = 0; i <= n - m; i ++) printf("%lld ", b[i] * fac[i + m] % mod * ifac[m] % mod);
    return 0;
}
posted @ 2021-08-23 01:31  lahlah  阅读(45)  评论(0编辑  收藏  举报