#include <bits/stdc++.h>
using namespace std;
#define rep(i,a,b) for (auto i=(a); i<(b); ++i)
#define iter(it,c) for (auto it = (c).begin(); it != (c).end(); ++it)
typedef pair<int, int> ii;
typedef vector<int> vi;
typedef vector<ii> vii;
typedef vector<vi> vvi;
typedef long long ll;
const int INF = 2147483647;
const ll DINF = 9000000000000000000LL;

const double pi = acos(-1.0);
const double EPS = 1e-9;

typedef complex<double> point;
double cross(point a, point b) { return imag(conj(a) * b); }
double ccw(point a, point b, point c) { return cross(b - a, c - b); }

#define MAXN 5000
point hull[MAXN];
bool cmp(const point &a, const point &b) {
  return abs(real(a) - real(b)) > EPS ?
    real(a) < real(b) : imag(a) < imag(b); }
int convex_hull(vector<point> p) {
  int n = size(p), l = 0;
  sort(p.begin(), p.end(), cmp);
  rep(i,0,n) {
    if (i > 0 && p[i] == p[i - 1]) continue;
    while (l >= 2 &&
        ccw(hull[l - 2], hull[l - 1], p[i]) >= 0) l--;
    hull[l++] = p[i]; }
  int r = l;
  for (int i = n - 2; i >= 0; i--) {
    if (p[i] == p[i + 1]) continue;
    while (r - l >= 1 &&
        ccw(hull[r - 2], hull[r - 1], p[i]) >= 0) r--;
    hull[r++] = p[i]; }
  return l == 1 ? 1 : r - 1; }

template <class K> bool eq(K a, K b) { return a == b; }
template <> bool eq<double>(double a, double b) {
    return abs(a - b) < EPS; }
template <class T> struct matrix {
  int rows, cols, cnt; vector<T> data;
  inline T& at(int i, int j) { return data[i * cols + j]; }
  matrix(int r, int c) : rows(r), cols(c), cnt(r * c) {
    data.assign(cnt, T(0)); }
  matrix(const matrix& other) : rows(other.rows),
    cols(other.cols), cnt(other.cnt), data(other.data) { }
  T& operator()(int i, int j) { return at(i, j); }
  matrix<T> rref(T &det, int &rank) {
    matrix<T> mat(*this); det = T(1), rank = 0;
    for (int r = 0, c = 0; c < cols; c++) {
      int k = r;
      rep(i,k+1,rows) if (abs(mat(i,c)) > abs(mat(k,c))) k = i;
      if (k >= rows || eq<T>(mat(k, c), T(0))) continue;
      if (k != r) {
        det *= T(-1);
        rep(i,0,cols) swap(mat.at(k, i), mat.at(r, i));
      } det *= mat(r, r); rank++;
      T d = mat(r,c);
      rep(i,0,cols) mat(r, i) /= d;
      rep(i,0,rows) {
        T m = mat(i, c);
        if (i != r && !eq<T>(m, T(0)))
          rep(j,0,cols) mat(i, j) -= m * mat(r, j);
      } r++;
    } return mat; }
};

#define P(p) const point3d &p
#define L(p0, p1) P(p0), P(p1)
#define PL(p0, p1, p2) P(p0), P(p1), P(p2)
struct point3d {
  double x, y, z;
  point3d() : x(0), y(0), z(0) {}
  point3d(double _x, double _y, double _z)
    : x(_x), y(_y), z(_z) {}
  point3d operator+(P(p)) const {
    return point3d(x + p.x, y + p.y, z + p.z); }
  point3d operator-(P(p)) const {
    return point3d(x - p.x, y - p.y, z - p.z); }
  point3d operator-() const {
    return point3d(-x, -y, -z); }
  point3d operator*(double k) const {
    return point3d(x * k, y * k, z * k); }
  point3d operator/(double k) const {
    return point3d(x / k, y / k, z / k); }
  double operator%(P(p)) const {
    return x * p.x + y * p.y + z * p.z; }
  point3d operator*(P(p)) const {
    return point3d(y*p.z - z*p.y,
                   z*p.x - x*p.z, x*p.y - y*p.x); }
  double length() const {
    return sqrt(*this % *this); }
  double distTo(P(p)) const {
    return (*this - p).length(); }
  double distTo(P(A), P(B)) const {
    // A and B must be two different points
    return ((*this - A) * (*this - B)).length() / A.distTo(B);}
  double signedDistTo(PL(A,B,C)) const {
    // A, B and C must not be collinear
    point3d N = (B-A)*(C-A); double D = A%N;
    return ((*this)%N - D)/N.length(); }
  bool isZero() const {
    return abs(x) < EPS && abs(y) < EPS && abs(z) < EPS; }
  bool isOnLine(L(A, B)) const {
    return ((A - *this) * (B - *this)).isZero(); }
};

long double volume(vector<point3d> pts) {
    int n = pts.size();
    point3d inside(0,0,0);
    iter(it,pts) {
        inside = inside + *it;
    }
    inside = inside / n;
    long double vol = 0;
    set<vi> seen;
    rep(i,0,n) {
        rep(j,i+1,n) {
            rep(k,j+1,n) {

                point3d p = pts[j]-pts[i];
                point3d q = pts[k]-pts[i];
                if (abs((p*q).length()) < EPS) {
                    // collinear
                    continue;
                }

                vi face;
                int left = 0,
                    right = 0;
                rep(l,0,n) {
                    double d = pts[l].signedDistTo(pts[i], pts[j], pts[k]);
                    if (abs(d) < EPS) {
                        face.push_back(l);
                    } else if (d < 0) {
                        left++;
                    } else {
                        right++;
                    }
                }

                if (seen.find(face) == seen.end() && (left == 0 || right == 0)) {
                    seen.insert(face);

                    vector<point> plane;
                    point3d mid(0,0,0);
                    iter(it,face) {
                        mid = mid + pts[*it];
                        point3d t = pts[*it] - pts[i];
                        matrix<double> mat(3, 3);
                        mat(0,0) = p.x;
                        mat(0,1) = q.x;
                        mat(0,2) = t.x;
                        mat(1,0) = p.y;
                        mat(1,1) = q.y;
                        mat(1,2) = t.y;
                        mat(2,0) = p.z;
                        mat(2,1) = q.z;
                        mat(2,2) = t.z;
                        double det;
                        int rank;
                        mat = mat.rref(det, rank);
                        assert(rank == 2);
                        plane.emplace_back(mat(0,2), mat(1,2));
                    }
                    mid = mid / face.size();

                    int m = face.size();
                    point3d normal;
                    rep(a,0,m) {
                        rep(b,a+1,m) {
                            rep(c,b+1,m) {
                                normal = (pts[face[b]]-pts[face[a]])*(pts[face[c]]-pts[face[a]]);
                                if (abs(normal.length()) > EPS) {
                                    normal = normal / normal.length();
                                    goto found_normal;
                                }
                            }
                        }
                    }
                    assert(false);
                    found_normal:
                    if ((mid + normal * 0.001 - inside).length() < (mid - normal * 0.001 - inside).length()) {
                        normal = -normal;
                    }

                    int cnt = convex_hull(plane);
                    long double area = 0;
                    rep(l,1,cnt-1) {
                        int r = l+1;
                        point3d a = p * real(hull[0]) + q * imag(hull[0]);
                        point3d b = p * real(hull[l]) + q * imag(hull[l]);
                        point3d c = p * real(hull[r]) + q * imag(hull[r]);
                        area += ((b-a)*(c-a)).length() / 2.0;
                    }

                    vol += (mid % normal) * area;
                }
            }
        }
    }

    return vol / 3.0;
}

ll a, b, c, v;

ll rand_int(ll a, ll b) {
    if (a > b) assert(false);
    return rand() % (b-a+1) + a;
}

tuple<ll,ll,ll> get_ear(ll target, ll mnx, ll mny, ll mnz) {
    ll x = -1, y = -1, z = -1;

    ll st = min((ll)sqrt(target)+10LL, min(mnx,min(mny,mnz)));
    while (x == -1 || y == -1 || z == -1) {
        if (target % st == 0) {
            if (st <= mnx && x == -1) { x = st; target /= st; continue; }
            if (st <= mny && y == -1) { y = st; target /= st; continue; }
            if (st <= mnz && z == -1) { z = st; target /= st; continue; }
        }
        st--;
    }

    if (target != 1) {
        return {0,0,0};
    }

    return {x,y,z};
}

int main() {
    srand(time(NULL));

    cin >> a >> b >> c >> v;

    if (v < 1e9) {
        auto [x,y,z] = get_ear(v, a,b,c);
        if (x != 0) {
            cout << 4 << endl;
            cout << 0 << " " << 0 << " " << 0 << endl;
            cout << x << " " << 0 << " " << 0 << endl;
            cout << 0 << " " << y << " " << 0 << endl;
            cout << 0 << " " << 0 << " " << z << endl;
            return 0;
        }
    }

    while (true) {
        ll mn = DINF;

        ll mnx = -1, mny = -1, mnz = -1;

        while (true) {
            ll x = rand_int(max(1LL,v/6/b/c), a);
            if (v/6/x/c > b) continue;
            ll y = rand_int(max(1LL,v/6/x/c), b);
            if (v/6/x/y > c) continue;
            ll z = rand_int(max(1LL,v/6/x/y), c);
            if (6*x*y*z < v) continue;

            if (x*y*z < mn) {
                mn = x*y*z;
                mnx = x;
                mny = y;
                mnz = z;

                if (6*mn - v < 100000000000L) {
                    break;
                }
            }
        }

        assert(mnx != -1);

        if (mnx*mny*mnz < 100) {
            rep(it,0,100) {
                vector<point3d> pts;
                rep(jt,0,100) {
                    point3d pt(rand_int(0,mnx), rand_int(0,mny), rand_int(0,mnz));
                    bool on_boundary = false;
                    rep(i,0,pts.size()) {
                        if ((pt - pts[i]).isZero()) on_boundary = true;
                        rep(j,i+1,pts.size()) {
                            if (pt.isOnLine(pts[i], pts[j])) on_boundary = true;
                        }
                    }
                    if (on_boundary) continue;
                    pts.push_back(pt);
                    if (pts.size() < 4) continue;

                    double v2 = volume(pts)*6;
                    if (abs(v2 - v) < 1e-6) {
                        cout << pts.size() << endl;
                        iter(it,pts) {
                            cout << setprecision(0) << fixed << (*it).x << " " << (*it).y << " " << (*it).z << endl;
                        }
                        return 0;
                    }
                    if (v2 > v) {
                        break;
                    }
                }
            }
        }

        ll target = 6*mnx*mny*mnz - v;

        if (target == 0) {
            cout << 8 << endl;
            cout << 0 << " " << 0 << " " << 0 << endl;
            cout << mnx << " " << mny << " " << mnz << endl;
            cout << mnx << " " << mny << " " << 0 << endl;
            cout << mnx << " " << 0 << " " << mnz << endl;
            cout << 0 << " " << mny << " " << mnz << endl;
            cout << 0 << " " << 0 << " " << mnz << endl;
            cout << 0 << " " << mny << " " << 0 << endl;
            cout << mnx << " " << 0 << " " << 0 << endl;

            break;
        }

        auto [x,y,z] = get_ear(target, mnx,mny,mnz);

        if (x != 0) {
            cout << 10 << endl;
            cout << x << " " << 0 << " " << 0 << endl;
            cout << 0 << " " << y << " " << 0 << endl;
            cout << 0 << " " << 0 << " " << z << endl;
            cout << mnx << " " << mny << " " << mnz << endl;
            cout << mnx << " " << mny << " " << 0 << endl;
            cout << mnx << " " << 0 << " " << mnz << endl;
            cout << 0 << " " << mny << " " << mnz << endl;
            cout << 0 << " " << 0 << " " << mnz << endl;
            cout << 0 << " " << mny << " " << 0 << endl;
            cout << mnx << " " << 0 << " " << 0 << endl;

            break;
        }
    }

    return 0;
}