Problem:
 Problem C: Unit Fraction Partition A fraction whose numerator is 1
and whose denominator is a positive integer
is called a unit fraction.
A representation of a positive rational number p / q as the sum of finitely many unit fractions is called a partition of p / q into unit fractions.
For example,
1/2 + 1/6
is a partition of
2/3
into unit fractions.
The difference in the order of addition is disregarded.
For example, we do not distinguish
1/6 + 1/2
from
1/2 + 1/6. For given four positive integers p , q , a , and n ,
count the number
of partitions of p / q into unit fractions satisfying the following two conditions.  The partition is the sum of at most n many unit fractions.  The product of the denominators of the unit fractions in the
partition is less than or equal to a .  For example, if
( p , q , a , n ) = (2,3,120,3),
you should report 4 since   enumerates all of the valid partitions. Input The input is a sequence of at most 1000 data sets followed by a terminator. A data set is a line containing four positive integers p , q , a , and n satisfying p , q <= 800, a <= 12000
and n <= 7.
The integers are separated by a space. The terminator is composed of just one line
which contains four zeros separated by a space.
It is not a part of the input data but a mark for the end of the input. Output The output should be composed of lines
each of which contains a single integer.
No other characters should appear in the output. The output integer corresponding to a data set p , q , a , n should be the number of all partitions of p / q into at most n many unit fractions such that the product of the denominators of the
unit fractions is less than or equal to a . Sample Input 2 3 120 3
2 3 300 3
2 3 299 3
2 3 12 3
2 3 12000 7
54 795 12000 7
2 3 300 1
2 1 200 5
2 4 54 2
0 0 0 0 Output for the Sample Input 4
7
6
2
42
1
0
9
3

Solution (Python):
import math,string,itertools,fractions,heapq,collections,re,array,bisect,sys,random,time,copy,functools

sys.setrecursionlimit(10**7)
inf = 10**20
eps = 1.0 / 10**10
mod = 10**9+7
dd = [(-1,0),(0,1),(1,0),(0,-1)]
ddn = [(-1,0),(-1,1),(0,1),(1,1),(1,0),(1,-1),(0,-1),(-1,-1)]

def LI(): return [int(x) for x in sys.stdin.readline().split()]
def LI_(): return [int(x)-1 for x in sys.stdin.readline().split()]
def LF(): return [float(x) for x in sys.stdin.readline().split()]
def LS(): return sys.stdin.readline().split()
def I(): return int(sys.stdin.readline())
def F(): return float(sys.stdin.readline())
def S(): return input()
def pf(s): return print(s, flush=True)

def main():
    rr = []

    def f(p, q, a, n):
        r = 0
        fr = fractions.Fraction(p,q)
        q = [(fr, 0, 1)]
        if fr.numerator == 1 and fr.denominator <= a:
            r += 1
        n1 = n-1
        for i in range(1, 70):
            nq = []
            ti = fractions.Fraction(1, i)
            i2 = (i+1) ** 2
            for t, c, m in q:
                mik = m
                tt = t
                for k in range(1, n-c):
                    mik *= i
                    tt -= ti
                    if tt <= 0 or a < mik * i:
                        break
                    if tt.numerator == 1 and tt.denominator >= i and mik * tt.denominator <= a:
                        r += 1
                    if c+k < n1 and mik * (i ** max(math.ceil(tt / ti), 2)) < a:
                        nq.append((tt, c+k, mik))
                if m * (i ** max(math.ceil(t / ti), 2)) < a:
                    nq.append((t,c,m))
            if not nq:
                break
            q = nq
        return r

    while True:
        p, q, a, n = LI()
        if p == 0 and q == 0:
            break
        rr.append(f(p, q, a, n))

    return '\n'.join(map(str, rr))



print(main())

Solution (Python):
from functools import lru_cache

@lru_cache(maxsize=1<<10)
def solve(p,q,a,n):
  def _solve(num,dem,d,m,s):
    if num==0:
      return 1
    if d == 0:
      return 0
    if num * a //m < dem:
      return 0
    return sum((_solve(num*i-dem, dem*i, d-1, m*i, i) for i in range(s, min(dem*n//num, a//m)+1)),0)
  return _solve(p,q,n,1,1)
ans=[]
while True:
  p,q,a,n=map(int,input().split())
  if p==0:
    break
  ans.append(solve(p,q,a,n))
print(*ans, sep="\n")

Solution (Python):
from fractions import gcd

def solve(p, q, a, n, l=1):
    ans = 1 if p==1 and q<=a and q>=l else 0
    denom = max(l, q//p)
    p_denom = denom*p
    while n*q >= p_denom and denom <= a:  
        p_, q_ = p_denom-q, q*denom
        if p_ <= 0:
            denom += 1
            p_denom += p
            continue
        gcd_ = gcd(p_, q_)
        p_ //= gcd_
        q_ //= gcd_
        if n==2 and p_==1 and q_*denom<=a and q_>=denom:
            ans += 1
        else:
            ans += solve(p_, q_, a//denom, n-1, denom)
        denom += 1
        p_denom += p
    
    return ans

while True:
    p, q, a, n = map(int, input().split())
    if p==q==a==n==0:
        break
    gcd_ = gcd(p, q)
    p //= gcd_
    q //= gcd_
    if n==1:
        print(1 if p==1 and q<=a else 0)
    else:
        print(solve(p, q, a, n))

Solution (Python):
from functools import lru_cache

@lru_cache(maxsize=1<<10)
def solve(p, q, a, n):
    def _solve(num, dem, d, m, s):
        if num == 0:
            return 1
        if d == 0:
            return 0
        if num * a // m < dem:
            return 0
        return sum((_solve(num*i-dem, dem*i, d-1, m*i, i)
                    for i in range(s, min(dem*n//num, a//m)+1)), 0)
    return _solve(p, q, n, 1, 1)

if __name__ == "__main__":
    ans = []
    while True:
        p, q, a, n = [int(i) for i in input().split()]
        if p == 0:
            break
        ans.append(solve(p, q, a, n))
    print(*ans, sep="\n")

Solution (Python):
from functools import lru_cache

@lru_cache(maxsize=1<<10)
def solve(p, q, a, n):
    def _solve(num, dem, d, m, s):
        if num == 0:
            return 1
        if d == 0:
            return 0
        if num * a // m < dem:
            return 0
        return sum((_solve(num*i-dem, dem*i, d-1, m*i, i)
                    for i in range(s, min(dem*n//num, a//m)+1)), 0)
    return _solve(p, q, n, 1, 1)

if __name__ == "__main__":
    ans = []
    while True:
        p, q, a, n = [int(i) for i in input().split()]
        if p == 0:
            break
        ans.append(solve(p, q, a, n))
    print(*ans, sep="\n")

Solution (Python):


def dfs(num, dem, k, cur_num, prd, p, q, a, n):

    
    
    if num == 0:
        
        return 1

    result = 0

    i = cur_num
    while prd * i <= a:
        
        if (n - k) * dem < num * i:
            break
        if num * i < dem:
            i += 1
            continue

        result += dfs((num * i) - dem, dem * i, k + 1, i, prd * i, p, q, a, n)
        i += 1

    return result

def func():
    while True:

        p, q, a, n = map(int, raw_input().split())

        if p == 0 and q == 0 and a == 0 and n == 0:
            break

        print dfs(p, q, 0, 1, 1, p, q, a, n)
        
    return None

if __name__ == '__main__':
    func()