行列式の値の求め方、逆行列の作り方の C 言語プログラム

なぜかこの時期に「行列式の値」と「逆行列」を作成するプログラムを毎年書いている気がして、その都度作るのが非常に面倒臭いので、今後も困らないようにアップしておく。
行列式の値の求め方
・2 次の場合
たすきがけで簡単に解ける。

double a[2][2]={{1,2},{4,-1}};
double det=0.0;

det=a[0][0]*a[1][1];
det-=a[0][1]*a[1][0];
printf("%f\n",det); // -> -9.000000

・3 次の場合
サラスの方法で簡単に解ける。

double a[3][3]={{1,8,9},{-3,2,1},{4,1,5}};
double det=0.0;

det=a[0][0]*a[1][1]*a[2][2];
det+=a[1][0]*a[2][1]*a[0][2];
det+=a[2][0]*a[0][1]*a[1][2];
det-=a[2][0]*a[1][1]*a[0][2];
det-=a[1][0]*a[0][1]*a[2][2];
det-=a[0][0]*a[2][1]*a[1][2];
printf("%f\n",det); // -> 62.000000

・4 次以上の場合
通常の計算では 3×3 以下の余因子行列を作成して、符号を掛けた上で足し合わせるやり方を使うと思う。
これをプログラムで実装するには、再帰または多段ループが必要になり計算量も多い。
特に 5 次以降になると面倒。
でも、上三角行列を作って、対角部分 (左上から右下) を掛け合わせれば、簡単に行列式の値を求めることができる。
このやり方だと非常に簡単。

double a[4][4]={{2,-2,4,2},{2,-1,6,3},{3,-2,12,12},{-1,3,-4,4}};
double det=1.0,buf;
int n=4;  //配列の次数
int i,j,k;

//三角行列を作成
for(i=0;i<n;i++){
 for(j=0;j<n;j++){
  if(i<j){
   buf=a[j][i]/a[i][i];
   for(k=0;k<n;k++){
   a[j][k]-=a[i][k]*buf;
   }
  }
 }
}

//対角部分の積
for(i=0;i<n;i++){
 det*=a[i][i];
}

printf("%f\n",det); // -> 120.000000

逆行列の求め方
何といってもやはり掃き出し法が簡単。2 次から n 次まで対応。
掃き出し法の他に、LU 分解で解く方法もあります。

double a[4][4]={{1,2,0,-1},{-1,1,2,0},{2,0,1,1},{1,-2,-1,1}}; //入力用の配列
double inv_a[4][4]; //ここに逆行列が入る
double buf; //一時的なデータを蓄える
int i,j,k; //カウンタ
int n=4;  //配列の次数

//単位行列を作る
for(i=0;i<n;i++){
 for(j=0;j<n;j++){
 inv_a[i][j]=(i==j)?1.0:0.0;
 }
}
//掃き出し法
for(i=0;i<n;i++){
 buf=1/a[i][i];
 for(j=0;j<n;j++){
 a[i][j]*=buf;
 inv_a[i][j]*=buf;
}
for(j=0;j<n;j++){
 if(i!=j){
  buf=a[j][i];
  for(k=0;k<n;k++){
   a[j][k]-=a[i][k]*buf;
   inv_a[j][k]-=inv_a[i][k]*buf;
  }
 }
}
}
//逆行列を出力
for(i=0;i<n;i++){
 for(j=0;j<n;j++){
  printf(" %f",inv_a[i][j]);
 }
 printf("\n");
}

/* 出力
2.000000 2.000000 -1.000000 3.000000
-4.000000 -5.000000 3.000000 -7.000000
3.000000 4.000000 -2.000000 5.000000
-7.000000 -8.000000 5.000000 -11.000000
*/

プログラムは　Microsoft Visual Studio 2008 および Fedora 9 の gcc で実行確認済み。
手元に Linux のマシンがなかったので、稼働中の Web サーバを使って動作確認した。

最強最速アルゴリズマー養成講座　プログラミングコンテストTopCoder攻略ガイド

posted with amazlet at 16.06.02

SBクリエイティブ (2013-08-14)
売り上げランキング: 79,620

Amazon.co.jpで詳細を見る