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

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

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

1 2 3 4 5 6 7 8 9 10 11 12 13

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가 되겠다.

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

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!

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 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235

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);