#include <bits/stdc++.h>
using namespace std;
using i64 = long long;

struct point { i64 x, y, z; };

point operator+(point a, point b) { return {a.x+b.x, a.y+b.y, a.z+b.z}; }
point operator-(point a, point b) { return {a.x-b.x, a.y-b.y, a.z-b.z}; }
i64 operator*(point a, point b) { return a.x*b.x + a.y*b.y + a.z*b.z; } // dot product
point operator%(point a, point b) { // cross product
	return {a.y*b.z-a.z*b.y, a.z*b.x-a.x*b.z, a.x*b.y-a.y*b.x}; 
}

bool operator==(const point &a, const point &b) {
	return make_tuple(a.x,a.y,a.z) == make_tuple(b.x,b.y,b.z);
}
bool operator<(const point &a, const point &b) {
	return make_tuple(a.x,a.y,a.z) < make_tuple(b.x,b.y,b.z);
}

point normalize(point a) {
	i64 g = abs(gcd(a.x, gcd(a.y, a.z)));
	return {a.x/g, a.y/g, a.z/g};
}

const int N = 100;

struct convex_hull {
	struct face {
		int a, b, c;
		point normal;
	};
	
	int n;
	i64 volume;
	vector<point> p;
	vector<face> faces;
	array<bitset<N>,N> mark;

	void add_face(int a, int b, int c) {
		faces.push_back({a, b, c, (p[b]-p[a])%(p[c]-p[a])});
		mark[a][b] = mark[b][c] = mark[c][a] = 1;
		volume += p[a] * (p[c] % p[b]);
	}

	void remove_face(int a, int b, int c) {
		mark[a][b] = mark[b][c] = mark[c][a] = 0;
		volume -= p[a] * (p[c] % p[b]);
	}

	bool add_point(point q) {
		vector<face> kept_faces;
		for (face &f: faces) {
			if ((q-p[f.a]) * f.normal > 0) {
				remove_face(f.a, f.b, f.c);
			} else {
				kept_faces.push_back(f);
			}
		}
		if (kept_faces.size() == faces.size()) return false;
		faces = kept_faces, p.push_back(q);
		for (face &f: kept_faces) {
			if (!mark[f.a][f.c]) add_face(f.a, f.c, n);
			if (!mark[f.b][f.a]) add_face(f.b, f.a, n);
			if (!mark[f.c][f.b]) add_face(f.c, f.b, n);
		}
		if (++n == N) remove_extra_points();
		return true;
	}
	
	void remove_extra_points() {
		set<point> q(begin(p), end(p));
		
		vector<vector<face>> point_faces(n);
		for (face &f: faces) {
			point_faces[f.a].push_back(f);
			point_faces[f.b].push_back(f);
			point_faces[f.c].push_back(f);
		}
		for (int i = 0; i < n; i++) {
			if (!q.contains(p[i])) continue;
			if (point_faces[i].empty()) {
				q.erase(p[i]);
			} else {
				bool flat = true;
				point normal = normalize(point_faces[i][0].normal);
				for (face &f: point_faces[i]) {
					if (normalize(f.normal) != normal) flat = false;
				}
				if (flat) q.erase(p[i]);
			}
		}
		set<point> r = q;
		for (point a: q) for (point b: q) if (a < b) for (point c: q) if (b < c) {
			if ((b-a)%(c-a) == point{0,0,0}) r.erase(b);
		}
		
		vector<point> s(begin(r), end(r));
		if (int(s.size()) == n) return;

		int a = 0, b = 1, c = 2, d = 3;
		while ((s[b]-s[a]) * ((s[c]-s[a])%(s[d]-s[a])) == 0) d++;
		swap(s[d],s[3]);
		
		n = 3, volume = 0, p = vector{s[a],s[b],s[c]};
		faces.clear();
		for (int i = 0; i < N; i++) mark[i].reset();
		add_face(0,1,2), add_face(0,2,1);
		for (int i = 3; i < int(s.size()); i++) {
			add_point(s[i]);
		}
	}

	convex_hull(point a, point b, point c, point d): n(3), volume(0), p({a,b,c}) {
		add_face(0,1,2), add_face(0,2,1);
		add_point(d);
	}
};

mt19937_64 gen(987);

vector<point> solve_random(int a, int b, int c, i64 v) {
	auto random_point = [&]() -> point {
		using uid = uniform_int_distribution<i64>;
		return {uid(0,a)(gen), uid(0,b)(gen), uid(0,c)(gen)};
	};
	while (true) {
		point p = random_point();
		point q = random_point();
		point r = random_point();
		point s = random_point();
		if ((q-p) * ((r-p) % (s-p)) == 0) continue;
		convex_hull C(p,q,r,s);
		int stuck = 0;
		while (C.n < N && abs(C.volume) < v) {
			convex_hull D = C;
			if (C.add_point(random_point())) {
				if (abs(C.volume) > v) {
					C = D;
					if (++stuck == 100) break;
				} else {
					stuck = 0;
				} 
			}
		}
		if (abs(C.volume) == v) return C.p;
	}
	return {};
}

int main() {
	i64 a, b, c, v;
	cin >> a >> b >> c >> v;
	
	a = (v-1)/(6LL*b*c) + 1;
	b = (v-1)/(6LL*a*c) + 1;
	c = (v-1)/(6LL*a*b) + 1;

	auto res = solve_random(a,b,c,v);

	cout << res.size() << endl;
	for (point p: res) {
		cout << p.x << " " << p.y << " " << p.z << endl;
	}
}