anarchy golf の inverse matrix
実は私が出題者です。*1ということで、感想とか問題の種明かし的なものを書いてみます。
問題文
問題文の英語が変だったらごめんなさい。というか、どこが変か教えてください。
sample input と sample output
まず、(Square root に対する感想を読んで)入出力に小数があると、小数の出力フォーマットで言語毎に有利不利がでそうだなと思ったので、行列の要素が全部整数になるようにしてみました。こういう条件をつけると多分行列式は 1 か -1 になってしまうはずです。なので、sample input は、単位行列からスタートして、
- 行の入れ替え
- 列の入れ替え
- ある行に別の行の整数倍をしたものを足す
- ある列に別の列の整数倍をしたものを足す
をランダムに適当な回数繰り返せばいいかな、と思い、下の方にあるプログラムを使って作りました。対応する output は Maxima を利用して作りました。
繰り返しの回数が多かったせいか、Gauss-Jordan 法あたりで解く場合、行列の基本変形のうち、行の入れ替えが 1 回しか必要がない問題になってしまいました。
どの言語が有利か
行列を扱うライブラリを持っている言語が有利です。出題した後に気づいたのですが、例えば Ruby に行列を扱うライブラリがあるようです。
あと、exec を禁止していませんので、C 言語なんかだと Ruby を呼び出した方が短くなったりするんじゃないでしょうか。
PostScript で解いた方法
最初 Gauss-Jordan 法を使って解こうとしましたが、どうも GhostScript の浮動小数点が単精度であるらしく、桁が足りずに変な答えが出てしまいました。そこで、Gauss-Jordan 法をちょっと変形して
- A. 第 n 行から第 9 行の中で、第 n 列の要素が非 0 で絶対値が最も小さい行を探し、第 n 行と交換
- B. 交換後の第 n 行 n 列の要素が負なら n 行目に -1 を掛ける
- C. 第 m 行 n 列 (m ≠ n、m=1..9) の要素を第 n 行 n 列の要素で割ったときの商を q とし、第 m 行から第 n 行の q 倍を引く。
- D. 第 n 行 n 列が 1 の場合は他の行の第 n 列目は全て 0 になっているはずなので、n に 1 足して A に戻る。
- E. 第 n 行 n 列が 1 で無い場合は、n をそのままにして A に戻る。他の行の第 n 列目に、絶対値が第 n 行 n 列の要素より小さい 0 でないものがあるはず (無かったら行列式が 1 か -1 にならない) なので、無限ループにはならないはず。
というような方法を使って解きました。
Sample input 作成プログラム
出力は、逆行列を計算する Maxima プログラムです。
繰り返しの回数が多すぎると int の範囲からあふれるので注意。逆行列に分数が必要になることがあります。
#include<stdio.h> #include<stdlib.h> #define MAX 9 /* 行 i と行 j を入れ替え */ void pijL(int mat[MAX][MAX],int i,int j) { int k,tmp; for(k=0;k<MAX;k++){ tmp=mat[i][k]; mat[i][k]=mat[j][k]; mat[j][k]=tmp; } } /* 列 i と列 j を入れ替え */ void pijR(int mat[MAX][MAX],int i,int j) { int k,tmp; for(k=0;k<MAX;k++){ tmp=mat[k][i]; mat[k][i]=mat[k][j]; mat[k][j]=tmp; } } /* 行 i に行 j の c 倍を足す */ void cijL(int mat[MAX][MAX],int i,int j,int c) { int k; fprintf(stderr,"%d\n",c); for(k=0;k<MAX;k++){ mat[i][k]+=c*mat[j][k]; } } /* 列 i に列 j の c 倍を足す */ void cijR(int mat[MAX][MAX],int i,int j,int c) { int k; fprintf(stderr,"%d\n",c); for(k=0;k<MAX;k++){ mat[k][j]+=c*mat[k][i]; } } int randn(int n) { return (int)(n*(double)rand()/RAND_MAX); } int main() { int matrix[MAX][MAX]; int i,j,k; srand(1); for(i=0;i<MAX;i++){ for(j=0;j<MAX;j++){ matrix[i][j]=(i==j); } } for(k=0;k<100;k++){ i=randn(MAX); do{ j=randn(MAX); }while(i==j); switch(randn(4)){ case 0: pijL(matrix,i,j);break; case 1: pijR(matrix,i,j);break; case 2: cijL(matrix,i,j,randn(8)+randn(8)-7);break; case 3: cijR(matrix,i,j,randn(8)+randn(8)-7);break; default: break; } } puts("matrix("); for(i=0;i<MAX;i++){ printf("["); for(j=0;j<MAX;j++){ printf("%d",matrix[i][j]); if(j<MAX-1)printf(","); } printf("]"); if(i<MAX-1)printf(","); puts(""); } puts(");"); puts("%^^-1;"); return 0; }