Solution for https://codeforces.com/contest/1151/problem/F

sonya_and_informatics.cc

#include <bits/stdc++.h>
using namespace std;
using ll = long long;

constexpr ll mod = 1000000007;

ll fast_pow_mod(ll a, ll b, ll m) {
    a %= m;
    if(a == 0ll) return a;

    ll ret = 1;
    while(b) {
        if(b&1ll == 1ll) ret = (ret * a) % m;
        a = (a * a) % m;
        b >>= 1;
    }

    return ret;
}

ll mod_inv_prime(ll a, ll m) {
    // m is assumed to be prime here.
    return fast_pow_mod(a, m - 2, m);  // by Fermat's little theorem
}

struct RatP {
    ll nd, p;

    RatP(ll nom, ll denom, ll _p) : p(_p) {
        nd = (nom * mod_inv_prime(denom, p)) % p;
    }

    RatP(ll _nd, ll _p) : nd(_nd), p(_p) {
        nd = nd % p;
    }

    RatP() {
        nd = 0ll;
        p = mod;
    }
};

string show(const RatP &r) {
    stringstream ss;
    ss << "(" << r.nd << "," << r.p << ")";
    return ss.str();
}

RatP unit() {
    return RatP(1ll, mod);
}

RatP operator*(const RatP &lhs, const RatP &rhs) {
    return RatP((lhs.nd * rhs.nd) % lhs.p, lhs.p);
}

RatP operator+(const RatP &lhs, const RatP &rhs) {
    return RatP((lhs.nd + rhs.nd) % lhs.p, lhs.p);
}


template<typename T>
ostream &operator<<(ostream &out, const vector<vector<T>> &mat) {
    for(const auto &row : mat) {
        for(const auto &el : row) {
            out << "\t\t" << show(el);  
        }
        out << endl;
    }

    return out;
}

template<typename T>
vector<vector<T>> fast_pow_mat(const vector<vector<T>> &mat, ll b) {
    ll n = mat.size();
    assert(n == mat[0].size());
    vector<vector<T>> ret(n, vector<T>(n));
    for (int i = 0; i < n; ++i) {
        ret[i][i] = unit();
    }
    vector<vector<T>> curr = mat;

    while(b) {
        if(b & 1ll == 1ll) {
            ret = ret * curr;
        }
        curr = curr * curr;
        b >>= 1;
    }

    return ret;
}

template<typename T>
vector<vector<T>> operator*(const vector<vector<T>> &lhs, const vector<vector<T>> &rhs) {
    assert(lhs[0].size() == rhs.size());
    ll n = lhs.size(), m = rhs[0].size(), o = rhs.size();
    vector<vector<T>> ret(n, vector<T>(m));

    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < m; ++j) {
            for (int k = 0; k < o; ++k) {
                ret[i][j] = ret[i][j] + lhs[i][k] * rhs[k][j];
            }
        }
    }

    return ret;
}

template<typename T>
vector<vector<T>> operator+(const vector<vector<T>> &lhs, const vector<vector<T>> &rhs) {
    assert(lhs.size() == rhs.size() && lhs[0].size() == rhs[0].size());
    ll n = lhs.size(), m = lhs[0].size();
    
    vector<vector<T>> ret(n, vector<T>(m));
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < m; ++j) {
            ret[i][j] = lhs[i][j] + rhs[i][j];
        }
    }

    return ret;
}

template<typename T>
vector<T> operator*(const vector<vector<T>> &mat, const vector<T> &vec) {
    assert(mat[0].size() == vec.size());
    ll n = mat.size(), m = vec.size();
    vector<T> ret(n);

    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < m; ++j) {
            ret[i] = ret[i] + mat[i][j] * vec[j];
        }
    }

    return ret;
}

ll fac(ll a, ll p) {
    assert(a >= 0);
    ll ret = 1ll;
    while(a > 0) {
        ret = (ret * a) % p;
        a--;
    }

    return ret;
}

ll choose(ll a, ll b, ll p) {
    if(b > a) return 0ll;
    ll ret = fac(a, p);
    ret = (ret * mod_inv_prime(fac(b, p), p)) % p;
    ret = (ret * mod_inv_prime(fac(a - b, p), p)) % p;
    return ret;
}

int main() {
    ll n, k, one = 0ll;
    cin >> n >> k;
    ll n2 = choose(n, 2, mod);
    vector<ll> vec(n);
    for(int i = 0; i < n; i++) {
        cin >> vec[i];
        if(vec[i] == 1ll) one++;
    }
    if(one == 0ll || one == n) {
        cout << 1 << endl;
        return 0;
    }

    vector<vector<RatP>> mat(one + 1, vector<RatP>(one + 1, RatP(0ll, mod)));

    ll mi = max(0ll, 2 * one - n);
    for (int i = mi; i < one + 1; ++i) {
        ll right_one = i, left_one = one - i;
        ll right_zero = one - right_one, left_zero = n - one - left_one;
        mat[i][i] = RatP(
            choose(one, 2, mod) +
            choose(n - one, 2, mod) +
            choose(left_zero, 1, mod) * choose(right_zero, 1, mod) +
            choose(left_one, 1, mod) * choose(right_one, 1, mod),
            n2, mod);
        if(i != mi) {
            mat[i - 1][i] = RatP(choose(right_one, 1, mod) * choose(left_zero, 1, mod), n2, mod);
        }
        if(i != one) {
            mat[i + 1][i] = RatP(choose(left_one, 1, mod) * choose(right_zero, 1, mod), n2, mod);
        }
    }

    ll right = 0ll;
    for(int i = 0; i < one; i++) {
        if(vec[n - 1 - i] == 1ll) right++;
    }

    vector<RatP> state(one + 1, RatP(0ll, mod));
    state[right] = RatP(1ll, mod);

    auto matk = fast_pow_mat(mat, k);
    auto statek = matk * state;

    cout << statek[one].nd << endl;

    return 0;
}