码迷,mamicode.com
首页 > 其他好文 > 详细

【GDOI2020模拟03.11】我的朋友们(多项式求逆+分治NTT)

时间:2020-03-12 23:39:15      阅读:79      评论:0      收藏:0      [点我收藏+]

标签:else   turn   我的朋友   要求   alt   ini   大于   题解   统计   

题目大意:


技术图片
\(1 \le l \le n \le 10^5\)

题解:


写完之后觉得我死了,这应该是我做过的最复杂的分治NTT。

\(p[i]\)表示第\(i\)个点概率,为了表达方便,先把\(p[i]\) reverse一下。

\(P(i)=p[i]·x+1-p[i]\)

考虑设\(f[i]\)表示从第\(i\)个点出发的期望步数。
\(i<l,f[i]=0\)

\(H_i(x)=\prod_{j=i-l+1}^{i}P(j)\)
\(i\ge l,f[i]=1+\sum_{j=0}^l f[i-j]*H_i[x^j]\)

\(j=0\)时,\(f[i-j]=f[i]\),把这一项移到等式右边,解个方程可得:
\(i \ge l,f[i]=(1+\sum_{j=1}^l f[i-j]*H_i[x^j])/(1-H_i[x^0])\)

得到了\(O(n^2)\)的做法。


\(F_i=\sum_{j=1}^i f[j]*x^j\)
\(f[i]=F_{i-1}*H_i[x^i]\)

考虑分治求\(f[]\),现在已经求出了\(f[1..x-1]\),要求\(f[x..y]\)

在这个分治的过程中,动态的维护两个函数\(F[x][y],G[x][y]\)

\(G[x][y]=\prod_{j=y-l+1}^{x}P(j)\)

它的含义即为\(f[x,y]\)的转移里都会用到\(P\)的乘积。

注意到\(y-l+1\)可能会大于\(x\),此时我们把意义推广,设\(SP\)\(P\)的前缀积,则:
\(G[x][y]=\prod_{j=y-l+1}^{x}P(j)=SP[x]/SP[y-l]\)

\(F[x][y]=F_{x-1}*G[x][y]\)

如果能一直维护这两个函数,那么当\(x=y\)时直接取\(F[x][y]\)\(x\)项的系数就好了。

\(m=(x+y)/2\)

1.由\([x,y]\)递归至\([x,m]\)

\(G[x][m]=G[x][y]*\prod_{i=m}^{r-1}P(i-l+1)=G[x][y]*\prod_{i=m+1}^{r}P(i-l)\)

\(F[x][m]=F[x][y]*\prod_{i=m+1}^{r}P(i-l)\)

2.由\([x,y]\)递归至\([m+1,y]\)

\(G[m+1][y]=G[x][y]*\prod_{i=x+1}^{m+1}P(i)=G[x][y]*\prod_{i=x}^{m}P(i+1)\)

\(F[m+1][y]=F[x][y]*\prod_{i=x}^{m}P(i+1)+(F_m-F_{x-1})*G[m+1][y]\)

发现每次乘上的\(\prod P\)都可以通过分治NTT预处理出来

当然这么做并没有使复杂度降低。

观察\(F[x][y],G[x][y]\),在递归的过程中一直乘\(\prod P\),也就是次数最多+\(y-x\)

最后要统计答案取的那一项的次数也在\([x,y]\)里。

那么只维护\({F[x][y]\over x^{max(0,x-(y-x))}}\)的前\(2(y-x)\)项和\(G[x][y]\)的前\(y-x\)项就可以了。

因为一切多项式的长度都是\(O(y-x)\),所以时间复杂度是\(O(n~log^2~n)\)

最后一个问题就是\(F[1][n]\)\(G[1][n]\)怎么求?

\(F[1][n]\)显然是0.

\(G[1][n]=SP[1]/SP[n-l]\)

分治NTT之后再多项式求逆就可以求出\(G[1][n]\)

Code:

#include<bits/stdc++.h>
#define fo(i, x, y) for(int i = x, _b = y; i <= _b; i ++)
#define ff(i, x, y) for(int i = x, _b = y; i <  _b; i ++)
#define fd(i, x, y) for(int i = x, _b = y; i >= _b; i --)
#define ll long long
#define pp printf
#define hh pp("\n")
using namespace std;

const int mo = 998244353;

ll ksm(ll x, ll y) {
    ll s = 1;
    for(; y; y /= 2, x = x * x % mo)
        if(y & 1) s = s * x % mo;
    return s;
}

#define V vector<ll>
#define re resize
#define si size()
const int nm = 1 << 18;
namespace ntt {
    ll w[nm], a[nm], b[nm]; int r[nm];
    void build() {
        for(int i = 1; i < nm; i *= 2) {
            ll v = ksm(3, (mo - 1) / 2 / i);
            w[i] = 1; ff(j, 1, i) w[i + j] = w[i + j - 1] * v % mo;
        }
    }
    void dft(ll *a, int n, int f) {
        ff(i, 0, n) {
            r[i] = r[i / 2] / 2 + (i & 1) * (n / 2);
            if(i < r[i]) swap(a[i], a[r[i]]);
        } ll b;
        for(int i = 1; i < n; i *= 2) for(int j = 0; j < n; j += 2 * i)
            ff(k, 0, i) b = a[i + j + k] * w[i + k], a[i + j + k] = (a[j + k] - b) % mo, a[j + k] = (a[j + k] + b) % mo;
        if(f == -1) {
            reverse(a + 1, a + n);
            b = ksm(n, mo - 2);
            ff(i, 0, n) a[i] = (a[i] + mo) * b % mo;
        }
    }
    void fft(V &p, V &q) {
        int p0 = p.si + q.si - 1, n = 1;
        for(; n < p0; n *= 2);
        ff(i, 0, n) a[i] = b[i] = 0;
        ff(i, 0, p.si) a[i] = p[i];
        ff(i, 0, q.si) b[i] = q[i];
        dft(a, n, 1); dft(b, n, 1);
        ff(i, 0, n) a[i] = a[i] * b[i] % mo;
        dft(a, n, -1);
        p.re(p0);
        ff(i, 0, p0) p[i] = a[i];
    }
}

V operator * (V p, V q) {
    ntt :: fft(p, q);
    return p;
}
void dft(V &p, int f) {
    ff(i, 0, p.si) ntt :: a[i] = p[i];
    ntt :: dft(ntt :: a, p.si, f);
    ff(i, 0, p.si) p[i] = ntt :: a[i];
}
V qni(V a) {
    int a0 = a.si, n0 = 1;
    while(n0 < a0) n0 *= 2;
    V b; b.re(1); b[0] = ksm(a[0], mo - 2);
    for(int n = 2; n <= n0; n *= 2) {
        V d = b; d.re(n); b.re(2 * n);
        V c = a; c.re(n); c.re(2 * n);
        dft(c, 1); dft(b, 1);
        ff(i, 0, b.si) b[i] = c[i] * b[i] % mo * b[i] % mo;
        dft(b, -1); b.re(n);
        ff(i, 0, b.si) b[i] = (2 * d[i] - b[i] + mo) % mo;
    }
    b.re(a0);
    return b;
}
V operator + (V p, V q) {
    p.re(max(p.si, q.si));
    ff(i, 0, q.si) p[i] = (p[i] + q[i]) % mo;
    return p;
}

const int N = 1e5 + 5;

int n, l;
ll x, y, p[N];
ll sp[N], np[N];

void Init() {
    scanf("%d %d", &n, &l);
    fo(i, 1, n) {
        scanf("%lld %lld", &x, &y);
        p[i] = x * ksm(y, mo - 2) % mo;
    }
    reverse(p + 1, p + n + 1);
    sp[0] = 1; fo(i, 1, n) sp[i] = sp[i - 1] * (1 - p[i]) % mo;
    np[n] = ksm(sp[n], mo - 2); fd(i, n, 1) np[i - 1] = np[i] * (1 - p[i]) % mo;
}

ll calc(int x) {
    return sp[x] * np[x - l] % mo;
}

V t1[N * 4], t2[N * 4];

#define i0 i + i
#define i1 i + i + 1
void dg(int i, int x, int y) {
    if(x > y) return;
    if(x == y) {
        t1[i].re(2);
        t1[i][0] = 1 - p[x + 1];
        t1[i][1] = p[x + 1];
        
        ll v = x >= l ? p[x - l] : 0;
        t2[i].re(2);
        t2[i][0] = 1 - v;
        t2[i][1] = v;
        return;
    }
    int m = x + y >> 1;
    dg(i0, x, m); dg(i1, m + 1, y);
    t1[i] = t1[i0] * t1[i1];
    t2[i] = t2[i0] * t2[i1];
}

V g;

V dfs(int x, int y) {
    V a;
    if(x > y) {
        a.re(1); a[0] = 1;
        return a;
    }
    if(x == y) {
        a.re(2); a[0] = 1 - p[x]; a[1] = p[x];
        return a;
    }
    int m = x + y >> 1;
    return dfs(x, m) * dfs(m + 1, y);
}

ll a[N];

void zy(V &a, int b) {
    ff(i, b, a.si) a[i - b] = a[i];
    a.re(a.si - b);
}

void yy(V &a, int b) {
    int sa = a.si;
    a.re(sa + b);
    fd(i, sa - 1, 0) a[i + b] = a[i], a[i] = 0;
}

void qz(V &a, int b) {
    if(a.si > b + 1) a.re(b + 1);
}

V fz(int i, int x, int y, V f, V g, int w) {
    V b; b.re(1); b[0] = 0;
    if(x > y) return b;
    
    int nw = max(0, x - (y - x));
    if(nw > w) zy(f, nw - w);
    qz(f, 2 * (y - x)); qz(g, (y - x));
    
    if(x == y) {
        if(x < l) {
            a[x] = 0;
        } else {
            a[x] = (f[0] + 1) * ksm(1 - calc(x), mo - 2) % mo;
        }
        b.re(1); b[0] = a[x];
        return b;
    }
    
    int m = x + y >> 1;
    
    V nf = f * t2[i1], ng = g * t2[i1];
    b = fz(i0, x, m, nf, ng, nw);
    
    nf = f * t1[i0], ng = g * t1[i0];
    qz(ng, y - x);
    V d = b * ng; yy(d, x - nw);
    nf = nf + d;
    V c = fz(i1, m + 1, y, nf, ng, nw);
    
    yy(c, m - x + 1);
    return b + c;
}

int main() {
    freopen("friends.in", "r", stdin);
    freopen("friends.out", "w", stdout);
    ntt :: build();
    Init();
    dg(1, 1, n);
    V g;
    if(l == n) g = dfs(1, 1); else {
        g = dfs(2, n - l);
        g.re(n + 1);
        g = qni(g);
    }
    V f; f.re(n + 1); fo(i, 0, n) f[i] = 0;
    fz(1, 1, n, f, g, 0);
    pp("%lld\n", (a[n] + mo) % mo);
}

【GDOI2020模拟03.11】我的朋友们(多项式求逆+分治NTT)

标签:else   turn   我的朋友   要求   alt   ini   大于   题解   统计   

原文地址:https://www.cnblogs.com/coldchair/p/12483517.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!