ただ逆行列を求めたい、、、

フォーラム(掲示板)ルール
フォーラム(掲示板)ルールはこちら  ※コードを貼り付ける場合は [code][/code] で囲って下さい。詳しくはこちら
naponapo
記事: 38
登録日時: 11年前

ただ逆行列を求めたい、、、

#1

投稿記事 by naponapo » 11年前

こんにちは、いつもお世話になってます!
n次の正則行列から、逆行列を求めるプログラムを作ってみたくなり、コードを書いたのですが、うまくいきません
テストとして、3次の単位行列の逆行列を求めようとしたプログラムを以下のように作りました
デバッガを使って、間違っているところの検討はつけているのですが、

168行目?

コード:

inv[s][k]=((-1)^(k+s+2))*(det[k][s])/(det1(origin)); //逆行列の各要素の計算
この部分のinv[0][0]がすでに予想と違う値(ー3)になっているのですが、
この右辺の各オペランドは正常な値をとっていることがわかりました
なぜ、値が異なってしまうのでしょうか?
原因がわからなくて困っています
力を貸してもらえませんか
問題は関数⑧です

コード:

#include <stdio.h>
#define NUM 3
int Judge1=0;
int Judge2=0;
void mat_express1(int a[NUM][NUM]){  //NUM次正方行列を表示する    関数①    

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

		printf("\n");
	};



}

void mat_tri1(int b[][NUM]){      //NUM次正方行列を行基本変形して、上三角行列にする   関数②
	int i,j,k;   int mul;  
	Judge1=0;
	for (i=0;i<NUM;i++)
	{
		mul=b[i][i];//対角要素mulをi列掃き出しの要とする
		if(mul==0)//もし、対角要素mulが0なら上三角化強制終了
		{
			Judge1=1;//行列式の値が0であることのフラグを立てる(関数④につなげる)
			break;
		};

		for(k=i+1;k<NUM;k++)
		{
			int c=b[k][i];
			for(j=0;j<NUM;j++)
			{
				double sub=b[k][j]-c*b[i][j]/mul;
				b[k][j]-=c*b[i][j]/mul;

				if(b[k][j]<sub)//もし要素が整数にならないなら、上三角化を諦める
				{
					printf("正しくない");
					break;
				};
			};
		};


	};

}   //終了

void mat_tri2(int b[][NUM-1]){//NUM-1次正方行列の上三角化(関数②と同じ)       関数③
	int i,j,k;   int mul;
	Judge2=0;
	for (i=0;i<NUM-1;i++)
	{
		mul=b[i][i];

		if(mul==0)
		{
			Judge2=1;
			break;
		};

		for(k=i+1;k<NUM-1;k++)
		{
			int c=b[k][i];
			for(j=0;j<NUM-1;j++)
			{
				double sub=b[k][j]-c*b[i][j]/mul;
				b[k][j]-=c*b[i][j]/mul;

				if(b[k][j]<sub)
				{
					printf("正しくない");
					break;
				};
			};
		};


	};

}   //終了

int det1(int va[][NUM]){   //NUM次正方行列の行列式の値を求める          関数④
	int seki=va[0][0];
	if(Judge1==0)
	{
		mat_tri1(va);
		int i;
		for(i=1;i<NUM;i++)
			seki*=va[i][i];//上三角化した行列の対角要素をすべてかける
	}
	if(Judge1==0)
		return(seki);
	else
		return(0);
}

int det2(int va[][NUM-1]){  //NUM-1次正方行列の行列式の値を求める(関数④と同じ)       関数⑤
	int seki=va[0][0];
	if(Judge2==0)
	{
		mat_tri2(va); 
		int i; 
		for(i=0;i<NUM-1;i++)
			seki*=va[i][i];

	};

	if(Judge2==0)
		return(seki);
	else
		return(0);
}

void kurinuki(int origin[][NUM],int small[][NUM-1],int i,int j){ //行列originの i行j列をくり抜いた行列smallを作る  関数⑥
	int n,m;
	for(n=0;n<i-1;n++)
	{
		for(m=0;m<j-1;m++)
		{
			small[n][m]=origin[n][m];
		};
		for(m=j;m<NUM;m++)
		{
			small[n][m-1]=origin[n][m];
		};

	};

	for(n=i;n<NUM;n++)
	{
		for(m=0;m<j-1;m++)
		{
			small[n-1][m]=origin[n][m];
		}
		for(m=j;m<NUM;m++)
		{
			small[n-1][m-1]=origin[n][m];
		}

	};

}

void mat_copy(int origin[][NUM],int copy[][NUM]){//行列originを行列copyに写す   関数⑦
	for(int i=0;i<NUM;i++)
		for(int j=0;j<NUM;j++)
			copy[i][j]=origin[i][j];
}

void inv_mat(int origin[][NUM],int inv[][NUM]){    //行列originの逆行列invを得る  関数⑧
	int small[NUM-1][NUM-1];
	int det[NUM][NUM];
	int k,s;
	int copy[NUM][NUM];
	mat_copy(origin,copy);
	for(k=0;k<NUM;k++)
	{
		for(s=0;s<NUM;s++)
		{   
			mat_copy(copy,origin); //上三角化されたoriginの復元

			kurinuki(origin,small,k+1,s+1); //originのk+1行s+1列をくり抜いてsmall行列を作る
			det[k][s]=det2(small);   //smallの行列式を行列detの(k,s)成分とする

			inv[s][k]=((-1)^(k+s+2))*(det[k][s])/(det1(origin)); //逆行列の各要素の計算
		};

	};


}


int main(void){
	int vc[NUM][NUM]={1,0,0,
		0,1,0,
		0,0,1};

	int gyaku[NUM][NUM];

	mat_express1(vc);

	printf("\n");

	inv_mat(vc,gyaku);
	mat_express1(gyaku);

	return 0;
}

アバター
みけCAT
記事: 6734
登録日時: 14年前
住所: 千葉県
連絡を取る:

Re: ただ逆行列を求めたい、、、

#2

投稿記事 by みけCAT » 11年前

「((-1)^(k+s+2))」は「((k+s+2)&1?-1:1)」の間違いではないですか?
複雑な問題?マシンの性能を上げてOpenMPで殴ればいい!(死亡フラグ)

naponapo
記事: 38
登録日時: 11年前

Re: ただ逆行列を求めたい、、、

#3

投稿記事 by naponapo » 11年前

「((-1)^(k+s+2))」は累乗にならないのですか?

アバター
みけCAT
記事: 6734
登録日時: 14年前
住所: 千葉県
連絡を取る:

Re: ただ逆行列を求めたい、、、

#4

投稿記事 by みけCAT » 11年前

naponapo さんが書きました:「((-1)^(k+s+2))」は累乗にならないのですか?
Haskellではなりますが、C言語ではなりません。
複雑な問題?マシンの性能を上げてOpenMPで殴ればいい!(死亡フラグ)

naponapo
記事: 38
登録日時: 11年前

Re: ただ逆行列を求めたい、、、

#5

投稿記事 by naponapo » 11年前

では、定数の変数乗を表すには、何か関数を使わなければ表せないのですか?

アバター
みけCAT
記事: 6734
登録日時: 14年前
住所: 千葉県
連絡を取る:

Re: ただ逆行列を求めたい、、、

#6

投稿記事 by みけCAT » 11年前

naponapo さんが書きました:では、定数の変数乗を表すには、何か関数を使わなければ表せないのですか?
特別な関数を使わないと表せない、ということはありません。
「変数」が整数なら、素直に計算するか、繰り返し二乗法で求めることができます。
「変数」が実数なら、expとlogをマクローリン展開/テイラー展開で計算すれば求めることができますが、
通常の環境では素直にpow関数を使ったほうがいいでしょう。
複雑な問題?マシンの性能を上げてOpenMPで殴ればいい!(死亡フラグ)

naponapo
記事: 38
登録日時: 11年前

Re: ただ逆行列を求めたい、、、

#7

投稿記事 by naponapo » 11年前

すこし修正してみます

naponapo
記事: 38
登録日時: 11年前

Re: ただ逆行列を求めたい、、、

#8

投稿記事 by naponapo » 11年前

ありがとうございました!
おかげさまで、うまく求められるようになりました
(要素に分数が出てきたら求められませんけど、、、)
ご指摘のあった問題の行以外にもバグが一箇所あったのですが、そこは自己解決でなんとかなりました!
C言語での累乗の表現の仕方はとても勉強になりました(3項演算子、&演算子、pow関数について)
これから気をつけたいと思います

naponapo
記事: 38
登録日時: 11年前

Re: ただ逆行列を求めたい、、、

#9

投稿記事 by naponapo » 11年前

解決しました

アバター
みけCAT
記事: 6734
登録日時: 14年前
住所: 千葉県
連絡を取る:

Re: ただ逆行列を求めたい、、、

#10

投稿記事 by みけCAT » 11年前

解決でしたら、解決チェックをお願いします。
完成したコードも貼っていただけるとありがたいです。
複雑な問題?マシンの性能を上げてOpenMPで殴ればいい!(死亡フラグ)

naponapo
記事: 38
登録日時: 11年前

Re: ただ逆行列を求めたい、、、

#11

投稿記事 by naponapo » 11年前

行列の要素に分数が出てくると、うまく求められないですけど、
こんな感じです!

コード:

#include <stdio.h>
#define NUM 3
int Judge1=0;
int Judge2=0;
void mat_express1(int a[NUM][NUM]){  //NUM次正方行列を表示する    関数①    

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

		printf("\n");
	};



}

void mat_tri1(int b[][NUM]){      //NUM次正方行列を行基本変形して、上三角行列にする   関数②
	int i,j,k;   int mul;  
	Judge1=0;
	for (i=0;i<NUM;i++)
	{
		mul=b[i][i];//対角要素mulをi列掃き出しの要とする
		if(mul==0)//もし、対角要素mulが0なら上三角化強制終了
		{
			Judge1=1;//行列式の値が0であることのフラグを立てる(関数④につなげる)
			break;
		};

		for(k=i+1;k<NUM;k++)
		{
			int c=b[k][i];
			for(j=0;j<NUM;j++)
			{
				double sub=b[k][j]-c*b[i][j]/mul;
				b[k][j]-=c*b[i][j]/mul;

				if(b[k][j]<sub)//もし要素が整数にならないなら、上三角化を諦める
				{
					printf("正しくない");
					break;
				};
			};
		};


	};

}   //終了

void mat_tri2(int b[][NUM-1]){//NUM-1次正方行列の上三角化(関数②と同じ)       関数③
	int i,j,k;   int mul;
	Judge2=0;
	for (i=0;i<NUM-1;i++)
	{
		mul=b[i][i];

		if(mul==0)
		{
			Judge2=1;
			break;
		};

		for(k=i+1;k<NUM-1;k++)
		{
			int c=b[k][i];
			for(j=0;j<NUM-1;j++)
			{
				double sub=b[k][j]-c*b[i][j]/mul;
				b[k][j]-=c*b[i][j]/mul;

				if(b[k][j]<sub)
				{
					printf("正しくない");
					break;
				};
			};
		};


	};

}   //終了

int det1(int va[][NUM]){   //NUM次正方行列の行列式の値を求める          関数④
	int seki=va[0][0];
	if(Judge1==0)
	{
		mat_tri1(va);
		int i;
		for(i=1;i<NUM;i++)
			seki*=va[i][i];//上三角化した行列の対角要素をすべてかける
	}
	if(Judge1==0)
		return(seki);
	else
		return(0);
}

int det2(int va[][NUM-1]){  //NUM-1次正方行列の行列式の値を求める(関数④と同じ)       関数⑤
	int seki=va[0][0];
	
	
		mat_tri2(va); 
		int i; 
		for(i=0;i<NUM-1;i++)
			seki*=va[i][i];



	if(Judge2==0)
		return(seki);
	else
		return(0);
}

void kurinuki(int origin[][NUM],int small[][NUM-1],int i,int j){ //行列originの i行j列をくり抜いた行列smallを作る  関数⑥
	int n,m;
	for(n=0;n<i-1;n++)
	{
		for(m=0;m<j-1;m++)
		{
			small[n][m]=origin[n][m];
		};
		for(m=j;m<NUM;m++)
		{
			small[n][m-1]=origin[n][m];
		};

	};

	for(n=i;n<NUM;n++)
	{
		for(m=0;m<j-1;m++)
		{
			small[n-1][m]=origin[n][m];
		}
		for(m=j;m<NUM;m++)
		{
			small[n-1][m-1]=origin[n][m];
		}

	};

}

void mat_copy(int origin[][NUM],int copy[][NUM]){//行列originを行列copyに写す   関数⑦
	for(int i=0;i<NUM;i++)
		for(int j=0;j<NUM;j++)
			copy[i][j]=origin[i][j];
}

void inv_mat(int origin[][NUM],int inv[][NUM]){    //行列originの逆行列invを得る  関数⑧
	int small[NUM-1][NUM-1];
	int det[NUM][NUM];
	int k,s;
	int copy[NUM][NUM];
	mat_copy(origin,copy);
	for(k=0;k<NUM;k++)
	{
		for(s=0;s<NUM;s++)
		{   
			mat_copy(copy,origin); //上三角化されたoriginの復元

			kurinuki(origin,small,k+1,s+1); //originのk+1行s+1列をくり抜いてsmall行列を作る
			det[k][s]=det2(small);   //smallの行列式を行列detの(k,s)成分とする

			inv[s][k]=((k+s+2)&1?-1:1)*(det[k][s])/(det1(origin)); //逆行列の各要素の計算
		};

	};


}


int main(void){
	int vc[NUM][NUM]={1,0,0,
		0,1,0,
		0,0,3};

	int gyaku[NUM][NUM];

	mat_express1(vc);

	printf("\n");

	inv_mat(vc,gyaku);
	mat_express1(gyaku);

	return 0;
}

閉鎖

“C言語何でも質問掲示板” へ戻る