p476.ans1
======================
10

=================
p476.ans2
======================
8

=================
p476.cpp
======================
#include <bits/stdc++.h>

using namespace std;

template<typename T1, typename T2>
ostream& operator<<(ostream& out, const pair<T1, T2>& x) {
    return out << x.first << ' ' << x.second;
}

template<typename T1, typename T2>
istream& operator>>(istream& in, pair<T1, T2>& x) {
    return in >> x.first >> x.second;
}

template<typename T>
istream& operator>>(istream& in, vector<T>& a) {
    for(auto& x: a) {
        in >> x;
    }
    return in;
};

template<typename T>
ostream& operator<<(ostream& out, const vector<T>& a) {
    for(auto x: a) {
        out << x << ' ';
    }
    return out;
};

// base and base_digits must be consistent
const int base = 1000000000;
const int base_digits = 9;

struct bigint {
    vector<int> z;
    int sign;

    bigint() : sign(1) {}

    bigint(long long v) { *this = v; }

    bigint(const string& s) { read(s); }

    void operator=(const bigint& v) {
        sign = v.sign;
        z = v.z;
    }

    void operator=(long long v) {
        sign = 1;
        if(v < 0) {
            sign = -1, v = -v;
        }
        z.clear();
        for(; v > 0; v = v / base) {
            z.push_back(v % base);
        }
    }

    bigint operator+(const bigint& v) const {
        if(sign == v.sign) {
            bigint res = v;

            for(int i = 0, carry = 0;
                i < (int)max(z.size(), v.z.size()) || carry; ++i) {
                if(i == (int)res.z.size()) {
                    res.z.push_back(0);
                }
                res.z[i] += carry + (i < (int)z.size() ? z[i] : 0);
                carry = res.z[i] >= base;
                if(carry) {
                    res.z[i] -= base;
                }
            }
            return res;
        }
        return *this - (-v);
    }

    bigint operator-(const bigint& v) const {
        if(sign == v.sign) {
            if(abs() >= v.abs()) {
                bigint res = *this;
                for(int i = 0, carry = 0; i < (int)v.z.size() || carry; ++i) {
                    res.z[i] -= carry + (i < (int)v.z.size() ? v.z[i] : 0);
                    carry = res.z[i] < 0;
                    if(carry) {
                        res.z[i] += base;
                    }
                }
                res.trim();
                return res;
            }
            return -(v - *this);
        }
        return *this + (-v);
    }

    void operator*=(int v) {
        if(v < 0) {
            sign = -sign, v = -v;
        }
        for(int i = 0, carry = 0; i < (int)z.size() || carry; ++i) {
            if(i == (int)z.size()) {
                z.push_back(0);
            }
            long long cur = z[i] * (long long)v + carry;
            carry = (int)(cur / base);
            z[i] = (int)(cur % base);
            // asm("divl %%ecx" : "=a"(carry), "=d"(a[i]) : "A"(cur),
            // "c"(base));
        }
        trim();
    }

    bigint operator*(int v) const {
        bigint res = *this;
        res *= v;
        return res;
    }

    friend pair<bigint, bigint> divmod(const bigint& a1, const bigint& b1) {
        int norm = base / (b1.z.back() + 1);
        bigint a = a1.abs() * norm;
        bigint b = b1.abs() * norm;
        bigint q, r;
        q.z.resize(a.z.size());

        for(int i = a.z.size() - 1; i >= 0; i--) {
            r *= base;
            r += a.z[i];
            int s1 = b.z.size() < r.z.size() ? r.z[b.z.size()] : 0;
            int s2 = b.z.size() - 1 < r.z.size() ? r.z[b.z.size() - 1] : 0;
            int d = ((long long)s1 * base + s2) / b.z.back();
            r -= b * d;
            while(r < 0) {
                r += b, --d;
            }
            q.z[i] = d;
        }

        q.sign = a1.sign * b1.sign;
        r.sign = a1.sign;
        q.trim();
        r.trim();
        return make_pair(q, r / norm);
    }

    friend bigint sqrt(const bigint& a1) {
        bigint a = a1;
        while(a.z.empty() || a.z.size() % 2 == 1) {
            a.z.push_back(0);
        }

        int n = a.z.size();

        int firstDigit = (int)sqrt((double)a.z[n - 1] * base + a.z[n - 2]);
        int norm = base / (firstDigit + 1);
        a *= norm;
        a *= norm;
        while(a.z.empty() || a.z.size() % 2 == 1) {
            a.z.push_back(0);
        }

        bigint r = (long long)a.z[n - 1] * base + a.z[n - 2];
        firstDigit = (int)sqrt((double)a.z[n - 1] * base + a.z[n - 2]);
        int q = firstDigit;
        bigint res;

        for(int j = n / 2 - 1; j >= 0; j--) {
            for(;; --q) {
                bigint r1 =
                    (r - (res * 2 * base + q) * q) * base * base +
                    (j > 0 ? (long long)a.z[2 * j - 1] * base + a.z[2 * j - 2]
                           : 0);
                if(r1 >= 0) {
                    r = r1;
                    break;
                }
            }
            res *= base;
            res += q;

            if(j > 0) {
                int d1 =
                    res.z.size() + 2 < r.z.size() ? r.z[res.z.size() + 2] : 0;
                int d2 =
                    res.z.size() + 1 < r.z.size() ? r.z[res.z.size() + 1] : 0;
                int d3 = res.z.size() < r.z.size() ? r.z[res.z.size()] : 0;
                q = ((long long)d1 * base * base + (long long)d2 * base + d3) /
                    (firstDigit * 2);
            }
        }

        res.trim();
        return res / norm;
    }

    bigint operator/(const bigint& v) const { return divmod(*this, v).first; }

    bigint operator%(const bigint& v) const { return divmod(*this, v).second; }

    void operator/=(int v) {
        if(v < 0) {
            sign = -sign, v = -v;
        }
        for(int i = (int)z.size() - 1, rem = 0; i >= 0; --i) {
            long long cur = z[i] + rem * (long long)base;
            z[i] = (int)(cur / v);
            rem = (int)(cur % v);
        }
        trim();
    }

    bigint operator/(int v) const {
        bigint res = *this;
        res /= v;
        return res;
    }

    int operator%(int v) const {
        if(v < 0) {
            v = -v;
        }
        int m = 0;
        for(int i = z.size() - 1; i >= 0; --i) {
            m = (z[i] + m * (long long)base) % v;
        }
        return m * sign;
    }

    void operator+=(const bigint& v) { *this = *this + v; }
    void operator-=(const bigint& v) { *this = *this - v; }
    void operator*=(const bigint& v) { *this = *this * v; }
    void operator/=(const bigint& v) { *this = *this / v; }

    bool operator<(const bigint& v) const {
        if(sign != v.sign) {
            return sign < v.sign;
        }
        if(z.size() != v.z.size()) {
            return z.size() * sign < v.z.size() * v.sign;
        }
        for(int i = z.size() - 1; i >= 0; i--) {
            if(z[i] != v.z[i]) {
                return z[i] * sign < v.z[i] * sign;
            }
        }
        return false;
    }

    bool operator>(const bigint& v) const { return v < *this; }
    bool operator<=(const bigint& v) const { return !(v < *this); }
    bool operator>=(const bigint& v) const { return !(*this < v); }
    bool operator==(const bigint& v) const {
        return !(*this < v) && !(v < *this);
    }
    bool operator!=(const bigint& v) const { return *this < v || v < *this; }

    void trim() {
        while(!z.empty() && z.back() == 0) {
            z.pop_back();
        }
        if(z.empty()) {
            sign = 1;
        }
    }

    bool isZero() const { return z.empty() || (z.size() == 1 && !z[0]); }

    bigint operator-() const {
        bigint res = *this;
        res.sign = -sign;
        return res;
    }

    bigint abs() const {
        bigint res = *this;
        res.sign *= res.sign;
        return res;
    }

    long long longValue() const {
        long long res = 0;
        for(int i = z.size() - 1; i >= 0; i--) {
            res = res * base + z[i];
        }
        return res * sign;
    }

    friend bigint gcd(const bigint& a, const bigint& b) {
        return b.isZero() ? a : gcd(b, a % b);
    }
    friend bigint lcm(const bigint& a, const bigint& b) {
        return a / gcd(a, b) * b;
    }

    void read(const string& s) {
        sign = 1;
        z.clear();
        int pos = 0;
        while(pos < (int)s.size() && (s[pos] == '-' || s[pos] == '+')) {
            if(s[pos] == '-') {
                sign = -sign;
            }
            ++pos;
        }
        for(int i = s.size() - 1; i >= pos; i -= base_digits) {
            int x = 0;
            for(int j = max(pos, i - base_digits + 1); j <= i; j++) {
                x = x * 10 + s[j] - '0';
            }
            z.push_back(x);
        }
        trim();
    }

    friend istream& operator>>(istream& stream, bigint& v) {
        string s;
        stream >> s;
        v.read(s);
        return stream;
    }

    friend ostream& operator<<(ostream& stream, const bigint& v) {
        if(v.sign == -1) {
            stream << '-';
        }
        stream << (v.z.empty() ? 0 : v.z.back());
        for(int i = (int)v.z.size() - 2; i >= 0; --i) {
            stream << setw(base_digits) << setfill('0') << v.z[i];
        }
        return stream;
    }

    static vector<int> convert_base(
        const vector<int>& a, int old_digits, int new_digits
    ) {
        vector<long long> p(max(old_digits, new_digits) + 1);
        p[0] = 1;
        for(int i = 1; i < (int)p.size(); i++) {
            p[i] = p[i - 1] * 10;
        }
        vector<int> res;
        long long cur = 0;
        int cur_digits = 0;
        for(int i = 0; i < (int)a.size(); i++) {
            cur += a[i] * p[cur_digits];
            cur_digits += old_digits;
            while(cur_digits >= new_digits) {
                res.push_back(int(cur % p[new_digits]));
                cur /= p[new_digits];
                cur_digits -= new_digits;
            }
        }
        res.push_back((int)cur);
        while(!res.empty() && res.back() == 0) {
            res.pop_back();
        }
        return res;
    }

    typedef vector<long long> vll;

    static vll karatsubaMultiply(const vll& a, const vll& b) {
        int n = a.size();
        vll res(n + n);
        if(n <= 32) {
            for(int i = 0; i < n; i++) {
                for(int j = 0; j < n; j++) {
                    res[i + j] += a[i] * b[j];
                }
            }
            return res;
        }

        int k = n >> 1;
        vll a1(a.begin(), a.begin() + k);
        vll a2(a.begin() + k, a.end());
        vll b1(b.begin(), b.begin() + k);
        vll b2(b.begin() + k, b.end());

        vll a1b1 = karatsubaMultiply(a1, b1);
        vll a2b2 = karatsubaMultiply(a2, b2);

        for(int i = 0; i < k; i++) {
            a2[i] += a1[i];
        }
        for(int i = 0; i < k; i++) {
            b2[i] += b1[i];
        }

        vll r = karatsubaMultiply(a2, b2);
        for(int i = 0; i < (int)a1b1.size(); i++) {
            r[i] -= a1b1[i];
        }
        for(int i = 0; i < (int)a2b2.size(); i++) {
            r[i] -= a2b2[i];
        }

        for(int i = 0; i < (int)r.size(); i++) {
            res[i + k] += r[i];
        }
        for(int i = 0; i < (int)a1b1.size(); i++) {
            res[i] += a1b1[i];
        }
        for(int i = 0; i < (int)a2b2.size(); i++) {
            res[i + n] += a2b2[i];
        }
        return res;
    }

    bigint operator*(const bigint& v) const {
        vector<int> a6 = convert_base(this->z, base_digits, 6);
        vector<int> b6 = convert_base(v.z, base_digits, 6);
        vll a(a6.begin(), a6.end());
        vll b(b6.begin(), b6.end());
        while(a.size() < b.size()) {
            a.push_back(0);
        }
        while(b.size() < a.size()) {
            b.push_back(0);
        }
        while(a.size() & (a.size() - 1)) {
            a.push_back(0), b.push_back(0);
        }
        vll c = karatsubaMultiply(a, b);
        bigint res;
        res.sign = sign * v.sign;
        for(int i = 0, carry = 0; i < (int)c.size(); i++) {
            long long cur = c[i] + carry;
            res.z.push_back((int)(cur % 1000000));
            carry = (int)(cur / 1000000);
        }
        res.z = convert_base(res.z, 6, base_digits);
        res.trim();
        return res;
    }
};

int n, k;
vector<array<int, 3>> triplets;

void read() {
    cin >> n >> k;
    triplets.resize(k);
    for(int i = 0; i < k; i++) {
        cin >> triplets[i][0] >> triplets[i][1] >> triplets[i][2];
    }
}

void solve() {
    // The number of ways to split 3n students into n unordered teams of three
    // with no constraints is (3n)! / (6^n * n!): order the students, then
    // divide out the 6 internal orderings of each team and the n! orderings of
    // the teams.
    //
    // The k forbidden triples are removed by inclusion-exclusion: for a chosen
    // subset S of forbidden triples that are pairwise disjoint (occupy 3*|S|
    // distinct students), the number of partitions that contain all of them is
    // the free count of the remaining 3n - 3|S| students, i.e.
    // (3n - cnt)! / (6^(n - cnt/3) * (n - cnt/3)!) with cnt = 3|S|, weighted by
    // (-1)^|S|. Subsets containing two triples sharing a student contribute 0.
    //
    // We aggregate the signs into coef[cnt] = sum of (-1)^|S| over all disjoint
    // subsets S with 3|S| = cnt, found by backtracking that marks used students
    // and skips conflicting triples. Then the answer is the sum over cnt of
    // coef[cnt] * (3n - cnt)! / (6^(n - cnt/3) * (n - cnt/3)!). All arithmetic
    // uses bigint since the factorials grow far beyond 64 bits.

    int total = 3 * n;
    vector<bigint> factorial(total + 1), power_6(total + 1);
    factorial[0] = 1;
    power_6[0] = 1;
    for(int i = 1; i <= total; i++) {
        factorial[i] = factorial[i - 1] * i;
        power_6[i] = power_6[i - 1] * 6;
    }

    vector<char> used(total, 0);
    vector<int64_t> coef(total + 1, 0);

    function<void(int, int, int64_t)> rec = [&](int pos, int cnt, int64_t p) {
        if(pos == k) {
            coef[cnt] += p;
            return;
        }

        rec(pos + 1, cnt, p);

        bool has_conflict = false;
        vector<int> to_undo;
        for(int j: triplets[pos]) {
            if(used[j - 1]) {
                has_conflict = true;
                break;
            }

            used[j - 1] = 1;
            to_undo.push_back(j - 1);
        }

        if(!has_conflict) {
            rec(pos + 1, cnt + 3, -p);
        }

        for(int j: to_undo) {
            used[j] = 0;
        }
    };

    rec(0, 0, 1);

    bigint ans = 0;
    for(int cnt = 0; cnt <= total; cnt++) {
        if(coef[cnt] != 0) {
            int rem = n - cnt / 3;
            bigint term = factorial[total - cnt] /
                          (power_6[rem] * factorial[rem]);
            ans += term * (int)coef[cnt];
        }
    }

    cout << ans << '\n';
}

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(nullptr);

    int T = 1;
    // cin >> T;
    for(int test = 1; test <= T; test++) {
        read();
        // cout << "Case #" << test << ": ";
        solve();
    }

    return 0;
}



=================
p476.in1
======================
2 0

=================
p476.in2
======================
2 3
1 2 3
4 5 6
1 4 5

=================
p476.py
======================
import sys

sys.set_int_max_str_digits(100000)


def main():
    # We solve this with inclusion-exclusion principle. If k = 0, the answer is
    # simply (3n)!/((3!)^n * n!) - think of all permutations and then there are
    # 3! ways to rearrange each triplet and n! ways to rearrange the triplets.
    # If k > 0, we can look at all the banned triplets and subtract the count
    # given that we surely include some triplet. However, we then have some overlap
    # so we subtract an answer more than once. This yields the inclusion-exclusion
    # solution. We should only be careful about the overlap between the triplets.
    # We implement this in python as otherwise we would need to use a big integers.
    #
    # Calculating this naively would take O(2^k * k) time, but due to the bit integers
    # it will actually be O(2^k * k * n) in practice. To speed it up, we calculate the
    # coefficients of the inclusion-exclusion and do the actual calculation for each cnt
    # only once as (here coef[cnt] is the sum of the inclusion-exclusion (-1)^|mask|):
    # coef[cnt] * (3n - cnt)! / ((3!)^(n - cnt // 3) * (n - cnt // 3)!) in the end.
    #
    # The complexity this way is O(2^k * k + n^2) which is feasible for k <= 20 and n <= 10^3.
    # Depending on implementation, bitmasks might be a bit slow, so we do this via a backtracking.

    n, k = map(int, input().strip().split())
    triplets = [list(map(int, input().strip().split())) for _ in range(k)]

    factorial = [1] * (3 * n + 1)
    power_6 = [1] * (3 * n + 1)
    for i in range(1, 3 * n + 1):
        factorial[i] = factorial[i - 1] * i
        power_6[i] = power_6[i - 1] * 6

    used = [False] * (3 * n)
    coef = [0] * (3 * n + 1)

    def rec(pos, cnt, p):
        if pos == k:
            coef[cnt] += p
            return

        # Don't include triplet at pos.
        rec(pos + 1, cnt, p)

        # Include triplet at pos.
        has_conflict = False
        to_undo = []
        for j in triplets[pos]:
            if used[j - 1]:
                has_conflict = True
                break

            used[j - 1] = True
            to_undo.append(j - 1)

        if not has_conflict:
            rec(pos + 1, cnt + 3, -p)

        for j in to_undo:
            used[j] = False


    # Faster than iterating all masks.
    rec(0, 0, 1)

    ans = 0
    for cnt in range(0, 3 * n + 1):
        if coef[cnt] != 0:
            assert cnt % 3 == 0, "cnt must be a multiple of 3"
            ans += (
                coef[cnt]
                * factorial[3 * n - cnt]
                // (power_6[n - cnt // 3] * factorial[n - cnt // 3])
            )

    print(ans)


if __name__ == "__main__":
    main()

=================
statement.txt
======================
476. Coach's Trouble
Time limit per test: 0.5 second(s)
Memory limit: 65536 kilobytes
input: standard
output: standard

Berland is a modern country, and one could hardly find a man, who would believe evil spirits. Everybody knows that celebrating Halloween is just tribute to old tradition. However, a coach of Berland University student programming teams wants to divide students (3N in total) into teams of three to have as many demonic teams as possible at the moment of future Halloween contest. But he knows for some triples of students that they can't make a demonic team, as they haven't performed well in action in previous contests. He has a list of these K forbidden triples. The coach supposes that any three students can make a demonic team unless they do not form a triple contained in the list. And now he wants to know the number of such partitions that all teams are demonic.

Input
The first line of the input contains two integers N and K, separated by one space (1 ≤ N ≤ 1000, 0 ≤ K ≤ 20). Next K lines contain three integers each ai, bi, ci (1 ≤ i ≤ K, 1 ≤ ai, bi, ci ≤ 3N). All triples are unique, that is they all are diffent as sets, and ai ≠q bi, ai ≠q ci, bi ≠q ci.

Output
The output should contain the only number without leading zeroes — the answer to the task.

Example(s)
sample input
sample output
2 0
10

sample input
sample output
2 3
1 2 3
4 5 6
1 4 5
8

=================
