stablegradients/pie-gem5-pairs huggingface dataset:
// Slower C++ solution:
#include <bits/stdc++.h>

using namespace std;

const int N = 1414 ;

int n , m , answer = 0 , matrix [N][N] , dp [N][N] ;

int isval (int i , int j , int len )

{

	if ( len == 0 )

	{

		return 1 ;

	}

	int z = len ;

	for (int q = 1 ; q <= len ; q ++)

	{

		if ( matrix [i - q ][j] == 1 || matrix[i][j - q ] == 1)

		{

			z = q - 1 ; 

			break ;

		}

	}

 	return z + 1;

}

int main ()

{

	ios::sync_with_stdio(false), cin.tie(0);

	cin >> n >> m ; 

	for (int i = 0 ; i < n ; i ++)

	{

		for (int j = 0 ; j < m ; j ++)

		{

			cin >> matrix[j][i] ;

			if ( matrix[j][i] == 0 )

			{

				dp[j][i] = 1 ;

				answer = 1 ;

			}

		}

	}

	for (int i = 1 ; i < n ; i ++)

	{

		for (int j = 1 ; j < m ; j ++)

		{

			if ( matrix[j][i] == 0)

			{

				dp[j][i] = isval ( j , i , dp[j - 1][i - 1 ] ) ;

			}

			answer = max ( answer , dp[j][i] ) ;

		}

	}

	answer *= answer ;

	cout << answer << endl ; 

	/*cerr << "*********************************" << endl ;

	for (int i = 0 ; i < n ; i ++)

	{

		for (int j = 0 ; j < m ; j ++)

		{

			cerr << dp[j][i] << " " ;

		}

		cerr << endl ;

	}*/

	return  0 ; 

}

// Faster C++ solution for the same problem (measured speedup: 5.2467358603577905x):
#include<cstdio>

#include<cstring>

#include<algorithm>

#define MAXN 1400

#define INF 0x3fffffff

using namespace std;



int n,m,ans=0;

int d[MAXN+5][MAXN+5];



int main()

{

    scanf("%d%d",&n,&m);

    for(int i=1;i<=n;i++)

        for(int j=1;j<=m;j++)

        {

            scanf("%d",&d[i][j]);

            d[i][j]=!d[i][j];

            ans=max(ans,d[i][j]);

        }

    for(int i=2;i<=n;i++)

        for(int j=2;j<=m;j++)

            if(d[i][j])

            {

                d[i][j]=min(d[i-1][j],min(d[i][j-1],d[i-1][j-1]))+1;

                ans=max(ans,d[i][j]);

            }

    /*for(int i=1;i<=n;i++)

    {

        for(int j=1;j<=m;j++)

            printf("%d  ",d[i][j]);

        printf("\n");

    }*/

    printf("%d\n",ans*ans);

}
// Slower C++ solution:
#include <bits/stdc++.h>

using namespace std;

const int N = 1414 ;

int n , m , answer = 0 , matrix [N][N] , dp [N][N] ;

int isval (int i , int j , int len )

{

	if ( len == 0 )

	{

		return 1 ;

	}

	int z = len ;

	for (int q = 1 ; q <= len ; q ++)

	{

		if ( matrix [i - q ][j] == 1 || matrix[i][j - q ] == 1)

		{

			z = q - 1 ; 

			break ;

		}

	}

 	return z + 1;

}

int main ()

{

	ios::sync_with_stdio(false), cin.tie(0);

	cin >> n >> m ; 

	for (int i = 0 ; i < n ; i ++)

	{

		for (int j = 0 ; j < m ; j ++)

		{

			cin >> matrix[j][i] ;

			if ( matrix[j][i] == 0 )

			{

				dp[j][i] = 1 ;

				answer = 1 ;

			}

		}

	}

	for (int i = 1 ; i < n ; i ++)

	{

		for (int j = 1 ; j < m ; j ++)

		{

			if ( matrix[j][i] == 0)

			{

				dp[j][i] = isval ( j , i , dp[j - 1][i - 1 ] ) ;

			}

			answer = max ( answer , dp[j][i] ) ;

		}

	}

	answer *= answer ;

	cout << answer << endl ; 

	/*cerr << "*********************************" << endl ;

	for (int i = 0 ; i < n ; i ++)

	{

		for (int j = 0 ; j < m ; j ++)

		{

			cerr << dp[j][i] << " " ;

		}

		cerr << endl ;

	}*/

	return  0 ; 

}

// Faster C++ solution for the same problem (measured speedup: 5.249226981472117x):
#include <stdio.h>

#include <algorithm>



int h, w;

int a[1401][1401];

int dp[1401][1401];



int main()

{

	scanf("%d%d", &h, &w);

	int res = 0;

	for (int i = 0; i < h; i++)

		for (int j = 0; j < w; j++)

		{

			scanf("%d", &a[i][j]);

			dp[i][j] = !a[i][j];

			res |= dp[i][j];

		}



	for (int i = 1; i < h; i++)

		for (int j = 1; j < w; j++)

		{

			if (a[i][j] == 0)

				dp[i][j] = std::min(std::min(dp[i - 1][j - 1], dp[i - 1][j]), dp[i][j - 1]) + 1;

			res = std::max(res, dp[i][j]);

		}



	printf("%d\n", res * res);



	return 0;

}
// Slower C++ solution:
#include <bits/stdc++.h>

using namespace std;

const int N = 1414 ;

int n , m , answer = 0 , matrix [N][N] , dp [N][N] ;

int isval (int i , int j , int len )

{

	if ( len == 0 )

	{

		return 1 ;

	}

	int z = len ;

	for (int q = 1 ; q <= len ; q ++)

	{

		if ( matrix [i - q ][j] == 1 || matrix[i][j - q ] == 1)

		{

			z = q - 1 ; 

			break ;

		}

	}

 	return z + 1;

}

int main ()

{

	ios::sync_with_stdio(false), cin.tie(0);

	cin >> n >> m ; 

	for (int i = 0 ; i < n ; i ++)

	{

		for (int j = 0 ; j < m ; j ++)

		{

			cin >> matrix[j][i] ;

			if ( matrix[j][i] == 0 )

			{

				dp[j][i] = 1 ;

				answer = 1 ;

			}

		}

	}

	for (int i = 1 ; i < n ; i ++)

	{

		for (int j = 1 ; j < m ; j ++)

		{

			if ( matrix[j][i] == 0)

			{

				dp[j][i] = isval ( j , i , dp[j - 1][i - 1 ] ) ;

			}

			answer = max ( answer , dp[j][i] ) ;

		}

	}

	answer *= answer ;

	cout << answer << endl ; 

	/*cerr << "*********************************" << endl ;

	for (int i = 0 ; i < n ; i ++)

	{

		for (int j = 0 ; j < m ; j ++)

		{

			cerr << dp[j][i] << " " ;

		}

		cerr << endl ;

	}*/

	return  0 ; 

}

// Faster C++ solution for the same problem (measured speedup: 5.245256733466173x):
#include <cstdio>

#include <cstring>

#include <cstdlib>

#include <cmath>

#include <algorithm>

#define inf 0x3f3f3f3f

#define inf64 0x3f3f3f3f3f3f3f3f

using namespace std;

typedef long long ll;

const int maxn = 2e3+10;

int dp[maxn][maxn],a[maxn][maxn];



int main(){

    int n,m,ans=0;

    scanf("%d%d",&n,&m);

    for(int i=1;i<=n;i++){

        for(int j=1;j<=m;j++) scanf("%d",&a[i][j]);

    }

    for(int i=1;i<=n;i++){

        for(int j=m;j>=1;j--){

            if(a[i][j]) dp[i][j]=0;

            else dp[i][j]=min(min(dp[i-1][j],dp[i][j+1]),dp[i-1][j+1])+1;

            ans=max(ans,dp[i][j]);

//            printf("dp[%d][%d]=%d\n",i,j,dp[i][j]);

        }

    }

    printf("%d\n",ans*ans);

    return 0;

}
// Slower C++ solution:
#include<bits/stdc++.h>

using namespace std;



#define ALL(a)  a.begin(),a.end()

#define SZ(a) int((a).size())

#define SORT(c) sort(ALL((c)))

#define UNIQ(c) (c).erase(unique(ALL((c))), end((c)))

#define FOR(i, s, e) for (int(i) = (s); (i) < (e); (i)++)

#define debug(x) cerr << #x << ": " << x << endl



using LL = long long; template<class T> using V = vector<T>;





int main() {

	cin.tie(0);

	ios_base::sync_with_stdio(false);



	int H, W;

	cin >> H >> W;

	V<V<int>>a(H + 1, V<int>(W + 1));

	FOR(i, 0, H) {

		FOR(j, 0, W) {

			cin >> a[i + 1][j + 1];

		}

	}

	LL ans = 0;

	V<V<int>>dp(H + 1, V<int>(W + 1, 0));

	FOR(i, 1, H + 1) {

		FOR(j, 1, W + 1) {

			if (a[i][j] == 1)dp[i][j] = 0;

			else dp[i][j] = min(dp[i - 1][j - 1], min(dp[i - 1][j], dp[i][j - 1])) + 1;

		}

	}



	FOR(i, 1, H + 1) {

		FOR(j, 1, W + 1) {

			ans = max(ans, (LL)dp[i][j]);

		}

	}

	cout << ans * ans << endl;



	return 0;

}


// Faster C++ solution for the same problem (measured speedup: 5.30061216463847x):
#include <cstdio>

#include <cstring>

#include <cstdlib>

#include <cmath>

#include <algorithm>

#define inf 0x3f3f3f3f

#define inf64 0x3f3f3f3f3f3f3f3f

using namespace std;

typedef long long ll;

const int maxn = 2e3+10;

int dp[maxn][maxn],a[maxn][maxn];



int main(){

    int n,m,ans=0;

    scanf("%d%d",&n,&m);

    for(int i=1;i<=n;i++){

        for(int j=1;j<=m;j++) scanf("%d",&a[i][j]);

    }

    for(int i=1;i<=n;i++){

        for(int j=m;j>=1;j--){

            if(a[i][j]) dp[i][j]=0;

            else dp[i][j]=min(min(dp[i-1][j],dp[i][j+1]),dp[i-1][j+1])+1;

            ans=max(ans,dp[i][j]);

//            printf("dp[%d][%d]=%d\n",i,j,dp[i][j]);

        }

    }

    printf("%d\n",ans*ans);

    return 0;

}
// Slower C++ solution:
#include <iostream>

#include <algorithm>

using namespace std;

int dp[1401][1401], G[1401][1401];



int main()

{

	int H, W;

	cin >> H >> W;

	for (int i = 0; i < H; ++i)

		for (int j = 0; j < W; ++j)

			cin >> G[i][j];

	int	maxWidth = 0;

	for (int i = 0; i < H; ++i)

		for (int j = 0; j < W; ++j)

		{

			dp[i][j] = (G[i][j] == 1 ? 0 : 1);

			if (maxWidth < dp[i][j])

				maxWidth = dp[i][j];

		}

		

	for (int i = 1; i < H; ++i)

		for (int j = 1; j < W; ++j)

		{

			if (G[i][j])

				dp[i][j] = 0;

			else

			{

				dp[i][j] = min(dp[i - 1][j - 1], min(dp[i - 1][j], dp[i][j - 1])) + 1;

				maxWidth = max(maxWidth, dp[i][j]);

			}

		}

	cout << maxWidth * maxWidth << endl;

	return 0;

}


// Faster C++ solution for the same problem (measured speedup: 5.194001145980882x):
#include <cstdio>

#include <cstring>

#include <cstdlib>

#include <cmath>

#include <algorithm>

#define inf 0x3f3f3f3f

#define inf64 0x3f3f3f3f3f3f3f3f

using namespace std;

typedef long long ll;

const int maxn = 2e3+10;

int dp[maxn][maxn],a[maxn][maxn];



int main(){

    int n,m,ans=0;

    scanf("%d%d",&n,&m);

    for(int i=1;i<=n;i++){

        for(int j=1;j<=m;j++) scanf("%d",&a[i][j]);

    }

    for(int i=1;i<=n;i++){

        for(int j=m;j>=1;j--){

            if(a[i][j]) dp[i][j]=0;

            else dp[i][j]=min(min(dp[i-1][j],dp[i][j+1]),dp[i-1][j+1])+1;

            ans=max(ans,dp[i][j]);

//            printf("dp[%d][%d]=%d\n",i,j,dp[i][j]);

        }

    }

    printf("%d\n",ans*ans);

    return 0;

}
// Slower C++ solution:
#include<iostream>

#include<algorithm>

using namespace std;



const int MAX = 1400;



int minSquare(int tile[MAX][MAX], int x, int y) {

	if (tile[y][x] == 1) {

		tile[y][x] = 0;

		return 0;

	}

	return min(tile[y - 1][x - 1], min(tile[y - 1][x], tile[y][x - 1])) + 1;

}



int main() {

	int H, W;

	int tile[MAX][MAX];

	int ans = 0;





	cin >> H >> W;

	for (int i = 0; i < H; i++) {

		for (int j = 0; j < W; j++) {

			cin >> tile[i][j];

		}

	}



	for (int i = 0; i < H; i++) {

		if (tile[i][0] == 1)tile[i][0] = 0;

		else {

			tile[i][0] = 1;

			ans = 1;

		}

	}

	for (int i = 1; i < W; i++) {

		if (tile[0][i] == 1)tile[0][i] = 0;

		else {

			tile[0][i] = 1;

			ans = 1;

		}

	}



	int x = 1;

	int y = 1;

	int n = min(H, W);

	for (int i = 1; i < n; i++) {

		for (int j = y; j < H; j++) {

			tile[j][x] = minSquare(tile, x, j);

			ans = max(tile[j][x], ans);

		}

		for (int j = x + 1; j < W; j++) {

			tile[y][j] = minSquare(tile, j, y);

			ans = max(tile[y][j], ans);

		}

		x++;

		y++;

	}



	cout << ans * ans << endl;

	return 0;

}

// Faster C++ solution for the same problem (measured speedup: 3.1572286564662737x):
#include<stdio.h>

#include<algorithm>

using namespace std;



const int MAX = 1400;



int main() {

	int H, W;

	int tile[MAX][MAX];

	int ans = 0;



	scanf("%d%d", &H, &W);

	for (int i = 0; i < H; i++) {

		for (int j = 0; j < W; j++) {

			scanf("%d", &tile[i][j]);

		}

	}



	for (int y = 0; y < H; y++) {

		for (int x = 0; x < W; x++) {

			if (tile[y][x] == 1) tile[y][x] = 0;

			else if (x < 1 || y < 1) tile[y][x] = 1;

			else tile[y][x] = min(tile[y - 1][x - 1], min(tile[y - 1][x], tile[y][x - 1])) + 1;

			ans = max(tile[y][x], ans);

		}

	}



	printf("%d\n", ans * ans);

	return 0;

}
// Slower C++ solution:
#include <iostream>

using namespace std;

void func(int,int,double,int);

int p,q,a,n,cnt;

double k;

 

int main(){

  while(1){

    cin>>p>>q>>a>>n;

    if(!p&&!q&&!a&&!n) break;

    k=1.0*p/q;

    cnt=0;

    func(0,1,0,1);

    cout<<cnt<<endl;

  }

  return 0;

}

 

void func(int x,int w,double y,int I){

  if(-0.00000001<=k-y&&k-y<=0.00000001) cnt++;

  if(x==n) return ;

  for(int i=I;w*i<=a;i++){

    if(y+1.0/i>k+0.1) continue;

    func(x+1,w*i,1.0/i+y,i);

  }

}

// Faster C++ solution for the same problem (measured speedup: 1.0799527325793312x):
#include <iostream>

using namespace std;

void func(int,int,double,int);

int p,q,a,n,cnt;

double k;

 

int main(){

  while(1){

    cin>>p>>q>>a>>n;

    if(!p&&!q&&!a&&!n) break;

    k=1.0*p/q;

    cnt=0;

    func(0,1,0,1);

    cout<<cnt<<endl;

  }

  return 0;

}

 

void func(int x,int w,double y,int I){

  if(-0.00000001<=k-y&&k-y<=0.00000001) cnt++;

  if(x==n) return ;

  for(int i=I;w*i<=a;i++){

    if(y+1.0/i>k+0.0000001) continue;

    func(x+1,w*i,1.0/i+y,i);

  }

}
// Slower C++ solution:
#include <bits/stdc++.h>

using namespace std;



#define rep(i,x,y) for(int i=(x);i<(y);++i)

#define debug(x) #x << "=" << (x)



#ifdef DEBUG

#define _GLIBCXX_DEBUG

#define dump(x) std::cerr << debug(x) << " (L:" << __LINE__ << ")" << std::endl

#else

#define dump(x)

#endif



typedef long long int ll;

typedef pair<int,int> pii;

//template<typename T> using vec=std::vector<T>;



const int inf=1<<30;

const long long int infll=1LL<<58;

const double eps=1e-9;

const int dx[]={1,0,-1,0},dy[]={0,1,0,-1};



template <typename T> ostream &operator<<(ostream &os, const vector<T> &vec){

	os << "[";

	for (const auto &v : vec) {

		os << v << ",";

	}

	os << "]";

	return os;

}



int p,q,a,n;



int gcd(int x,int y){

    if(y==0) return x;

    return gcd(y,x%y);

}



int dfs(int denom,int numer,int prev,int product,int num){

    int res=0;

    if((double)p/q-eps<=(double)numer/denom  and (double)numer/denom<=(double)p/q+eps) ++res;

    if(num==0) return res;



    {

        int tmp1=denom*prev,tmp2=numer*prev+denom*num;

        if((double)tmp2/tmp1<(double)p/q) return 0;

    }



    if((double)numer/denom>(double)p/q) return 0;



    int x=prev;

    while(true){

        if(product*x>a) break;

        res+=dfs(denom*x,numer*x+denom,x,product*x,num-1);

        ++x;

    }

    return res;

}



void solve(){

    while(cin >> p >> q >> a >> n){

        if(p==0 and q==0 and a==0 and n==0) break;

        cout << dfs(1,0,1,1,n) << endl;

    }

}



int main(){

	std::ios::sync_with_stdio(false);

	std::cin.tie(0);

	solve();

	return 0;

}

// Faster C++ solution for the same problem (measured speedup: 1.3135775431752694x):
#include <bits/stdc++.h>

using namespace std;



#define rep(i,x,y) for(int i=(x);i<(y);++i)

#define debug(x) #x << "=" << (x)



#ifdef DEBUG

#define _GLIBCXX_DEBUG

#define dump(x) std::cerr << debug(x) << " (L:" << __LINE__ << ")" << std::endl

#else

#define dump(x)

#endif



typedef long long int ll;

typedef pair<int,int> pii;

//template<typename T> using vec=std::vector<T>;



const int inf=1<<30;

const long long int infll=1LL<<58;

const double eps=1e-9;

const int dx[]={1,0,-1,0},dy[]={0,1,0,-1};



template <typename T> ostream &operator<<(ostream &os, const vector<T> &vec){

	os << "[";

	for (const auto &v : vec) {

		os << v << ",";

	}

	os << "]";

	return os;

}



int p,q,a,n;



int dfs(int denom,int numer,int prev,int product,int num){

    {

        int tmp1=denom*prev,tmp2=numer*prev+denom*num;

        if(tmp2*q<p*tmp1) return 0;

    }

    if(numer*q>p*denom) return 0;



    int res=0;

    if(numer*q==p*denom) ++res;

    if(num==0) return res;



    int x=prev;

    while(true){

        if(product*x>a) break;

        res+=dfs(denom*x,numer*x+denom,x,product*x,num-1);

        ++x;

    }

    return res;

}



void solve(){

    while(cin >> p >> q >> a >> n){

        if(p==0 and q==0 and a==0 and n==0) break;

        cout << dfs(1,0,1,1,n) << endl;

    }

}



int main(){

	std::ios::sync_with_stdio(false);

	std::cin.tie(0);

	solve();

	return 0;

}
// Slower C++ solution:
#include <bits/stdc++.h>

using namespace std;



#define rep(i,x,y) for(int i=(x);i<(y);++i)

#define debug(x) #x << "=" << (x)



#ifdef DEBUG

#define _GLIBCXX_DEBUG

#define dump(x) std::cerr << debug(x) << " (L:" << __LINE__ << ")" << std::endl

#else

#define dump(x)

#endif



typedef long long int ll;

typedef pair<int,int> pii;

//template<typename T> using vec=std::vector<T>;



const int inf=1<<30;

const long long int infll=1LL<<58;

const double eps=1e-9;

const int dx[]={1,0,-1,0},dy[]={0,1,0,-1};



template <typename T> ostream &operator<<(ostream &os, const vector<T> &vec){

	os << "[";

	for (const auto &v : vec) {

		os << v << ",";

	}

	os << "]";

	return os;

}



int p,q,a,n;



int gcd(int x,int y){

    if(y==0) return x;

    return gcd(y,x%y);

}



int dfs(int denom,int numer,int prev,int product,int num){

    {

        int tmp1=denom*prev,tmp2=numer*prev+denom*num;

        if((double)tmp2/tmp1+eps<(double)p/q) return 0;

    }

    if((double)numer/denom-eps>(double)p/q) return 0;



    int res=0;

    if((double)p/q-eps<=(double)numer/denom  and (double)numer/denom<=(double)p/q+eps) ++res;

    if(num==0) return res;



    int x=prev;

    while(true){

        if(product*x>a) break;

        res+=dfs(denom*x,numer*x+denom,x,product*x,num-1);

        ++x;

    }

    return res;

}



void solve(){

    while(cin >> p >> q >> a >> n){

        if(p==0 and q==0 and a==0 and n==0) break;

        cout << dfs(1,0,1,1,n) << endl;

    }

}



int main(){

	std::ios::sync_with_stdio(false);

	std::cin.tie(0);

	solve();

	return 0;

}

// Faster C++ solution for the same problem (measured speedup: 1.208915788019769x):
#include <bits/stdc++.h>

using namespace std;



#define rep(i,x,y) for(int i=(x);i<(y);++i)

#define debug(x) #x << "=" << (x)



#ifdef DEBUG

#define _GLIBCXX_DEBUG

#define dump(x) std::cerr << debug(x) << " (L:" << __LINE__ << ")" << std::endl

#else

#define dump(x)

#endif



typedef long long int ll;

typedef pair<int,int> pii;

//template<typename T> using vec=std::vector<T>;



const int inf=1<<30;

const long long int infll=1LL<<58;

const double eps=1e-9;

const int dx[]={1,0,-1,0},dy[]={0,1,0,-1};



template <typename T> ostream &operator<<(ostream &os, const vector<T> &vec){

	os << "[";

	for (const auto &v : vec) {

		os << v << ",";

	}

	os << "]";

	return os;

}



int p,q,a,n;



int dfs(int denom,int numer,int prev,int product,int num){

    {

        int tmp1=denom*prev,tmp2=numer*prev+denom*num;

        if(tmp2*q<p*tmp1) return 0;

    }

    if(numer*q>p*denom) return 0;



    int res=0;

    if(numer*q==p*denom) ++res;

    if(num==0) return res;



    int x=prev;

    while(true){

        if(product*x>a) break;

        res+=dfs(denom*x,numer*x+denom,x,product*x,num-1);

        ++x;

    }

    return res;

}



void solve(){

    while(cin >> p >> q >> a >> n){

        if(p==0 and q==0 and a==0 and n==0) break;

        cout << dfs(1,0,1,1,n) << endl;

    }

}



int main(){

	std::ios::sync_with_stdio(false);

	std::cin.tie(0);

	solve();

	return 0;

}
// Slower C++ solution:
#include <bits/stdc++.h>

using namespace std;



#define rep(i,x,y) for(int i=(x);i<(y);++i)

#define debug(x) #x << "=" << (x)



#ifdef DEBUG

#define _GLIBCXX_DEBUG

#define dump(x) std::cerr << debug(x) << " (L:" << __LINE__ << ")" << std::endl

#else

#define dump(x)

#endif



typedef long long int ll;

typedef pair<int,int> pii;

//template<typename T> using vec=std::vector<T>;



const int inf=1<<30;

const long long int infll=1LL<<58;

const double eps=1e-9;

const int dx[]={1,0,-1,0},dy[]={0,1,0,-1};



template <typename T> ostream &operator<<(ostream &os, const vector<T> &vec){

	os << "[";

	for (const auto &v : vec) {

		os << v << ",";

	}

	os << "]";

	return os;

}



int p,q,a,n;



int dfs(int denom,int numer,int prev,int product,int num){

    {

        int tmp1=denom*prev,tmp2=numer*prev+denom*num;

        if((double)tmp2/tmp1+eps<(double)p/q) return 0;

    }

    if((double)numer/denom-eps>(double)p/q) return 0;



    int res=0;

    if((double)p/q-eps<=(double)numer/denom  and (double)numer/denom<=(double)p/q+eps) ++res;

    if(num==0) return res;



    int x=prev;

    while(true){

        if(product*x>a) break;

        res+=dfs(denom*x,numer*x+denom,x,product*x,num-1);

        ++x;

    }

    return res;

}



void solve(){

    while(cin >> p >> q >> a >> n){

        if(p==0 and q==0 and a==0 and n==0) break;

        cout << dfs(1,0,1,1,n) << endl;

    }

}



int main(){

	std::ios::sync_with_stdio(false);

	std::cin.tie(0);

	solve();

	return 0;

}

// Faster C++ solution for the same problem (measured speedup: 1.2087300995920114x):
#include <bits/stdc++.h>

using namespace std;



#define rep(i,x,y) for(int i=(x);i<(y);++i)

#define debug(x) #x << "=" << (x)



#ifdef DEBUG

#define _GLIBCXX_DEBUG

#define dump(x) std::cerr << debug(x) << " (L:" << __LINE__ << ")" << std::endl

#else

#define dump(x)

#endif



typedef long long int ll;

typedef pair<int,int> pii;

//template<typename T> using vec=std::vector<T>;



const int inf=1<<30;

const long long int infll=1LL<<58;

const double eps=1e-9;

const int dx[]={1,0,-1,0},dy[]={0,1,0,-1};



template <typename T> ostream &operator<<(ostream &os, const vector<T> &vec){

	os << "[";

	for (const auto &v : vec) {

		os << v << ",";

	}

	os << "]";

	return os;

}



int p,q,a,n;



int dfs(int denom,int numer,int prev,int product,int num){

    {

        int tmp1=denom*prev,tmp2=numer*prev+denom*num;

        if(tmp2*q<p*tmp1) return 0;

    }

    if(numer*q>p*denom) return 0;



    int res=0;

    if(numer*q==p*denom) ++res;

    if(num==0) return res;



    int x=prev;

    while(true){

        if(product*x>a) break;

        res+=dfs(denom*x,numer*x+denom,x,product*x,num-1);

        ++x;

    }

    return res;

}



void solve(){

    while(cin >> p >> q >> a >> n){

        if(p==0 and q==0 and a==0 and n==0) break;

        cout << dfs(1,0,1,1,n) << endl;

    }

}



int main(){

	std::ios::sync_with_stdio(false);

	std::cin.tie(0);

	solve();

	return 0;

}
// Slower C++ solution:
#include <iostream>

using namespace std;

int N, M, L, R, s[50003], t[50003], e[50003], c[10003]; long long d[10003];

int main() {

	cin >> N >> M;

	while(true) {

		cin >> s[R] >> t[R] >> e[R]; s[R]--, t[R]--;

		if (!e[R]) break; R++;

	}

	cin >> L;

	for (int i = 0; i < L; i++) {

		for (int j = 0; j < M; j++) cin >> c[j];

		for (int j = 0; j < N; j++) d[j] = 0;

		for (int j = 0; j < R; j++) d[s[j]] += c[t[j]] * e[j];

		for (int j = 0; j < N; j++) cout << (j ? " " : "") << d[j];

		cout << endl;

	}

}

// Faster C++ solution for the same problem (measured speedup: 5.07500886190697x):
#include <stdio.h>

int N, M, L, R, s[50003], t[50003], e[50003], c[10003]; long long d[10003];

int main() {

	scanf("%d%d", &N, &M);

	while(true) {

		scanf("%d%d%d", s + R, t + R, e + R); s[R]--, t[R]--;

		if (!e[R]) break; R++;

	}

	scanf("%d", &L);

	for (int i = 0; i < L; i++) {

		for (int j = 0; j < M; j++) scanf("%d", &c[j]);

		for (int j = 0; j < N; j++) d[j] = 0;

		for (int j = 0; j < R; j++) d[s[j]] += c[t[j]] * e[j];

		for (int j = 0; j < N; j++) printf(j ? " %d" : "%d", d[j]);

		printf("\n");

	}

}