템플릿과 역행렬 3. 2^n차가 아닌 행렬

5차 정사각행렬 같은 경우는 어떻게 역행렬을 구할수 있을까??? 간단한 방법은 5차를 억지로 8차로 늘려버리는 것이다.

남는 자리에 0 (대각선은 1)을 채우면 된다. 그리고 얘의 역행렬을 구하고 왼쪽위에서 5x5개만 찝어내면 그게 원래의 역행렬이다. 그러면 이제 할일은 2^n꼴이 아닌 matrix를 2^n꼴로 확장시키는 것.

template<int n> class roundbit { public: enum{value=roundbit<(n+1)/2>::value*2}; }; template<> class roundbit<2> { public: enum{value=2}; };
2^m꼴이 아닌 수 n을 n이상의 2^m꼴로 바꿔주기 위해서 만든 클래스이다. roundbit<7>::value하면 8이 되고, roundbit<33>::value하면 64가 되겠다.

template<int dim> class matrix { protected: enum{halfdim=roundbit<dim>::value/2, round=roundbit<dim>::value,}; matrix<halfdim> a,b,c,d; public: static matrix<dim> zero; matrix<dim>() { if(round>dim) { for(int i=0;i<dim;i++) { for(int j=dim;j<round;j++) { get(i,j)=0; get(j,i)=0; } } for(int i=dim;i<round;i++) { for(int j=dim;j<round;j++) { get(i,j)=(double)(i==j); } } } } matrix<dim>(matrix<halfdim> _a, matrix<halfdim> _b, matrix<halfdim> _c, matrix<halfdim> _d) { a=_a; b=_b; c=_c; d=_d; matrix<dim>(); } };
halfdim의 정의도 약간 손보고, round라는 놈도 새로 만들었다. dim이 2^n이 아닐때, 즉 round>dim일때 나머지 공간에 0과 1을 채워준다.

이렇게 해서 2^n꼴이 아닌 임의의 정사각행렬의 역행렬도 구할수 있게되었다. 템플릿 프로그래밍으로 간단하게(? 좀 복잡헌듸) 역행렬 구하는 프로그램을 짜봤다. (근데 이게 과연 효율적인지는 모르겠다. 가우스 조던 소거법이 더 좋을지도.)

이것은 서비스정신에서 비롯한 FullCode!
template<int n> class roundbit { public: enum{value=roundbit<(n+1)/2>::value*2}; }; template<> class roundbit<2> { public: enum{value=2}; }; template<int null=0> class nullvalue { public: static double v; }; template<int null> double nullvalue<null>::v; template<int dim> class matrix { protected: enum{halfdim=roundbit<dim>::value/2, round=roundbit<dim>::value,}; matrix<halfdim> a,b,c,d; public: static matrix<dim> zero; matrix<dim>() { if(round>dim) { for(int i=0;i<dim;i++) { for(int j=dim;j<round;j++) { get(i,j)=0; get(j,i)=0; } } for(int i=dim;i<round;i++) { for(int j=dim;j<round;j++) { get(i,j)=(double)(i==j); } } } } matrix<dim>(matrix<halfdim> _a, matrix<halfdim> _b, matrix<halfdim> _c, matrix<halfdim> _d) { a=_a; b=_b; c=_c; d=_d; matrix<dim>(); } matrix<dim> operator+(matrix<dim> &t) { return matrix<dim>(a+t.a, b+t.b, c+t.c, d+t.d); } matrix<dim> operator-(matrix<dim> &t) { return matrix<dim>(a-t.a, b-t.b, c-t.c, d-t.d); } matrix<dim> operator-() { return matrix<dim>(-a, -b, -c, -d); } matrix<dim> operator*(matrix<dim> &t) { return matrix<dim>(a*t.a + b*t.c, a*t.b + b*t.d, c*t.a + d*t.c, c*t.b + d*t.d); } double det() { if(a.has_inverse()) { return (a*(d-c*a.inverse()*b)).det(); } if(d.has_inverse()) { return (d*(a-b*d.inverse()*c)).det(); } return 0.0; } matrix<dim> inverse() { matrix<dim> r; matrix<halfdim> q; if(a.has_inverse()) { matrix<halfdim> ai=a.inverse(); if((q=(d-c*ai*b)).has_inverse()) { r.d=q.inverse(); r.c=r.d*(-c*ai); r.b=-ai*b*r.d; r.a=ai-ai*b*r.c; return r; } } if(d.has_inverse()) { matrix<halfdim> di=d.inverse(); if((q=(a-b*di*c)).has_inverse()) { r.a=q.inverse(); r.b=r.a*(-b*di); r.c=-di*c*r.a; r.d=di-di*c*r.b; return r; } } return zero; } bool has_inverse() { if(a.has_inverse() && (d-c*a.inverse()*b).has_inverse()) { return true; } if(d.has_inverse() && (a-b*d.inverse()*c).has_inverse()) { return true; } return false; } double& get(int row, int col) { switch(row/halfdim) { case 0: switch(col/halfdim) { case 0: return a.get(row, col); case 1: return b.get(row, col-halfdim); } case 1: switch(col/halfdim) { case 0: return c.get(row-halfdim, col); case 1: return d.get(row-halfdim, col-halfdim); } } return nullvalue<>::v; } }; template<int dim> matrix<dim> matrix<dim>::zero=matrix<dim> (matrix<halfdim>::zero, matrix<halfdim>::zero, matrix<halfdim>::zero, matrix<halfdim>::zero); template<> class matrix<2> { protected: double a,b,c,d; public: static matrix<2> zero; matrix<2>() { } matrix<2>(double _a, double _b, double _c, double _d) { a=_a; b=_b; c=_c; d=_d; } matrix<2> operator+(matrix<2> &t) { return matrix<2>(a+t.a, b+t.b, c+t.c, d+t.d); } matrix<2> operator-(matrix<2> &t) { return matrix<2>(a-t.a, b-t.b, c-t.c, d-t.d); } matrix<2> operator-() { return matrix<2>(-a, -b, -c, -d); } matrix<2> operator*(matrix<2> &t) { return matrix<2>(a*t.a + b*t.c, a*t.b + b*t.d, c*t.a + d*t.c, c*t.b + d*t.d); } double det() { return a*d-b*c; } matrix<2> inverse() { double dt=det(); if(dt!=0) { return matrix<2>(d/dt, -b/dt, -c/dt, a/dt); } return zero; } bool has_inverse() { return a*d!=b*c; } double& get(int row, int col) { switch(row) { case 0: switch(col) { case 0: return a; case 1: return b; } case 1: switch(col) { case 0: return c; case 1: return d; } } return nullvalue<>::v; } }; matrix<2> matrix<2>::zero=matrix<2>(0.0, 0.0, 0.0, 0.0);