行列式のプログラムでのループがわかりません。

質問

行列式を求めるプログラムを考えていて、一応はちゃんと答えがでるのですが、
ループを使って計算させようとしたのですがイマイチわかりません。

#include <stdio.h>

main()
{
　　　　int a[3][3]={{1,2,2},{-1,5,3},{2,0,1}};
　　　　int i,j,n;

　　　　for(i=0;i<3;i++)
　　　　{
　　　　　　　　　for(j=0;j<3;j++)
　　　　　　　　　　　　　printf("%3d",a[i][j]);
　　　　　　　　　　　　　printf("\n");
　　　　}

　　　　n=a[0][0]*(a[1][1]*a[2][2]-a[1][2]*a[2][1])
　　　　　-a[1][0]*(a[0][1]*a[2][2]-a[0][2]*a[2][1])
　　　　　+a[2][0]*(a[0][1]*a[1][2]-a[0][2]*a[1][1]);

　　　　printf("この行列の行列式は \n");
　　　　printf("%d \n",n);
}

1 2 2
-1 5 3
2 0 1
という行列式です。
nの処理をループを使うとどう書いたらいいでしょうか？
これはたまたま３行３列だからひとつひとつ指定して計算できますが
これが４行４列とかだったら絶対ループ使わなきゃ面倒ですからね～。
とりあえず３行３列のが理解できたらいいとは思うのですがなかなか複雑で難しいですね。

回答No.8

graphaffine · ベストアンサー

今日は、usui323さん。
私の方法は、#7で二つ目に紹介されているサイトと本質的に同じですね。
定義から求めているので当たり前ですが。
ですので、概略に留めさせてもらいます。

　nは、行列の次数、fact_nはｎの階乗とする。
　junは順列を格納する符号無し文字型配列
　GetPermutate　i番目の順列を返す関数。iは１以上fact_n以下
　Sign　順列の符号を返す関数
　A(i,j) 行列の(i,j)成分
　全角空白でインデントしてます

　det=0;
　for (cnt1=1;cnt1<=fact_n;cnt1++)
　{
　　　GetPermutate(cnt1,jun);
　　　term = Sign(jun);
　　　for (cnt2=0;cnt2<n;cnt2++)
　　　　　termproduct *= A(cnt2+1,jun(cnt2));
　　　det += termproduct;
　}

お礼

回答ありがとうございました。

あまりプログラムに慣れていないので理解するのが大変でしたが、大変参考になりました。

回答No.9

graphaffine

#8の訂正

誤　termproduct *= A(cnt2+1,jun(cnt2));

正　termproduct *= A(cnt2+1,jun[cnt2]);

回答No.7

hofuhofu

#2で書いた場所はここです
http://mukun_mmg.tripod.co.jp/mmg/bncpp/al047.html

ちなみにこのままでは1行1列からn行n列への対角線上に0がでてきた場合はエラーが発生します。

再帰を使ったプログラムがあったので載せておきます。
ここの「おまけプログラム」にある「行列式を定義に従って求める」です。
http://www.altum.jp/math/exp1/

お礼

回答ありがとうございました。

大変参考になりました。m(-_-)m

回答No.6

absurd0rt

下のコードをエクセルに貼り付け、
"A1"=1,"B1"=2,"C1"=2
"A2"=-1,"B2"=5,"C2"=3
"A3"=2,"B3"=0,"C3"=1
と打ち下のコードを実行すれば"行列式"がでます。
ただVBですし、例外処理をすっぽかして作ってあります。
作ってて正攻法の方が絶対にいいかなと思いました。
参考程度にお願いします。

Private Sub CommandButton1_Click()
　　'*GYAKU
　　Dim ArrayA() As Double
　　Dim ArrayL() As Double
　　Dim ArrayR() As Double
　　Dim InvMatrix() As Double
　　Dim ArrayB() As Double
　　Dim TempB() As Double
　　Dim lngCnt1 As Long
　　Dim lngCnt2 As Long
　　Dim lngCnt3 As Long
　　Dim lngCnt4 As Long
　　Dim lngNCount As Long
　　Dim lngGyoShiki As Long

　　'初期設定
　　For lngCnt1 = 1 To 100
　　　　If IsEmpty(Worksheets("Sheet1").Cells(lngCnt1, 1)) = True Then
　　　　　　Exit For
　　　　End If
　　Next lngCnt1
　　lngNCount = lngCnt1 - 1
　　If lngNCount < 1 Then
　　　　Exit Sub
　　End If

　　ReDim ArrayA(lngNCount, lngNCount)
　　ReDim ArrayL(lngNCount, lngNCount)
　　ReDim ArrayR(lngNCount, lngNCount)
　　ReDim InvMatrix(lngNCount, lngNCount)
　　ReDim ArrayB(lngNCount)
　　ReDim TempB(lngNCount)

　　'初期値
　　For lngCnt1 = 1 To lngNCount
　　　　For lngCnt2 = 1 To lngNCount
　　　　　　ArrayA(lngCnt1, lngCnt2) = Worksheets("Sheet1").Cells(lngCnt1, lngCnt2)
　　　　Next lngCnt2
　　Next lngCnt1

　　'Main
　　For lngCnt1 = 1 To lngNCount
　　　　For lngCnt2 = 1 To lngNCount
　　　　　　If lngCnt2 = lngCnt1 Then
　　　　　　　　ArrayB(lngCnt2) = 1
　　　　　　Else
　　　　　　　　ArrayB(lngCnt2) = 0
　　　　　　End If
　　　　Next lngCnt2

　　　　'STEP1 行列Ｒを作る
　　　　For lngCnt2 = 1 To lngNCount
　　　　　　For lngCnt3 = 1 To lngNCount
　　　　　　　　ArrayR(lngCnt2, lngCnt3) = ArrayA(lngCnt2, lngCnt3)
　　　　　　　　TempB(lngCnt2) = ArrayB(lngCnt2)
　　　　　　Next lngCnt3
　　　　Next lngCnt2
　　　　For lngCnt2 = 1 To lngNCount
　　　　　　For lngCnt3 = 1 To lngNCount
　　　　　　　　If lngCnt3 <> lngCnt2 Then
　　　　　　　　　　ArrayL(lngCnt3, lngCnt2) = ArrayR(lngCnt3, lngCnt2) / ArrayR(lngCnt2, lngCnt2)
　　　　　　　　　　For lngCnt4 = lngCnt2 To lngNCount
　　　　　　　　　　　　ArrayR(lngCnt3, lngCnt4) = ArrayR(lngCnt3, lngCnt4) - ArrayL(lngCnt3, lngCnt2) * ArrayR(lngCnt2, lngCnt4)
　　　　　　　　　　Next lngCnt4
　　　　　　　　　　TempB(lngCnt3) = TempB(lngCnt3) - ArrayL(lngCnt3, lngCnt2) * TempB(lngCnt2)
　　　　　　　　End If
　　　　　　Next lngCnt3
　　　　ArrayL(lngCnt2, lngCnt2) = 1
　　　　Next lngCnt2

　　　　'STEP2 逆行列ＡＡを求める
　　　　For lngCnt2 = 1 To lngNCount
　　　　　　InvMatrix(lngCnt2, lngCnt1) = TempB(lngCnt2) / ArrayR(lngCnt2, lngCnt2)
　　　　Next lngCnt2
　　Next lngCnt1

　　'STEP3 行列式
　　lngGyoShiki = 1
　　For lngCnt1 = 1 To lngNCount
　　　　lngGyoShiki = lngGyoShiki * ArrayR(lngCnt1, lngCnt1)
　　Next lngCnt1

　　'表示
　　Worksheets("Sheet1").Cells(lngNCount * 1 + 2, 1) = "逆行列"
　　For lngCnt1 = 1 To lngNCount
　　　　For lngCnt2 = 1 To lngNCount
　　　　　　Worksheets("Sheet1").Cells(lngCnt1 + (lngNCount * 1 + 2), lngCnt2) = InvMatrix(lngCnt1, lngCnt2)
　　　　Next lngCnt2
　　Next lngCnt1
　　Worksheets("Sheet1").Cells(lngNCount * 2 + 4, 1) = "固有値"
　　For lngCnt1 = 1 To lngNCount
　　　　Worksheets("Sheet1").Cells((lngNCount * 2 + 4) + 1, lngCnt1) = ArrayR(lngCnt1, lngCnt1)
　　Next lngCnt1
　　Worksheets("Sheet1").Cells(lngNCount * 2 + 7, 1) = "行列式"
　　Worksheets("Sheet1").Cells(lngNCount * 2 + 8, 1) = lngGyoShiki

　　'後始末
　　Erase ArrayA
　　Erase ArrayL
　　Erase ArrayR
　　Erase InvMatrix
　　Erase ArrayB
　　Erase TempB
End Sub

お礼

回答ありがとうございます。

VBですか。VBは授業でとりあつかっていないので
全然わからないのですが、こんなに複雑になるんですね。

詳細なプログラムまで書いていただきありがとうございました。

回答No.5

graphaffine

私は１０次までの行列式計算を行うプログラムを作っています。別にそれ以上ができないというわけではないですが、項数が多い（１０次の場合で約３６０万項）ので１０次で切っているだけです。

で、プログラムはＣ＋＋でテンプレートを使っていますし、処理内容は行列式の定義にきちんと則った（順列計算を使った）やり方です。

このへんが、理解できるという事なら、プログラムの所在をお知らせします。

お礼

回答ありがとうございます。
一応問題解決しました。

10次ですか、すごいですね。

100%理解できるというわけではありませんが、参考の為教えてくれるとありがたいです。

回答No.4

kary

No.2の方のご回答のように、ガウスの消去法を使う方法があると思います。

#include <stdio.h>
#include <math.h>

#define N 5

void main()
{
double X[N][N] = {
{ 1.0, 6.0, 4.0, 3.0, 9.0},
{ 2.0, 1.0, 8.0, 6.0, 7.0},
{ 3.0, 7.0, 1.0, 9.0, 5.0},
{ 4.0, 2.0, 5.0, 1.0, 3.0},
{ 5.0, 8.0, 9.0, 4.0, 1.0},
};
double A = 1.0;
int i, j, k;

for (i=N-1; i>=0; i--) {
for (j=i-1; j>=0; j--) {
for (k=0; k<i; k++) {
X[j][k] -= X[i][k] * X[j][i] / X[i][i];
}
}
A *= X[i][i] * pow(-1, i+i);
}

printf("%f\n", A);
}
こんな感じになると思います。(とりあえず、この例では正解が得られるようですが、間違いがあるかもしれません。間違っていたらごめんなさい。)

回答No.3

absurd0rt

行列式にしろ、固有値にしろ対角行列に変換さえしてしまえば楽に求められるものです。そうすれば４行４列など問題ではありません。
プログラムは気が向いたらヒントくらい書くかもしれませんがあまり期待しないでください。

お礼

回答ありがとうございました。

問題解決しました。

回答No.2

hofuhofu

とあるメールマガジンのバックナンバーにn行n列の行列式を求めるサンプルプログラムがありました。
Gaussの消去法を使って出す方法ですが、割と機械的に処理できているようです。
計算はdoubleで行っていました。
まったく誤差を出さないようにと、すべてintでしようとすると、かなり手間が増えそうな感じです。

お礼

回答ありがとうございます。

問題解決しました。

回答No.1

jmh

こんな公式があったような気がします:
det(Ａ)＝Σ(-1)^(i-1)×(Ａの(ｉ,１)成分)
　　　×det(Ａから第ｉ行と第１列を取り除いた行列)
　（Σはｉ＝１～Ａの次数）
これを使って、行列の次数を下げていく再帰を書いたらできるかもしれないです。

お礼

回答ありがとうございました。

問題解決しました。