BNUOJ 51637 Quailty's Life | BZOJ 4587 推箱子

题目

https://acm.bnu.edu.cn/v3/problem_show.php?pid=51637
http://www.lydsy.com/JudgeOnline/problem.php?id=4587

题意

有一个(n+1)×(m+1)的网格，初始在点(0,0)，每次可以从(x,y)走到(x+1,y)或(x,y+1)或(x+1,y+1)。有k个点(ai,bi)(1≤i≤k)不能到达，问从初始点走到横坐标为n或纵坐标为m的方案数模p，p不一定是质数。注意当横坐标为n或纵坐标为m时不能再走动。
1≤n≤109,1≤m≤104,0≤k≤10,1≤p<231,∑m≤105

题解

1.一个暴力

定义dp[i][j]表示从(0,0)走到(i,j)的方案数，定义横坐标为n或纵坐标为m的(i,j)为终止状态。
初始状态dp[0][0]=1。
如果(i,j)不能到达则dp[i][j]=0。
如果(i,j)能到达，dp[i][j]包含dp[i−1][j−1]，如果(i−1,j)不是终止状态则dp[i][j]也包含dp[i−1][j]，如果(i,j−1)不是终止状态则dp[i][j]也包含dp[i][j−1]。
答案是所有终止状态的dp[i][j]之和。
状态转移，只涉及加减，模运算很好套用，时间复杂度O(nm)。

int ans = 0, f[n + 1][m + 1];for(int i = 0; i <= n; ++i)    for(int j = 0; j <= m; ++j)    {        if(canreach[i][j] == -1)        {            f[i][j] = 0;            continue;        }        if(i == 0 && j == 0)        {            f[i][j] = 1;            continue;        }        long long tmp = 0;        if(i > 0 && j < m)            tmp += f[i - 1][j];        if(i < n && j > 0)            tmp += f[i][j - 1];        if(i > 0 && j > 0)            tmp += f[i - 1][j - 1];        f[i][j] = tmp % mod;        if(i == n || j == m)            ans += f[i][j];    }

用于对拍

2. k=0的情况

考虑从(0,0)走到(i,j)的方案数，然后直接计算i=n和j=m的情况。
设从(0,0)走到(i,j)的一种方案使用了a次x+1，b次y+1，c次x+1同时y+1，则有

⎧⎩⎨a+c=ib+c=ja,b,c≥0⇒⎧⎩⎨a=i−cb=j−c0≤c≤min(i,j)

对应的方案数是

(i+j−c)!(i−c)!(j−c)!c!=(i+j−ci)(ic)=(j+i−cj)(jc)

考虑i=n的情况，则0≤j≤m的情况为

∑j=0m∑c=0min(n,j)(n+j−cn)(nc)=∑c=0min(n,m)(nc)∑j=cm(n+j−cn)=∑c=0min(n,m)(nc)(n+1+j−cn+1)

那么枚举c就可以了，不用枚举j。对于0≤i≤n,j=m的情况也这么分析一下。
不过这个题不是求无阻碍到(n,∗)或(∗,m)的方案数，因为(n,i)走不到(n,i+1)。
但是到终止状态前的一个点是无阻碍的，枚举前一个点算贡献即可。
考虑怎么计算(nc)modp。一种方法是将p分解成质因子的幂次pk，O(c)枚举分子分母，通过消p因子的方法算出模值，再用中国剩余定理合并。还有一种方法是不分解成质因子的幂次，直接在模p意义下做，记录与p互质的部分的值，记录不与p互质的部分每个质因子的幂次，两者复杂度相近，但编程复杂度低。
由于还要枚举c，看上去做法是O(m2logp)的，实际上可以由(nc)在O(logp)的时间内推到(nc+1)（乘一个数除一个数），所以复杂度是O(mlogp)的。

3. k>0的情况

3.1 我会容斥？

”不经过任意一个坏点的方案数“可以通过计算”至少经过其中k个已知坏点的方案数“容斥得到。
枚举坏点集合S的子集T，计算从(0,0)走过T中每个坏点后走到终止状态的方案数f(T)，该方案数对答案的贡献是(−1)|T|⋅f(T)，这个方案不为0的情况一定是坏点可以按坐标排序后顺次经过，需要计算从一个点(i,j)走到(i′,j′)的方案数，以及从最后一个坏点走到终止状态的方案数，前者依然可以枚举c来算。
看上去时间复杂度是O(2kmlogp)的，容易构造数据使其超时。

3.2 我会容斥！

对于每个坏点i，定义f[i]表示走到这个点且之前没碰到坏点的方案数，这个就好容斥了，枚举在这个坏点左上角的坏点j，不合法的方案是f[j]乘上从(aj,bj)到(ai,bi)的方案数，这些方案里遇到的第一个坏点是j，所以不重不漏，最后用从(0,0)到(ai,bi)的方案数减一下即可。
答案也可以这么统计，算从(0,0)到终止状态的方案，减去所有f[i]乘上从(ai,bi)到终止状态的方案即可。
看上去时间复杂度是O(k2mlogp)的，好好写就不会超时了。
注意坏点在终止状态上的情况。写个暴力对拍一下就可以发现问题

代码

#include <cstdio>#include <algorithm>typedef long long LL;const int maxp = 46341, maxc = 11, maxk = 11, maxe = 51;void exgcd(int a, int b, int &x, int &y){    if(!b)    {        x = 1;        y = 0;        return;    }    exgcd(b, a % b, y, x);    y -= a / b * x;}int mod_inv(int x, int p){    int s, t;    exgcd(x, p, s, t);    return s < 0 ? s + p : s;}int tot, prime[maxp], fir[maxp], mod, cnt, fact[maxc], Exp[maxc], Coeff, Lim[maxc], Pw[maxc][maxe];inline void mod_dec(int &x, int y, int c = 1){    while(c--)        if((x -= y) < 0)            x += mod;}void init(){    Coeff = 1;    for(int i = 0; i < cnt; ++i)    {        Exp[i] = Lim[i] = 0;        Pw[i][0] = 1;    }}void update(int val, int flag){    for(int i = 0; i < cnt; ++i)        for( ; val % fact[i] == 0; Exp[i] += flag, val /= fact[i]);    Coeff = (LL)Coeff * (flag == 1 ? val : mod_inv(val, mod)) % mod;}int query(){    int ret = Coeff;    for(int j = 0; j < cnt && ret; ++j)    {        if(!Exp[j])            continue;        for( ; Lim[j] < Exp[j]; ++Lim[j])            Pw[j][Lim[j] + 1] = (LL)Pw[j][Lim[j]] * fact[j] % mod;        ret = (LL)ret * Pw[j][Exp[j]] % mod;    }    return ret;}int calc_1(int n, int m){    if(n > m)        std::swap(n, m);    if(n < 0)        return 0;    else if(!n)        return 1;    int ret = 0;    init();    for(int i = 1; i <= n; ++i)    {        update(m + i, 1);        update(i, -1);    }    for(int i = 0; i <= n; ++i)    {        mod_dec(ret, mod - query());        if(i == n)            break;        update(n - i, 1);        update(m - i, 1);        update(n + m - i, -1);        update(i + 1, -1);    }    return ret;}int calc_2(int n, int m){    if(n > m)        std::swap(n, m);    if(n < 0)        return 0;    else if(!n)        return 1;    int ret = 0;    --n;    --m;    init();    for(int i = 1; i <= n; ++i)    {        update(m + i, 1);        update(i, -1);    }    update(n + m + 1, 1);    for(int i = 0; i <= n; ++i)    {        update(n + 1, -1);        mod_dec(ret, mod - query(), 2);        update(n + 1, 1);        update(m + 1, -1);        mod_dec(ret, mod - query(), 2);        update(m + 1, 1);        update(n + m + 1 - i, -1);        mod_dec(ret, query());        if(i == n)            break;        update(n - i, 1);        update(m - i, 1);        update(i + 1, -1);    }    return ret;}int t, n, m, k, f[maxk], ans;std::pair<int, int> lim[maxk];int main(){    for(int i = 2; i < maxp; ++i)    {        if(!fir[i])            prime[tot++] = fir[i] = i;        for(int j = 0, k; (k = i * prime[j]) < maxp; ++j)        {            fir[k] = prime[j];            if(fir[i] == prime[j])                break;        }    }    scanf("%d", &t);    while(t--)    {        scanf("%d%d%d%d", &n, &m, &k, &mod);        int tmp = mod;        cnt = 0;        for(int i = 0; i < tot && prime[i] * prime[i] <= tmp; ++i)            if(tmp % prime[i] == 0)                for(fact[cnt++] = prime[i]; tmp % prime[i] == 0; tmp /= prime[i]);        if(tmp > 1)            fact[cnt++] = tmp;        for(int i = 0; i < k; ++i)            scanf("%d%d", &lim[i].first, &lim[i].second);        std::sort(lim, lim + k);        k = std::unique(lim, lim + k) - lim;        for(int i = 0; i < k; ++i)        {            int &ai = lim[i].first, &bi = lim[i].second;            if(ai < n && bi < m)            {                f[i] = calc_1(ai, bi);                for(int j = 0; j < i; ++j)                {                    int &aj = lim[j].first, &bj = lim[j].second;                    int coeff = calc_1(ai - aj, bi - bj);                    if(coeff)                        f[i] = (f[i] - (LL)coeff * f[j]) % mod;                }            }            else            {                f[i] = calc_1(ai - 1, bi - 1);                if(ai < n)                    mod_dec(f[i], mod - calc_1(ai, bi - 1));                if(bi < m)                    mod_dec(f[i], mod - calc_1(ai - 1, bi));                for(int j = 0; j < i; ++j)                {                    int &aj = lim[j].first, &bj = lim[j].second;                    int coeff = calc_1(ai - 1 - aj, bi - 1 - bj);                    if(ai < n)                        mod_dec(coeff, mod - calc_1(ai - aj, bi - 1 - bj));                    if(bi < m)                        mod_dec(coeff, mod - calc_1(ai - 1 - aj, bi - bj));                    if(coeff)                        f[i] = (f[i] - (LL)coeff * f[j]) % mod;                }            }        }        ans = calc_2(n, m);        for(int i = 0; i < k; ++i)        {            int coeff = calc_2(n - lim[i].first, m - lim[i].second);            if(coeff)                ans = (ans - (LL)coeff * f[i]) % mod;        }        if(ans < 0)            ans += mod;        printf("%d\n", ans);    }    return 0;}