include iostream include vector include algorithm using namespace std

#include <iostream>
#include <vector>
#include <algorithm>

using namespace std;

#define mp make_pair
#define pb push_back

#define all(x) x.begin(),x.end()

int m, n, d;

struct v2 {
    int x, y;
	v2(int x_, int y_): x(x_), y(y_) {}

	double fi() const { return atan2(y, x); }
};

inline bool operator < (const v2 _v1, const v2 _v2) { return _v1.fi() < _v2.fi(); }


typedef unsigned int uint;
inline vector<uint> mult(vector<uint> l, vector<uint> r){
    vector<uint> ans(4, 0);
    ans[0] = l[0] * r[0] + l[1] * r[2];
    ans[1] = l[0] * r[1] + l[1] * r[3];
    ans[2] = l[2] * r[0] + l[3] * r[2];
    ans[3] = l[2] * r[1] + l[3] * r[3];
    return ans;
}


inline vector<uint> binpow(vector<uint> a, int n){
    vector<uint> ans(4, 0);
    ans[0] = ans[3] = 1;
    while (n){
        if (n & 1) ans = mult(ans, a);
        a = mult(a, a);
        n >>= 1;
	}
    return ans;
}

uint s0;
inline uint gets(int i, int j) {
    int a = i * n + j + 1;
    vector<uint> ar (4, 0);
    ar[0] = 1664525;
    ar[1] = 1013904223;
    ar[3] = 1;
	ar = binpow(ar, a);
    return (ar[0] * s0 + ar[1]) % 3;
}

vector<v2> v2s;

const double pi = 2 * acos(0);
const double fieps = pi / 5;

const pair<int, int> undef = mp(-2, -2);
int d2, EPS2;

inline long long sqr(long long a) {
    return a * a;
}

void init() {
    for(int x = -d; x <= d; x++) {
        for(int y = -d; y <= d; y++) {
            int R2 = sqr(x) + sqr(y);
            if(R2 >= d2 - EPS2 && R2 <= d2) v2s.pb(v2(x, y));
		}
	}
	sort(all(v2s));
}

inline long long dist2(int i1, int j1, int i2, int j2) {
    return sqr(1LL * i1 - i2) + sqr(1LL * j1 - j2);
}

int main() {
    cin >>m >>n >>d >>s0;
    d2 = sqr(d), EPS2 = d2 / 2;
    init();
    vector< pair<int, int> > ans;
    int i = 0, j = 0;
    while(dist2(i, j, m - 1, n - 1) > d2) {
        ans.pb(mp(i + 1, j + 1));
        double curfi = atan2(n - j, m - i);
        uint sij = gets(i, j);
        bool found = false;
        v2 toFind(cos(curfi - fieps) * 1000000, sin(curfi - fieps) * 1000000);
        for(auto it = lower_bound(all(v2s), toFind); it != v2s.end() && it->fi() < curfi + fieps; it++) {
            int i_ = i + it->x, j_ = j + it->y;
            if(i_ >= 0 && i_ < m && j_ >= 0 && j_ < n && gets(i_, j_) == (sij + 1) % 3) {
                i = i_, j = j_;
                found = true;
                break;
			}
		}
        if(!found) {
            cout <<"ERROR POINT NOT FOUND!" <<endl;
            exit(1);
		}
	}
    ans.pb(mp(i + 1, j + 1));

    uint stateDist = (1LL * gets(m - 1, n - 1) - gets(i, j) + 3) % 3;

    if(stateDist == 0) {
        uint s1 = (1LL * gets(i, j) + 1) % 3;
        pair<int, int> p1 = undef, p2 = undef;
        uint s2 = (1LL * gets(i, j) + 2) % 3;
        for(int x = m - d / 2; x < m; x++) {
            for(int y = n - d / 2; y < n; y++) {
                if(p1 == undef && gets(x, y) == s1) p1 = mp(x, y);
                if(p2 == undef && gets(x, y) == s2) p2 = mp(x, y);
			}
		}
        if(p1 == undef || p2 == undef) exit(1);
        ans.pb(mp(p1.first + 1, p1.second + 1));
        ans.pb(mp(p2.first + 1, p2.second + 1));
	} else if(stateDist == 2) {
        uint s1 = (1LL * gets(i, j) + 1) % 3;
        pair<int, int> p1 = undef;
        for(int x = m - d / 2; x < m; x++) {
            for(int y = n - d / 2; y < n; y++) {
                if(p1 == undef && gets(x, y) == s1) p1 = mp(x, y);
			}
		}

        if(p1 == undef) exit(1);
        ans.pb(mp(p1.first + 1, p1.second + 1));
	}

    ans.pb(mp(m, n));

    cout <<ans.size() - 1 <<endl;

    for(int i = 0; i < ans.size(); i++) {
        cout <<ans[i].first <<" " <<ans[i].second <<endl;
	}

    double mxmoves = 2 * sqrt(sqr(n-1)+sqr(m-1)) / d + 3;
    if(ans.size() - 1 > mxmoves) {
        exit(1);
	}

    return 0;
}