#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}; 
}

const int N = 100;

struct convex_hull {
	struct face {
		int a, b, c;
		point normal;
	};
	
	int n;
	i64 volume = 0;
	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);
		}
		n++;
		return true;
	}

	convex_hull(point a, point b, point c, point d): n(3), 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);
		while (C.n < N && abs(C.volume) < v) {
			C.add_point(random_point());
		}
		if (abs(C.volume) == v) return C.p;
	}
	return {};
}

tuple<i64,i64,i64> largest_cut(i64 a, i64 b, i64 c, i64 r) {
	if (a*b*c <= r) return {a,b,c};
	if (a < b) { auto [y,x,z] = largest_cut(b,a,c,r); return {x,y,z}; }
	if (a < c) { auto [z,y,x] = largest_cut(c,b,a,r); return {x,y,z}; }
	if (b < c) { auto [x,z,y] = largest_cut(a,c,b,r); return {x,y,z}; }
	for(; ; r--) {
		for (i64 x = a; x*x*x >= r; x--) if (r%x == 0) {
			for (i64 y = min(b,x); x*y*y >= r; y--) if ((r/x)%y == 0) {
				i64 z = r/x/y;
				if (z <= c) return {x,y,z};
			}
		}
	}
}

vector<point> solve(int a, int b, int c, i64 v) {
	a = (v-1)/(6LL*b*c) + 1;
	b = (v-1)/(6LL*a*c) + 1;
	c = (v-1)/(6LL*a*b) + 1;

	i64 r = 6LL*a*b*c - v;

	vector<point> p;
	auto corner = [&](int x, int y, int z, int u, int v, int w) {
		p.emplace_back(x > 0 ? x-u : u, y, z);
		p.emplace_back(x, y > 0 ? y-v : v, z);
		p.emplace_back(x, y, z > 0 ? z-w : w);
	};

	for (int mask = 0; mask < 8; mask++) {
		int x = mask & 1 ? a : 0;
		int y = mask & 2 ? b : 0;
		int z = mask & 4 ? c : 0;
		int aa = mask & 1 ? a/2 : (a+1)/2;
		int bb = mask & 2 ? b/2 : (b+1)/2;
		int cc = mask & 4 ? c/2 : (c+1)/2;
		auto [u,v,w] = largest_cut(aa,bb,cc,r);
		r -= u*v*w;
		corner(x,y,z,u,v,w);
	}
	
	if (r == 0) return p;
	return solve_random(a,b,c,v);
}

int main() {
	i64 a, b, c, v;
	cin >> a >> b >> c >> v;
	auto res = solve(a,b,c,v);

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