容斥原理 训练笔记

​

容斥原理

设S是一个有限集，A_1,A_2…A_n是S的n个子集，则

$|S-{\bigcup }_{i=1}^n{A}_i|={\sum }_{i=0}^n(-1{)}^i{\sum }_{1\le j_1<j_2...<j_i\le n}|{\bigcap }_{k=1}^i{A}_{j_k}|$

基本应用：

m件不同的物品，分给n个人，要求每一个人至少分得一件物品，求不同的分配方案数

令$A_i$表示第i个人没有物品，$S$表示$m$个物品分给$n$个人的总方案数

则$|S-{\bigcup }_{i=1}^n{A}_i|={\sum }_{i=0}^n(-1{)}^i{\sum }_{1\le j_1<j_2...<j_i\le n}|{\bigcap }_{k=1}^i{A}_{j_k}|$

$={\sum }_{i=0}^n(-1{)}^i(\genfrac{}{}{0ex}{}{n}{i})(n-i{)}^m$

有 $2n$个元素 $a_1,a_2,...,a_n$ 和$b_1,b_2,...,b_n$，求有多少个它们的全排列，满足任意的

$1\le i\le n$,$a_i$ 和$b_i$ 都不相邻。

同样的，令$A_i$表示$a_i$和$b_i$相邻，

则$|S-{\bigcup }_{i=1}^n{A}_i|={\sum }_{i=0}^n(-1{)}^i{\sum }_{1\le j_1<j_2...j_i\le n}|{\bigcap }_{k=1}^i{A}_{j_k}|$
$={\sum }_{i=0}^n(-1{)}^i(\genfrac{}{}{0ex}{}{n}{i})|有i对相邻的方案数|$
$={\sum }_{i=0}^n(-1{)}^i(\genfrac{}{}{0ex}{}{n}{i})2^i\ast (2n-i)!$

有了以上基本常识就可以上大招了

例题

CF449D

大意：
给出一个长度为n的序列$a_1,a_2...a_n$。求从中选择一个非空子集使得他们的按位与之和等于0的方案数

思路：
不考虑复杂度的话我们有一个非常套路的容斥做法。考虑性质Ai表示子集与之后第i位为1，那么我们的答案其实就是
$|\Omega -A_1\bigcup A_2...\bigcup A_{20}|={\sum }_{i=0}^{20}(-1{)}^i{\sum }_{1\le j_1<j_2...<j_i\le 20}|A_{j_1}\bigcup A_{j_2}...A_{j_i}|$
其中||符号就表示集合的大小

显然就可以状压枚举，这样的时间复杂度是$O(n\ast 1e6)$，考虑优化。

注意到对于$|A_{j_1}\bigcup A_{j_2}...A_{j_i}|$，我们记满足对应所有性质的元素的个数为k，则该集合的大小就是$2^k-1$,那么什么元素会满足这些性质呢？就是二进制为其超集的元素呗，其价值就是1.

所以我们只要做一遍超集后缀和即可，时间复杂度来到$O(20\ast 1e6)$

code

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define endl '\n'
const ll N=1e5+10;
const ll mod=1e9+7;
ll n,cnt=0,a;
ll mas[N];
ll up=20;
ll vis[30];
ll dp[(1<<20)+10];
ll ksm(ll x,ll y)
{
    ll ans=1;
    while(y)
    {
        if(y&1) ans=ans*x%mod;
        x=x*x%mod;
        y>>=1;
    }
    return ans;
}
ll gt()
{
    ll tot=0;
    ll fl;
    for(int i=1;i<=n;++i)
    {
        fl=1;
        for(int j=0;j<up;++j)
        {
            if(!vis[j]) continue;
            if((mas[i]&(1<<j))==0)
            {
                fl=0;
                break;
            }
        }
        if(fl) tot++;
    }

    return ((ksm(2,tot)-1)%mod+mod)%mod;
}
void solve()
{
    cin>>n;
    for(int i=1;i<=n;++i) cin>>a,dp[a]++;
    for(int j=0;j<up;++j)
    {
        for(int i=0;i<(1<<up);++i)
        {
            if((i&(1<<j))==0) dp[i]+=dp[i^(1<<j)];
        }
    }
    ll ans=0;
    for(int s=0;s<(1<<up);++s)
    {
        cnt=0;
        for(int i=0;i<up;++i) if(s&(1<<i)) cnt++;
        if(cnt%2) ans=((ans-ksm(2,dp[s])+1)%mod+mod)%mod;
        else ans=((ans+ksm(2,dp[s])-1)%mod+mod)%mod;
    }
    cout<<ans<<endl;
}
int main()
{
    ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);
    solve();
    return 0;
}

CF1799G

大意：

思路：
个人认为这道题还是很顶的

题目有两个限制：每个人要得到$c_i$票，每个人不能投给自己组

然后这里就有一个很巧妙的转化：我们可以先将同一个组的人放在一起来看，这样就可以将第一个限制稍微弱化一些，最后我们再想办法处理组内的方案数

那么此时的问题如下：每一个组$i$总共有$nu{m}_i$人,总共需要$f_i$票，组内的人不能投给自己组，求合法方案数

考虑生成函数处理

我们用$a_i$来表示第$i$组的变元。那么对于第$i$组的每一个人来说，它能投给除了自己组外的任河组，所以它的贡献是$a_1+a_2+...+a_{i-1}+a_{i+1}+...+a_n$,我们可以简单记为$s-a_i$,其中s就表示${\sum }_{i=1}^n{a}_i$

那么我们最后得到式子为$(s-a_1{)}^{nu{m}_1}(s-a_2{)}^{nu{m}_2}...(s-a_n{)}^{nu{m}_n}$,记为$S$,而我们的答案显然就是$[a_1^{f_1}{a}_2^{f_2}...a_n^{f_n}](S)$.此时不难发现，对于$a_i$的指数，每一个s都可以提供1，而$(s-a_i{)}^{nu{m}_i}$中的$-a_i$也可以提供不多于$nu{m}_i$的指数。

所以我们考虑每一个$a_i$的次数${\det }_i(de{t}_i\le nu{m}_i,de{t}_i\le f_i)$,则每一个i有一个从$nu{m}_i$中选择$de{t}_i$的方案数$\left(\genfrac{}{}{0ex}{}{nu{m}_i}{de{t}_i}\right)$,并且$de{t}_i$会提供$(-1{)}^{de{t}_i}$的系数，然后剩下的所有指数都由$s$提供，总共是$\sum f_i-\sum de{t}_i=n-\sum de{t}_i$,要分成若干堆，第i堆有$f_i-de{t}_i$个，所以其实是一个可重集，不同组之间是相乘的关系，那么最终的答案其实就是
$$ans=\sum _{de{t}_i\le f_i}(-1{)}^{\sum de{t}_i}\prod (\left(\genfrac{}{}{0ex}{}{nu{m}_i}{de{t}_i}\right))(\genfrac{}{}{0ex}{}{(n-\sum de{t}_i)!}{(f_1-de{t}_1)!(f_2-de{t}_2)!...(f_n-de{t}_n)!})$$
$$=\prod nu{m}_i\ast \sum _{de{t}_i\le f_i}(-1{)}^{\sum de{t}_i}\frac{(n-\sum de{t}_i)!}{\prod de{t}_i!(nu{m}_i-de{t}_i)!(f_i-de{t}_i)!}$$

那么其实n只有200，所以我们可以比较轻松地来得到这个式子的结果

考虑枚举$\sum de{t}_i$,我们记为$d$,那么再记一个$d{p}_{i,j}$表示当前处理到了第i组，${\sum }_{x\le i}de{t}_i=j$,显然只要在过程中枚举合法的$de{t}_i$就能完成这个dp了，外层枚举d，内层枚举i,j,最内层枚举$de{t}_i$,乍一看时间复杂度好像是$O(n^4)$,但是因为$de{t}_i\le nu{m}_i$,而$\sum nu{m}_i=n$,所以实际复杂度只有$O(n^3)$

这样我们就完成了组与组之间的关系。考虑组内部的票数分配，$nu{m}_i$个人，每一个人需要$c_i$票，总共$\sum c_i=f_i$票，所以这其实还是一个多重集，那么我们只要乘上系数$\frac{f_i!}{\prod c_i!}$即可

推推式子还是有点累的~

code

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define endl '\n'
const ll N=210;
const ll mod=998244353;
ll n;
ll f[N],t[N],c[N],num[N];//每一组的总期望票数，组别，个人期望票数,每组人数
ll dp[N][N];
ll ksm(ll x,ll y)
{
    ll ans=1;
    while(y)
    {
        if(y&1) ans=ans*x%mod;
        x=x*x%mod;
        y>>=1;
    }
    return ans;
}
ll inv(ll x)
{
    return ksm(x,mod-2);
}
ll p[N],pp[N];
void init(ll n)
{
    p[0]=1;
    for(ll i=1;i<=n;++i) p[i]=p[i-1]*i%mod;
    pp[n]=inv(p[n]);
    for(int i=n-1;i>=0;--i) pp[i]=pp[i+1]*(i+1)%mod;

}
void solve()
{
    init(200);
    // for(int i=1;i<=10;++i) cout<<p[i]<<" "<<pp[i]<<endl;
    cin>>n;
    for(int i=1;i<=n;++i) cin>>c[i];//期望票数
    for(int i=1;i<=n;++i) cin>>t[i],f[t[i]]+=c[i];
    for(int i=1;i<=n;++i) num[t[i]]++;
    ll ans=1,sum=0;
    for(int i=1;i<=n;++i) ans=ans*p[num[i]]%mod;
    
    for(int d=0;d<=n;++d)//det_i之和
    {
        for(int i=1;i<=n;++i) for(int j=0;j<=n;++j) dp[i][j]=0;
        dp[0][0]=p[n-d];

        for(int i=1;i<=n;++i)
        {
            //当前枚举到第i组
            for(int j=0;j<=d;++j)//到第i个人位置det_i的总和为j
            {
                for(int k=0;k<=min(f[i],num[i])&&k<=j;++k)//det_i=k
                {
                    ll det=pp[k]*pp[num[i]-k]%mod*pp[f[i]-k]%mod;
                    dp[i][j]=(dp[i][j]+det*dp[i-1][j-k]%mod)%mod;
                }
            }
        }
        if((d%2)==0) sum=(sum+ans*dp[n][d]%mod)%mod;
        else sum=((sum-ans*dp[n][d]%mod)%mod+mod)%mod;
    }
    for(int i=1;i<=n;++i)//组内多重集
    {
        sum=sum*p[f[i]]%mod*pp[c[i]]%mod;
    }
    cout<<sum<<endl;
}
int main()
{
    ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);
    solve();
    return 0;
}

http://acm.hdu.edu.cn/contest/problem?cid=1102&pid=1011
23HDU 09K Cargo
大意：

n家店，每家店卖一种商品，总共有m种商品在卖。k次操作，每次随机选择一家店并且买下对应种类的商品。对于每一种商品i，记ci为售卖该种商品的店数，如果k次操作手里刚好有ci种i商品，且分别来自不同的店，则该商品合法。问k次操作之后所有种类的商品都不合法的方案数%998244353

$n,m\le 1e5,k<998244353$

思路：

可能数据范围并不是很支持直接容斥，但是可以发现基本上就是容斥的样子，所以试试看优化容斥。

令$A_i$表示第种商品合法的方案数，

则$|S-{\bigcup }_{i=1}^m{A}_i|={\sum }_{i=0}^m(-1{)}^i{\sum }_{1\le j_1<j_2...j_i}|{\bigcap }_{k=1}^i{A}_{j_k}|$

只考虑$A_i$，方案数为$A_n^{c_i}(n-c_i{)}^{k-c_i}$,意义显然。多个$A_i$取交的话，方案数就是$A_n^{\sum c_i}(n-\sum c_i{)}^{k-\sum c_i}$

发现这个方案数只与$\sum c_i$有关，而与顺序，某一个具体值无关，所以我们可以改变一下思路，转为枚举$\sum c_i$

具体来说，对于每一个$\sum c_i=j$,它的贡献就是用偶数个$c_i$累加得到j的方案数减去用奇数个$c_i$累加得到j的方案数，（这一步其实就是容斥原理公式的转化）

既然如此，我们考虑生成函数$1-x^{c_i}$，表示对于每一个$c_i$，取一个数得到和为$c_i$的方案数为-1，这是因为1是一个奇数，那么显然，在$\prod (1-x^{c_i})$中，如果某一个指数是由偶数个数累加得到的，其系数自然会是正数，奇数同理。由此该生成函数就满足了我们的条件，我们只要对每一项$[x_j](\prod (1-x^{c_i}))$乘上对应的方案数系数$A_n^{\sum c_i}(n-\sum c_i{)}^{k-\sum c_i}$即可.

个人感觉这种题目应该还是有点套路可循的，因为如果数据范围小一点的话，就是明显显的容斥，然后找到性质用生成函数来优化就不是那么难想到的了

code

#include <map>
#include <set>
#include <array>
#include <cmath>
#include <queue>
#include <stack>
#include <tuple>
#include <bitset>
#include <cctype>
#include <cstdio>
#include <random>
#include <string>
#include <vector>
#include <cassert>
#include <cstring>
#include <iomanip>
#include <limits.h>
#include <iostream>
#include <algorithm>
#include <functional>
#include <unordered_map>
using namespace std;

#define N maxn
#define db double
#define il inline
#define fir first
#define sec second
#define eps (1e-8)
#define pb push_back
#define ll long long
#define mkp make_pair
#define eb emplace_back
#define pii pair<int, int>
#define lowbit(a) (a & (-a))
#define SZ(a) ((int)a.size())
#define ull unsigned long long
#define all(a) a.begin(), a.end()
#define split cout << "=========\n";
#define GG { cout << "NO\n"; return; }
#define pll pair<long long, long long>
#define equals(a, b) (fabs((a) - (b)) < eps)

constexpr int ON = 0;
constexpr int CW = -1;
constexpr int CCW = 1;
constexpr int BACK = 2;
constexpr int FRONT = -2;
const db pi = acos(-1.000);
constexpr int maxn = 2e5 + 50;
constexpr int INF = 0x3f3f3f3f;
constexpr ll LINF =  0x3f3f3f3f3f3f3f3f;
constexpr int mod = 998244353; /* 1e9 + 7 */
constexpr int dir[8][2] = {-1, 0, -1, 1, 0, 1, 1, 1, 1, 0, 1, -1, 0, -1, -1, -1};

mt19937_64 rnd(random_device {}());
uniform_int_distribution<ull> dist(0, ULLONG_MAX);//use dist(rnd)

bool BEGIN;
ll qpow(ll x,ll y)
{
    ll ans=1;
    while(y)
    {
        if(y&1) ans=ans*x%mod;
        x=x*x%mod;
        y>>=1;
    }
    return ans;
}

namespace Poly
{
    #define mul(x, y) (1ll * x * y >= mod ? 1ll * x * y % mod : 1ll * x * y)
    #define minus(x, y) (1ll * x - y < 0 ? 1ll * x - y + mod : 1ll * x - y)
    #define plus(x, y) (1ll * x + y >= mod ? 1ll * x + y - mod : 1ll * x + y)//上面其实没用到
    #define ck(x) (x >= mod ? x - mod : x)//取模运算太慢了

    typedef vector<int> poly;
    const int G = 3;//根据具体的模数而定，原根可不一定不一样！！！
    //一般模数的原根为 2 3 5 7 10 6
    const int inv_G = qpow(G, mod - 2),tt = 22;
    int deer[2][tt][(1 << tt)];
    vector<int>RR(1 << (tt + 1), 0),inv(1 << tt, 0);

    void init(const int t) {//预处理出来NTT里需要的w和wn，砍掉了一个log的时间
        assert(t < tt);//一定要注意！！
        for(int p = 1; p <= t; ++ p) {
            int buf1 = qpow(G, (mod - 1) / (1 << p));
            int buf0 = qpow(inv_G, (mod - 1) / (1 << p));
            deer[0][p][0] = deer[1][p][0] = 1;
            for(int i = 1; i < (1 << p); ++ i) {
                deer[0][p][i] = 1ll * deer[0][p][i - 1] * buf0 % mod;//逆
                deer[1][p][i] = 1ll * deer[1][p][i - 1] * buf1 % mod;
            }
        }
        inv[1] = 1;
        for(int i = 2; i <= (1 << t); ++ i)
            inv[i] = 1ll * inv[mod % i] * (mod - mod / i) % mod;
    }

    int NTT_init(int n) {//快速数论变换预处理
        int limit = 1, L = 0;
        while(limit <= n) limit <<= 1, L ++ ;
        assert(L < tt);
        assert(limit < 1 << (tt + 1));
        for(int i = 0; i < limit; ++ i)
            RR[i] = (RR[i >> 1] >> 1) | ((i & 1) << (L - 1));
        return limit;
    }

    void NTT(poly &A, bool type, int limit) {//快速数论变换
        A.resize(limit);
        for(int i = 0; i < limit; ++ i)
            if(i < RR[i])
                swap(A[i], A[RR[i]]);
        for(int mid = 2, j = 1; mid <= limit; mid <<= 1, ++ j) {
            int len = mid >> 1;
            for(int pos = 0; pos < limit; pos += mid) {
//                auto wn = deer[type][j].begin();
                for(int i = pos, p = 0; i < pos + len; ++ i, ++ p) {
                    int tmp = 1ll * deer[type][j][p] * A[i + len] % mod;
                    A[i + len] = ck(A[i] - tmp + mod);
                    A[i] = ck(A[i] + tmp);
                }
            }
        }
        if(type == 0) {
            for(int i = 0; i < limit; ++ i)
                A[i] = 1ll * A[i] * inv[limit] % mod;
        }
    }

    poly poly_mul(poly A, poly B) {//多项式乘法
        int deg = A.size() + B.size() - 1;
        int limit = NTT_init(deg);
        poly C(limit);
        NTT(A, 1, limit);
        NTT(B, 1, limit);
        for(int i = 0; i < limit; ++ i)
            C[i] = 1ll * A[i] * B[i] % mod;
        NTT(C, 0, limit);
        C.resize(deg);
        return C;
    }

    poly poly_inv(poly &f, int deg) {//多项式求逆 deg<f.szie()
        if(deg == 1)
            return poly(1, qpow(f[0], mod - 2));

        poly A(f.begin(), f.begin() + deg);
        poly B = poly_inv(f, (deg + 1) >> 1);
        int limit = NTT_init(deg << 1);
        NTT(A, 1, limit), NTT(B, 1, limit);
        for(int i = 0; i < limit; ++ i)
            A[i] = B[i] * (2 - 1ll * A[i] * B[i] % mod + mod) % mod;
        NTT(A, 0, limit);
        A.resize(deg);
        return A;
    }

    poly poly_dev(poly f) {//多项式求导
        int n = f.size();
        for(int i = 1; i < n; ++ i) f[i - 1] = 1ll * f[i] * i % mod;
        if(n > 1)f.resize(n - 1);
        else f[0] = 0;
        return f.resize(n - 1), f;//求导整体左移，第0项不要
    }

    poly poly_idev(poly f) {//多项式求积分
        int n = f.size();
        for(int i = n - 1; i ; -- i) f[i] = 1ll * f[i - 1] * inv[i] % mod;
        return f[0] = 0, f;//积分整体右移，第0项默认为0
    }

    poly poly_ln(poly f, int deg) {//多项式求对数，第一项为1
        poly A = poly_idev(poly_mul(poly_dev(f), poly_inv(f, deg)));
        return A.resize(deg), A;
    }

    poly poly_exp(poly &f, int deg) {//多项式求指数，第一项为0
        if(deg == 1)
            return poly(1, 1);

        poly B = poly_exp(f, (deg + 1) >> 1);
        B.resize(deg);
        poly lnB = poly_ln(B, deg);
        for(int i = 0; i < deg; ++ i)
            lnB[i] = ck(f[i] - lnB[i] + mod);

        int limit = NTT_init(deg << 1);//n -> n^2
        NTT(B, 1, limit), NTT(lnB, 1, limit);
        for(int i = 0; i < limit; ++ i)
            B[i] = 1ll * B[i] * (1 + lnB[i]) % mod;
        NTT(B, 0, limit);
        B.resize(deg);
        return B;
    }

    poly poly_pow(poly f, int k) {//多项式快速幂,第一项得是1
        f = poly_ln(f, f.size());
        for(auto &x : f) x = 1ll * x * k % mod;
        return poly_exp(f, f.size());
    }
     poly poly_ksm(poly f, int k,int m) {//多项式快速幂,适用于初始只有几项，同时所有项都需要的情况,会比上面那个快一点
        poly res(1, 1);
        while(k){
            if(k & 1)
                {
                    res = poly_mul(res, f);
                    res.resize(m,0);
                }
            f = poly_mul(f, f);
            f.resize(m,0);
            k >>= 1;
        }
        return res;
    }



}

using Poly::poly;
using Poly::poly_pow;
using Poly::poly_ksm;
using Poly::poly_mul;
using Poly::poly_inv;
vector<ll> vt;
poly poly_f(int l,int r,int num)
    {
        if(l>r)
        {
            return poly(1,1);
        }
        if(l==r)
        {
            poly ans(vt[l]+1,0);
            ans[0]=1;
            ans[vt[l]]=(mod-1)%mod;
            return ans;
        }
        int mid=(l+r)>>1;
        poly f = poly_f(l,mid,num);
        poly g = poly_f(mid+1,r,num);
        f = poly_mul(f,g);
        // f.resize(num+1);
        return f;
    }

ll p[N],pp[N],C[N];
void init(ll n)
{
    p[0]=1;
    for(ll i=1;i<=n;++i) p[i]=p[i-1]*i%mod;
    pp[1]=1;
    for(ll i=2;i<=n;++i)
    {
        pp[i]=(mod-mod/i)*pp[mod%i]%mod;
    }
    // for(ll i=1;i<=n;++i) C[i]=C[i-1]*qpow(i,mod-2)%mod*(k-i+1)%mod;
}
void init2(ll n,ll k)
{
    C[0]=1;
    for(ll i=1;i<=n;++i) C[i]=C[i-1]*pp[i]%mod*(k-i+1)%mod;
    for(int i=1;i<=n;++i) C[i]=C[i]*p[i]%mod;
}

ll n,m,k;
map<ll,ll> num1;
void solve()
{
    num1.clear();
    cin>>n>>m>>k;
    init2(max(n,m),k);
    // for(int i=0;i<=10;++i) cout<<C[i]<<' ';
    //     cout<<endl;
    for(int i=1;i<=n;++i)
    {
        ll a;cin>>a;
        num1[a]++;
    }

    vt.clear();
    for(auto i:num1) vt.push_back(i.second);
    Poly::init(20);
    poly f=poly_f(0,vt.size()-1,n);


    ll ans=0;
    for(ll i=0;i<=min(n,k);++i)
    {
        // cout<<f[i]<<' '<<C[i]<<" "<<qpow(n-i,k-i)<<endl;
        ans=((ans+f[i]*C[i]%mod*qpow(n-i,k-i)%mod)%mod+mod)%mod;
    }
    cout<<ans*qpow(qpow(n,k),mod-2)%mod<<endl;
}

bool END;
signed main() {
    // cout << fixed << setprecision(10);
    ios::sync_with_stdio(false); cin.tie(nullptr);
    init(200000);
    int T; cin >> T; while (T--)
    solve();
    // cout << ((&END - & BEGIN) >> 21) << '\n';
    return 0;
}

18青岛G

大意：
给定一个由0,1，2组成的长度为n的数列，每次操作可以选择一个不含2的连续区间，将其中所有元素置为0.求恰好操作m次并且所有2都被消除的方案数,$n\le 100,m\le 10^9$
思路：

不妨记一个不含1的连续区间是一个合法区间，那么仔细观察一下，对于一个长度为$len$的合法区间，其操作的合法区间数是$\left(\genfrac{}{}{0ex}{}{len}{2}\right)+len=\frac{len(len+1)}{2}$，所以其答案就是$(\frac{len(len+1)}{2}{)}^m$,是可以$O(1)$计算得到的。这启示我们使用容斥来处理该问题 ，令$A_i$表示从左往右数第$i$个2没有消除的方案数，要求$ans=|S-A_1\bigcup A_2...\bigcup A_x|={\sum }_{i=0}^x(-1{)}^i|A_{j_1}\bigcap A_{j_2}...\bigcap A_{j_i}|$
但是显然数据范围并不允许我们直接暴力枚举，所以考虑dp优化（很常见的思路）。

显然$|A_{j_1}\bigcap A_{j_2}...\bigcap A_{j_i}|$前面的系数由$i$，也就是不消除的2的个数决定,如果是奇数，系数就是-1，并且此时其他的2是不能消除的。所以不妨记$d{p}_{i,k,0/1}$表示前i个位置，合法操作区间数为k，总共有奇数/偶数个不能消除的2

为了方便计算，令$a_0=a_{n+1}=1$,这不会影响结果

那么最终我们的答案就是${\sum }_{j=1}^{n(n+1)/2}(d{p}_{n+1,j,0}-d{p}_{n+1,j,1})j^m$,也就是枚举合法操作方案数为j的情况数

考虑转移，对于当前位置i，我们枚举上一个不消除的位置j，$pr{e}_i\le j<i,a_j=2$,其中$pr{e}_i$表示i之前的第一个1的位置，因为1是一定不能消除的。此时$d{p}_{i,k+\frac{(i-j)(i-j+1)}{2},0/1}={\sum }_{j=pr{e}_i}^{i-1}d{p}_{j,k,1/0}$

最后总体复杂度$O(n^4+log(m))$

最后要特判一下1的个数，如果是奇数的话，显然改变了实际2的个数，所以答案要取负

ll dp[N][N*N][3];
ll n,m;
ll mas[N];
ll lim[N];//前i个点的最多的合法区间数
void solve()
{
    cin>>n>>m;
    for(int i=1;i<=n;++i) cin>>mas[i];
    mas[0]=mas[n+1]=1;
    dp[0][0][1]=1;dp[0][0][0]=0;lim[0]=0;
    ll pre=0;//上一个1的位置
    for(int i=1;i<=n+1;++i)
    {
        if(mas[i]==0) continue;
        lim[i]=lim[pre]+(i-pre)*(i-pre-1)/2;
        for (int j=0;j<=lim[i];++j) dp[i][j][0]=dp[i][j][1]=0;
        for(int j=pre;j<i;++j)
        {
            if(mas[j]==0) continue;
            ll add=(i-j)*(i-j-1)/2;
            for(int k=0;k<=lim[j];++k)
            {
                dp[i][k+add][0]=(dp[i][k+add][0]+dp[j][k][1])%mod;
                dp[i][k+add][1]=(dp[i][k+add][1]+dp[j][k][0])%mod;
            }
        }
        if(mas[i]==1) pre=i;
    }
    ll ans=0;
    for(int i=0;i<=lim[n+1];++i)
    {
        ll det=((dp[n+1][i][0]-dp[n+1][i][1])%mod+mod)%mod;
        ans=(ans+det*ksm(i,m)%mod)%mod;
    }

    ll cn=0;
    for(int i=1;i<=n;++i) if(mas[i]==1) cn++;
    if(cn%2) ans=(mod-ans)%mod;

    cout<<ans<<endl;
}

未完待续

​