RT

#include<bits/stdc++.h>
#define int long long
using namespace std;

int n,r[3001],q[3001],v[3001],f[3001],g[3001],fac[3001],inv[3001],Mod=1e9+7,now;
vector<int>Map[3001];

int dfs(int u,int lenth)
{
    if(v[u])
        return 0;
    v[u]=1;
    for(int i:Map[u])
        lenth=max(lenth,dfs(i,lenth+1));
    return lenth;
}

int Pow(int a,int b)
{
    int ans=1;
    while(b)
    {
        if(b&1)
            (ans*=a)%=Mod;
        (a*=a)%=Mod;
        b>>=1;
    }
    return ans;
}

int C(int a,int b)
{return fac[a]*inv[b]%Mod*inv[a-b]%Mod;}

signed main()
{
    cin>>n;
    fac[0]=1;
    for(int i=1;i<=n<<1;i++)
        fac[i]=fac[i-1]*i%Mod;
    inv[n<<1]=Pow(fac[n<<1],Mod-2);
    for(int i=n<<1;i>=1;i--)
        inv[i-1]=inv[i]*i%Mod;
    for(int i=1;i<=n;i++)
        cin>>r[i];
    for(int i=1;i<=n;i++)
        cin>>q[i];
    for(int i=1;i<=n;i++)
        Map[r[i]].push_back(q[i]),
        Map[q[i]].push_back(r[i]);
    f[0]=1;
    for(int i=1;i<=n;i++)
        if(!v[i])
        {
            int lenth=dfs(i,1);
            now+=lenth;
            if(lenth==1)
            {
                for(int i=n;i;i--)
                    f[i]+=f[i-1];
                continue;
            }
            for(int i=now;i;i--)
            {
                int g=0;
                for(int j=0;j<=lenth&&j<=i;j++)
                    (g+=f[i-j]*(C(2*lenth-j,j)+C(2*lenth-j-1,j-1)))%=Mod;
                f[i]=g;
            }
        }
    int ans=0;
    for(int i=0;i<=n;i++)
        (ans+=(f[i]*fac[n-i]%Mod)*(i&1?-1:1))%=Mod;
    cout<<(ans+Mod)%Mod<<"\n";
}

他甚至说我样例挂了