[to record the solution for my own record for easy understanding]

Problem statement: http://usaco.org/index.php?page=viewproblem2&cpid=1068

dp[i, j, k] where

• i-th barn, where the corresponding cows that fits is match[i]

• out of match[i] cows, j of them are assigned to barn i+1 ... N

• and flag = 0 means no mis-match cows / all cows are in barns, 1 = allow some cows to left behind

so, each time we add a new barn (barn and cow are sorted low to high), we have an additional match[i] ... match[i-1] cow to assign. Let's take j from dp[i-1] and k from those now cows

dp[i, j+k, f] += dp[i, j, f] * (probability of picking k out of match[i] ... match[i-1], which is a simple c(n,m) )

note that barn i is the largest, so it must be occupied unless all cows are assigned (i.e. f = 0). So, in above case

• if we pick j from dp[i-1], and one of them landed at barn i, it is good, but the number of cow after barn i+1 is actually j+k-1

• if it does not land at barn i, for the k cow we picked, one of them must be used for barn i. And again then, the correct index should be j+k-1.

So, after the above dp, before we finish barn i, we need to do j+1 -> j index shift for f = 1. Note that we can pick any of (j+1) cow to fit at barn i, so times (j+1).

if f=0, all the cows are assigned so that we can leave barn i open (not used). So, we need the original count (not assign to barn i) and add the j+1 -> j shift

code as below

/*
ID: USACO_template
LANG: C++
PROG: http://usaco.org/index.php?page=viewproblem2&cpid=1068
*/
#include <iostream>  //cin , cout
#include <fstream>   //fin, fout
#include <stdio.h>   // scanf , pringf
#include <cstdio>
#include <algorithm> // sort , stuff
#include <stack>     // stacks
#include <queue>     // queues
#include <map>
#include <string>

using namespace std;

typedef pair<int, int>          pii;
typedef vector<int>             vi;     /// adjlist without weight
typedef vector<pii>             vii;    /// adjlist with weight
typedef vector<pair<int,pii>>   vpip;   /// edge with weight
typedef long long               ll;

#define mp  make_pair
#define fst first
#define snd second
#define pb  push_back
#define sz(x) (int)(x).size()
#define trav(u, adj_v) for (auto& u: adj_v)

const int MOD = 1e9+7;  // 998244353;
const int MX  = 2e5+5;   //
const ll  INF = 1e18;    //

#define MAXV 3007
#define MAXE 100007

bool debug = false, submit = true;

int N;
int s[MAXV], t[MAXV];

int match[MAXV];
ll dp[MAXV][MAXV][2];
ll c[MAXV][MAXV];

int main() {
    debug = true;
    ios_base::sync_with_stdio(0); cin.tie(0);

    int i, j, k;

    cin >> N ;
    for(i=1;i<=N;i++) cin >> s[i];
    for(i=1;i<=N;i++) cin >> t[i];
    sort(s+1, s+1+N); sort(t+1, t+1+N);

    /// match[j] means the j-th barn can be house upto match[j]-th cow (after the small to large sort is completed
    i=1; j = 1;
    while(i<=N && j<=N) {
        match[j] = match[j-1];
        while(i<=N && t[j] >= s[i]) {
            match[j] = i; i++;
        }
        j++;
    }
    while(i>=N && j<=N) { match[j] = match[j-1]; j++;}

    /// pre caculate probability of pick m from n, C(n, m), which = c(n-1, m) + c(n-1, m-1)
    c[0][0] = 1LL;
    for(i=1; i<=N; i++) {
        for(j=0; j<=N; j++) {
            c[i][j] = c[i-1][j] + ( (j-1>=0) ? c[i-1][j-1] : 0 );
            c[i][j] %= MOD;
        }
    }

    /// dp(i, j, f) is use upto i-th barn, and from 1 ... match[i] cow, at least j cows are assigned to bar i+1  ... N,
    /// f=0 is all cows are used, f=1 is not all cows are used
    dp[0][0][0] = 1LL;
    for(i=1; i<=N; i++) {
        for(j=0; j<=match[i-1]; j++) {
            for(int f=0; f<=1; f++) {
                int klen = match[i] - match[i-1];
                for(k=0; k <= klen; k++) {
                    int nf = 1;
                    if(f==0 && k == klen) nf = 0; /// all cows used
                    dp[i][j+k][nf] += (dp[i-1][j][f] * c[klen][k])%MOD;   
                    dp[i][j+k][nf] %= MOD;
                }
            }
        }
        /// then we can not let j cow sit at barn i, but barn i must be used unless it is all assigned.
        for(j=0; j<match[i];j++) {
            dp[i][j][0] += ( dp[i][j+1][0] * (ll)(j+1) )%MOD; dp[i][j][0] %=MOD;
            dp[i][j][1] = ( dp[i][j+1][1] * (ll)(j+1) )%MOD; dp[i][j][1] %=MOD;     ///
        }
        dp[i][match[i]][1] = 0LL;
    }

    cout << (dp[N][0][0]%MOD + dp[N][0][1]%MOD) % MOD << endl;

}