Language: C++

#include <bits/stdc++.h>
#include <boost/multiprecision/cpp_int.hpp>
using namespace std;

//#define int long long
//using ll = boost::multiprecision::cpp_int;
using ll = long long;
//constexpr int INF = 1e9;//INT_MAX=(1<<31)-1=2147483647
constexpr ll INF = (ll)1e18;//(1LL<<63)-1=9223372036854775807
constexpr ll MOD = (ll)1e9 + 7;
constexpr double EPS = 1e-9;
constexpr int dx[4]={1,0,-1,0};
constexpr int dy[4]={0,1,0,-1};
#ifdef LOCAL_ENV
	#define debug(var) std::cout<<#var" = "<<var<<std::endl
#else
	#define debug(var)
#endif
#define p(var) std::cout<<var<<std::endl
#define rep(i, n)   for(ll i=0, i##_length=(n); i< i##_length; ++i)
#define repeq(i, n) for(ll i=1, i##_length=(n); i<=i##_length; ++i)
#define all(a) (a).begin(), (a).end()
#define PI (acos(-1))
template<typename T> inline T power(T x, T n){ return n ? power(x*x%MOD,n/2)*(n%2?x:1)%MOD : 1; }
template<typename T> inline T gcd(T a, T b) { return b ? gcd(b,a%b) : a; }
template<typename T> inline T lcm(T a, T b) { return a / gcd(a,  b) * b; }
template<typename T1, typename T2>inline bool chmax(T1 &a, T2 b) { return a < b && (a = b, true); }
template<typename T1, typename T2>inline bool chmin(T1 &a, T2 b) { return a > b && (a = b, true); }
template<typename T1, typename T2> ostream& operator<<(ostream& s, const pair<T1, T2>& p) {return s << "(" << p.first << ", " << p.second << ")";}
template<typename T> ostream& operator<<(ostream& s, const vector<T>& v) {
	for (int i = 0, len = v.size(); i < len; ++i){
		s << v[i]; if (i < len - 1) s << "\t";
	}
	return s;
}
template<typename T> ostream& operator<<(ostream& s, const vector< vector<T> >& vv) {
	for (int i = 0, len = vv.size(); i < len; ++i){
		s << vv[i] << endl;
	}
	return s;
}
template<typename T1, typename T2> ostream& operator<<(ostream& s, const map<T1, T2>& m) {
	s << "{" << endl;
	for (typeof(m.begin()) itr = m.begin(); itr != m.end(); ++itr){
		s << "\t" << (*itr).first << " : " << (*itr).second << endl;
	}
	s << "}" << endl;
	return s;
}

/*-----8<-----8<-----*/


vector<double> simequ(const vector<vector<double>> &A, const vector<double> &b){
	ll N = A.size();
	vector<vector<double>> B(N, vector<double>(N + 1));
	for(ll i = 0; i < N; i++)
		for(ll j = 0; j < N; j++)
			B[i][j] = A[i][j];
	//行列Aの後ろにbを並べて同時に処理する
	for(ll i = 0; i < N; i++)
		B[i][N] = b[i];

	// 前進消去
	for (int i = 0; i < N - 1; i++) {
		//注目している変数の係数の絶対値が大きい式をi番目に持ってくる
		ll pivot = i;
		for(ll j = i + 1; j < n; j++){
			if(abs(B[j][i]) > abs(B[pivot][i])) pivot = j;
		}
		if(pivot != i)swap(B[pivot], B[i]);//???

		for (int j = i + 1; j < N; j++) {
			double d = B[j][i] / B[i][i];
			for (int k = i + 1; k <= N; k++)
				B[j][k] -= B[i][k] * d;
		}
	}

	// 後退代入
	for (int i = N - 1; i >= 0; i--) {
		double d = B[i][N];
		for (int j = i + 1; j < N; j++)
			d -= B[i][j] * B[j][N];
		B[i][N] = d / B[i][i];
	}

	// 結果出力
	vector<double> result(N);
	for (int i = 0; i < N; i++)
		result[i] = B[i][N];

	return result;
}

//Ax=bを解く Aは正方行列
//解がない、もしくは一意でない場合は長さ0の配列を返す
vector<double> gauss_jordan(const vector<vector<double>> &A, const vector<double> &b){
	ll n = A.size();
	vector<vector<double>> B(n, vector<double>(n + 1));
	for(ll i = 0; i < n; i++)
		for(ll j = 0; j < n; j++)
			B[i][j] = A[i][j];
	//行列Aの後ろにbを並べて同時に処理する
	for(ll i = 0; i < n; i++)
		B[i][n] = b[i];

	for(ll i = 0; i < n; i++){
		//注目している変数の係数の絶対値が大きい式をi番目に持ってくる
		ll pivot = i;
		for(ll j = i + 1; j < n; j++){
			if(abs(B[j][i]) > abs(B[pivot][i])) pivot = j;
		}
		//if(pivot != i)swap(B[pivot], B[i]);//???
		for(ll j = i; j <= n; j++){
			double t=B[i][j];
			B[i][j]=B[pivot][j];
			B[pivot][j]=t;
		}

		//解がない、もしくは一意でない
		if(abs(B[i][i]) < EPS) return vector<double>();

		//注目している変数の係数を1にする
		for(ll j = i + 1; j <= n; j++)
			B[i][j] /= B[i][i];
		
		for(ll j = 0; j < n ; j++){
			if(i != j){
				//j番目の式からi番目の変数を消去
				for(ll k = i + 1 ; k <= n; k++)
					B[j][k] -= B[j][i] * B[i][k];
			}
		}
	}
	vector<double> x(n);
	//後ろに並べたbが解になる
	for(ll i = 0; i < n; i++)
		x[i] = B[i][n];

	return x;
}

inline bool isnatural(double x) { return x > -EPS && ceil(x) == floor(x); }

signed main(){
	double N, Y;
	cin >> N >> Y;
	rep(i, N+1){
		double t = Y - 10000.0 * i;
		if(t < 0)break;
		vector<vector<double>> a;
		vector<double> b;
		a.push_back({    1.0,     1.0 });	b.push_back(N - i);
		a.push_back({ 5000.0,  1000.0 });	b.push_back(t);
		vector<double> x = simequ(a,b);
		if(x.size()>0 && isnatural(x[0]) && isnatural(x[1])){
			cout << i << " " << (int)round(x[0]) << " " << (int)round(x[1]) << endl;
			return 0;
		}
	}
	cout << "-1 -1 -1" << endl;
	return 0;
}


signed main2(){
	//連立方程式を解く 1.0w + -2.0x + 3.0y + -4.0z = 5.0 ...の(w,x,y,z)が解ける
	vector<vector<double>> a;
	vector<double> b;
	a.push_back({ 1.0, -2.0,  3.0, -4.0}); b.push_back( 5.0);
	a.push_back({-2.0,  5.0,  8.0, -3.0}); b.push_back( 9.0);
	a.push_back({ 5.0,  4.0,  7.0,  1.0}); b.push_back(-1.0);
	a.push_back({ 9.0,  7.0,  3.0,  5.0}); b.push_back( 4.0);

	vector<double> c=gauss_jordan(a,b);

	for(int i=0;i<c.size();i++){
		cout << c[i] << endl;
	}

	return 0;
}


Compiler output:
a.cc:2:10: fatal error: boost/multiprecision/cpp_int.hpp: No such file or directory
    2 | #include <boost/multiprecision/cpp_int.hpp>
      |          ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
compilation terminated.