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