#include <bits/stdc++.h>
                                        #include<iostream>
                                        #include<cstdio>
                                        #include<vector>
                                        #include<queue>
                                        #include<map>
                                        #include<cstring>
                                        #include<string>
                                        #include <math.h>
                                        #include<algorithm>
                                    //    #include <boost/multiprecision/cpp_int.hpp>
                                        #include<functional>
                                #define int long long
                                        #define inf  1000000007
                                        #define pa pair<int,int>
                                        #define ll long long
                                        #define pal pair<double,pa>
                                        #define ppa pair<string,int>
                                        #define ppap pair<int,pa>
                                        #define ssa pair<string,int>
                                        #define  mp make_pair
                                        #define  pb push_back
                                        #define EPS (1e-12)
                                        #define equals(a,b) (fabs((a)-(b))<EPS)
                                  #define VI vector<int>
                                        using namespace std;
                                        
                                  class pas{
                                    	public:
                                    	int x,y;
                                    	pas(int x=0,int y=0):x(x),y(y) {}
                                    	pas operator + (pas p) {return pas(x+p.x,y+p.y);}
                                    	pas operator - (pas p) {return pas(x-p.x,y-p.y);}
                                    	pas operator * (int a) {return pas(x*a,y*a);}
                                    	pas operator / (int a) {return pas(x/a,y/a);}
                                    //	double absv() {return sqrt(norm());}
                                    	int norm() {return x*x+y*y;}
                                    	bool operator < (const pas &p) const{
                                   		return x != p.x ? x<p.x: y<p.y;
                                   	}
                                  //		bool operator < (const pas &p) const{
                                    //		return y != p.y ? y<p.y: x<p.x;
                                    //	}
                                    	bool operator == (const pas &p) const{
                                    		return x==p.x && y==p.y;
                                    	}
                                    };

                                        class Point{
                                        	public:
                                        	double x,y;
                                        	Point(double x=0,double y=0):x(x),y(y) {}
                                        	Point operator + (Point p) {return Point(x+p.x,y+p.y);}
                                        	Point operator - (Point p) {return Point(x-p.x,y-p.y);}
                                        	Point operator * (double a) {return Point(x*a,y*a);}
                                        	Point operator / (double a) {return Point(x/a,y/a);}
                                        	double absv() {return sqrt(norm());}
                                        	double norm() {return x*x+y*y;}
                                        	bool operator < (const Point &p) const{
                                        		return x != p.x ? x<p.x: y<p.y;
                                        	}
                                        	bool operator == (const Point &p) const{
                                        		return fabs(x-p.x)<EPS && fabs(y-p.y)<EPS;
                                        	}
                                        };
                                        typedef Point Vector;
                                 #define pl pair<int,pas>
                                        struct Segment{
                                        Point p1,p2;
                                        };
                                        
struct star{
Segment se[5];
};
                                         double dot(Vector a,Vector b){
                                        	return a.x*b.x+a.y*b.y;
                                        }
                                        double cross(Vector a,Vector b){
                                        	return a.x*b.y-a.y*b.x;
                                        }
                                    
            bool parareru(Point a,Point b,Point c,Point d){
            //	if(abs(cross(a-b,d-c))<EPS)cout<<"dd "<<cross(a-b,d-c)<<endl;
            	return abs(cross(a-b,d-c))<EPS;
            }
            double distance_ls_p(Point a, Point b, Point c) {
              if ( dot(b-a, c-a) < EPS ) return (c-a).absv();
              if ( dot(a-b, c-b) < EPS ) return (c-b).absv();
              return abs(cross(b-a, c-a)) / (b-a).absv();
            }
            bool is_intersected_ls(Segment a,Segment b) {
            	if(a.p1==b.p1||a.p2==b.p1||a.p1==b.p2||a.p2==b.p2) return false;
            	if(parareru((a.p2),(a.p1),(b.p1),(b.p2))&&parareru((a.p2),(a.p1),(b.p1),(b.p1))){
            //		cout<<"sss"<<endl;
            		if(dot(a.p1-b.p1,a.p1-b.p2)<EPS) return true;
            		if(dot(a.p2-b.p1,a.p2-b.p2)<EPS) return true;
            		if(dot(a.p1-b.p1,a.p2-b.p1)<EPS) return true;
            		if(dot(a.p1-b.p2,a.p2-b.p2)<EPS) return true;
            		return false;
            	}
              else return ( cross(a.p2-a.p1, b.p1-a.p1) * cross(a.p2-a.p1, b.p2-a.p1) < EPS ) && ( cross(b.p2-b.p1, a.p1-b.p1) * cross(b.p2-b.p1, a.p2-b.p1) < EPS );
            }
            double segment_len(Segment a){
            	
            	return (a.p1-a.p2).absv();
            }
            double segment_dis(Segment a,Segment b){
            	if(is_intersected_ls(a,b))return 0;
            	double r=distance_ls_p(a.p1, a.p2, b.p1);
            	r=min(r,distance_ls_p(a.p1, a.p2, b.p2));
            	r=min(r,distance_ls_p(b.p1, b.p2, a.p2));
            	r=min(r,distance_ls_p(b.p1, b.p2, a.p1));
            	return r;
            }
            Point intersection_ls(Segment a, Segment b) {
              Point ba = b.p2-b.p1;
              double d1 = abs(cross(ba, a.p1-b.p1));
              double d2 = abs(cross(ba, a.p2-b.p1));
              double t = d1 / (d1 + d2);
             
              return a.p1 + (a.p2-a.p1) * t;
            }
			pair<Point,Point> circle_intersection(Point c1,double r1,Point c2,double r2){
				double d=(c1-c2).absv();
				double h=(r1*r1-r2*r2+d*d)/2.0/d;
				double l=sqrt(r1*r1-h*h);
			//	cout<<d<<" "<<h<<" "<<l<<endl;
				Point asi=c1+(c2-c1)*h/((c2-c1).absv());
				Vector r1r2=(c2-c1)/((c2-c1).absv());
				Vector sui={r1r2.y,-r1r2.x};
			//	cout<<sui.x<<" "<<sui.y<<endl;
				pair<Point,Point> z=mp(asi+sui*l,asi-sui*l);
				if(z.first.x>z.second.x) swap(z.first,z.second);
				return z;
			}

int gcd(int x,int y){
	if(x<y) return gcd(y,x);
	if(x==y) return x;
	if(x%y==0) return y;
	return gcd(y,x%y);
}
 class pa2{
                                            public:
                                            int x,y;
                                            pa2(int x=0,int y=0):x(x),y(y) {}
                                            pa2 operator + (pa2 p) {return pa2(x+p.x,y+p.y);}
                                            pa2 operator - (pa2 p) {return pa2(x-p.x,y-p.y);}
                                            bool operator < (const pa2 &p) const{
                                                return x != p.x ? x<p.x: y<p.y;
                                            }
										 	 bool operator > (const pa2 &p) const{
                                                return x != p.x ? x>p.x: y>p.y;
                                            }
                                            bool operator == (const pa2 &p) const{
                                                return abs(x-p.x)==0 && abs(y-p.y)==0;
                                            }
                                            bool operator != (const pa2 &p) const{
                                                return !(abs(x-p.x)==0 && abs(y-p.y)==0);
                                            }
                                                 
                                         
                                        };

/*　　　---memo---
add_edge(始点,終点,youryou ,コスト)で枝追加
min_cost_flow(始点,終点,ryou)でhiyouをintで返す
*/
int p,q,ans,a,n;

pa bun(int si1,int bo1,int si2,int bo2){
	int g=gcd(bo1,bo2);
	int bo3=bo1*bo2/g;
	int si3=si1*bo2/g+si2*bo1/g;
	g=gcd(si3,bo3);
	return mp(si3/g,bo3/g);
}


void dfs(int si,int bo,int seki,int kosuu,int saigo){
	
	if(p*bo==si*q){
	ans++;
		return;
	}
	if(p*bo<si*q) return;
	if(kosuu==n) return;
	//cout<<si<<" "<<bo<<" "<<seki<<" "<<kosuu<<" "<<saigo<<endl;
	for(int i=saigo;i*seki<=a;i++){
	//	cout<<i<<endl;
		pa z=bun(si,bo,n-kosuu,i);
		if(z.second*p>z.first*q) return;
		
		int g=gcd(bo,i);
		int s2=bo/g;
		int bo2=bo*(i/g);
		int si2=si*(i/g)+s2;
		g=gcd(si2,bo2);
		if(g>1){
		si2/=g;
			bo2/=g;
		}
		
		dfs(si2,bo2,i*seki,kosuu+1,i);
	}
}

     signed main(){
     	while(1){
     	cin>>p>>q>>a>>n;
     		if(p==0 && q==0) return 0;
     		ans=0;
     		dfs(0,1,1,0,1);
     		cout<<ans<<endl;
     	
     	}
    
           	return 0;
           }


Measured memory_bytes: 3032