3元1次連立方程式の解を求めるコードを書いた

ゼミが終わって寝て起きたら「3元1次連立方程式の解を求めるプログラムを実装してください！C言語で！」という地獄のようなDMが来てました。DM内では全て係数が既知である場合の解を求めよというのだったので、C言語滅べと思いながら実装しました。

しかし、係数がすべて既知のコードを書いても一切面白くなかったので、係数が任意の場合に解を求めるコードを書きました。流石にこっちはC++で書きました。

連立方程式の解を求める方法は色々あります。掃き出し法だったり、余因子を用いて逆行列を出したり。今回は3元1次連立方程式なので、Cramerの公式を用いました。一番実装がわかりやすいので。

書いたコードは次の通り。

	#include <bits/stdc++.h>
	#define rep(i, n) for (int i = 0; i < (n); i++)
	#define rep2(i,s,n) for(int i = (s); i < (n); i++)
	using namespace std;
	using ll = long long;
	using P = pair<int, int>;

	double detA(vector<double> &a) {
	double d = 0;
	d += a[0] * a[4] * a[8];
	d += a[2] * a[3] * a[7];
	d += a[1] * a[5] * a[6];
	d -= a[2] * a[4] * a[6];
	d -= a[1] * a[3] * a[8];
	d -= a[0] * a[5] * a[7];
	return d;
	}

	int main(){
	vector<double> a;
	vector<double> b;
	cout << "左辺の各係数を入力" << endl;
	rep(i, 9) {
	double tmp;
	cin >> tmp;
	a.push_back(tmp);
	}
	cout << "右辺の各値を入力" << endl;
	rep(i, 3) {
	double tmp;
	cin >> tmp;
	b.push_back(tmp);
	}

	double d = detA(a);

	if (d == 0) {
	cout << "｡оО(｡´•ㅅ•｡)Оо。<無理" << endl;
	return 0;
	}

	vector<double> copy_a;
	copy(a.begin(), a.end(),back_inserter(copy_a));

	//１列目にbを代入
	copy_a[0] = b[0], copy_a[3] = b[1], copy_a[6] = b[2];
	//det計算
	double detx = detA(copy_a);
	// vector内をもとに戻す
	copy_a = a;
	//2列目にbを代入
	copy_a[1] = b[0], copy_a[4] = b[1], copy_a[7] = b[2];
	//det計算
	double dety = detA(copy_a);
	cout << dety << endl;
	//もとに戻す
	copy_a = a;
	//3列目にbを代入
	copy_a[2] = b[0], copy_a[5] = b[1], copy_a[8] = b[2];
	//det計算
	double detz = detA(copy_a);
	cout << detz << endl;

	cout << "x=" << detx / d << endl;
	cout << "y=" << dety / d << endl;
	cout << "z=" << detz / d << endl;
	return 0;
	}

view raw cramer.cpp hosted with

by GitHub

ぶっちゃけ愚直に実装しただけなので、もっとエレガントに実装できそうです。

vectorのコピーをするときは、operator=が使えますが、43行目にcopy関数を使っています。バグってたときにここだろと思って変えた部分を戻し忘れただけなのですが、vectorのコピーもいろいろな方法があるようです。

コードの確認は、適当にサンプルケースを作って確認しました。多分合ってるはず。(ケース漏れありそう)