結果

問題 No.310 2文字しりとり
ユーザー anta
提出日時 2015-12-02 23:34:23
言語 C++11（廃止可能性あり）
(gcc 13.3.0)
結果
AC  
実行時間 630 ms / 6,000 ms
コード長 27,316 bytes
コンパイル時間 2,183 ms
コンパイル使用メモリ 116,088 KB
実行使用メモリ 131,072 KB
最終ジャッジ日時 2024-09-14 08:09:48
合計ジャッジ時間 6,830 ms
ジャッジサーバーID
（参考情報） 		judge1 / judge5
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
other AC * 28
権限があれば一括ダウンロードができます
コンパイルメッセージ
main.cpp: In function ‘int main()’:
main.cpp:634:26: warning: format ‘%d’ expects argument of type ‘int’, but argument 2 has type ‘mint’ {aka ‘ModInt<1000000007>’} [-Wformat=]
  634 |                 printf("%d\n", ans);
      |                         ~^     ~~~
      |                          |     |
      |                          int   mint {aka ModInt<1000000007>}

ソースコード

diff #
プレゼンテーションモードにする

////////////////////////////////////////////////////////
//テストとかジェネレーターとか全部残してある
#include <iostream>
#include <cstdio>
#include <algorithm>
#include <numeric>
#include <vector>
#include <random>
#ifdef NDEBUG
#undef NDEBUG
#endif
#include <cassert>
#define rep(i,n) for(int (i)=0;(i)<(int)(n);++(i))
#define rer(i,l,u) for(int (i)=(int)(l);(i)<=(int)(u);++(i))
#define reu(i,l,u) for(int (i)=(int)(l);(i)<(int)(u);++(i))
using namespace std;
////////////////////////////////////////////////////////
//ModInt
template<int MOD>
struct ModInt {
    static const int Mod = MOD;
    unsigned x;
    ModInt() : x(0) {}
    ModInt(signed sig) { int sigt = sig % MOD; if(sigt < 0) sigt += MOD; x = sigt; }
    ModInt(signed long long sig) { int sigt = sig % MOD; if(sigt < 0) sigt += MOD; x = sigt; }
    int get() const { return (int)x; }
    ModInt &operator+=(ModInt that) { if((x += that.x) >= MOD) x -= MOD; return *this; }
    ModInt &operator-=(ModInt that) { if((x += MOD - that.x) >= MOD) x -= MOD; return *this; }
    ModInt &operator*=(ModInt that) { x = (unsigned long long)x * that.x % MOD; return *this; }
    ModInt &operator/=(ModInt that) { return *this *= that.inverse(); }
    ModInt operator+(ModInt that) const { return ModInt(*this) += that; }
    ModInt operator-(ModInt that) const { return ModInt(*this) -= that; }
    ModInt operator*(ModInt that) const { return ModInt(*this) *= that; }
    ModInt operator/(ModInt that) const { return ModInt(*this) /= that; }
    ModInt inverse() const {
        signed a = x, b = MOD, u = 1, v = 0;
        while(b) {
            signed t = a / b;
            a -= t * b; std::swap(a, b);
            u -= t * v; std::swap(u, v);
        }
        if(u < 0) u += Mod;
        ModInt res; res.x = (unsigned)u;
        return res;
    }
    bool operator==(ModInt that) const { return x == that.x; }
    bool operator!=(ModInt that) const { return x != that.x; }
    ModInt operator-() const { ModInt t; t.x = x == 0 ? 0 : Mod - x; return t; }
};
typedef ModInt<1000000007> mint;
////////////////////////////////////////////////////////
//Black box linear algebra
void truncateLeadingZerosOfPolynomial(const vector<mint> &p, vector<mint> &res) {
    int n = (int)p.size();
    while(n > 0 && p[n-1].x == 0) -- n;
    res.assign(p.begin(), p.begin() + n);
}
int extendedEuclideanAlgorithm(const vector<mint> &f, const vector<mint> &g, int k, vector<mint> &out_r, vector<mint> &out_t, bool traditional = 
    false) {
    if((!f.empty() && f.back().x == 0) || (!g.empty() && g.back().x == 0)) {
        vector<mint> nf, ng;
        truncateLeadingZerosOfPolynomial(f, nf);
        truncateLeadingZerosOfPolynomial(g, ng);
        return extendedEuclideanAlgorithm(nf, ng, k, out_r, out_t, traditional);
    }
    if(f.size() < g.size())
        return extendedEuclideanAlgorithm(g, f, k, out_r, out_t, traditional);
    assert(k >= 0);
    int n[2];
    vector<mint> r[2], t[2];//, s[2];
    int nf = (int)f.size() - 1, ng = (int)g.size() - 1;
    n[0] = nf, n[1] = ng;
    mint ip0 = traditional || f.empty() ? mint(1) : f.back().inverse();
    mint ip1 = traditional || g.empty() ? mint(1) : g.back().inverse();
    int bufsize = n[0] == -1 ? 1 : n[0] + 1;
    rep(i, 2) {
        r[i].resize(bufsize);
//      s[i].resize(bufsize);
        t[i].resize(bufsize);
    }
    rer(j, 0, n[0]) r[0][j] = f[j] * ip0;
    rer(j, 0, n[1]) r[1][j] = g[j] * ip1;
//  s[0][0] = ip0;
    t[1][0] = ip1;
    if(n[0] < k) {
        out_r.swap(r[0]); out_r.resize(n[0] + 1);
        out_t.swap(t[0]); out_t.resize(1);
        return 0;
    }else if(n[1] < k) {
        out_r.swap(r[1]); out_r.resize(n[1] + 1);
        out_t.swap(t[1]); out_t.resize(1);
        return 1;
    }
    for(int i = 1; n[1] >= 0; ++ i) {
        swap(n[0], n[1]);
        r[0].swap(r[1]);
//      s[0].swap(s[1]);
        t[0].swap(t[1]);
        //r[0] = r_i; r[1] = r_{i-1}, r_{i+1}
        mint ilc = r[0][n[0]].inverse();
        assert(traditional || ilc.x == 1);
        //q_i     <- r_{i-1} quo r_i
        //r_{i+1} <- r_{i-1} - q_i r_i
        //s_{i+1} <- s_{i-1} - q_i s_i
        //t_{i+1} <- t_{i-1} - q_i t_i
        int nr = n[0], ns = ng - n[1], nt = nf - n[1];
        for(int j = n[1] - nr; j >= 0; -- j) {
            mint c = r[1][nr + j] * ilc;
            for(int l = 0; l <= nr; ++ l)
                r[1][j + l] -= c * r[0][l];
//          for(int l = 0; l <= ns; ++ l)
//              s[1][j + l] -= c * s[0][l];
            for(int l = 0; l <= nt; ++ l)
                t[1][j + l] -= c * t[0][l];
        }
        while(n[1] >= 0 && r[1][n[1]].x == 0)
            -- n[1];
        mint ip = n[1] == -1 ? mint(1) : r[1][n[1]].inverse();
        if(!traditional) {
            rer(j, 0, n[1]) r[1][j] *= ip;
//          rer(j, 0, ng - n[0]) s[1][j] *= ip;
            rer(j, 0, nf - n[0]) t[1][j] *= ip;
        }
        if(n[1] < k) {
            out_r.swap(r[1]); out_r.resize(n[1] + 1);
            out_t.swap(t[1]); out_t.resize(nf - n[0] + 1);
            return i+1;
        }
    }
    assert(false);
    return -1;
}
//体上の数列のminimum polynomialを計算する。
//入力がlinearly recurrentでない場合は、res.back() == 0 になる。
void computeMinimumPolynomialForLinearlyRecurrentSequence(const vector<mint> &a, vector<mint> &res) {
    int n2 = (int)a.size(), n = n2 / 2;
    assert(n2 % 2 == 0);
    vector<mint> x(n2 + 1), s, t;
    x[n2] = 1;
    int i = extendedEuclideanAlgorithm(x, a, n, s, t);
    if(s.size() == n2 + 1) {    //a is the zero sequence
        res.assign(1, mint());
        res[0] = mint(1);
        return;
    }
    
    int d = max((int)s.size(), (int)t.size() - 1);
    int e = 0;
    while(t[e].x == 0) ++ e;
    mint c = t[e].inverse();
    res.assign(d + 1, mint());
    reu(i, e, t.size())
        res[d - i] = t[i] * c;
}
struct RandomModInt {
    default_random_engine re;
    uniform_int_distribution<int> dist;
#ifndef _DEBUG
    RandomModInt() : re(random_device{}()), dist(1, mint::Mod - 1) { }
#else
    RandomModInt() : re(), dist(1, mint::Mod - 1) { }
#endif
    mint operator()() {
        mint r;
        r.x = dist(re);
        return r;
    }
} randomModInt;
void randomModIntVector(vector<mint> &v) {
    int n = (int)v.size();
    for(int i = 0; i < n; ++ i)
        v[i] = randomModInt();
}
bool verifyMinimalPolynomialForKrylovSubspace(const vector<vector<mint> > &Aib, const vector<mint> &m) {
    int n = (int)Aib[0].size();
    assert(!m.empty());
    int d = (int)m.size() - 1;
    assert(d < (int)Aib.size());
    vector<mint> sum(n);
    for(int i = 0; i <= d; ++ i) {
        mint coeff = m[i];
        const vector<mint> &v = Aib[i];
        for(int j = 0; j < n; ++ j)
            sum[j] += coeff * v[j];
    }
    for(int j = 0; j < n; ++ j)
        if(sum[j] != mint())
            return false;
    return true;
}
//誤った答えを返すことがある (ただし minimal polynomial の約数である)
//確実にしたい場合、verifyMinimalPolynomialForKrylovSubspaceでverifyされるまでイテレーションすればよい。
void minimalPolynomialForKrylovSubspace(const vector<vector<mint> > &Aib, vector<mint> &res) {
    int n = (int)Aib[0].size();
    assert(Aib.size() == n * 2);
    vector<mint> u(n);
    randomModIntVector(u);
    vector<mint> uTAib(n * 2);
    for(int i = 0; i < n * 2; ++ i)
        uTAib[i] = inner_product(u.begin(), u.end(), Aib[i].begin(), mint());
    computeMinimumPolynomialForLinearlyRecurrentSequence(uTAib, res);
}
//誤った答えを返すことがある (ただし minimal polynomial の約数である)
//A(b,res): res += Ab
template<typename BlackBox>
void blackBoxMinimalPolynomial(const BlackBox &A, int n, vector<mint> &m) {
    if(n == 0) {
        m.assign(1, mint(1));
        return;
    }
    vector<vector<mint> > Aib(n * 2, vector<mint>(n));
    vector<mint> &b = Aib[0];
    randomModIntVector(b);
    for(int i = 0; i < n * 2 - 1; ++ i)
        A(Aib[i], Aib[i + 1]);
    minimalPolynomialForKrylovSubspace(Aib, m);
}
template<typename BlackBox>
struct DiagonalPreconditionBox {
    const BlackBox &A;
    const vector<mint> D;
    DiagonalPreconditionBox(const BlackBox &A, const vector<mint> &D): A(A), D(D) { }
    //res = ADb
    void operator()(const vector<mint> &b, vector<mint> &res) const {
        int n = (int)b.size();
        assert(D.size() == n);
        vector<mint> Db(n);
        for(int i = 0; i < n; ++ i)
            Db[i] = D[i] * b[i];
        A(Db, res);
    }
};
//常に正しい答えを返す
template<typename BlackBox>
mint blackBoxDeterminant(const BlackBox &A, int n) {
    vector<mint> D(n);
    vector<mint> m;
    do {
        randomModIntVector(D);
        DiagonalPreconditionBox<BlackBox> AD(A, D);
        blackBoxMinimalPolynomial(AD, n, m);
        assert(!m.empty());
        if(m[0] == mint()) return mint();
    } while((int)m.size() < n + 1);
    //m = char(AD)
    assert(m.size() == n + 1);
    assert(m.back() == mint(1));
    mint detD = 1;
    for(int i = 0; i < n; ++ i)
        detD *= D[i];
    mint invdetD = detD.inverse();
    mint detA = m[0] * invdetD;
    if(n % 2 == 1)
        detA = mint() - detA;
    return detA;
}
////////////////////////////////////////////////////////
//ユーティリティー
#ifdef NDEBUG
#error assert is disabled!
#endif
template<typename T>
T readNatural(T lo, T up) {
    assert(0 <= lo && lo <= up);
    T x = 0;
    while(1) {
        int d = getchar();
        if(!('0' <= d && d <= '9')) {
            ungetc(d, stdin);
            break;
        }
        d -= '0';
        assert(d <= up && x <= (up - d) / 10);
        x = x * 10 + d;
    }
    assert(lo <= x && x <= up);
    return x;
}
void readSpace() { int c = getchar(); assert(c == ' '); }
static bool read_eof = false;
void readEOL() {
    int c = getchar();
    if(c == EOF) {
        assert(!read_eof);
        read_eof = true;
    }else {
        assert(c == '\n');
    }
}
void readEOF() {
    assert(!read_eof);
    int c = getchar();
    assert(c == EOF);
    read_eof = true;
}
struct UnionFind {
    vector<int> data;
    void init(int n) { data.assign(n, -1); }
    bool unionSet(int x, int y) {
        x = root(x); y = root(y);
        if(x != y) {
            if(data[y] < data[x]) swap(x, y);
            data[x] += data[y]; data[y] = x;
        }
        return x != y;
    }
    bool findSet(int x, int y) { return root(x) == root(y); }
    int root(int x) { return data[x] < 0 ? x : data[x] = root(data[x]); }
    int size(int x) { return -data[root(x)]; }
};
////////////////////////////////////////////////////////
//解答メイン部分
#ifdef MY_LOCAL_RUN
//#define RANDOM_TEST
//#define TIME_CHECK
//#define NAIVE_CHECK
#endif
mint countEulerianCyclesOnComplementGraph(int N, const vector<pair<int,int> > &edges, int src, int dst) {
    assert(N > 0);
    vector<int> loops(N, 1);
    vector<int> outDeg(N, N), inDeg(N, N);
    //サイクルにならない場合は辺を1つ足してサイクルにするようにする
    if(src != -1) {
        assert(src != dst);
        ++ outDeg[dst];
        ++ inDeg[src];
    }
    for(auto e : edges) {
        -- outDeg[e.first];
        -- inDeg[e.second];
        if(e.first == e.second)
            -- loops[e.first];
    }
    rep(i, N)
        assert(inDeg[i] == outDeg[i]);
    //BEST theoremによりオイラー閉路を数えられる
    //<https://en.wikipedia.org/wiki/BEST_theorem>
    vector<mint> fact(N + 1);
    fact[0] = 1;
    rer(n, 1, N)
        fact[n] = fact[n - 1] * n;
    //matrix tree theorem によって有向木を数える
    //<https://en.wikipedia.org/wiki/Kirchhoff%27s_theorem#Kirchhoff.27s_theorem_for_directed_multigraphs>
    //頂点 N - 1 を根とする。
    //1行削除されるため、それに接する辺も一緒に削除しておく
    //また、ループも削除しておく
    vector<pair<int, int> > validEdges;
    for(auto e : edges) if(e.first < N - 1 && e.second < N - 1 && e.first != e.second)
        validEdges.push_back(e);
    //なんとなくソートしておく
    sort(validEdges.begin(), validEdges.end());
    //行列の対角要素
    vector<mint> diag(N - 1);
    rep(i, N - 1)
        diag[i] = outDeg[i] - loops[i];
    mint res = blackBoxDeterminant([N, &diag, &validEdges, src, dst](const vector<mint> &b, vector<mint> &res) {
        assert(b.size() == N - 1 && res.size() == N - 1);
        //A[i][j] where i != j:
        //= -(1 + (i == dst && j == src) - ((i,j) in edges))
        mint sumb;
        rep(i, N - 1)
            sumb += b[i];
        //自分以外の場所から1引かれる
        rep(i, N - 1)
            res[i] -= sumb - b[i];
        if(src != -1) {
            //dst -> src の辺
            if(dst < N - 1 && src < N - 1)
                res[dst] -= b[src];
        }
        for(auto e : validEdges)
            res[e.first] += b[e.second];
        //A[i][i] = diag[i]
        rep(i, N - 1)
            res[i] += diag[i] * b[i];
    }, N - 1);
    rep(i, N) {
        int deg = outDeg[i];
        assert(deg > 0);
        res *= fact[deg - 1];
    }
    return res;
}
mint naiveSolve(const vector<vector<bool> > &g) {
    int N = (int)g.size();
    vector<pair<int, int> > aliveEdges;
    rep(i, N) rep(j, N) if(g[i][j])
        aliveEdges.push_back({i, j});
    int E = (int)aliveEdges.size();
    assert(E <= 25);
    if(E == 0)
        return 1;
    vector<mint> dp(N << E);
    rep(i, N)
        dp[0 * N + i] += 1;
    rep(s, 1 << E) rep(i, N) {
        mint x = dp[s * N + i];
        if(x.x == 0) continue;
        rep(j, E) if((~s >> j & 1) && aliveEdges[j].first == i)
            dp[(s | 1 << j) * N + aliveEdges[j].second] += x;
    }
    mint res;
    rep(i, N)
        res += dp[((1 << E)-1) * N + i];
    return res;
}
void naiveCheck(const vector<vector<bool> > &g, mint ans) {
#ifdef NAIVE_CHECK
    mint ansnaive = naiveSolve(g);
    if(ansnaive.get() != ans.get()) {
        cerr << ansnaive.get() << " != " << ans.get() << endl;
        abort();
    }
#endif
}
#ifdef TIME_CHECK
struct FILETIME {
    unsigned dwLowDateTime, dwHighDateTime;
};
extern "C" void* __stdcall GetCurrentProcess(void);
extern "C" int __stdcall GetProcessTimes(void *hProcess, FILETIME *lpCreationTime, FILETIME *lpExitTime, FILETIME *lpKernelTime, FILETIME 
    *lpUserTime);
void getCPUTime(double &userTime, double &systemTime) {
    void *handle = GetCurrentProcess();
    FILETIME dummy1, dummy2, kernel, user;
    GetProcessTimes(handle, &dummy1, &dummy2, &kernel, &user);
    userTime = user.dwHighDateTime * 429.4967296 + user.dwLowDateTime * 1e-7;
    systemTime = kernel.dwHighDateTime * 429.4967296 + kernel.dwLowDateTime * 1e-7;
}
double getCPUTime() {
    double user, sys;
    getCPUTime(user, sys);
    return user + sys;
}
struct CPUTimeIt {
    double user, sys;
    const char *msg;
    CPUTimeIt(const char *msg_): msg(msg_) { getCPUTime(user, sys); }
    ~CPUTimeIt() {
        double userEnd, sysEnd;
        getCPUTime(userEnd, sysEnd);
        fprintf(stderr, "%s: user %.6fs / sys %.6fs\n", msg, userEnd - user, sysEnd - sys);
    }
    operator bool() { return false; }
};
#define CPUTIMEIT(s) if(CPUTimeIt cputimeit_##__LINE__ = s); else
#endif
namespace test_det {
void test();
void test_classical_time();
}
int main() {
//  freopen("C:/test/sagyo/874/test_0.txt", "r", stdin);
#ifdef RANDOM_TEST
    default_random_engine re(random_device{}());
#endif
#ifdef TIME_CHECK
    test_det::test_classical_time();
#endif
    do
#ifdef TIME_CHECK
    CPUTIMEIT("main")
#endif
    {
        const int MaxN = 4000, MaxM = 4000;
#ifndef RANDOM_TEST
        int N = readNatural(1, MaxN); readSpace();
        int M = readNatural(0, MaxM); readEOL();
#else
#ifndef TIME_CHECK
        int N = uniform_int_distribution<int>(3, 3)(re);
        int M = uniform_int_distribution<int>(0, N * N)(re);
#else
        int N = MaxN;
        int M = MaxM;
#endif
#endif
        vector<vector<bool> > g(N, vector<bool>(N, true));
        vector<int> outDeg(N, N), inDeg(N, N);
        vector<pair<int, int> > edges(M);
#if defined(RANDOM_TEST) && defined(TIME_CHECK)
        vector<int> randPerm(N);
        iota(randPerm.begin(), randPerm.end(), 0);
        shuffle(randPerm.begin(), randPerm.end(), re);
#endif
        rep(i, M) {
#ifndef RANDOM_TEST
            int A = readNatural(1, N) - 1; readSpace();
            int B = readNatural(1, N) - 1; readEOL();
#else
#ifndef TIME_CHECK
            int A, B;
            do {
                A = uniform_int_distribution<int>(0, N - 1)(re);
                B = uniform_int_distribution<int>(0, N - 1)(re);
            } while(!g[A][B]);
#else
            int A = randPerm[i % N], B = randPerm[(i + 1) % N];
#endif
#endif
            assert(g[A][B]);
            g[A][B] = false;
            edges[i] = {A, B};
            -- outDeg[A];
            -- inDeg[B];
        }
        //空グラフは場合分けしておく
        if(M == N * N) {
            naiveCheck(g, 1);
            puts("1");
            continue;
        }
        //ちなみに、ある補グラフ上でのDFSやBFSは線形時間で可能である
        UnionFind uf; uf.init(N);
        rep(i, N) rep(j, N) if(g[i][j])
            uf.unionSet(i, j);
        //まず、eulerianかどうか判定する
        int mainCC = -1;
        int src = -1, dst = -1;
        bool eulerian = true;
        rep(i, N) {
            if(outDeg[i] == inDeg[i]) {
            } else if(outDeg[i] == inDeg[i] + 1) {
                eulerian &= src == -1;
                src = i;
            } else if(outDeg[i] == inDeg[i] - 1) {
                eulerian &= dst == -1;
                dst = i;
            } else {
                eulerian = false;
            }
            if(mainCC == -1 && outDeg[i] > 0)
                mainCC = i;
        }
        assert(mainCC != -1);
        vector<int> vertices;
        rep(i, N) {
            if(outDeg[i] > 0 || inDeg[i] > 0)
                eulerian &= uf.findSet(mainCC, i);
        }
        if(src != -1 || dst != -1)
            eulerian &= src != -1 && dst != -1;
        if(!eulerian) {
            naiveCheck(g, 0);
            puts("0");
            continue;
        }
        //孤立点を取り除いたグラフを作る
        vector<int> vertexIndex(N, -1);
        int V = 0;
        rep(i, N)
            if(outDeg[i] > 0 || inDeg[i] > 0)
                vertexIndex[i] = V ++;
        if(src != -1) {
            src = vertexIndex[src];
            dst = vertexIndex[dst];
            assert(src != -1 && dst != -1);
        }
        for(auto &e : edges) {
            e.first = vertexIndex[e.first];
            e.second = vertexIndex[e.second];
            if(e.first == -1 || e.second == -1)
                e = {-1, -1};
        }
        edges.erase(remove(edges.begin(), edges.end(), make_pair(-1, -1)), edges.end());
        mint ans = countEulerianCyclesOnComplementGraph(V, edges, src, dst);
        //サイクルの場合、どの辺から始めるかの分をかける
        if(src == -1)
            ans *= N * N - M;
        naiveCheck(g, ans);
        printf("%d\n", ans);
#ifdef MY_LOCAL_RUN
    } while(1);
#else
    } while(0);
#endif
    readEOF();
    return 0;
}
////////////////////////////////////////////////////////
//テスト
namespace test_det {
namespace classical {
typedef mint Num;
typedef vector<Num> Vec;
typedef vector<Vec> Mat;
 
Num determinant(Mat A) {
    int n = (int)A.size();
    if(!A.empty() && n != A[0].size()) return Num();
    Num det = 1;
    for(int i = 0; i < n; ++ i) {
        int pi = i;
        while(pi < n && A[pi][i] == Num()) ++ pi;
        if(pi == n) return Num();
        A[pi].swap(A[i]);
        det *= A[i][i];
        if(pi != i) det = -det;
        Num inv = Num(1) / A[i][i];
        for(int j = i + 1; j < n; ++ j) {
            Num t = A[j][i] * inv;
            if(t == Num()) continue;
            for(int k = i; k < n; ++ k)
                A[j][k] -= t * A[i][k];
        }
    }
    return det;
}
}
void test() {
    default_random_engine re;
    uniform_int_distribution<int> dist(1, mint::Mod);
    rep(ii, 100000) {
        const int X = mint::Mod - 1;
        int n = uniform_int_distribution<int>(1, 10)(re);
        classical::Mat A(n, classical::Vec(n));
        rep(i, n) rep(j, n)
            A[i][j].x = uniform_int_distribution<int>(0, X)(re);
        mint det1 = classical::determinant(A);
        mint det2 = blackBoxDeterminant([n, &A](const vector<mint> &b, vector<mint> &res) {
            assert(res.size() == n && count(res.begin(), res.end(), mint()) == n);
            rep(i, n)
                res[i] = inner_product(b.begin(), b.end(), A[i].begin(), res[i]);
        }, n);
        if(det1 != det2) {
            cerr << det1.get() << " != " << det2.get() << endl;
            rep(i, n) {
                rep(j, n)
                    cerr << A[i][j].get() << " ";
                cerr << endl;
            }
        }
    }
}
void test_classical_time() {
#ifdef TIME_CHECK
    default_random_engine re{777};
    int N = 1000;
    cerr << "N = " << N << endl;
    vector<vector<mint> > A(N, vector<mint>(N, mint(1)));
    rep(i, N)
        A[i][i] = mint();
    rep(k, N) {
        int i = uniform_int_distribution<int>(0, N - 1)(re);
        int j = uniform_int_distribution<int>(0, N - 1)(re);
        A[i][j] = randomModInt();
    }
    mint det;
    CPUTIMEIT("det") {
        det = classical::determinant(A);
    }
    cerr << "det = " << det.get() << endl;
#else
    abort();
#endif
}
}
////////////////////////////////////////////////////////
//テストケースジェネレーター
#if 0
//テストケース生成のためのもの
#include <string>
#include <vector>
#include <algorithm>
#include <numeric>
#include <set>
#include <map>
#include <iostream>
#include <sstream>
#include <fstream>
#include <limits>
#include <functional>
#include <random>
#include <tuple>
#include <memory>
#include <cstdio>
#include <cassert>
using namespace std;
#ifdef NDEBUG
#undef NDEBUG
#endif
#define rep(i,n) for(int (i)=0;(i)<(int)(n);++(i))
#define rer(i,l,u) for(int (i)=(int)(l);(i)<=(int)(u);++(i))
#define reu(i,l,u) for(int (i)=(int)(l);(i)<(int)(u);++(i))
#include "C:\Dropbox\backup\implements\Util\TimeIt_Simple.hpp"
//問題の識別子
const string problem_name = "874";
//テストケースを生成するディレクトリ
const string testcase_dir = "C:\\test\\sagyo\\" + problem_name;
//テストケースの拡張子。この拡張子のファイルはなんでも自動的に消されるので注意
const string testcase_postfix = ".txt";
struct RNG {
    default_random_engine re;
    RNG(): re() {}
    void addSeed(unsigned x) {
        re.seed(re() ^ x);
    }
    unsigned operator()() {
        return re();
    }
    int operator()(int N) {
        return uniform_int_distribution<int>(0, N - 1)(re);
    }
    int operator()(int min, int max) {
        return uniform_int_distribution<int>(min, max)(re);
    }
};
class TestCase {
public:
    TestCase(int no = -1) : no(no) {}
    typedef vector<pair<int, int> > Edges;
    static const int MaxN = 4000, MaxM = 4000;
    virtual void generate(RNG &rng) = 0;
    virtual string getName() const {
        return "test_" + to_string(no);
    }
    void write(ostream &os) const {
        int N;
        vector<pair<int, int> > edges;
        tie(N, edges) = get();
        int M = (int)edges.size();
        assert(1 <= N && N <= MaxN);
        assert(0 <= M && M <= MaxM);
        os << N << ' ' << M << '\n';
        vector<vector<bool> > g(N, vector<bool>(N, true));
        for(auto p : edges) {
            assert(0 <= p.first && p.first < N);
            assert(0 <= p.second && p.second < N);
            assert(g[p.first][p.second]);
            g[p.first][p.second] = false;
            os << p.first + 1 << ' ' << p.second + 1 << '\n';
        }
        os.flush();
    }
protected:
    virtual tuple<int,Edges> get() const = 0;
private:
    int no;
};
class Manual : public TestCase {
protected:
    int N;
    Edges edges;
public:
    Manual(int no, int N, Edges edges) : TestCase(no), N(N), edges(edges) {}
    void generate(RNG &rng) override {
        //nope
    }
protected:
    tuple<int,Edges> get() const override { return make_tuple(N, edges); }
};
class RandomBalanced : public TestCase {
protected:
    int N, M;
    bool notCycle;
    Edges edges;
public:
    RandomBalanced(int no, int N, int M, bool notCycle = false)
        : TestCase(no), N(N), M(M), notCycle(notCycle) {}
    void generate(RNG &rng) override {
        edges.clear();
        int E = 0;
        vector<int> v;
        vector<vector<bool> > g(N, vector<bool>(N, true));
        if(notCycle) ++ M;
        while(E < M) {
            int len, s;
            while(1) {
                len = rng(1, M - E);
                s = rng(N);
                v.clear();
                if(len == 1) {
                    if(!g[s][s]) {
                        continue;
                    } else {
                        v.push_back(s);
                        break;
                    }
                }
                v.push_back(s);
                bool ok = gen_dfs(1, rng, g, len, v);
                if(!ok)
                    continue;
                break;
            }
            for(int k = 0; k < (int)v.size() - 1; ++ k)
                edges.push_back({v[k], v[k + 1]});
            g[v.back()][v.front()] = false;
            edges.push_back({v.back(), v.front()});
            E = (int)edges.size();
        }
        if(notCycle) {
            -- M;
            edges.erase(edges.begin() + rng((int)edges.size()));
        }
    }
private:
    bool gen_dfs(int k, RNG &rng, vector<vector<bool> > &g, int len, vector<int> &v) {
        if(k == len)
            return true;
        for(int trys = 0; trys < 100; ++ trys) {
            int j = rng(N);
            if(!g[v.back()][j])
                continue;
            g[v.back()][j] = false;
            if(k == len - 1 && !g[j][v.front()]) {
                g[v.back()][j] = true;
                continue;
            }
            v.push_back(j);
            if(gen_dfs(k + 1, rng, g, len, v))
                return true;
            v.pop_back();
            g[v.back()][j] = true;
        }
        return false;
    }
protected:
    tuple<int,Edges> get() const override { return make_tuple(N, edges); }
};
class RandomBigCycle : public TestCase {
protected:
    int N;
    Edges edges;
public:
    RandomBigCycle(int no, int N) : TestCase(no), N(N) {}
    void generate(RNG &rng) override {
        vector<int> perm(N);
        iota(perm.begin(), perm.end(), 0);
        shuffle(perm.begin(), perm.end(), rng.re);
        edges.clear();
        rep(i, N)
            edges.push_back({perm[i], perm[(i + 1) % N]});
    }
protected:
    tuple<int,Edges> get() const override { return make_tuple(N, edges); }
};
class RandomBidirected : public TestCase {
protected:
    int N;
    int halfM;
    Edges edges;
public:
    RandomBidirected(int no, int N, int halfM) : TestCase(no), N(N), halfM(halfM) {}
    void generate(RNG &rng) override {
        edges.clear();
        vector<vector<bool> > g(N, vector<bool>(N, true));
        for(int i = 0; i < halfM; ++ i) {
            int a, b;
            do {
                a = rng(N), b = rng(N);
            }while(a == b || !g[a][b]);
            edges.push_back({a, b});
            edges.push_back({b, a});
            g[a][b] = g[b][a] = false;
        }
    }
protected:
    tuple<int,Edges> get() const override { return make_tuple(N, edges); }
};
int main() {
    vector<unique_ptr<TestCase> > testcases;
    auto A = [&](TestCase *t) { testcases.emplace_back(t); };
    const int MaxN = TestCase::MaxN, MaxM = TestCase::MaxM;
    int no = 0;
    //サンプル
    A(new Manual(no ++, 2, {{0, 1}}));
    A(new Manual(no ++, 2, {}));
    A(new Manual(no ++, 2, {{1, 1}, {0, 1}, {0, 0}, {1, 0}}));
    A(new Manual(no ++, 3, {{0, 2}, {2, 0}, {2, 1}}));
    A(new Manual(no ++, 3, {}));
    A(new Manual(no ++, 10, {}));
    //コーナーケースたち
    A(new Manual(no ++, 1, {}));
    A(new Manual(no ++, 1, {{0, 0}}));
    A(new Manual(no ++, 2, {{0, 1}, {1, 0}}));  //未連結
    A(new Manual(no ++, 2, {{0, 0}, {0, 1}, {1, 0}}));  //孤立点がある
    A(new Manual(no ++, 4, {{0, 1}, {0, 2}}));  //unbalanced
    //ランダム
    A(new RandomBidirected(no ++, 5, 5));
    A(new RandomBalanced(no ++, 5, 10));
    A(new RandomBalanced(no ++, 5, 10, true));
    A(new RandomBalanced(no ++, 5, 15));
    A(new RandomBalanced(no ++, 5, 15, true));
    A(new RandomBidirected(no ++, 20, 50));
    A(new RandomBalanced(no ++, 20, 100));
    A(new RandomBalanced(no ++, 100, 200));
    A(new RandomBalanced(no ++, 70, MaxM));
    A(new RandomBalanced(no ++, 100, MaxM));
    //最大ケース1
    A(new Manual(no ++, MaxN, {}));
    //最大ケース2
    A(new RandomBidirected(no ++, MaxN, MaxM / 2));
    A(new RandomBalanced(no ++, MaxN, MaxM));
    A(new RandomBalanced(no ++, MaxN, MaxM, true));
    //最大ケース3
    A(new RandomBigCycle(no ++, MaxN));
    
    if(testcase_dir[0] != 'C') abort();
    if(testcase_postfix[0] != '.') abort();
    if(0) {
        cerr << "compiling answer..." << endl;
        {
            stringstream ss;
            ss << "cd " << testcase_dir << " && g++ answer.cpp -std=c++11 -O2 -o answer.exe";
            string cmd = ss.str();
            if(system(cmd.c_str()) != 0) {
                cerr << "an error occurred" << endl;
                abort();
            }
        }
    }
    cerr << "deleting files..." << endl;
    {
        string cmd = string("del \"") + testcase_dir + "\\*" + testcase_postfix + "\"";
        if(system(cmd.c_str()) != 0) {
            cerr << "an error occurred" << endl;
            abort();
        }
    }
    int numcases = (int)testcases.size();
    cerr << "number of test cases: " << numcases << endl;
    set<string> testnames;
    for(int ii = 0; ii < numcases; ii ++) {
        TestCase &t = *testcases[ii];
        const string testname = t.getName();
        cerr << "generating " << testname << "..." << endl;
        if(testnames.count(testname)) {
            cerr << "testname duplication: " << testname << endl;
            abort();
        }
        testnames.count(testname);
        //テストケース生成で使う乱数の初期化
        RNG rng;
        rep(i, problem_name.size())
            rng.addSeed(problem_name[i]);
        rep(i, testname.size())
            rng.addSeed(testname[i]);
        rep(i, 100) rng();
        //テストケース生成
        t.generate(rng);
        //ファイルへ書き込む
        string filename = string(testcase_dir) + "\\" + testname + testcase_postfix;
        fstream f(filename.c_str(), ios_base::out);
        if(!f) {
            cerr << "can't open the file" << endl;
            abort();
        }
        t.write(f);
        f.close();
        if(1) {
            cerr << "run " << testname << "..." << endl;
            string run_cmd = string(testcase_dir) + "\\answer.exe < \"" + filename + "\"";
            TIMEIT("run total") {
                if(system(run_cmd.c_str()) != 0) {
                    cerr << "an error occurred" << endl;
                    abort();
                }
            }
        }
        cerr << endl;
    }
    return 0;
}
#endif
הההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההה
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
0