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

*1:実は Tak y も私が出題者です。