Cod sursă (job #786108)

// Author: Tanasescu Andrei-Rares
#include <iostream>
#include <fstream>
#include <algorithm>
#include <cmath>
#include <map>
#include <unordered_map>
#include <set>
#include <unordered_set>
#include <queue>
#include <stack>
#include <deque>
#include <iomanip>
#include <vector>
#include <cassert>

#pragma GCC optimize("O3")

#define fi first
#define se second
#define pb push_back
#define pf push_front

using namespace std;

ifstream fin ("dejavu.in");
ofstream fout ("dejavu.out");

typedef long long ll;
typedef pair<int, int> pii;
typedef pair<ll, ll> pll;

const ll Nmax=5e3+5, inf=1e9+5, MOD=1e9+7;

int n, q;

// dp[i][j] = modalitatea de a continua un prefix ce are i litere singuratice
//            si j perechi de litere
// lungimea prefixului = i+j*2
ll dp[Nmax][Nmax/2];
int fr[Nmax];

struct AIB{
    int v[Nmax];
    inline void update(int pos, int val){
        while (pos<Nmax){
            v[pos]+=val;
            pos+=pos&-pos;
        }
    }
    inline int query(int pos){
        int sum=0;
        while (pos>0){
            sum+=v[pos];
            pos-=pos&-pos;
        }
        return sum;
    }
    inline void clr(){
        for (int i=1; i<Nmax; i++)
            v[i]=0;
    }

}aib1, aib2;

int main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);

    fin>>n>>q;

    // initializare dp
    for (int i=0; i<=n; i++){
        int r=n-i;
        if (r%2==0)
            dp[i][r/2]=1;
    }
    // build dp
    for (int l=n-1; l>=0; l--)
        for (int i=0; i<=l; i++){
            int r=l-i;
            if (r%2==0)
                dp[i][r/2]=dp[i+1][r/2]*(n-i-r/2);
            if (i!=0)
                dp[i][r/2]=(dp[i][r/2]+dp[i-1][r/2+1]*i)%MOD;
        }

    // answer queries
    while (q--){
        int nr, cnt1=0, cnt2=0;
        ll sol=0;
        for (int i=0; i<n; i++){
            fin>>nr;
            int crt1=aib1.query(nr-1);
            int crt0=nr-1-crt1-aib2.query(nr-1);

            sol=(sol+crt0*dp[cnt1+1][cnt2])%MOD;
            if (cnt1!=0)
                sol=(sol+crt1*dp[cnt1-1][cnt2+1])%MOD;

            fr[nr]++;
            if (fr[nr]==1){
                cnt1++;
                aib1.update(nr, 1);
            }
            else{
                cnt2++;
                cnt1--;
                aib1.update(nr, -1);
                aib2.update(nr, 1);
            }
        }
        fout<<sol<<'\n';

        for (int i=1; i<=n; i++)
            fr[i]=0;
        aib1.clr();
        aib2.clr();

    }

    return 0;
}