CodeEff/ECCO huggingface dataset:
# Original solution
import sys
input = sys.stdin.readline
n,Q=list(map(int,input().split()))
b_num = 2**(n-1).bit_length()
mx = 10**9
segl=[mx]*2*b_num 
def update(k, x):
    k += b_num-1
    segl[k] = x
    while k+1:
        k = (k-1)//2
        segl[k] = min(segl[k*2+1],segl[k*2+2])
if __name__ == '__main__':
    xs = []
    for i in range(n):
        s,t,x=list(map(int,input().split()))
        xs += [(s-x,True,x,i), (t-x,False,x,i)]
    qs = []
    for _ in range(Q):
        qs.append(int(eval(input())))
    xs.sort(key=lambda x:x[0])
    j = 0
    for x in xs:
        while Q-j:
            if qs[j] < x[0]:
                print((segl[0] if segl[0]!= mx else -1))
                j += 1
                continue
            break
        update(x[3], x[2] if x[1] else mx)
    while Q-j:
        print((-1))
        j+=1


# Optimized solution
import sys
input = sys.stdin.readline
n,Q=list(map(int,input().split()))
b_num = 2**(n-1).bit_length()
mx = 10**9
segl=[mx]*2*b_num 
def update(k, x):
    k += b_num-1
    segl[k] = x
    while k+1:
        k = (k-1)//2
        segl[k] = min(segl[k*2+1],segl[k*2+2])
if __name__ == '__main__':
    xs = []
    for i in range(n):
        s,t,x=list(map(int,input().split()))
        xs += [(s-x,True,x,i), (t-x,False,x,i)]
    qs = []
    for _ in range(Q):
        qs.append(int(eval(input())))
    xs.sort(key=lambda x:x[0])
    j = 0
    for x in xs:
        while Q-j and qs[j] < x[0]:
            print((segl[0] if segl[0]!= mx else -1))
            j += 1
        update(x[3], x[2] if x[1] else mx)
    while Q-j:
        print((-1))
        j+=1
# Original solution
import sys;input=sys.stdin.readline
N, Q = list(map(int, input().split()))
X = []
for _ in range(N):
    s, t, x = list(map(int, input().split()))
    X.append((2*s-2*x-1, 0, x))
    X.append((2*t-2*x-1, 1, x))

inf = 10**18
X.append((-inf, 0, inf))
X.append((inf, 1, inf))

X.sort(key=lambda x: x[0])

from heapq import *
from collections import defaultdict
#print(X[1:-1])
j = 0
ll = []
vs = set()
d = defaultdict(int)
for _ in range(Q):
    q = int(eval(input()))
    while 1:
        u, f, x = X[j]
        if u > 2*q:
            break
        if f:
            d[x] += 1
        else:
            heappush(ll, x)
        j += 1
    while 1:
        x = heappop(ll)
        if d[x] > 0:
            d[x] -= 1
            continue
        heappush(ll, x)
        if x == inf:
            print((-1))
        else:
            print(x)
        break



# Optimized solution

from heapq import *
def main():
    import sys;input=sys.stdin.readline
    N, Q = list(map(int, input().split()))
    X = []
    for _ in range(N):
        s, t, x = list(map(int, input().split()))
        X.append((2*s-2*x-1, 0, x))
        X.append((2*t-2*x-1, 1, x))
    inf = 10**18
    X.append((-inf, 0, inf))
    X.append((inf, 1, inf))
    X.sort(key=lambda x: x[0])

    j = 0
    ll = []
    vs = set()
    d = dict()
    for _ in range(Q):
        q = int(eval(input()))
        while 1:
            u, f, x = X[j]
            if u > 2*q:
                break
            if f:
                if x not in d:
                    d[x] = 0
                d[x] += 1
            else:
                heappush(ll, x)
            j += 1
        while 1:
            x = heappop(ll)
            if x in d and d[x] > 0:
                d[x] -= 1
                continue
            heappush(ll, x)
            if x == inf:
                print((-1))
            else:
                print(x)
            break

if __name__ == '__main__':
    main()
# Original solution
from bisect import bisect_left
inf = float('inf')

class SegmentTree:
    def __init__(self, N):
        tmp = 1
        while tmp < N:
            tmp <<= 1
        self.N = tmp
        self.data = [inf for _ in range(2*self.N-1)]
    
    def update(self, l, r, x):
        L = l + self.N
        R = r + self.N

        while L < R:
            if R & 1:
                R -= 1
                self.data[R] = min(self.data[R], x)
            if L & 1:
                self.data[L] = min(self.data[L], x)
                L += 1
            L >>= 1; R >>= 1
    
    def value(self, k):
        k += self.N
        res = self.data[k]
        while k != 0:
            k >>= 1
            res = min(res, self.data[k])
        return res

N, Q = list(map(int,input().split()))
STX = [list(map(int,input().split())) for _ in range(N)]
D = [int(eval(input())) for _ in range(Q)]

st = SegmentTree(Q)
for S, T, X in STX:
    l = bisect_left(D, S-X)
    r = bisect_left(D, T-X)
    st.update(l, r, X)

for k in range(Q):
    ans = st.value(k)
    if ans == inf:
        ans = -1
    print(ans)

# Optimized solution
import sys
input = sys.stdin.readline
from bisect import bisect_left
inf = float('inf')

N, Q = list(map(int,input().split()))
STX = [list(map(int,input().split())) for _ in range(N)]
D = [int(eval(input())) for _ in range(Q)]

res = 1
while res < Q:
    res <<= 1
N0 = res
data = [inf for _ in range(2*N0-1)]

def update(l, r, x):
    L = l + N0
    R = r + N0
    while L < R:
        if R & 1:
            R -= 1
            data[R] = min(data[R], x)
        if L & 1:
            data[L] = min(data[L], x)
            L += 1
        L >>= 1; R >>= 1

def value(k):
    k += N0
    res = data[k]
    while k != 0:
        k >>= 1
        res = min(res, data[k])
    return res


for S, T, X in STX:
    l = bisect_left(D, S-X)
    r = bisect_left(D, T-X)
    update(l, r, X)

for k in range(Q):
    ans = value(k)
    if ans == inf:
        ans = -1
    print(ans)
# Original solution
from heapq import heappush, heappop

n, q = map(int, input().split())
STX = [list(map(int, input().split())) for _ in range(n)]
task = []
for s, t, x in STX:
    task.append((t - x, 0, x))
    task.append((s - x, 1, x))
for i in range(q):
    d = int(input())
    task.append((d, 2, i))

task.sort()
kouho = set()
kouho_hp = []
ans = [-1] * q
for a, b, c in task:
    if b == 0:
        kouho.remove(c)
    if b == 1:
        kouho.add(c)
        heappush(kouho_hp, c)
    if b == 2:
        while kouho_hp and kouho_hp[0] not in kouho:
            heappop(kouho_hp)
        if not kouho_hp:
            continue
        else:
            ans[c] = kouho_hp[0]
        # heapは0-indexが最小値になっているからheapified_list[0]でもよい
print(*ans, sep="\n")


# Optimized solution
from heapq import heappush, heappop
import sys

input = sys.stdin.buffer.readline


n, q = map(int, input().split())
STX = [list(map(int, input().split())) for _ in range(n)]
task = []
for s, t, x in STX:
    task.append((t - x, 0, x))
    task.append((s - x, 1, x))
for i in range(q):
    d = int(input())
    task.append((d, 2, i))

task.sort()
kouho = set()
kouho_hp = []
ans = [-1] * q
for a, b, c in task:
    if b == 0:
        kouho.remove(c)
    if b == 1:
        kouho.add(c)
        heappush(kouho_hp, c)
    if b == 2:
        while kouho_hp and kouho_hp[0] not in kouho:
            heappop(kouho_hp)
        if not kouho_hp:
            continue
        else:
            ans[c] = kouho_hp[0]
        # heapは0-indexが最小値になっているからheapified_list[0]でもよい
print(*ans, sep="\n")
# Original solution
#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""
input:
4 5
0 0 1 0 0
1 0 0 0 0
0 0 0 1 0
0 0 0 1 0

output:
4
"""
import sys


def solve():
    dp = [[0] * W for _ in range(H)]
    max_width = 0
    for m in range(H):
        for n in range(W):
            dp[m][n] = (int(carpet_info[m][n]) + 1) % 2
            max_width |= dp[m][n]

    for i in range(1, H):
        for j in range(1, W):
            if int(carpet_info[i][j]):
                dp[i][j] = 0
            else:
                dp[i][j] = min(dp[i - 1][j - 1], dp[i - 1][j], dp[i][j - 1]) + 1
                max_width = max(max_width, dp[i][j])

    return pow(max_width, 2)


if __name__ == '__main__':
    _input = sys.stdin.readlines()
    H, W = list(map(int, _input[0].split()))
    carpet_info = list([x.split() for x in _input[1:]])
    print((solve()))

# Optimized solution
#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""
input:
4 5
0 0 1 0 0
1 0 0 0 0
0 0 0 1 0
0 0 0 1 0

output:
4
"""
import sys


def solve():
    dp = [[0] * (W + 1) for _ in range(H + 1)]

    for i in range(H):
        for j in range(W):
            if not int(carpet_info[i][j]):
                dp[i + 1][j + 1] = min(dp[i][j], dp[i][j + 1], dp[i + 1][j]) + 1

    cand = (max(each) for each in dp)
    return pow(max(cand), 2)


if __name__ == '__main__':
    _input = sys.stdin.readlines()
    H, W = list(map(int, _input[0].split()))
    carpet_info = list([x.split() for x in _input[1:]])
    print((solve()))
# Original solution
import sys
h, w = list(map(int, sys.stdin.readline().split()))
dp = [[0 for j in range(w)] for i in range(h)]
G = [[int(j) for j in sys.stdin.readline().split()] for i in range(h)]
maxWidth = 0
for i in range(0, h):
    for j in range(0, w):
        dp[i][j] = (G[i][j] + 1) % 2
        maxWidth |= dp[i][j]
for i in range(1, h):
    for j in range(1, w):
        if(G[i][j] == 0):
            dp[i][j] = min(dp[i - 1][j - 1], dp[i - 1][j], dp[i][j - 1]) + 1
            if maxWidth < dp[i][j]:
                maxWidth = dp[i][j]
print((maxWidth * maxWidth))

# Optimized solution
import sys
h, w = list(map(int, sys.stdin.readline().split()))
dp = [[0] * w for _ in range(h)]
G = [[int(j) for j in sys.stdin.readline().split()] for _ in range(h)]
for i in range(0, h):
    for j in range(0, w):
        dp[i][j] = (G[i][j] + 1) % 2
for i in range(1, h):
    for j in range(1, w):
        if(G[i][j] == 0):
            dp[i][j] = min(dp[i - 1][j - 1], dp[i - 1][j], dp[i][j - 1]) + 1
print((max([max(item) for item in dp])**2))
# Original solution
import sys
h, w = list(map(int, sys.stdin.readline().split()))
dp = [[0] * w for _ in range(h)]
G = [[int(j) for j in sys.stdin.readline().split()] for _ in range(h)]
for i in range(0, h):
    for j in range(0, w):
        dp[i][j] = (G[i][j] + 1) % 2
for i in range(1, h):
    for j in range(1, w):
        if(G[i][j] == 0):
            dp[i][j] = min(dp[i - 1][j - 1], dp[i - 1][j], dp[i][j - 1]) + 1
print((max([max(item) for item in dp])**2))

# Optimized solution
import sys
h, w = list(map(int, sys.stdin.readline().split()))
dp = [[0] * w for _ in range(h)]
G = [[int(j) for j in sys.stdin.readline().split()] for _ in range(h)]
for x in range(h):
    dp[x][0] = 1 if G[x][0] == 0 else 0
for y in range(w):
    dp[0][y] = 1 if G[0][y] == 0 else 0
for i in range(1, h):
    for j in range(1, w):
        if(G[i][j] == 0):
            dp[i][j] = min(dp[i - 1][j - 1], dp[i - 1][j], dp[i][j - 1]) + 1
print((max([max(item) for item in dp])**2))
# Original solution
import sys
import itertools
h, w = list(map(int, sys.stdin.readline().split()))
dp = [[0] * w for _ in range(h)]
G = [[int(j) for j in sys.stdin.readline().split()] for _ in range(h)]
for x in range(h):
    dp[x][0] = 1 if G[x][0] == 0 else 0
for y in range(w):
    dp[0][y] = 1 if G[0][y] == 0 else 0
for i in range(1, h):
    for j in range(1, w):
        if(G[i][j] == 0):
            dp[i][j] = min(dp[i - 1][j - 1], dp[i - 1][j], dp[i][j - 1]) + 1
print((max(itertools.chain(*dp)) ** 2))

# Optimized solution
import sys
import itertools
h, w = list(map(int, sys.stdin.readline().split()))
dp = [[0] * w for _ in range(h)]
G = [[int(j) for j in sys.stdin.readline().split()] for _ in range(h)]
for x in range(h):
    dp[x][0] = 1 if G[x][0] == 0 else 0
for y in range(w):
    dp[0][y] = 1 if G[0][y] == 0 else 0
for i in range(1, h):
    for j in range(1, w):
        if(G[i][j] == 0):
            dp[i][j] = min(dp[i - 1][j - 1], dp[i - 1][j], dp[i][j - 1]) + 1
print((max(max(i) for i in dp) ** 2))
# Original solution
H, W = list(map(int, input().split()))
c = []
for i in range(H):
    c.append([int(i) for i in input().split()])
dp = [[0 for i in range(W)] for i in range(H)]
maxwidth = 0
for i in range(H):
    for j in range(W):
        dp[i][j] = (c[i][j] + 1) % 2
        maxwidth |= dp[i][j]

for i in range(1, H):
    for j in range(1, W):
        if c[i][j]:
            dp[i][j] = 0
        else:
            dp[i][j] = min(dp[i-1][j], dp[i-1][j-1], dp[i][j-1]) + 1
            maxwidth = max(maxwidth,dp[i][j])

print((maxwidth**2))


# Optimized solution
H, W = list(map(int, input().split()))
c = []

for i in range(H):
    c.append([int(i) for i in input().split()])
dp = [[0 for i in range(W)] for i in range(H)]

maxwidth = 0

for i in range(H):
    for j in range(W):
        dp[i][j] = (c[i][j] + 1) % 2
        maxwidth |= dp[i][j]

for i in range(1, H):
    for j in range(1, W):
        if c[i][j]:
            continue
        else:
            dp[i][j] = min(dp[i-1][j], dp[i-1][j-1], dp[i][j-1]) + 1
            maxwidth = max(maxwidth, dp[i][j])

print((maxwidth**2))
# Original solution
def counting(k, v, square):
    square += 1
    for _ in range(square):

        for r in range(k, k + square):
            if r > len(matrix) - 1:
                return False
            for c in range(v, v + square):
                if c > len(matrix[k]) - 1:
                    return False
                if matrix[r][c] == '1':
                    return False
        return True


line = input().split(' ')

answer = [[0 for _ in range(int(line[1]))] for _ in range(int(line[0]))]
matrix = []
maximum = 0
is_over = False
ggg = 0
for _ in range(int(line[0])):
    row = input().split(' ')
    matrix.append(row)

for i in range(len(matrix)):
    for j in range(len(matrix[i])):
        ggg += 1
        for m in range(min(int(line[0]), int(line[1]))):
            res = counting(i, j, maximum)
            if res:
                maximum += 1
            else:
                break

print((maximum * maximum))

# Optimized solution
def solve2():
    r, c = list(map(int, input().split(' ')))

    ans_matrix = [[0 for _ in range(c + 1)] for _ in range(r + 1)]
    input_matrix = []

    for _ in range(r):
        input_matrix.append(list(map(int, input().split(' '))))

    largest_square = 0

    for i in range(len(input_matrix)):
        for j in range(len(input_matrix[i])):
            new_i, new_j = i + 1, j + 1
            if input_matrix[i][j] == 0:
                ans_matrix[new_i][new_j] = min(ans_matrix[new_i - 1][new_j], ans_matrix[new_i][new_j - 1], ans_matrix[new_i - 1][new_j - 1]) + 1
                largest_square = ans_matrix[new_i][new_j] if ans_matrix[new_i][new_j] > largest_square else largest_square
            else:
                ans_matrix[new_i][new_j] = 0

    print((largest_square ** 2))


def main():
    solve2()


if __name__ == '__main__':
    main()
# Original solution
# -*- coding: utf-8 -*-

result = 0

def dfs(num, dem, k, cur_num, p, q, a, n):
    global result
    #print num, dem, k, cur_num
    #time.sleep(0.5)
    if (num * q) == (dem * p):
        result += 1
        #print "result: " + str(result)
        return

    if k == n:
        return

    if (num * q) > (dem * p):
        return

    i = cur_num
    while i * dem <= a:
        dfs((num * i) + dem, dem * i, k + 1, i, p, q, a, n)
        i += 1

    return

def func():
    while True:
        global result
        p, q, a, n = list(map(int, input().split()))

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

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

if __name__ == '__main__':
    func()
    

# Optimized solution
# -*- coding: utf-8 -*-

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

    #print num, dem, k, cur_num, prd
    #time.sleep(0.5)
    if num == 0:
        #print "result: " + str(result)
        return 1

    result = 0

    i = cur_num
    while prd * i <= a:
        #print i, num, dem, k, cur_num, prd
        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 = list(map(int, 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()
# Original solution
from functools import lru_cache
from fractions import Fraction

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

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

# Optimized solution
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")
# Original solution
r ,d, x = list(map(int,input().split()))
xi = 0
for i in range(10):
  if i ==0:
    xi = r * x - d
  else:
    xi = r * xi - d
  print(xi)
  

# Optimized solution
r ,d, x = list(map(int,input().split()))
for i in range(10):
  x = r * x - d
  print(x)
# Original solution
r, d, x = list(map(int, input().split()))

for i in range(10):
    x= r*x -d
    print(x)

# Optimized solution
r, d, x = list(map(int, input().split()))

for _ in range(10):
    x= r*x -d
    print(x)
# Original solution
#!/usr/bin/env python3
# -*- coding: utf-8 -*-

r, D, x = list(map(int, input().split()))

for _ in range(10):
    x = r * x - D
    print(x)


# Optimized solution
r, d, x = list(map(int, input().strip().split()))

for i in range(10):
    x = r * x - d
    print(x)
# Original solution
#!/usr/bin/env python3



# Generated by 1.1.4 https://github.com/kyuridenamida/atcoder-tools  (tips: You use the default template now. You can remove this line by using your custom template)
def main():
    # Failed to predict input format
    r,d,x = list(map(int, input().split()))

    for i in range(10):
        x = r*x-d
        print(x)
    

if __name__ == '__main__':
    main()


# Optimized solution
#!/usr/bin/env python3




def main():
    r, D, x2000 = list(map(int, input().split()))

    x = x2000
    for i in range(1,11):
        x = r*x-D
        print(x)
        
    

if __name__ == '__main__':
    main()
# Original solution
r, D, x_2000 = list(map(int, input().split()))
x = x_2000
for i in range(10):
    x = r * x - D
    print(x)

# Optimized solution
r, D, x = list(map(int, input().split()))
for i in range(10):
    x = r * x - D
    print(x)
# Original solution
r,d,x = list(map(int,input().split()))
for i in range(10):
  temp = r*x - d
  print(temp)
  x = temp

# Optimized solution
r,d,x=list(map(int,input().split()))
for i in range(10):x=r*x-d;print(x)
# Original solution
R_up, Down, First_number = list(map(int, input().split()))
Tmp_number = First_number
for _ in range(10):
    Tmp_number = R_up*Tmp_number-Down
    print(Tmp_number)

# Optimized solution
#B
r, D, x = list(map(int, input().split()))
for i in range(10):
    xnew = r*x - D
    print(xnew)
    x = xnew
# Original solution
r, d, x = list(map(int, input().split()))

tmp = x
for i in range(10):
    ans = r * tmp - d
    print(ans)
    tmp = ans

# Optimized solution
r, D, x = list(map(int, input().split()))

for i in range(10):
    ans = r * x - D
    print(ans)
    x = ans
# Original solution
r, D, x_2000 = list(map(int, input().split()))

def x(i):
  if i == 2000:
    return x_2000
  elif i >= 2001:
    return r * x(i - 1) - D
  
for i in range(2001, 2011):
  print((x(i)))

# Optimized solution
r, D, x2000 = list(map(int, input().split()))
x = [x2000]
for i in range(10):
  x.append(x[i] * r - D)
  print((x[i+1]))
# Original solution
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 = []
    M = 7368792

    def f(m, n):
        g = [None] * M
        i = m
        while True:
            if g[i] is None:
                if n == 0:
                    return i
                else:
                    n -= 1
                for j in range(i,M,i):
                    g[j] = 1
            i += 1

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

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




print(main())



# Optimized solution
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)


class Prime():
    def __init__(self, m):
        self.M = m + 1
        self.A = a = [True] * self.M
        a[0] = a[1] = 0
        self.T = t = []
        sq = int(math.sqrt(self.M))
        q = 0
        for i in range(2, sq + 1):
            if not a[i]:
                a[i] = q
                continue
            t.append(i)
            a[i] = q
            q += 1
            for j in range(i*i,m,i):
                a[j] = False

        for i in range(sq + 1, self.M):
            if a[i]:
                a[i] = q
                q += 1
                t.append(i)
            else:
                a[i] = q

    def is_prime(self, n):
        return self.A[n]


def main():
    rr = []

    M = 7368792
    pr = Prime(M)

    def f(m, n):
        nn = n
        g = [None] * (m*m)
        for i in range(m, m*m):
            if not g[i]:
                if n == 0:
                    return i
                n -= 1
                for j in range(i, m*m, i):
                    g[j] = 1
        d = pr.A[m*m]
        return pr.T[d+n]

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

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




print(main())