1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103 | namespace Matrix
{
static const double eps=1e-9;
template<typename T=ModInt<Mod>> struct matrix
{
vector<vector<T>> a;
inline matrix(int n=0):a(n) {}
inline matrix(const vector<vector<T>> &b):a(b) {}
inline matrix(const matrix &b) {a=b.a;}
inline vector<T> &operator [] (int i) {return a[i];}
inline matrix operator = (const matrix &b) {this->a=b.a; return *this;}
inline int size() {return a.size();}
inline void resize(int n) {a.resize(n,vector<T>(n));}
inline friend matrix operator + (matrix a,matrix b)
{
int N=a.size();
for(int i=0;i<N;++i) for(int j=0;j<N;++j) a[i][j]+=b[i][j];
return a;
}
inline friend matrix operator - (matrix a,matrix b)
{
int N=a.size();
for(int i=0;i<N;++i) for(int j=0;j<N;++j) a[i][j]-=b[i][j];
return a;
}
template<typename W> inline friend matrix operator * (matrix a,W b)
{
int N=a.size();
for(int i=0;i<N;++i) for(int j=0;j<N;++j) a[i][j]=(T)(a[i][j]*b);
return a;
}
inline friend matrix operator * (matrix a,matrix b)
{
int N=a.size(); matrix res; res.resize(N);
for(int i=0;i<N;++i) for(int j=0;j<N;++j) for(int k=0;k<N;++k) res[i][j]+=a[i][k]*b[k][j];
return res;
}
inline friend matrix operator ^ (matrix a,int b)
{
int N=a.size(); matrix res; res.resize(N);
for(int i=0;i<N;++i) for(int j=0;j<N;++j) res[i][i]=((T)1);
while(b) {if(b&1) res=res*a; a=a*a; b>>=1;}
return res;
}
inline pair<vector<T>,int> Gauss(vector<T> coef)
{
int N=a.size(),x=0,y=0;
for(;y<N;++y)
{
int now=x;
T mix=fabs(a[x][y]);
for(int i=x+1;i<N;++i) if(fabs(a[i][y])-mix>eps) now=i,mix=fabs(a[i][y]);
if(mix<eps) continue;
swap(coef[x],coef[now]),swap(a[x],a[now]);
for(int i=0;i<N;++i)
{
if(i==x) continue;
for(int j=y+1;j<N;++j) a[i][j]-=a[i][y]/a[x][y]*a[x][j];
coef[i]-=a[i][y]/a[x][y]*coef[x];
}
++x;
}
--x;
vector<T> ZERO(0),ans(N);
if(x<N-1)
{
for(int i=x+1;i<N;++i) if(fabs(coef[i])>eps) return make_pair(ZERO,-1);
return make_pair(ZERO,0);
}
for(int i=0;i<N;++i)
{
ans[i]=coef[i]/a[i][i];
if(fabs(ans[i])<eps) ans[i]=0;
}
return make_pair(ans,1);
}
inline pair<matrix,int> Inv()
{
int N=a.size();
matrix A(a),B,ZERO(0);
B.resize(N); for(int i=0;i<N;++i) B[i][i]=1;
for(int i=0;i<N;++i)
{
if(!A[i][i]) for(int j=i+1;j<N;++j) if(A[j][i]!=0) {swap(A[i],A[j]),swap(B[i],B[j]); break;}
if(!A[i][i]) return make_pair(ZERO,0);
T inv=1/A[i][i];
for(int j=0;j<N;++j) A[i][j]*=inv,B[i][j]*=inv;
for(int j=i+1;j<N;++j)
{
T ret=-A[j][i];
for(int k=0;k<N;++k) A[j][k]+=ret*A[i][k],B[j][k]+=ret*B[i][k];
}
}
for(int i=N-1;~i;--i) for(int j=i+1;j<N;++j)
{
T ret=-A[i][j];
for(int k=0;k<N;++k) A[i][k]+=ret*A[j][k],B[i][k]+=ret*B[j][k];
}
return make_pair(B,1);
}
};
}
|