#920337 (C++14) No.1907 DETERMINATION

• 問題: No.1907 DETERMINATION
• ユーザー: Noszály Áron
• 提出日時: 2023-10-16 01:17:37
• 言語: C++14 (gcc 13.3.0 + boost 1.87.0)
• 結果: AC
• 実行時間: 1,264 ms / 4,000 ms
• コード長: 7,917 bytes
• コンパイル時間: 3,424 ms
• コンパイル使用メモリ: 168,136 KB
• 実行使用メモリ: 9,088 KB
• 最終ジャッジ日時: 2024-09-16 22:09:28
• 合計ジャッジ時間: 39,575 ms
• ジャッジサーバーID （参考情報）: judge6 / judge1

ファイルパターン	結果
sample	AC * 4
other	AC * 63

ソースコード

#include<iostream>
#include<vector>
#include<map>
#include<set>
#include<cassert>
#include<unordered_map>
#include<unordered_set>
#include<functional>
#include<queue>
#include<stack>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<sstream>
#include<iomanip>
#include<cstdio>
#include<cstdlib>
#include<numeric>
#include<random>
#include<chrono>
#include<bitset>
#include<ext/pb_ds/assoc_container.hpp>
#include<ext/pb_ds/tree_policy.hpp>
using namespace std;
using namespace __gnu_pbds;

#define all(x) (x).begin(), (x).end()
#define pb push_back
#define eb emplace_back
#define xx first
#define yy second
#define sz(x) (int)(x).size()
#define gc getchar
#define IO ios_base::sync_with_stdio(false);cin.tie(0);cout.tie(0)
#define mp make_pair
#define ins insert

#ifndef ONLINE_JUDGE
#  define LOG(x) (cerr << #x << " = " << (x) << endl)
#else
#  define LOG(x) ((void)0)
#endif

using ll = long long;
using ull = unsigned long long ;
using ld = long double ;
using str = string;
using ordered_set=tree<pair<int,int>, null_type, less<pair<int,int>>, rb_tree_tag, tree_order_statistics_node_update>;

const double PI=acos(-1);
const ll INF = 1LL<<62;
const ll MINF = -(1LL<<62);

template<typename T> T getint() {
    T val=0;
    char c;
    
    bool neg=false;
    while((c=gc()) && !(c>='0' && c<='9')) {
        neg|=c=='-';
    }

    do {
        val=(val*10)+c-'0';
    } while((c=gc()) && (c>='0' && c<='9'));

    return val*(neg?-1:1);
}

//mt19937 rng(chrono::steady_clock::now().time_since_epoch().count()); uniform_int_distribution<int>(0, n-1)(rng)

const ll mod=998244353;

template<int m>
struct modular {
    int64_t r;
    modular() : r(0) {}
    modular(int64_t rr) : r(rr) {if(abs(r) >= m) r %= m; if(r < 0) r += m;}
    modular inv() const {return bpow(*this, m - 2);}
    modular operator * (const modular &t) const {return (r * t.r) % m;}
    modular operator / (const modular &t) const {return *this * t.inv();}
    modular operator += (const modular &t) {r += t.r; if(r >= m) r -= m; return *this;}
    modular operator -= (const modular &t) {r -= t.r; if(r < 0) r += m; return *this;}
    modular operator + (const modular &t) const {return modular(*this) += t;}
    modular operator - (const modular &t) const {return modular(*this) -= t;}
    modular operator *= (const modular &t) {return *this = *this * t;}
    modular operator /= (const modular &t) {return *this = *this / t;}
    
    bool operator == (const modular &t) const {return r == t.r;}
    bool operator != (const modular &t) const {return r != t.r;}
    
    operator int64_t() const {return r;}

     
    int64_t bpow(int64_t a, int64_t b) const {
		int64_t res=1;
		while(b) {
			if(b&1) res=(res*a)%m;
			a=(a*a)%m;
			b/=2;
		}
		return res;
	}
};

using Mint = modular<mod>;
using matrix = vector<vector<Mint>>;

ostream& operator<<(ostream& os, matrix M) {
    for(auto& i:M) {
        for(auto& j:i) {
            os<<j<<" ";
        }
        os<<"\n";
    }
    return os;
}
 

void swap_row(matrix& M, int i, int j) {
    assert(i!=j);
    for(int k=0;k<sz(M);++k) {
        swap(M[i][k], M[j][k]);
    }
} 

void swap_col(matrix& M, int i, int j) {
    assert(i!=j);
    for(int k=0;k<sz(M);++k) {
        swap(M[k][i], M[k][j]);
    }
}

void row_op(matrix& M, int i, int j, Mint c) {
    for(int k=0;k<sz(M);++k) {
        M[i][k]+=M[j][k]*c;
    }
}

void col_op(matrix& M, int i, int j, Mint c) {
    for(int k=0;k<sz(M);++k) {
        M[k][i]+=M[k][j]*c;
    }
}

void make_upper_heisenberg(matrix& M) {
    const int n=M.size();
    for(int col=0;col<n-2;++col) {
        int piv=-1;
        for(int row=col+1;row<n;++row) {
            if(M[row][col]!=Mint(0)) {
                piv=row;
                break;
            }
        }
        
        if(piv<0) continue ;
        if(col+1!=piv) {
            swap_row(M, col+1, piv);
            swap_col(M, col+1, piv);
        }
    
        auto inv=M[col+1][col].inv();
        for(int i=col+2;i<n;++i) {
            auto coeff=-M[i][col]*inv;
            row_op(M, i, col+1, coeff); //unitér legyen
            col_op(M, col+1, i, -coeff);
        }
    }
}

vector<Mint> char_poly(matrix& M) {
    const int n=sz(M);
    
    make_upper_heisenberg(M);
    
    vector<vector<Mint>> p(n+1);
    p[0]={1};

    //det(Ix-M)
    //kifejtés miatt: P[i+1]=1*(-P[i]*M(i,i)+xP[i])-\sum j=0..i-1 (P[j]*M(j,i)* prod M(j+1..i,j..i-1))
    for(int i=0;i<n;++i) { 
        p[i+1].assign(i+2, 0);
        for(int j=0;j<sz(p[i]);++j) p[i+1][j]-=p[i][j]*M[i][i];
        for(int j=0;j<sz(p[i]);++j) p[i+1][j+1]+=p[i][j];
        
        Mint mul=1;
        for(int j=i-1;j>=0;j--) {
            mul*=M[j+1][j];
            Mint hb=-M[j][i]*mul;
            for(int k=0;k<sz(p[j]);++k) p[i+1][k]+=hb*p[j][k];
        }
    }
    
    return p[n];
}

//det(M0+M1x) kiszámítása
vector<Mint> det_linear(matrix M0, matrix M1) {
    //sor és oszlop műveletekkel M1-et I-vé transzformáljuk det(M0+M1x)=1/(det(A)*det(B))*det(Ix+AM0B) alakra hozzuk => -AM0B karakterisztikus polinomja a megoldás
    //ha nem reguláris M1 "M0-ból hozunk át oszlopokat"
    const int N=M0.size();
    Mint invdetAB=1;
    int powx=0;
    for(int p=0;p<N;++p) {
        int piv=-1;
        for(int row=p;row<N;++row) {
            if(M1[row][p]!=Mint(0)) {
                piv=row;
                break ;
            }
        }
        if(piv<0) {
            powx++;
            if(powx>N) return vector<Mint>(N+1);
            
            for(int i=0;i<p;++i) {
                Mint val=-M1[i][p];
                col_op(M0, p, i, val);
                col_op(M1, p, i, val); //felesleges valójában
            }
            for(int i=0;i<N;++i) swap(M0[i][p], M1[i][p]);
            p--;
            
            continue ;
        }
        
        if(p!=piv) {
            swap_row(M0, piv, p);
            swap_row(M1, piv, p);
            invdetAB*=-1;
        }
        
        Mint v=M1[p][p];
        Mint vinv=v.inv();
        
        invdetAB*=v;
        for(int col=0;col<N;++col) {
            M0[p][col]*=vinv;
            M1[p][col]*=vinv;
        }
        
        for(int row=0;row<N;++row) {
            if(row!=p) {
                Mint val=-M1[row][p];
                row_op(M1, row, p, val);
                row_op(M0, row, p, val);
            }
        }
    }
    
    for(auto& i:M0) 
        for(auto& j:i) 
            j=-j;
    
    auto res=char_poly(M0);
    for(auto& i:res) i*=invdetAB;
    
    res.erase(res.begin(), res.begin()+powx);
    res.resize(N+1);
    
    return res;
}
 

void test_make_upper_hessenberg() {
    matrix m;
    
    m={{1,1},{1,1}};
    make_upper_heisenberg(m);
    assert((m==matrix{{1,1},{1,1}}));
    
    m={{1,2,3},{3,4,5},{0,9,9}};
    make_upper_heisenberg(m);
    assert((m==matrix{{{1,2,3},{3,4,5},{0,9,9}}}));
    
    m={{1,2,3},{4,5,6},{7,8,9}};
    make_upper_heisenberg(m);
    
    m={{1,2},{3,4}};
    assert((char_poly(m)==vector<Mint>{998244351,998244348,1}));
    //~ for(auto i:char_poly(m)) cerr<<i<<" ";cerr<<"\n";
}


int main() {
    IO;
    //~ test_make_upper_hessenberg();
    //~ int n;
    //~ cin>>n;
    //~ matrix m(n);
    //~ for(int i=0;i<n;++i) {
        //~ m[i].resize(n);
        //~ for(int j=0;j<n;++j) {
            //~ ll x;
            //~ cin>>x;
            //~ m[i][j]=x;
        //~ }
    //~ }
    
    //~ for(auto i:char_poly(m)) {
        //~ cout<<i<<" ";
    //~ }cout<<"\n";
    
    int n;
    cin>>n;
    matrix M0(n), M1(n);
    for(int i=0;i<n;++i) {
        M0[i].resize(n);
        for(int j=0;j<n;++j) {
            ll x;
            cin>>x;
            M0[i][j]=x;
        }
    }
    
    for(int i=0;i<n;++i) {
        M1[i].resize(n);
        for(int j=0;j<n;++j) {
            ll x;
            cin>>x;
            M1[i][j]=x;
        }
    }
    
    for(auto i:det_linear(M0, M1)) cout<<i<<"\n";
    
    
    return 0;
}