Cod sursă (job #786354)

#include <fstream>

using namespace std;

ifstream f("dejavu.in");
ofstream g("dejavu.out");

using ll = long long;

const int MAXN = 5000;
const int MAXFACT = 5000;
const long long MOD = 1000000007LL;

void euclid(ll a, ll b, ll &x, ll &y)
{
    if (!b) {
        x = 1;
        y = 0;
        return;
    }
    euclid(b, a%b, y, x);
    y -= a/b * x;
}
ll invMod(ll nr)
{
    ll x, y;
    euclid(nr, MOD, x, y);
    return x % MOD;
}

struct Fenwick
{
    int v[MAXN+2];
    int n;
    void init(int n, int val)
    {
        this->n = n;
        for (int i = 1; i <= n; i++) {
            v[i] = val * (i & (-i));
        }
    }
    void update(int poz, int val)
    {
        while (poz <= n) {
            v[poz] += val;
            poz += poz & (-poz);
        }
    }
    int query(int poz)
    {
        int ans = 0;
        while (poz) {
            ans += v[poz];
            poz &= poz-1;
        }
        return ans;
    }
};

int n, q, v[MAXN+2], fr[MAXN+2];
Fenwick aib1, aib2;
ll fact[MAXFACT+2], invfact[MAXFACT+2];
ll dp[MAXN+2][MAXN+2];
/// dp[i][j] = numarul de moduri de a crea un sir deja vu
///            daca avem i numere ce pot aparea
///            si j numere ce pot aparea de exact 2 ori
///            (i >= j)
///            (lungimea sirului: i+j-n)

void calcFact()
{
    fact[0] = 1;
    for (int i = 1; i <= MAXFACT; i++) {
        fact[i] = (fact[i-1] * i) % MOD;
    }
    invfact[MAXFACT] = invMod(fact[MAXFACT]);
    for (int i = MAXFACT-1; i >= 0; i--) {
        invfact[i] = (invfact[i+1] * (i+1)) % MOD;
    }
}

void calcDp()
{
    for (int i = (n+1)/2; i <= n; i++) {
        dp[i][n-i] = 1;
        for (int j = n-i+1; j <= i; j++) {
            dp[i][j] = ((j > 0 ? dp[i][j-1] * j : 0) + (i > 0 ? dp[i-1][j] * (i-j) : 0)) % MOD;
        }
    }
}

void solveQuery()
{
    ll ans = 0;
    for (int i = 1; i <= n; i++) {
        f >> v[i];
        fr[i] = 2;
    }
    int d1 = n, d2 = n;
    aib1.init(n, 0);
    aib2.init(n, 1);
    for (int i = 1; i <= n; i++) {
        ans = (ans + dp[d1-1][d2] * aib1.query(v[i]-1) + dp[d1][d2-1] * aib2.query(v[i]-1)) % MOD;
        fr[v[i]]--;
        if (fr[v[i]] == 1) {
            d2--;
            aib2.update(v[i], -1);
            aib1.update(v[i], +1);
        }
        else {
            d1--;
            aib1.update(v[i], -1);
        }
    }
    g << (ans % MOD + MOD) % MOD << '\n';
}

int main()
{
    f >> n >> q;
    calcFact();
    calcDp();
    for (int i = 1; i <= q; i++) {
        solveQuery();
    }
    return 0;
}