【算法】生成函数

推荐博客组合数学之三 —— 生成函数
生成函数是一种解决计数问题的好办法,一般多项式幂次表示方案的样子,系数表示有多少种方式。
比如这个1 1 1 1 1……的数列,就可以表示成f(x) = \sum_{i=0}^{+\infty}x^i
然后首先第一个应用最简单的应用,就是根据这个性质加速背包问题,我们把单个物品的体积看作幂次,然后如果有这种体积的物品就在这次上+1,用fft/ntt进行卷积,可以把一个n变成logn
比如这个题目luoguP3779[SDOI2017]龙与地下城
这个题目就是一个典型,我们要分情况考虑,对于较小的数据,我们可以直接列出一个初始多项式f(x)=\sum_{i=0}^{n-1}\frac{1}{n}x^i
然后我们对于这个多项式进行快速幂,求出它y次方。
那么对大数据的范围呢?
这个东西就要考虑正态分布等一堆东西就是了。
题解 P3779 【[SDOI2017]龙与地下城】
题目里面给你了这个函数的方差和期望公式,而你会发现60分做法和这两个东西毫无关系。他们真正的目的是想让你想到正态分布上面去。
事实上你感性想一下就会知道,这个概率会在中间那段特别大,会画出一个波峰,那个东西就叫做正态分布。
正态分布的密度函数是f(x)=\frac{1}{\sqrt{2\pi\sigma^2}}e^{\frac{-(x-\mu)^2}{2\sigma^2}}
看起来可以用c++自带的exp计算了,但是遗憾的是精度达不到。
但是根据中心极限定理
Y_n=\frac{\sum_{i=1}^{n}x_i-n\mu}{\sqrt{n\sigma^2}}
这个东西当n足够大是服从正态分布N(0,1)的。
也就说f(x)=\frac{1}{\sqrt{2\pi}}e^{\frac{-x^2}{2}}
然后我们就得到了一个可以大力辛普森积分的东西。
辛普森积分这个东西思想挺简单的,我们考虑大力用抛物线拟合它。
抛物线的面积公式是
\frac{(r-l)[f(l)+4f(mid)+f(r)]}{6}
然后你每次计算这个区间拆成两半的答案是不是和不拆基本一样,一样就返回,否则递归。

// luogu-judger-enable-o2
#include <bits/stdc++.h>
#include <algorithm>
#include <iostream>
#include <cstring>
#include <cstdlib>
#include <cstdio>
#include <vector>
#include <bitset>
#include <queue>
#include <cmath>
#include <map>
#include <set>
#define ri register int

using namespace std;

inline char gch()
{
    static char buf[100010], *h = buf, *t = buf;
    return h == t && (t = (h = buf) + fread(buf, 1, 100000, stdin), h == t) ? EOF : *h ++;
}

typedef unsigned int uint;

typedef long long lo;

template <typename int_qwq>

inline void re(int_qwq &x)
{
    x = 0;
    char a; bool b = 0;
    while(!isdigit(a = getchar()))
        b = a == '-';
    while(isdigit(a))
        x = x * 10 + a - '0', a = getchar();
    if(b)
        x = -x;
}

const int ms = 8e5 + 10;

const double pi = acos(-1.0);

int lx, ly, T, x, y, pos[ms];

struct node
{
    int a, b;
}poi[20];

struct in
{
    double r, i;
}tmpa[ms], tmpb[ms];

struct Poly
{
    int len;

    in a[ms];
}a, ans;


inline in operator + (in a, in b)
{
    return (in){a.r + b.r, a.i + b.i};
}

inline in operator - (in a, in b)
{
    return (in){a.r - b.r, a.i - b.i};
}

inline in operator * (in a, in b)
{
    return (in){a.r * b.r - a.i * b.i, a.r * b.i + a.i * b.r};
}

inline void fft(in *A, int len, int tp)
{
    for(ri i = 1; i < len; i ++)
        if(i < pos[i])
            swap(A[i], A[pos[i]]);
    in la, lb;
    for(ri i = 1; i < len; i <<= 1)
    {
        in wn = (in){cos(pi / i), tp * sin(pi / i)};
        for(ri j = 0; j < len; j += (i << 1))
        {
            in w = (in){1, 0};
            for(ri k = j; k < j + i; k ++)
            {
                la = A[k], lb = w * A[k + i], w = w * wn;
                A[k] = la + lb, A[k + i] = la - lb;
            }
        }
    }
}

inline void Mul(Poly &A, Poly B)
{
    int len = 1, num = 0;
    while(len < A.len + B.len)
        len <<= 1, num ++;
    for(ri i = 0; i < A.len; i ++)
        tmpa[i] = A.a[i];
    for(ri i = A.len; i < len; i ++)
        tmpa[i] = (in){0, 0};
    for(ri i = 0; i < B.len; i ++)
        tmpb[i] = B.a[i];
    for(ri i = B.len; i < len; i ++)
        tmpb[i] = (in){0, 0};
    for(ri i = 1; i < len; i ++)
        pos[i] = ((pos[i >> 1] >> 1) | ((i & 1) << (num - 1)));
    fft(tmpa, len, 1), fft(tmpb, len, 1);
    for(ri i = 0; i < len; i ++)
        tmpa[i] = tmpa[i] * tmpb[i];
    fft(tmpa, len, -1);
    for(ri i = 0; i < len; i ++)
        tmpa[i] = (in){tmpa[i].r / len, 0};
    while(len - 1 >= 0 && tmpa[len - 1].r == 0.0)
        len --;
    A.len = len;
    for(ri i = 0; i < len; i ++)
        A.a[i] = tmpa[i];
}

const double eps = 1.0e-12, K = 1 / sqrt(2 * pi);

/*inline double f(double x) { return K * exp(x * x / -2.0); }  //正态分布的概率密度函数
inline double psp(double l, double r) {
    double mid = (l + r) / 2.0;
    return (r - l) * (f(l) + 4.0 * f(mid) + f(r)) / 6.0;
}
inline double ask(double l, double r)  //然后这是Simpson积分的板子
{
    double mid = (l + r) / 2.0;
    double f1 = psp(l, r);
    double f2 = psp(l, mid) + psp(mid, r);
    double res = (-eps < f1 - f2 && f1 - f2 < eps) ? f1 : ask(l, mid) + ask(mid, r);
    return res;
}*/

inline double f(double x)
{
    return K * exp(x * x / -2.0);
}

inline double calc(double l, double r)
{
    double mid = (l + r) / 2;
    return (r - l) * (f(l) + 4.0 * f(mid) + f(r)) / 6.0;
}

double ask(double l, double r)
{
    double mid = (l + r) / 2;
    double f1 = calc(l, r), f2 = calc(l, mid) + calc(mid, r);
    double rt = (fabs(f1 - f2) < eps) ? f1 : ask(l, mid) + ask(mid, r);
    return rt;  
}

int main()
{
    re(T);
    while(T --)
    {
        re(x), re(y);
        int mx = 0;
        /**/if(x * y <= 524288)
        {
            for(ri i = 1; i <= 10; i ++)
                re(poi[i].a), re(poi[i].b), mx = max(mx, poi[i].b);
            a.len = x;
            double TT = x;
            for(ri i = 0; i < x; i ++)
                a.a[i].r = 1 / TT;
            ans.len = 1, ans.a[0].r = 1;
            int tmpy = y;
            while(tmpy)
            {
                if(tmpy & 1)
                    Mul(ans, a);
                Mul(a, a), tmpy >>= 1;
            }
            for(ri i = 1; i <= 10; i ++)
            {
                double pri = 0;
                for(ri j = poi[i].a; j <= poi[i].b; j ++)
                    pri += ans.a[j].r;
                printf("%.6lf\n", pri);
            }
        }
        else
        {
            double mu = (double)(x - 1) / 2, simga = (double)(x * x - 1) / 12;
            for(ri i = 1; i <= 10; i ++)
            {
                double aa, bb;
                scanf("%lf%lf", &aa, &bb);
                aa = (aa - y * mu) / sqrt(y * simga);
                bb = (bb - y * mu) / sqrt(y * simga);
                //cout << aa << ' ' << bb << '\n';
                printf("%.6lf\n", ask(0, bb) - ask(0, aa));
            }
        }
    }
}

但生成函数的形式还很多样,比如最开始的这个生成函数,我们现在想对它求个通项……至于求通项何用,后面会有例题的。
我们考虑利用高中的等比数列求和知识,错位相减
S_n=\sum_{i=0}^{n}x^i
xS_n=\sum_{i=0}^{n}x^{i+1}
(1-x)S_n=1-x^{n+1}
S_n = \frac{1-x^{n+1}}{1-x}
注意到在生成函数里面,我们并不关心x具体是几,而比较关心它的幂次和系数,所以这个x的大小就任我们所取。
我们考虑当n趋近于无限,x大于0小于1的时候,x^{n+1}可以近似看作0,那么
S_n=\frac{1}{1-x}
其他的,取奇数/偶数/k的倍数其实也是一个道理,只不过有的时候你需要提出一个x而已
我们还需要了解一个东西——负数情况下的二项式定理
(1+x)^{-n}=\sum_{i=0}^{\infty}=(-1)^iC_{n+i-1}^{i}x^i
(1-x)^{-n}=\sum_{i=0}^{\infty}=(-1)^iC_{n+i-1}^{i}(-x)^i=\sum_{i=0}^{\infty}=C_{n+i-1}^{i}x^i
那就上例题吧BZOJ3028食物
这就是我们刚才所写的一个集合
(1+x)^2(1+x+x^2)(1+x+x^2+x^3)(x+x^3+x^5+……)(1+x^4+x^8+……)(1+x^3+x^6+……)=\frac{x}{(1-x)^4}
\frac{x}{(1-x)^4}=\sum_{i=0}^{n}C_{n+3}^3 x^{i+1}
对于x^{n}来说,它的系数就是C_{n+2}^3
考虑卢卡斯,将n分解成kp+x,求C_{x+2}^3就是了

#include <bits/stdc++.h>
#include <algorithm>
#include <iostream>
#include <cstring>
#include <cstdlib>
#include <cstdio>
#include <vector>
#include <bitset>
#include <queue>
#include <cmath>
#include <map>
#include <set>
#define ri register int

using namespace std;

inline char gch()
{
    static char buf[100010], *h = buf, *t = buf;
    return h == t && (t = (h = buf) + fread(buf, 1, 100000, stdin), h == t) ? EOF : *h ++;
}

typedef unsigned int uint;

typedef long long lo;

template <typename int_qwq>

inline void re(int_qwq &x)
{
    x = 0;
    char a; bool b = 0;
    while(!isdigit(a = getchar()))
        b = a == '-';
    while(isdigit(a))
        x = x * 10 + a - '0', a = getchar();
    if(b)
        x = -x;
}

char s[1010];

const int mo = 10007;

int n, fac[11000], inv[11000];

inline int ksm(int x, int k)
{
    int rt = 1, a = x;
    for(; k; k >>= 1, a = a * a % mo)
        if(k & 1)
            rt = rt * a % mo;
    return rt;
}

int main()
{
    scanf("%s", s + 1);
    int len = strlen(s + 1);
    fac[1] = 1, fac[0] = 1, inv[1] = 1, inv[0] = 1;
    for(ri i = 2; i < mo; i ++)
        fac[i] = fac[i - 1] * i % mo;
    inv[mo - 1] = ksm(fac[mo - 1], mo - 2);
    for(ri i = mo - 2; i >= 2; i --)
        inv[i] = inv[i + 1] * (i + 1) % mo;
    for(ri i = 1; i <= len; i ++)
        n = (n * 10 + s[i] - '0') % mo;
    printf("%d", fac[(n + 2) % mo] * inv[3] % mo * inv[(n + mo - 1) % mo] % mo);
}

丢个大招luoguP3784[SDOI2017]遗忘的集合
这个题目先考虑S的生成函数,对于S内的某个元素i,其生成函数为[i\in S]\sum_{j=0}^{\infty}x^{ij}=(\frac{1}{1-x^i})^{a_i}
其中a_i=0/1
F(x)=\prod_{i=1}^{\infty}(\frac{1}{1-x^i})^{a_i}
现在我们已知F(x)的系数,想通过它来得到a_i
首先,多项式两边同时取-ln
-ln(F(x))=\sum_{i=1}^{\infty}a_i ln(1-x^i)
考虑对ln(1-x^i)进行先求导再积分(想办法换个形式表达)
ln(1-x^i)=\frac{-ix^{i-1}}{1-x^i}
这个时候我们再把\frac{1}{1-x^i}给还原回去。
-ix^{i-1}\sum_{j=0}^{\infty}x^{ij}=-\sum_{j=0}^{\infty}ix^{i(j+1)-1}
然后再利用高中数学积分回去
-\sum_{j=0}^{\infty}\frac{ix^{i(j+1)}}{i(j+1)}=-\sum_{j=1}^{\infty}\frac{x^{ij}}{j}
将这个结果带入到刚才的式子里面
lnF(x)=\sum_{i=1}^{\infty}a_i\sum_{j=1}^{\infty}\frac{x^{ij}}{j}
然后利用莫反的套路,设T = ij
lnF(x)=\sum_{T=1}^{\infty}x^{ij}\sum_{d|T}a_d\frac{d}{T}
然后f_x=\sum_{d|x}a_d\frac{d}{x}
xf_x=\sum_{d|x}a_x x
如果a_x x>0说明,x在S这个集合中。
然后调和级数的复杂度取算这个系数就行了。但是你要写mtt
loj上疯狂在tle的边缘旋转跳跃。真的考虑去背个新的mtt了。

#include <bits/stdc++.h>
#include <algorithm>
#include <iostream>
#include <cstring>
#include <cstdlib>
#include <cstdio>
#include <vector>
#include <bitset>
#include <queue>
#include <cmath>
#include <map>
#include <set>
#define ri register int

using namespace std;

inline char gch()
{
    static char buf[100010], *h = buf, *t = buf;
    return h == t && (t = (h = buf) + fread(buf, 1, 100000, stdin), h == t) ? EOF : *h ++;
}

typedef unsigned int uint;

typedef long long lo;

template <typename int_qwq>

inline void re(int_qwq &x)
{
    x = 0;
    char a; bool b = 0;
    while(!isdigit(a = getchar()))
        b = a == '-';
    while(isdigit(a))
        x = x * 10 + a - '0', a = getchar();
    if(b)
        x = -x;
}

const double pi = acos(-1.0);

const int ms = 1 << 20 | 1, Len = 1 << 19, mo = 32767;

int n, pos[ms];

lo p, F[ms], G[ms], C[ms], D[ms];

inline lo ksm(lo x, lo k)
{
    lo rt = 1, a = x;
    for(; k; k >>= 1, a = a * a % p)
        if(k & 1)
            rt = rt * a % p;
    return rt;
}

struct in
{
    double r, i;

    in(double rr = 0.0, double ii = 0.0)
    {
        r = rr, i = ii;
    }

    inline in conj()
    {
        return (in){r, -i};
    }
}ap[ms], P[ms], Q[ms], f[ms], g[ms], apf[ms], apg[ms];

inline in operator + (in a, in b)
{
    return (in){a.r + b.r, a.i + b.i};
}

inline in operator - (in a, in b)
{
    return (in){a.r - b.r, a.i - b.i};
}

inline in operator * (in a, in b)
{
    return (in){a.r * b.r - a.i * b.i, a.i * b.r + a.r * b.i};
}

inline in operator / (in a, in b)
{
    return (in){(a.r * b.r + a.i * b.i) / (b.r * b.r + b.i * b.i), (a.i * b.r - a.r * b.i) / (b.r * b.r + b.i * b.i)};
}

inline void fft(in *a, lo len, lo tp)
{
    for(ri i = 1; i < len; i ++)
        if(i < pos[i])
            swap(a[i], a[pos[i]]);
    in la, lb;
    for(ri i = 1; i < len; i <<= 1)
        for(ri j = 0; j < len; j += (i << 1))
            for(ri k = j; k < j + i; k ++)
            {
                in w(cos((k - j) * pi / i), tp * sin((k - j) * pi / i));
                la = a[k], lb = w * a[k + i];
                a[k] = la + lb, a[k + i] = la - lb;
            }
    if(tp == -1)
        for(ri i = 0; i < len; i ++)
            a[i].r /= 1.0 * len, a[i].i /= 1.0 * len;
}

inline void mtt(lo *a, lo *b, lo *c, lo len)
{
    int x = 1, num = 0;
    while(x < len + len)
        x <<= 1, num ++;
    for(ri i = 0; i < len; i ++)
    {
        a[i] %= p, b[i] %= p;
        ap[i] = (in){a[i] & mo, a[i] >> 15}, P[i] = (in){b[i] & mo, b[i] >> 15};
    }
    for(ri i = len; i < x; i ++)
        P[i] = (in){0, 0}, ap[i] = P[i];
    for(ri i = 0; i < x; i ++)
        pos[i] = ((pos[i >> 1] >> 1) | ((i & 1) << (num - 1)));
    fft(ap, x, 1), fft(P, x, 1);
    for(ri i = 0; i < x; i ++)
        Q[i] = P[(x - i) & (x - 1)].conj();
    //in g1(0, 2), g2(2, 0);
    in g1(0, -0.5), g2(0.5, 0);
    for(ri i = 0; i < x; i ++)
        f[i] = (P[i] + Q[i]) * g2, g[i] = (P[i] - Q[i]) * g1;
    for(ri i = 0; i < x; i ++)
        apf[i] = ap[i] * f[i], apg[i] = ap[i] * g[i];
    fft(apf, x, -1), fft(apg, x, -1);
    for(ri i = 0; i < len; i ++)
    {
        lo tmpa = (lo)(apf[i].r + 0.5) % p;
        lo tmpb = (((lo)(apf[i].i + 0.5) % p) << 15) % p;
        lo tmpc = (((lo)(apg[i].r + 0.5) % p) << 15) % p;
        lo tmpd = (((lo)(apg[i].i + 0.5) % p) << 30ll) % p;
        //cout << tmpa << ' ' << tmpb << ' ' << tmpc << ' ' << tmpd << '\n';
        c[i] = (tmpa + tmpb + tmpc + tmpd) % p;
    }
    for(ri i = len; i < x; i ++)
        c[i] = 0;
}

void merge(lo *A, lo *B, lo len)
{
    if(len == 1)
    {
        A[0] = ksm(B[0], p - 2); return;
    }
    merge(A, B, (len + 1) >> 1);
    mtt(A, B, C, len), mtt(A, C, D, len);
    for(ri i = 0; i < len; i ++)
        A[i] = (A[i] + A[i]) % p;
    for(ri i = 0; i < len; i ++)
        A[i] = ((A[i] - D[i]) % p + p) % p;
}

inline void askln(lo *A, lo n, lo *B)
{
    lo tmp[ms], inv[ms];
    for(ri i = 1; i <= n; i ++)
        tmp[i - 1] = A[i] * i % p;
    memset(inv, 0, sizeof(inv));
    merge(inv, A, n + 1);
    mtt(inv, tmp, B, n);
    //for(ri i = 0; i <= n; i ++)
    //  cout << B[i] << ' ';
    //cout << '\n';
    for(ri i = n; i; i --)
        B[i] = B[i - 1] * ksm(i, p - 2) % p;
}

int main()
{
    re(n), re(p), F[0] = 1;
    memset(G, 0, sizeof(G));
    for(ri i = 1; i <= n; re(F[i ++]));
    askln(F, n, G);
    for(ri i = 1; i <= n; i ++)
        G[i] = G[i] * i % p;
    for(ri i = 1; i <= n; i ++)
        for(ri j = (i << 1); j <= n; j += i)
            G[j] += p - G[i], G[j] -= (G[j] >= p) ? p : 0;
    int ans = 0;
    for(ri i = 1; i <= n; i ++)
        if(G[i])
            ans ++;
    printf("%d\n", ans);
    for(ri i = 1; i <= n; i ++)
        if(G[i])
            printf("%d ", i);
}

发表评论

电子邮件地址不会被公开。 必填项已用*标注