と定義すると、$(V(G), \leq)$ は半順序集合である。$(V(G), \leq)$ に対して、共通部分のない鎖に分解したときの最小サイズは最小パス被覆のサイズを意味する。
#line 1 "include/emthrm/graph/flow/matching/maximum_matching.hpp"
#include <random>
#include <vector>
#line 1 "include/emthrm/math/matrix/gauss_jordan.hpp"
#include <utility>
#line 1 "include/emthrm/math/matrix/matrix.hpp"
#line 5 "include/emthrm/math/matrix/matrix.hpp"
namespace emthrm {
template < typename T >
struct Matrix {
explicit Matrix ( const int m , const int n , const T def = 0 )
: data ( m , std :: vector < T > ( n , def )) {}
int nrow () const { return data . size (); }
int ncol () const { return data . empty () ? 0 : data . front (). size (); }
Matrix pow ( long long exponent ) const {
const int n = nrow ();
Matrix < T > res ( n , n , 0 ), tmp = * this ;
for ( int i = 0 ; i < n ; ++ i ) {
res [ i ][ i ] = 1 ;
}
for (; exponent > 0 ; exponent >>= 1 ) {
if ( exponent & 1 ) res *= tmp ;
tmp *= tmp ;
}
return res ;
}
inline const std :: vector < T >& operator []( const int i ) const { return data [ i ]; }
inline std :: vector < T >& operator []( const int i ) { return data [ i ]; }
Matrix & operator = ( const Matrix & x ) = default ;
Matrix & operator += ( const Matrix & x ) {
const int m = nrow (), n = ncol ();
for ( int i = 0 ; i < m ; ++ i ) {
for ( int j = 0 ; j < n ; ++ j ) {
data [ i ][ j ] += x [ i ][ j ];
}
}
return * this ;
}
Matrix & operator -= ( const Matrix & x ) {
const int m = nrow (), n = ncol ();
for ( int i = 0 ; i < m ; ++ i ) {
for ( int j = 0 ; j < n ; ++ j ) {
data [ i ][ j ] -= x [ i ][ j ];
}
}
return * this ;
}
Matrix & operator *= ( const Matrix & x ) {
const int m = nrow (), l = ncol (), n = x . ncol ();
std :: vector < std :: vector < T >> res ( m , std :: vector < T > ( n , 0 ));
for ( int i = 0 ; i < m ; ++ i ) {
for ( int k = 0 ; k < l ; ++ k ) {
for ( int j = 0 ; j < n ; ++ j ) {
res [ i ][ j ] += data [ i ][ k ] * x [ k ][ j ];
}
}
}
data . swap ( res );
return * this ;
}
Matrix operator + ( const Matrix & x ) const { return Matrix ( * this ) += x ; }
Matrix operator - ( const Matrix & x ) const { return Matrix ( * this ) -= x ; }
Matrix operator * ( const Matrix & x ) const { return Matrix ( * this ) *= x ; }
private:
std :: vector < std :: vector < T >> data ;
};
} // namespace emthrm
#line 7 "include/emthrm/math/matrix/gauss_jordan.hpp"
namespace emthrm {
template < bool IS_EXTENDED = false , typename T >
int gauss_jordan ( Matrix < T >* a , const T eps = 1e-8 ) {
const int m = a -> nrow (), n = a -> ncol ();
int rank = 0 ;
for ( int col = 0 ; col < ( IS_EXTENDED ? n - 1 : n ); ++ col ) {
int pivot = - 1 ;
T mx = eps ;
for ( int row = rank ; row < m ; ++ row ) {
const T abs = (( * a )[ row ][ col ] < 0 ? - ( * a )[ row ][ col ] : ( * a )[ row ][ col ]);
if ( abs > mx ) {
pivot = row ;
mx = abs ;
}
}
if ( pivot == - 1 ) continue ;
std :: swap (( * a )[ rank ], ( * a )[ pivot ]);
T tmp = ( * a )[ rank ][ col ];
for ( int col2 = 0 ; col2 < n ; ++ col2 ) {
( * a )[ rank ][ col2 ] /= tmp ;
}
for ( int row = 0 ; row < m ; ++ row ) {
if ( row != rank &&
(( * a )[ row ][ col ] < 0 ? - ( * a )[ row ][ col ] : ( * a )[ row ][ col ]) > eps ) {
tmp = ( * a )[ row ][ col ];
for ( int col2 = 0 ; col2 < n ; ++ col2 ) {
( * a )[ row ][ col2 ] -= ( * a )[ rank ][ col2 ] * tmp ;
}
}
}
++ rank ;
}
return rank ;
}
} // namespace emthrm
#line 1 "include/emthrm/math/modint.hpp"
#ifndef ARBITRARY_MODINT
# include <cassert>
#endif
#include <compare>
#include <iostream>
// #include <numeric>
#line 12 "include/emthrm/math/modint.hpp"
namespace emthrm {
#ifndef ARBITRARY_MODINT
template < unsigned int M >
struct MInt {
unsigned int v ;
constexpr MInt () : v ( 0 ) {}
constexpr MInt ( const long long x ) : v ( x >= 0 ? x % M : x % M + M ) {}
static constexpr MInt raw ( const int x ) {
MInt x_ ;
x_ . v = x ;
return x_ ;
}
static constexpr int get_mod () { return M ; }
static constexpr void set_mod ( const int divisor ) {
assert ( std :: cmp_equal ( divisor , M ));
}
static void init ( const int x ) {
inv < true > ( x );
fact ( x );
fact_inv ( x );
}
template < bool MEMOIZES = false >
static MInt inv ( const int n ) {
// assert(0 <= n && n < M && std::gcd(n, M) == 1);
static std :: vector < MInt > inverse { 0 , 1 };
const int prev = inverse . size ();
if ( n < prev ) return inverse [ n ];
if constexpr ( MEMOIZES ) {
// "n!" and "M" must be disjoint.
inverse . resize ( n + 1 );
for ( int i = prev ; i <= n ; ++ i ) {
inverse [ i ] = - inverse [ M % i ] * raw ( M / i );
}
return inverse [ n ];
}
int u = 1 , v = 0 ;
for ( unsigned int a = n , b = M ; b ;) {
const unsigned int q = a / b ;
std :: swap ( a -= q * b , b );
std :: swap ( u -= q * v , v );
}
return u ;
}
static MInt fact ( const int n ) {
static std :: vector < MInt > factorial { 1 };
if ( const int prev = factorial . size (); n >= prev ) {
factorial . resize ( n + 1 );
for ( int i = prev ; i <= n ; ++ i ) {
factorial [ i ] = factorial [ i - 1 ] * i ;
}
}
return factorial [ n ];
}
static MInt fact_inv ( const int n ) {
static std :: vector < MInt > f_inv { 1 };
if ( const int prev = f_inv . size (); n >= prev ) {
f_inv . resize ( n + 1 );
f_inv [ n ] = inv ( fact ( n ). v );
for ( int i = n ; i > prev ; -- i ) {
f_inv [ i - 1 ] = f_inv [ i ] * i ;
}
}
return f_inv [ n ];
}
static MInt nCk ( const int n , const int k ) {
if ( n < 0 || n < k || k < 0 ) [[ unlikely ]] return MInt ();
return fact ( n ) * ( n - k < k ? fact_inv ( k ) * fact_inv ( n - k ) :
fact_inv ( n - k ) * fact_inv ( k ));
}
static MInt nPk ( const int n , const int k ) {
return n < 0 || n < k || k < 0 ? MInt () : fact ( n ) * fact_inv ( n - k );
}
static MInt nHk ( const int n , const int k ) {
return n < 0 || k < 0 ? MInt () : ( k == 0 ? 1 : nCk ( n + k - 1 , k ));
}
static MInt large_nCk ( long long n , const int k ) {
if ( n < 0 || n < k || k < 0 ) [[ unlikely ]] return MInt ();
inv < true > ( k );
MInt res = 1 ;
for ( int i = 1 ; i <= k ; ++ i ) {
res *= inv ( i ) * n -- ;
}
return res ;
}
constexpr MInt pow ( long long exponent ) const {
MInt res = 1 , tmp = * this ;
for (; exponent > 0 ; exponent >>= 1 ) {
if ( exponent & 1 ) res *= tmp ;
tmp *= tmp ;
}
return res ;
}
constexpr MInt & operator += ( const MInt & x ) {
if (( v += x . v ) >= M ) v -= M ;
return * this ;
}
constexpr MInt & operator -= ( const MInt & x ) {
if (( v += M - x . v ) >= M ) v -= M ;
return * this ;
}
constexpr MInt & operator *= ( const MInt & x ) {
v = ( unsigned long long ){ v } * x . v % M ;
return * this ;
}
MInt & operator /= ( const MInt & x ) { return * this *= inv ( x . v ); }
constexpr auto operator <=> ( const MInt & x ) const = default ;
constexpr MInt & operator ++ () {
if ( ++ v == M ) [[ unlikely ]] v = 0 ;
return * this ;
}
constexpr MInt operator ++ ( int ) {
const MInt res = * this ;
++* this ;
return res ;
}
constexpr MInt & operator -- () {
v = ( v == 0 ? M - 1 : v - 1 );
return * this ;
}
constexpr MInt operator -- ( int ) {
const MInt res = * this ;
--* this ;
return res ;
}
constexpr MInt operator + () const { return * this ; }
constexpr MInt operator - () const { return raw ( v ? M - v : 0 ); }
constexpr MInt operator + ( const MInt & x ) const { return MInt ( * this ) += x ; }
constexpr MInt operator - ( const MInt & x ) const { return MInt ( * this ) -= x ; }
constexpr MInt operator * ( const MInt & x ) const { return MInt ( * this ) *= x ; }
MInt operator / ( const MInt & x ) const { return MInt ( * this ) /= x ; }
friend std :: ostream & operator << ( std :: ostream & os , const MInt & x ) {
return os << x . v ;
}
friend std :: istream & operator >> ( std :: istream & is , MInt & x ) {
long long v ;
is >> v ;
x = MInt ( v );
return is ;
}
};
#else // ARBITRARY_MODINT
template < int ID >
struct MInt {
unsigned int v ;
constexpr MInt () : v ( 0 ) {}
MInt ( const long long x ) : v ( x >= 0 ? x % mod () : x % mod () + mod ()) {}
static constexpr MInt raw ( const int x ) {
MInt x_ ;
x_ . v = x ;
return x_ ;
}
static int get_mod () { return mod (); }
static void set_mod ( const unsigned int divisor ) { mod () = divisor ; }
static void init ( const int x ) {
inv < true > ( x );
fact ( x );
fact_inv ( x );
}
template < bool MEMOIZES = false >
static MInt inv ( const int n ) {
// assert(0 <= n && n < mod() && std::gcd(x, mod()) == 1);
static std :: vector < MInt > inverse { 0 , 1 };
const int prev = inverse . size ();
if ( n < prev ) return inverse [ n ];
if constexpr ( MEMOIZES ) {
// "n!" and "M" must be disjoint.
inverse . resize ( n + 1 );
for ( int i = prev ; i <= n ; ++ i ) {
inverse [ i ] = - inverse [ mod () % i ] * raw ( mod () / i );
}
return inverse [ n ];
}
int u = 1 , v = 0 ;
for ( unsigned int a = n , b = mod (); b ;) {
const unsigned int q = a / b ;
std :: swap ( a -= q * b , b );
std :: swap ( u -= q * v , v );
}
return u ;
}
static MInt fact ( const int n ) {
static std :: vector < MInt > factorial { 1 };
if ( const int prev = factorial . size (); n >= prev ) {
factorial . resize ( n + 1 );
for ( int i = prev ; i <= n ; ++ i ) {
factorial [ i ] = factorial [ i - 1 ] * i ;
}
}
return factorial [ n ];
}
static MInt fact_inv ( const int n ) {
static std :: vector < MInt > f_inv { 1 };
if ( const int prev = f_inv . size (); n >= prev ) {
f_inv . resize ( n + 1 );
f_inv [ n ] = inv ( fact ( n ). v );
for ( int i = n ; i > prev ; -- i ) {
f_inv [ i - 1 ] = f_inv [ i ] * i ;
}
}
return f_inv [ n ];
}
static MInt nCk ( const int n , const int k ) {
if ( n < 0 || n < k || k < 0 ) [[ unlikely ]] return MInt ();
return fact ( n ) * ( n - k < k ? fact_inv ( k ) * fact_inv ( n - k ) :
fact_inv ( n - k ) * fact_inv ( k ));
}
static MInt nPk ( const int n , const int k ) {
return n < 0 || n < k || k < 0 ? MInt () : fact ( n ) * fact_inv ( n - k );
}
static MInt nHk ( const int n , const int k ) {
return n < 0 || k < 0 ? MInt () : ( k == 0 ? 1 : nCk ( n + k - 1 , k ));
}
static MInt large_nCk ( long long n , const int k ) {
if ( n < 0 || n < k || k < 0 ) [[ unlikely ]] return MInt ();
inv < true > ( k );
MInt res = 1 ;
for ( int i = 1 ; i <= k ; ++ i ) {
res *= inv ( i ) * n -- ;
}
return res ;
}
MInt pow ( long long exponent ) const {
MInt res = 1 , tmp = * this ;
for (; exponent > 0 ; exponent >>= 1 ) {
if ( exponent & 1 ) res *= tmp ;
tmp *= tmp ;
}
return res ;
}
MInt & operator += ( const MInt & x ) {
if (( v += x . v ) >= mod ()) v -= mod ();
return * this ;
}
MInt & operator -= ( const MInt & x ) {
if (( v += mod () - x . v ) >= mod ()) v -= mod ();
return * this ;
}
MInt & operator *= ( const MInt & x ) {
v = ( unsigned long long ){ v } * x . v % mod ();
return * this ;
}
MInt & operator /= ( const MInt & x ) { return * this *= inv ( x . v ); }
auto operator <=> ( const MInt & x ) const = default ;
MInt & operator ++ () {
if ( ++ v == mod ()) [[ unlikely ]] v = 0 ;
return * this ;
}
MInt operator ++ ( int ) {
const MInt res = * this ;
++* this ;
return res ;
}
MInt & operator -- () {
v = ( v == 0 ? mod () - 1 : v - 1 );
return * this ;
}
MInt operator -- ( int ) {
const MInt res = * this ;
--* this ;
return res ;
}
MInt operator + () const { return * this ; }
MInt operator - () const { return raw ( v ? mod () - v : 0 ); }
MInt operator + ( const MInt & x ) const { return MInt ( * this ) += x ; }
MInt operator - ( const MInt & x ) const { return MInt ( * this ) -= x ; }
MInt operator * ( const MInt & x ) const { return MInt ( * this ) *= x ; }
MInt operator / ( const MInt & x ) const { return MInt ( * this ) /= x ; }
friend std :: ostream & operator << ( std :: ostream & os , const MInt & x ) {
return os << x . v ;
}
friend std :: istream & operator >> ( std :: istream & is , MInt & x ) {
long long v ;
is >> v ;
x = MInt ( v );
return is ;
}
private:
static unsigned int & mod () {
static unsigned int divisor = 0 ;
return divisor ;
}
};
#endif // ARBITRARY_MODINT
} // namespace emthrm
#line 10 "include/emthrm/graph/flow/matching/maximum_matching.hpp"
namespace emthrm {
int maximum_matching ( const std :: vector < std :: vector < int >>& graph ) {
constexpr unsigned int P = 1000000007 ;
using ModInt = MInt < P > ;
ModInt :: set_mod ( P );
static std :: mt19937_64 engine ( std :: random_device {} ());
static std :: uniform_int_distribution <> dist ( 1 , P - 1 );
const int n = graph . size ();
Matrix < ModInt > tutte_matrix ( n , n , 0 );
for ( int i = 0 ; i < n ; ++ i ) {
for ( const int j : graph [ i ]) {
if ( j > i ) {
const ModInt x = ModInt :: raw ( dist ( engine ));
tutte_matrix [ i ][ j ] = x ;
tutte_matrix [ j ][ i ] = - x ;
}
}
}
return gauss_jordan ( & tutte_matrix , ModInt ( 0 )) / 2 ;
}
} // namespace emthrm