CodeEff/ECCO huggingface dataset:
# Original solution
N, A, B, K = list(map(int, input().split()))

mod = 998244353
F = [1]
RF = []


def inv(x):
    res = 1
    k = mod - 2
    y = x
    while k:
        if k & 1:
            res = (res * y) % mod
        y = (y * y) % mod
        k //= 2
    return res


def init():
    for i in range(N):
        F.append((F[i] * (i + 1)) % mod)
    for i in F:
        RF.append(inv(i))


def comb(x, y):
    if y < 0 or x < y:
        return 0
    return (F[x] * ((inv(F[x - y]) * inv(F[y])) % mod)) % mod


ans = 0
init()
for a in range(N):
    b = (K - a * A) / B
    if b < 0:
        break
    if b.is_integer():
        ans += (comb(N, a) * (comb(N, int(b)))) % mod
print((ans % mod))


# Optimized solution
N, A, B, K = list(map(int, input().split()))

mod = 998244353
fact = [1] * (N + 1)
ifact = [1] * (N + 1)


def init():
    for i in range(N):
        fact[i + 1] = fact[i] * (i + 1) % mod
        ifact[i + 1] = pow(fact[i + 1], mod - 2, mod)


def comb(x, y):
    if y < 0 or y > x:
        return 0
    return (fact[x] * ifact[x - y] * ifact[y]) % mod


ans = 0
init()
for a in range(N + 1):
    n = K - a * A
    if n < 0:
        break
    if n % B == 0:
        ans += (comb(N, a) * (comb(N, n // B))) % mod
print((ans % mod))
# Original solution
def cmb(n, r, mod):
    if ( r<0 or r>n ):
        return 0
    r = min(r, n-r)
    return g1[n] * g2[r] * g2[n-r] % mod

mod = 998244353 #出力の制限
N = 10**6
g1 = [1, 1] # 元テーブル
g2 = [1, 1] #逆元テーブル
inverse = [0, 1] #逆元テーブル計算用テーブル

for i in range( 2, N + 1 ):
    g1.append( ( g1[-1] * i ) % mod )
    inverse.append( ( -inverse[mod % i] * (mod//i) ) % mod )
    g2.append( (g2[-1] * inverse[-1]) % mod )

N, A, B, K = list(map(int, input().split()))
ans = 0
for i in range(N+1):
  m = K-A*i
  if m%B!=0:
    continue
  ans += cmb(N,i,mod)*cmb(N,m//B,mod)
  ans %= mod

print(ans)

# Optimized solution
def cmb(n, r, mod):
    if ( r<0 or r>n ):
        return 0
    r = min(r, n-r)
    return g1[n] * g2[r] * g2[n-r] % mod

mod = 998244353 #出力の制限
N = 10**6
g1 = [1, 1] # 元テーブル
g2 = [1, 1] #逆元テーブル
inverse = [0, 1] #逆元テーブル計算用テーブル

for i in range( 2, N + 1 ):
    g1.append( ( g1[-1] * i ) % mod )
    inverse.append( ( -inverse[mod % i] * (mod//i) ) % mod )
    g2.append( (g2[-1] * inverse[-1]) % mod )

N,A,B,K = list(map(int, input().split()))
ans = 0
for i in range(N+1):
  if (K-A*i)%B!=0:
    continue
  j = (K-A*i)//B
  m = cmb(N,i,mod)*cmb(N,j,mod)%mod
  ans += m
  ans %= mod
print(ans)
# Original solution
def read():
  return int(eval(input()))

def reads():
  return [int(x) for x in input().split()]

P = 998244353

N, A, B, K = reads()

fact = [1]
for i in range(1, N+1):
  fact.append(fact[i-1] * i % P)

def egcd(a, b):
  if b == 0:
    return a, 1, 0
  else:
    q, r = divmod(a, b)
    g, x, y = egcd(b, r)
    return g, y, x - q * y

def inv(n):
  g, x, _ = egcd(n, P)
  # assert g == 1
  return x

invfact = [inv(fact[i]) for i in range(N+1)]

def comb(n, k):
  # assert 0 <= k <= n
  return fact[n] * invfact[k] * invfact[n-k] % P

result = 0
for x in range(N+1):
  y, r = divmod(K - A * x, B)
  if r == 0 and 0 <= y <= N:
    result = result + comb(N, x) * comb(N, y) % P

print((result % P))

# Optimized solution
def read():
  return int(eval(input()))

def reads():
  return [int(x) for x in input().split()]

P = 998244353

N, A, B, K = reads()

def inv(n):
  return pow(n, P-2, P)

fact = [1] * (N+1)
for i in range(1, N+1):
  fact[i] = fact[i-1] * i % P

invfact = [inv(fact[N])] * (N+1)
for i in range(N, 0, -1):
  invfact[i-1] = invfact[i] * i % P

def comb(n, k):
  # assert 0 <= k <= n
  return fact[n] * invfact[k] * invfact[n-k] % P

result = 0
for x in range(N+1):
  y, r = divmod(K - A * x, B)
  if r == 0 and 0 <= y <= N:
    result = result + comb(N, x) * comb(N, y) % P

print((result % P))
# Original solution
import math
L = input().split()
sp = 998244353
# approx 10**9
N = int(L[0])
A = int(L[1])
B = int(L[2])
K = int(L[3])

def binom(n,m):
	if n == m:
		return 1
	elif m == 1:
		return n
	elif m > n:
		return 0
	else:
		a = math.factorial(n)
		b = math.factorial(m)
		c = math.factorial(n - m)
		return a // (b*c)

# max is A*N + B*N
sum = 0
for i in range(0,N):
	if (K - A*i) % B == 0:
		j = (K-A*i)//B
		sum += (binom(N,i) * binom(N,j)) % sp

print(sum)



# Optimized solution
N,A,B,K = list(map(int,input().split()))
sprime = 998244353 # approx 10**9
A,B = sorted((A,B))
if (A+B) * N < K:
	print((0))
	exit()
 
def egcd(a, b):
    if a == 0:
        return (b, 0, 1)
    else:
        g, y, x = egcd(b % a, a)
        return (g, x - (b // a) * y, y)
 
def modinv(a, m):
    g, x, y = egcd(a, m)
    if g != 1:
        raise Exception('modular inverse does not exist')
    else:
        return x % m
 
memo = [0] * (N+1)
memo[0] = memo[N] = 1
for i in range(1,N//2+1):
	v = memo[i-1] * (N-i+1) * modinv(i,sprime)
	v = v % sprime
	memo[i] = memo[N-i] = v
 
# max is A*N + B*N
sum = 0
for b in range(N+1):
  a = (K-b*B)//A
  if 0 <= a <= N and a*A+b*B == K:
    sum += memo[a]*memo[b]
    sum %= sprime

print(sum)
# Original solution
N, A, B, K = list(map(int, input().split()))
mod = 998244353

fac = [1]
inv = [1]
for i in range(1, N + 1) :
    fac.append(fac[-1] * i % mod)
    inv.append(pow(fac[-1], mod - 2, mod))

ret = 0
for na in range(N + 1) :
    if (K - na * A) % B != 0 :
        continue
        
    nb = (K - na * A) // B
    for ng in range(max(0, na + nb - N), min(na, nb) + 1) :
        ret += fac[N] * inv[na-ng] * inv[nb-ng] * inv[ng] * inv[N-na-nb+ng] % mod
        ret %= mod
        
print(ret)


# Optimized solution
N, A, B, K = list(map(int, input().split()))
mod = 998244353

fac = [1]
inv = [1]
for i in range(1, N + 1) :
    fac.append(fac[-1] * i % mod)
    inv.append(pow(fac[-1], mod - 2, mod))

ret = 0
for na in range(N + 1) :
    if (K - na * A) % B != 0 :
        continue
        
    nb = (K - na * A) // B
    if 0 <= na <= N and 0 <= nb <= N :
        ret += fac[N] * inv[N-na] * inv[na] * fac[N] * inv[N-nb] * inv[nb]
        ret %= mod
        
print(ret)
# Original solution
from math import factorial as fac
def cc(n,r):
  def fac2(n,r):
    m = 1
    for i in range(n-r+1,n+1):
      m*=i
    return m
  r = min(r,n-r)
  return fac2(n,r)//fac(r)
n, a, b, k = list(map(int, input().split()))
mod = 998244353
ans = 0
for i in range(min(k//a+1,n+1)):
  jof = k-a*i
  if jof%b==0 and jof<=n*b:
    j = jof//b
    ans += cc(n,i)*cc(n,j)
    ans %= mod
print(ans)

# Optimized solution
n,a,b,k = list(map(int,input().split()))
mod = 998244353
fun = [1]*(n+1)
for i in range(1,n+1):
  fun[i] = fun[i-1]*i%mod
rev = [1]*(n+1)
rev[n] = pow(fun[n],mod-2,mod)
for i in range(n-1,0,-1):
  rev[i] = rev[i+1]*(i+1)%mod
def nCr(n,r):
  if r > n:
    return 0
  return fun[n]*rev[r]%mod*rev[n-r]%mod
ans = 0
for i in range(min(n,k//a)+1):
  j = (k-a*i)//b
  if i*a+j*b == k and 0<=j<=n:
    ans += nCr(n,i)*nCr(n,j)
    ans %= mod
print(ans)
# Original solution
n, a, b, k = list(map(int, input().split()))
MOD = 998244353

a_max =  min(k // a, n) + 1
pair = []
for i in range(a_max):
    if (k - a * i) % b == 0:
        j = (k - a * i) // b
        if j <= n:
            pair.append([i % MOD,j % MOD])

inv = [0] * (n + 1)
inv[1] = 1
for i in range(2, n + 1):
    inv[i] = (-(MOD // i) * inv[MOD % i]) % MOD

nCr = [0] * (n + 1)
nCr[0] = 1
for r in range(1, n + 1):
    nCr[r] = (nCr[r - 1] * (n - r + 1) * inv[r]) % MOD

result = 0
for p in pair:
    result += (nCr[p[0]] % MOD) * (nCr[p[1]] % MOD)

if result == 0:
    print((1))
else:
    print((result % MOD))

# Optimized solution
# coding: utf-8
# Your code here!

n, a, b, k = list(map(int, input().split()))
MOD = 998244353

inv = [0] * ( n + 1)
inv[1] = 1
for i in range(2, n + 1):
    inv[i] = (-(MOD // i) * inv[MOD % i]) % MOD

nCr = [0] * ( n + 1)
nCr[0] = 1
for i in range(1, n+1):
    nCr[i] = (nCr[i-1] * (n - i + 1) * inv[i]) % MOD

result = 0
maxA = min(n, k//a)
for i in range(maxA + 1):
    if ( k - a * i) % b == 0:
        j = ( k - a * i) // b
        if j <= n:
            result += (nCr[i] * nCr[j]) % MOD
if result == 0:
    print((1))
else:
    print((result % MOD))
# Original solution
N,A,B,K=list(map(int,input().split()))
P=998244353
def egcd(a, b):
    (x, lastx) = (0, 1)
    (y, lasty) = (1, 0)
    while b != 0:
        q = a // b
        (a, b) = (b, a % b)
        (x, lastx) = (lastx - q * x, x)
        (y, lasty) = (lasty - q * y, y)
    return (lastx, lasty, a)
def inv(x):
    return egcd(x,P)[0]
Fact=[0 for i in range(N+1)]
Finv=[0 for i in range(N+1)]
Fact[0]=1
Finv[0]=1
for i in range(1,N+1):
    Fact[i]=(i*Fact[i-1])%P
    Finv[i]=(Finv[i-1]*inv(i))%P
def NC(k):
    tmp=(Finv[k]*Finv[N-k])%P
    return (Fact[N]*tmp)%P
ans=0
for s in range(N+1):
    t=(K-s*A)//B
    if s*A+t*B==K and 0<=t<=N:
        ans+=(NC(s)*NC(t))%P
print((ans%P))



# Optimized solution
N,A,B,K=list(map(int,input().split()))
P=998244353
def egcd(a, b):
    (x, lastx) = (0, 1)
    (y, lasty) = (1, 0)
    while b != 0:
        q = a // b
        (a, b) = (b, a % b)
        (x, lastx) = (lastx - q * x, x)
        (y, lasty) = (lasty - q * y, y)
    return (lastx, lasty, a)
def inv(x):
    return egcd(x,P)[0]
Fact=[0 for i in range(N+1)]
Finv=[0 for i in range(N+1)]
Fact[0]=1
Finv[0]=1
for i in range(1,N+1):
    Fact[i]=(i*Fact[i-1])%P
    Finv[i]=(Finv[i-1]*inv(i))%P
ans=0
def NC(k):
    tmp=(Finv[k]*Finv[N-k])%P
    return (Fact[N]*tmp)%P
for x in range(N+1):
    y=(K-x*A)//B
    if x*A+y*B==K and 0<=y<=N:
        ans+=(NC(x)*NC(y))%P
print((ans%P))
# Original solution
P=998244353
MAX_N=3*10**5+1
def egcd(a, b):
    (x, lastx) = (0, 1)
    (y, lasty) = (1, 0)
    while b != 0:
        q = a // b
        (a, b) = (b, a % b)
        (x, lastx) = (lastx - q * x, x)
        (y, lasty) = (lasty - q * y, y)
    return (lastx, lasty, a)
def inv(x):
    return egcd(x,P)[0]
Fact=[0 for i in range(MAX_N+1)]
Finv=[0 for i in range(MAX_N+1)]
Fact[0]=1
Finv[0]=1
for i in range(MAX_N):
    Fact[i+1]=(Fact[i]*(i+1))%P
    Finv[i+1]=inv(Fact[i+1])%P
def C(n,k):
    return (Fact[n]*(Finv[k]*Finv[n-k])%P)%P
N,A,B,K=list(map(int,input().split()))
L=[]
for i in range(K//A+1):
    if (K-A*i)%B==0:
        j=(K-A*i)//B
        if 0<=i<=N and 0<=j<=N:
            L.append((i,j))
ans=0
for x,y in L:
    ans+=C(N,x)*C(N,y)
    ans%=P
print(ans)



# Optimized solution
N,A,B,K=list(map(int,input().split()))
P=998244353
def egcd(a, b):
    (x, lastx) = (0, 1)
    (y, lasty) = (1, 0)
    while b != 0:
        q = a // b
        (a, b) = (b, a % b)
        (x, lastx) = (lastx - q * x, x)
        (y, lasty) = (lasty - q * y, y)
    return (lastx, lasty, a)
def inv(x):
    return egcd(x,P)[0]
Fact=[0 for i in range(N+1)]
Finv=[0 for i in range(N+1)]
Fact[0]=1
Finv[0]=1
for i in range(1,N+1):
    Fact[i]=(i*Fact[i-1])%P
    Finv[i]=(Finv[i-1]*inv(i))%P
def NC(k):
    tmp=(Finv[k]*Finv[N-k])%P
    return (Fact[N]*tmp)%P
ans=0
for s in range(N+1):
    t=(K-s*A)//B
    if s*A+t*B==K and 0<=t<=N:
        ans+=(NC(s)*NC(t))%P
print((ans%P))
# Original solution
import sys
input = sys.stdin.buffer.readline
import copy
 
def main():
    N,A,B,K = list(map(int,input().split()))
    
    MOD = 998244353
    fac = [0 for _ in range(N+1)]
    fac[0],fac[1] = 1,1
    inv = copy.deepcopy(fac)
    invfac = copy.deepcopy(fac)
    
    for i in range(2,N+1):
        fac[i] = (fac[i-1]*i)%MOD
        inv[i] = MOD-(MOD//i)*inv[MOD%i]%MOD
        invfac[i] = (invfac[i-1]*inv[i])%MOD
 
    def coef(x,y):
        num = (((fac[x]*invfac[y])%MOD)*invfac[x-y]%MOD)
        return num
        
    ans = 0
    for i in range(min(K//max(A,B),N)+1):
        rest = K-i*max(A,B)
        if rest%min(A,B) == 0:
            j = rest//min(A,B)
            if j <= N:
                ans += coef(N,i)*coef(N,j)
                ans %= MOD
            
    print(ans)
    
if __name__ == "__main__":
    main()


# Optimized solution
import sys
input = sys.stdin.buffer.readline
import copy
 
def main():
    N,A,B,K = list(map(int,input().split()))
    
    MOD = 998244353
    fac = [0 for _ in range(N+1)]
    fac[0],fac[1] = 1,1
    invfac = copy.deepcopy(fac)
    
    for i in range(2,N+1):
        fac[i] = (fac[i-1]*i)%MOD
    
    invfac[-1] = pow(fac[-1],MOD-2,MOD)
    for i in range(N,0,-1):
        invfac[i-1] = (invfac[i]*i)%MOD
 
    def coef(x,y):
        num = (((fac[x]*invfac[y])%MOD)*invfac[x-y]%MOD)
        return num

    ans = 0
    for i in range(min(K//max(A,B),N)+1):
        rest = K-i*max(A,B)
        if rest%min(A,B) == 0:
            j = rest//min(A,B)
            if j <= N:
                ans += coef(N,i)*coef(N,j)
                ans %= MOD
            
    print(ans)
    
if __name__ == "__main__":
    main()
# Original solution
from math import factorial as fa
f=lambda a,b:fa(a)//(fa(a-b)*fa(b))
n,a,b,k=list(map(int,input().split()))
v=998244353
l=[]

for i in range(n+1):
  if not (k-a*i)%b:
    j=(k-a*i)//b
    if 0<=j<=n+1:
      l.extend([(i,j)])

s=0
for (i,j) in l:
  s+=(f(n,i)*f(n,j))%v

print(s)

# Optimized solution
n,a,b,k=list(map(int,input().split()))
mod=998244353
s=0
fac=[1]*(n+1)
rev=[1]*(n+1)

for i in range(1,n+1):
  fac[i]=i*fac[i-1]%mod
  rev[i]=pow(fac[i],mod-2,mod)

f=lambda a,b:(fac[a]*rev[a-b]*rev[b])%mod

for i in range(n+1):
  j,r=divmod((k-a*i),b)
  if r==0 and 0<=j<=n:
    s+=(f(n,i)*f(n,j))%mod

print((s%mod))
# Original solution
def comb(n,k,p):
	from math import factorial
	if n<0 or k<0 or n<k: return 0
	if n==0 or k==0: return 1
	a,b,c = factorial(n)%p,factorial(k)%p,factorial(n-k)%p
	return (a*pf(b,p-2,p)*pf(c,p-2,p))%p
def pf(a,b,p):
	if b==0: return 1
	if b%2==0:
		d = pf(a,b//2,p)
		return d**2%p
	if b%2==1: return (a*pf(a,b-1,p))%p
n,a,b,k = (int(i) for i in input().split())
mod,ans = 998244353,0
for i in range(n+1):
	if (k-a*i)%b==0:
		j = (k-a*i)//b
		ans = (ans+comb(n,i,mod)*comb(n,j,mod))%mod
print(ans)

# Optimized solution
n,a,b,k = (int(i) for i in input().split())
mod,ans = 998244353,0
fn,fk = [1]*n,[1]*n
for i in range(n-1): fn[i+1] = (fn[i]*(i+2))%mod
def power(n,k):
	if k==1: return n
	elif k%2==0: return power((n**2)%mod,k//2)
	else: return (n*power(n,k-1))%mod
def comb(n,k):
	if n==0 or k==0 or n==k: return 1
	else: return (((fn[n-1]*fk[n-k-1])%mod)*fk[k-1])%mod
fk[-1] = power(fn[-1],mod-2)
for i in range(2,n+1): fk[-i] = (fk[-i+1]*(n+2-i))%mod
for i in range(n+1):
	if (k-a*i)%b==0 and k-a*i>=0 and (k-a*i)//b<=n:
		j = (k-a*i)//b
		ans = (ans+(comb(n,i)*comb(n,j))%mod)%mod
print(ans)
# Original solution
M=998244353;n,A,B,K=list(map(int,input().split()));F=[1];a=0
for i in range(n):F+=[F[-1]*(n-i)*pow(i+1,M-2,M)%M]
for i in range(n+1):
 r=K-i*A;q=r//B
 if(r%B==0)&(0<=q<=n):a+=F[i]*F[q]
print((a%M))

# Optimized solution
M=998244353;n,A,B,K=list(map(int,input().split()));F=[1];a=0
for i in range(n):F+=[F[-1]*(n-i)*pow(i+1,M-2,M)%M]
for i in range(n+1):
 r=K-i*A
 if r%B==0and 0<=r//B<=n:a+=F[i]*F[r//B]
print((a%M))
# Original solution
N,A,B,K = list(map(int,input().split()))
if A<B: A,B = B,A
MOD = 998244353

def mul(a,b):
    return (a*b) % MOD

def pow(a,n): # a^n
    ret = 1
    mag = a
    while n > 0:
        if n & 1:
            ret = mul(ret, mag)
        mag = mul(mag, mag)
        n >>= 1
    return ret

def inv(a):
    return pow(a, MOD-2)

fac = [1]
#fac_inv = [1]
for n in range (1,N+10):
    f = mul(fac[n-1], n)
    fac.append(f)
#    fac_inv.append(inv(f))

mem = {}
def ncr(n,r):
    if (n,r) in mem:
        return mem[(n,r)]
    else:
        ret = mul(mul(fac[n], inv(fac[n-r])), inv(fac[r]))
        if len(mem) <= 10**5:
            mem[(n,r)] = ret
        return ret

ans = 0
for a in range(N+1):
    if a*A > K: break
    rem = K-a*A
    if rem%B: continue
    b = rem//B
    for green in range(min(a,b) + 1):
        red = a - green
        blue = b - green
        if red + green + blue > N: continue
        tmp = 1
        tmp = mul(tmp, ncr(N, red+green+blue))
        tmp = mul(tmp, ncr(red+green+blue, green+blue))
        tmp = mul(tmp, ncr(green+blue, blue))
        ans += tmp
print((ans % MOD))


# Optimized solution
N,A,B,K = list(map(int,input().split()))
if A<B: A,B = B,A
MOD = 998244353

def mul(a,b):
    return (a*b) % MOD

def pow(a,n): # a^n
    ret = 1
    mag = a
    while n > 0:
        if n & 1:
            ret = mul(ret, mag)
        mag = mul(mag, mag)
        n >>= 1
    return ret

def inv(a):
    return pow(a, MOD-2)

fac = [1]
fac_inv = [1]
for n in range (1,N+10):
    f = mul(fac[n-1], n)
    fac.append(f)
    fac_inv.append(inv(f))

mem = [0] * (N+1)
def ncr(n,r):
    if mem[r] > 0:
        return mem[r]
    else:
        ret = mul(mul(fac[n], fac_inv[n-r]), fac_inv[r])
        mem[r] = ret
        return ret

ans = 0
for a in range(N+1):
    if a*A > K: break
    rem = K-a*A
    if rem%B: continue
    b = rem//B
    if b > N: continue
    ans += mul(ncr(N,a), ncr(N,b))
print((ans % MOD))
# Original solution
N,A,B,K = list(map(int,input().split()))
MOD = 998244353

fac = [1,1] + [0]*N
finv = [1,1] + [0]*N
inv = [0,1] + [0]*N
for i in range(2,N+2):
    fac[i] = fac[i-1] * i % MOD
    inv[i] = -inv[MOD%i] * (MOD // i) % MOD
    finv[i] = finv[i-1] * inv[i] % MOD

def ncr(n,r):
    if n < r: return 0
    if n < 0 or r < 0: return 0
    return fac[n] * (finv[r] * finv[n-r] % MOD) % MOD

ans = 0
for a in range(N+1):
    if a*A > K: break
    rem = K-a*A
    if rem%B: continue
    b = rem//B
    for green in range(min(a,b) + 1):
        red = a - green
        blue = b - green
        if red + green + blue > N: continue
        tmp = 1
        tmp *= ncr(N, red+green+blue)
        tmp *= ncr(red+green+blue, green+blue)
        tmp *= ncr(green+blue, blue)
        ans += tmp % MOD
print((ans % MOD))

# Optimized solution
N,A,B,K = list(map(int,input().split()))
MOD = 998244353

if (A+B)*N < K:
    print((0))
    exit()

M = N+5
fac = [1,1] + [0]*M
finv = [1,1] + [0]*M
inv = [0,1] + [0]*M
for i in range(2,M+2):
    fac[i] = fac[i-1] * i % MOD
    inv[i] = -inv[MOD%i] * (MOD // i) % MOD
    finv[i] = finv[i-1] * inv[i] % MOD

def ncr(n,r):
    if n < r: return 0
    if n < 0 or r < 0: return 0
    return fac[n] * (finv[r] * finv[n-r] % MOD) % MOD

ans = 0
for a in range(N+1):
    rem = K - a*A
    if rem < 0: break
    b,m = divmod(rem, B)
    if m: continue
    ans += ncr(N,a) * ncr(N,b)
    ans %= MOD
print(ans)
# Original solution
N,A,B,K = list(map(int,input().split()))
MOD = 998244353

if (A+B)*N < K:
    print((0))
    exit()

M = N+5
fac = [1,1] + [0]*M
finv = [1,1] + [0]*M
inv = [0,1] + [0]*M
for i in range(2,M+2):
    fac[i] = fac[i-1] * i % MOD
    inv[i] = -inv[MOD%i] * (MOD // i) % MOD
    finv[i] = finv[i-1] * inv[i] % MOD

def ncr(n,r):
    if n < r: return 0
    if n < 0 or r < 0: return 0
    return fac[n] * (finv[r] * finv[n-r] % MOD) % MOD

ans = 0
for a in range(N+1):
    rem = K - a*A
    if rem < 0: break
    b,m = divmod(rem, B)
    if m: continue
    ans += ncr(N,a) * ncr(N,b)
    ans %= MOD
print(ans)

# Optimized solution
N,A,B,K = list(map(int,input().split()))
MOD = 998244353

MAXN = N+5
fac = [1,1] + [0]*MAXN
finv = [1,1] + [0]*MAXN
inv = [0,1] + [0]*MAXN
for i in range(2,MAXN+2):
    fac[i] = fac[i-1] * i % MOD
    inv[i] = -inv[MOD%i] * (MOD // i) % MOD
    finv[i] = finv[i-1] * inv[i] % MOD

def comb(n,r):
    if n < r: return 0
    if n < 0 or r < 0: return 0
    return fac[n] * (finv[r] * finv[n-r] % MOD) % MOD

ans = 0
for a in range(N+1):
    if (K-a*A)%B: continue
    b = (K-a*A)//B
    ans += comb(N,a) * comb(N,b)
    ans %= MOD
print(ans)
# Original solution
N, A, B, K = [int(_) for _ in input().split()]

def gcd(a, b):
    if a < b:
        a, b = b, a
    while b:
        a, b = b, a % b
    return a

L = 300000+1
M = 998244353

Fm = {}
inverseFm = {}
x = 1
for i in range(L):
    Fm[i] = x
    x = x * (i + 1) % M

def inverseFm(x, cache={}):
    if x in cache:
        return cache[x]
    result = pow(Fm[x], M - 2, M)
    cache[x] = result
    return result

def C(n, r):
    result = Fm[n] * inverseFm(r) * inverseFm(n - r) % M
    return result

def solve(N, A, B, K):
    g = gcd(gcd(A, B), K)
    A, B, K = A//g, B//g, K//g
    if A > B:
        A, B = B, A
    l = A * B // gcd(A, B)
    x = K
    for i in range(min(N + 1, l)):
        if x % B == 0:
            break
        x -= A
    else:
        return 0
    result = 0
    for a in range(i, N + 1, B):
        b = (K - A * a) // B
        for x in range(max(0, a-b), min(N-b+1, a+1)):
            y = a - x
            z = b - y
            d = Fm[N] * inverseFm(x) * inverseFm(y) * inverseFm(z) * inverseFm(N - x - y - z) % M
  #          print("s", x, y, z)
            result = (result + d) % M
    return result 

def solve0(N, A, B, K):
    result = 0
    for x in range(0, N + 1):
        for y in range(0, N + 1 - x):
            for z in range(0, N + 1 - x - y):
                if A * x + (A + B) * y + B * z == K:
 #                   print("s0", x, y, z)
                    d = Fm[N] * inverseFm(x) * inverseFm(y) * inverseFm(z) * inverseFm(N - x - y - z) % M
                    result = (result + d) % M
    return result                   


def check():
    import random
    for i in range(100):
        N = random.randint(1, 20)
        A = random.randint(1, 20)
        B = random.randint(1, 20)
        K = random.randint(1, 20)
        r = solve(N, A, B, K)
        r0 = solve0(N, A, B, K)
        if r != r0:
            print((N, A, B, K, r, r0))

#check()
print((solve(N, A, B, K)))
#print(solve0(N, A, B, K))


# Optimized solution
N, A, B, K = [int(_) for _ in input().split()]

L = 300000+1
M = 998244353

Fm = {}
inverseFm = {}
x = 1
for i in range(L):
    Fm[i] = x
    x = x * (i + 1) % M

def inverseFm(x, cache={}):
    if x in cache:
        return cache[x]
    result = pow(Fm[x], M - 2, M)
    cache[x] = result
    return result

def C(n, r):
    result = Fm[n] * inverseFm(r) * inverseFm(n - r) % M
    return result

def solve(N, A, B, K):
    result = 0
    for a in range(0, N+1):
        b, m = divmod(K - A * a, B)
        if m == 0 and 0 <= b <= N:
            result += C(N, a) * C(N, b) % M
    return result % M

print((solve(N, A, B, K)))