$\det(A+Bx)$
(Math/detaplusbx.hpp)

View this file on GitHub
Last update: 2024-04-26 03:18:17+09:00
Include: #include "Math/detaplusbx.hpp"

Depends on

Code

#pragma once
#include "Math/matrix.hpp"
#include "Math/charpoly.hpp"
#include "Utility/random.hpp"

template<typename T>vector<T> detApBx(Matrix<T> a,Matrix<T> b){
    Random gen;
    int n=a.h;
    vector<T> f(n+1);
    T ran=gen.get();
    rep(i,0,n)rep(j,0,n)a[i][j]-=b[i][j]*ran;
    auto ainv=a.inv();
    if(a.det==0)return f;
    b*=ainv;
    rep(i,0,n)rep(j,0,n)b[i][j]=-b[i][j];
    f=CharPoly(b);
    reverse(ALL(f));
    for(auto& x:f)x*=a.det;
    vector C(n+1,vector<T>(n+1));
    vector<T> pw(n+1,1);
    rep(i,0,n+1){
        if(i)pw[i]=pw[i-1]*ran;
        C[i][0]=C[i][i]=1;
        rep(j,1,i)C[i][j]=C[i-1][j-1]+C[i-1][j];
    }
    vector<T> ret(n+1);
    rep(i,0,n+1)rep(j,0,i+1)ret[j]+=f[i]*C[i][j]*pw[i-j];
    return ret;
}

/**
 * @brief $\det(A+Bx)$
*/

#line 2 "Math/matrix.hpp"

template<class T>struct Matrix{
    int h,w; vector<vector<T>> val; T det;
    Matrix(){}
    Matrix(int n):h(n),w(n),val(vector<vector<T>>(n,vector<T>(n))){}
    Matrix(int n,int m):h(n),w(m),val(vector<vector<T>>(n,vector<T>(m))){}
    vector<T>& operator[](const int i){return val[i];}
    Matrix& operator+=(const Matrix& m){
        assert(h==m.h and w==m.w);
        rep(i,0,h)rep(j,0,w)val[i][j]+=m.val[i][j];
        return *this;
    }
    Matrix& operator-=(const Matrix& m){
        assert(h==m.h and w==m.w);
        rep(i,0,h)rep(j,0,w)val[i][j]-=m.val[i][j];
        return *this;
    }
    Matrix& operator*=(const Matrix& m){
        assert(w==m.h);
        Matrix<T> res(h,m.w);
        rep(i,0,h)rep(j,0,m.w)rep(k,0,w)res.val[i][j]+=val[i][k]*m.val[k][j];
        *this=res; return *this;
    }
    Matrix operator+(const Matrix& m)const{return Matrix(*this)+=m;}
    Matrix operator-(const Matrix& m)const{return Matrix(*this)-=m;}
    Matrix operator*(const Matrix& m)const{return Matrix(*this)*=m;}
    Matrix pow(ll k){
        Matrix<T> res(h,h),c=*this; rep(i,0,h)res.val[i][i]=1;
        while(k){if(k&1)res*=c; c*=c; k>>=1;} return res;
    }
    vector<int> gauss(int c=-1){
        if(val.empty())return {};
        if(c==-1)c=w;
        int cur=0; vector<int> res; det=1;
        rep(i,0,c){
            if(cur==h)break;
            rep(j,cur,h)if(val[j][i]!=0){
                swap(val[cur],val[j]);
                if(cur!=j)det*=-1;
                break;
            }
            det*=val[cur][i];
            if(val[cur][i]==0)continue;
            rep(j,0,h)if(j!=cur){
                T z=val[j][i]/val[cur][i];
                rep(k,i,w)val[j][k]-=val[cur][k]*z;
            }
            res.push_back(i);
            cur++;
        }
        return res;
    }
    Matrix inv(){
        assert(h==w);
        Matrix base(h,h*2),res(h,h);
        rep(i,0,h)rep(j,0,h)base[i][j]=val[i][j];
        rep(i,0,h)base[i][h+i]=1;
        base.gauss(h);
        det=base.det;
        rep(i,0,h)rep(j,0,h)res[i][j]=base[i][h+j]/base[i][i];
        return res;
    }
    bool operator==(const Matrix& m){
        assert(h==m.h and w==m.w);
        rep(i,0,h)rep(j,0,w)if(val[i][j]!=m.val[i][j])return false;
        return true;
    }
    bool operator!=(const Matrix& m){
        assert(h==m.h and w==m.w);
        rep(i,0,h)rep(j,0,w)if(val[i][j]==m.val[i][j])return false;
        return true;
    }
    friend istream& operator>>(istream& is,Matrix& m){
        rep(i,0,m.h)rep(j,0,m.w)is>>m[i][j];
        return is;
    }
    friend ostream& operator<<(ostream& os,Matrix& m){
        rep(i,0,m.h){
            rep(j,0,m.w)os<<m[i][j]<<(j==m.w-1 and i!=m.h-1?'\n':' ');
        }
        return os;
    }
};

/**
 * @brief Matrix
 */
#line 3 "Math/charpoly.hpp"

template<typename T>vector<T> CharPoly(Matrix<T> a){
    // to Hessenberg
    //reference:http://www.oishi.info.waseda.ac.jp/~samukawa/eigvieta.pdf
    int n=a.h;
    rep(s,0,n-2){
        rep(j,s+1,n)if(a[j][s]!=0){
            swap(a[j],a[s+1]);
            rep(i,0,n)swap(a[i][j],a[i][s+1]);
            break;
        }
        if(a[s+1][s]==0)continue;
        T X=T(1)/a[s+1][s];
        rep(i,s+2,n){
            T base=a[i][s]*X;
            rep(j,0,n)a[i][j]-=a[s+1][j]*base;
            rep(j,0,n)a[j][s+1]+=a[j][i]*base;
        }
    }
    vector dp(n+1,vector<T>(n+1));
    dp[0][0]=1;
    rep(i,0,n){
        rep(k,0,i+1){
            dp[i+1][k+1]+=dp[i][k];
            dp[i+1][k]-=dp[i][k]*a[i][i];
        }
        T prod=1;
        for(int j=i-1;j>=0;j--){
            prod*=a[j+1][j];
            T base=prod*a[j][i];
            rep(k,0,i+1)dp[i+1][k]-=dp[j][k]*base;
        }
    }
    return dp[n];
}

/**
 * @brief Characteristic Polynomial
*/
#line 2 "Utility/random.hpp"

namespace Random {
mt19937_64 randgen(chrono::steady_clock::now().time_since_epoch().count());
using u64 = unsigned long long;
u64 get() {
    return randgen();
}
template <typename T> T get(T L) { // [0,L]

    return get() % (L + 1);
}
template <typename T> T get(T L, T R) { // [L,R]

    return get(R - L) + L;
}
double uniform() {
    return double(get(1000000000)) / 1000000000;
}
string str(int n) {
    string ret;
    rep(i, 0, n) ret += get('a', 'z');
    return ret;
}
template <typename Iter> void shuffle(Iter first, Iter last) {
    if (first == last)
        return;
    int len = 1;
    for (auto it = first + 1; it != last; it++) {
        len++;
        int j = get(0, len - 1);
        if (j != len - 1)
            iter_swap(it, first + j);
    }
}
template <typename T> vector<T> select(int n, T L, T R) { // [L,R]

    if (n * 2 >= R - L + 1) {
        vector<T> ret(R - L + 1);
        iota(ALL(ret), L);
        shuffle(ALL(ret));
        ret.resize(n);
        return ret;
    } else {
        unordered_set<T> used;
        vector<T> ret;
        while (SZ(used) < n) {
            T x = get(L, R);
            if (!used.count(x)) {
                used.insert(x);
                ret.push_back(x);
            }
        }
        return ret;
    }
}

void relabel(int n, vector<pair<int, int>> &es) {
    shuffle(ALL(es));
    vector<int> ord(n);
    iota(ALL(ord), 0);
    shuffle(ALL(ord));
    for (auto &[u, v] : es)
        u = ord[u], v = ord[v];
}
template <bool directed, bool simple> vector<pair<int, int>> genGraph(int n) {
    vector<pair<int, int>> cand, es;
    rep(u, 0, n) rep(v, 0, n) {
        if (simple and u == v)
            continue;
        if (!directed and u > v)
            continue;
        cand.push_back({u, v});
    }
    int m = get(SZ(cand));
    vector<int> ord;
    if (simple)
        ord = select(m, 0, SZ(cand) - 1);
    else {
        rep(_, 0, m) ord.push_back(get(SZ(cand) - 1));
    }
    for (auto &i : ord)
        es.push_back(cand[i]);
    relabel(n, es);
    return es;
}
vector<pair<int, int>> genTree(int n) {
    vector<pair<int, int>> es;
    rep(i, 1, n) es.push_back({get(i - 1), i});
    relabel(n, es);
    return es;
}
}; // namespace Random


/**
 * @brief Random
 */
#line 5 "Math/detaplusbx.hpp"

template<typename T>vector<T> detApBx(Matrix<T> a,Matrix<T> b){
    Random gen;
    int n=a.h;
    vector<T> f(n+1);
    T ran=gen.get();
    rep(i,0,n)rep(j,0,n)a[i][j]-=b[i][j]*ran;
    auto ainv=a.inv();
    if(a.det==0)return f;
    b*=ainv;
    rep(i,0,n)rep(j,0,n)b[i][j]=-b[i][j];
    f=CharPoly(b);
    reverse(ALL(f));
    for(auto& x:f)x*=a.det;
    vector C(n+1,vector<T>(n+1));
    vector<T> pw(n+1,1);
    rep(i,0,n+1){
        if(i)pw[i]=pw[i-1]*ran;
        C[i][0]=C[i][i]=1;
        rep(j,1,i)C[i][j]=C[i-1][j-1]+C[i-1][j];
    }
    vector<T> ret(n+1);
    rep(i,0,n+1)rep(j,0,i+1)ret[j]+=f[i]*C[i][j]*pw[i-j];
    return ret;
}

/**
 * @brief $\det(A+Bx)$
*/

$\det(A+Bx)$ (Math/detaplusbx.hpp)

Depends on

Code

$\det(A+Bx)$
(Math/detaplusbx.hpp)