템플릿과 역행렬 2. 본격 역행렬 구하기

행렬을 작은 크기의 행렬 4개로 나눠어서 계산하는 것에 대해서는 이미 수학자들이 많이 연구해 놓았다. 그러나 글쓴이는 그걸 모르고 공책 한 바닥을 써가며 역행렬 구하는 공식을 구해는 삽질을 했다. (아이구 시간 아까워라.) 무궁한 지식의 원천 위키백과님에게 도움을 요청합니다. 그 공식은 위키백과에 상세히 나와있다.

이제 그럼 본격적으로 역행렬을 구하는 함수를 짜보자.
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; }
역시나 template<int dim> class matrix 안에 들어가야 할 코드이다. 먼저 조심스럽게 역함수가 존재하는지 물어보는 함수부터 짜봤다. 위 수식에서 보이다시피 첫번째 경우는 A의 역행렬이 존재하고, D-CA^(-1)B의 역행렬이 존재해야만 역행렬이 존재하고, 두번째 경우는 D의 역행렬이 존재하고, A-BD^(-1)C의 역행렬이 존재해야만 역행렬이 존재한다. 그걸 코드로 그대로 옮긴것이다.

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; }
역시 template<int dim> class matrix 안에 들어가야할 함수로, 위의 수식에 따라 그대로 역행렬을 구현했다. 역행렬이 존재하지 않을 경우에는 영행렬을 리턴하도록 했다.

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; }
역시 template<int dim> class matrix 안에 들어가야할 함수다. 행렬식을 구해주는 함수로, 역시 2가지 경우로 나누어서 구현하면 장땡.

이제 m행n열의 값을 얻어오는 함수를 만들어보자. 간단한건데 앞에서 까먹고 안 만들었네...
template<int null=0> class nullvalue { public: static double v; }; template<int null> double nullvalue<null>::v; 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> class matrix안에 들어가야한다.*/ 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; } /*위 함수는 class matrix<2>에 들어가야한다.*/
좀 지저분하게 길어졌다. 앞부분의 nullvalue는 일단 무시하고, 핵심적인거는 row, col값을 확인해봐서 a,b,c,d 중 어디에 포함되는 성분인지 알아보고 재귀시키는 것이다. nullvalue는 범위를 벗어났을때 리턴하는 쓸모없는 변수를 만드려고 만든 클래스이다. get함수가 double&를 리턴하므로, NULL값을 줄수 없으므로 nullvalue::v를 리턴한다. (난 늘 이렇게 쓰는데 이게 좋은 코딩 습관인지는 모르겠다.)
그런데 nullvalue를 template<int null=0>으로 선언한 이유는, 위와 같이 템플릿으로 선언하면 변수정의를 헤더에 두어도 정적변수 static double v가 중복정의되지 않기 때문이다. 저게 마음에 안 들면 그냥 템플릿을 없애고 헤더 파일말고 소스파일에 정적변수를 정의하면 된다.

이게 끝이다. 이렇게 함으로써 2^n의 행렬의 역행렬을 구하는 코드가 완성된 것이다. 메인함수에 matrix<64> 짜리 변수를 만들어서 랜덤하게 성분값을 줘보고, 역행렬을 구해보시죠. 만족스러운 결과입니다. 단지 아쉬운 점이라면 템플릿클래스라서 임의 크기의 행렬을 지원하지 않는다는 것!