题目
建立一个二分图，左右各n个点，在左边的第x个点和右边的第y个点之间连一条权值为\(a_{x,y}\)的边。根据“积和式”的定义，我们是要在矩阵中选择n个位置，满足任意两个位置不共行、不共列，并计算所有选择方案的对应位置的值的积之和。显然，每一种选位置的方案对应这个二分图的一个完美匹配，我们现在要计算这个二分图每一个完美匹配的匹配边的边权乘积之和。左边和右边都有1e5个点，一共1e10条边，这是很多的；但是边权不是1的边只有50条，这就很好，考虑怎么略去边权为1的边的计算。对于每条\(a_{x,y}>1\)的边，令其中\(a_{x,y}-1\)条为特殊边，剩下1条为普通边。枚举所有选了特殊边的匹配点对，方案数\(\prod (a_{x,y}-1)\)；剩下的每队匹配点就只有1种选择了，也就是选普通边，方案数\((n-左边的选择了特殊边的匹配点个数)\)!,因为左边剩下的点都可以在右边的未匹配点中任选一个，且每对只有1种方案。考虑仅由\(a_{x,y}>1\)的边构成的图，注意到每个连通块的方案数是不互相影响的，对于一个连通块，可以计算其中有\(i\)对点选择用特殊边匹配的方案数为\(f_i\)，然后把每个连通块的f数组手动暴力卷积起来，最后乘上剩下节点数量的阶乘。
接下来对每个连通块分别计算。注意到连通块内点数是不超过边数+1的，所以一个连通块内最多51个点，记连通块内总点数为k。那么，左边和右边一定有一边的点数是\(\leq \frac k2\)的，假设左边是点少的那边。那么可以令\(dp_{i,msk}\)表示当前遍历到右边的第i个点，左边已经被匹配掉的点的集合是msk，的所有方案的贡献之和。转移就枚举右边第i个点和左边哪个点匹配，或者干脆不匹配。最后\(dp_{sz(l),msk}\)中msk二进制里1的个数就是匹配掉的点对数。由于一共50条边，所以i那维和转移的枚举合起来最多50次，这部分的总复杂度也就是\(50\cdot 2^{25}\)，但是25还是比较大的(虽然这样好像能过？),我们可以在\(sz(l)<20\)的时候这么干，对其他情况采用另一种办法。
在\(sz(l)>19\)的时候，总点数也就至少是38，由于边数最多50，所以如果我们求出这个图的一个生成树，非树边最多13条。接下来就可以\(2^{13}\)次枚举非树边的选法，然后在生成树上进行一个树形dp。具体来说，\(dp_{i,j,0/1}\)表示以i为根的子树，子树内有j对点选了特殊边匹配，以及i是否被特殊边匹配占用。转移比较简单，需要在i的儿子中再跑一个背包。这部分复杂度是\(2^{13}\cdot 50^2\)，50是平方不是3次方的原因是转移时另外跑一个背包的时候，枚举了i的一个儿子子树中选择的匹配点对数和其他前缀儿子中匹配点对数总和，但是发现树中每2个节点只会在他们的LCA处产生一次枚举。
这样把这两个做法结合一下就可以轻松通过了~
下面的代码有小错，在CF上是WA的，但是模拟赛第三方数据可以通过，大家先将就着看吧（
点击查看代码
#include <bits/stdc++.h>
#define rep(i,n) for(int i=0;i<n;++i)
#define repn(i,n) for(int i=1;i<=n;++i)
#define LL long long
#define pii pair <LL,LL>
#define pb push_back
#define fi first
#define se second
#define mpr make_pair
using namespace std;
const LL MOD=1e9+7;
LL n,k,fac[100010],val[100],dp[60][530000],mark[100],dp2[60][60][2],dp3[60][60][2],sz[60];
map <pii,LL> mp,te;
map <LL,LL> vs;
vector <pii> g[200010],gg[100];
vector <LL> res,l,r,all;
bool vis[200010];
vector <pair <pii,LL> > e,nte;
void dfs(LL pos)
{
if(pos<n) l.pb(pos);else r.pb(pos);
vis[pos]=true;
rep(i,g[pos].size())
{
e.pb(mpr(mpr(min(pos,g[pos][i].fi),max(pos,g[pos][i].fi)),g[pos][i].se));
if(!vis[g[pos][i].fi]) dfs(g[pos][i].fi),te[mpr(min(pos,g[pos][i].fi),max(pos,g[pos][i].fi))]=1;
}
}
vector <LL> convolution(vector <LL> a,vector <LL> b)
{
vector <LL> c(a.size()+b.size()-1,0);
rep(i,a.size()) rep(j,b.size()) (c[i+j]+=a[i]*b[j])%=MOD;
return c;
}
LL getid(LL x){return lower_bound(all.begin(),all.end(),x)-all.begin();}
void dfs2(LL pos,LL par)
{
vector <pii> son;
sz[pos]=1;
rep(i,gg[pos].size()) if(gg[pos][i].fi!=par)
{
dfs2(gg[pos][i].fi,pos);
son.pb(gg[pos][i]);sz[pos]+=sz[gg[pos][i].fi];
}
rep(i,son.size()+3) rep(j,all.size()+3) rep(k,2) dp3[i][j][k]=0;
dp3[0][0][0]=1;
rep(i,son.size()) rep(j,all.size()+1) rep(k,2) if(dp3[i][j][k]>0)
{
rep(jj,sz[son[i].fi]) rep(kk,2) if(dp2[son[i].fi][jj][kk]>0)
{
(dp3[i+1][j+jj][k]+=dp3[i][j][k]*dp2[son[i].fi][jj][kk])%=MOD;
if(k==0&&kk==0&&mark[pos]==0&&mark[son[i].fi]==0) (dp3[i+1][j+jj+1][1]+=dp3[i][j][k]*dp2[son[i].fi][jj][kk]%MOD*son[i].se)%=MOD;
}
}
rep(j,all.size()+1) rep(k,2) dp2[pos][j][k]=dp3[son.size()][j][k];
}
int main()
{
//freopen("c.in","r",stdin);
//freopen("c.out","w",stdout);
fac[0]=1;repn(i,100005) fac[i]=fac[i-1]*i%MOD;
cin>>n>>k;
LL x,y,z;
rep(i,k)
{
cin>>x>>y>>z;--x;--y;
mp[mpr(x,y)]=z;vs[x]=vs[y+n]=1;
}
for(auto it:mp)
{
g[it.fi.fi].pb(mpr(it.fi.se+n,it.se-1));
g[it.fi.se+n].pb(mpr(it.fi.fi,it.se-1));
}
res.pb(1);
for(auto it:vs) if(!vis[it.fi])
{
rep(i,99) val[i]=0;
e.clear();l.clear();r.clear();te.clear();
dfs(it.fi);
sort(e.begin(),e.end());e.erase(unique(e.begin(),e.end()),e.end());
if(min(l.size(),r.size())<=19)
{
if(r.size()<l.size()) swap(l,r);
rep(i,99) gg[i].clear();
rep(i,l.size())
{
rep(j,e.size()) if(e[j].fi.fi==l[i]||e[j].fi.se==l[i])
{
int v=e[j].fi.fi+e[j].fi.se-l[i],id;
rep(j,r.size()) if(r[j]==v) id=j;
gg[id].pb(mpr(i,e[j].se));
}
}
rep(j,59) rep(k,1<<l.size()) dp[j][k]=0;
dp[0][0]=1;
rep(i,r.size()) rep(k,1<<l.size()) if(dp[i][k]>0)
{
(dp[i+1][k]+=dp[i][k])%=MOD;
rep(p,gg[i].size()) if((k&(1<<gg[i][p].fi))==0) (dp[i+1][k|(1<<gg[i][p].fi)]+=dp[i][k]*gg[i][p].se)%=MOD;
}
rep(k,1<<l.size()) (val[__builtin_popcount(k)]+=dp[r.size()][k])%=MOD;
}
else
{
rep(i,99) gg[i].clear();
nte.clear();
all.clear();rep(i,l.size()) all.pb(l[i]);rep(i,r.size()) all.pb(r[i]);
sort(all.begin(),all.end());
rep(i,e.size())
{
if(te.find(e[i].fi)!=te.end())
{
LL u=getid(e[i].fi.fi),v=getid(e[i].fi.se);
gg[u].pb(mpr(v,e[i].se));gg[v].pb(mpr(u,e[i].se));
}
else nte.pb(mpr(mpr(getid(e[i].fi.fi),getid(e[i].fi.se)),e[i].se));
}
rep(i,all.size()) mark[i]=0;
rep(i,1<<nte.size())
{
LL addn=__builtin_popcount(i),mulv=1,bad=0;
rep(j,nte.size()) if((i&(1<<j))>0)
{
(mulv*=nte[j].se)%=MOD;
if(mark[nte[j].fi.fi]||mark[nte[j].fi.se]) bad=1;
mark[nte[j].fi.fi]=mark[nte[j].fi.se]=1;
}
if(bad) continue;
rep(j,all.size()) rep(k,all.size()) rep(p,2) dp2[j][k][p]=0;
dfs2(0,-1);
rep(j,all.size()+1) rep(k,2) if(dp2[0][j][k]>0) (val[j+addn]+=dp2[0][j][k]*mulv)%=MOD;
}
}
int lst=-1;rep(i,99) if(val[i]>0) lst=i;
if(lst==-1) continue;
vector <LL> tmp;rep(i,lst+1) tmp.pb(val[i]);
res=convolution(res,tmp);
}
LL ans=0;
rep(i,res.size()) (ans+=res[i]*fac[n-i])%=MOD;
cout<<ans<<endl;
return 0;
}
![P1941 [NOIP2014 提高组] 飞扬的小鸟 题解](https://www.kunjuke.com/file/css/thumb/1.jpg/280-180.jpg)
![P1880 [NOI1995] 石子合并 题解](https://www.kunjuke.com/file/css/thumb/7.jpg/280-180.jpg)
