Problem:
 Score : 700 points    Problem Statement Takahashi has a tower which is divided into N layers.
Initially, all the layers are uncolored. Takahashi is going to paint some of the layers in red, green or blue to make a beautiful tower.
He defines the beauty of the tower as follows:   The beauty of the tower is the sum of the scores of the N layers, where the score of a layer is A if the layer is painted red, A+B if the layer is painted green, B if the layer is painted blue, and 0 if the layer is uncolored.   Here, A and B are positive integer constants given beforehand. Also note that a layer may not be painted in two or more colors.  Takahashi is planning to paint the tower so that the beauty of the tower becomes exactly K .
How many such ways are there to paint the tower? Find the count modulo 998244353 .
Two ways to paint the tower are considered different when there exists a layer that is painted in different colors, or a layer that is painted in some color in one of the ways and not in the other.      Constraints  1 ≤ N ≤ 3×10^5  1 ≤ A,B ≤ 3×10^5  0 ≤ K ≤ 18×10^{10}  All values in the input are integers.         Input Input is given from Standard Input in the following format:  N  A  B  K       Output Print the number of the ways to paint tiles, modulo 998244353 .        Sample Input 1 4 1 2 5      Sample Output 1 40  In this case, a red layer worth 1 points, a green layer worth 3 points and the blue layer worth 2 points. The beauty of the tower is 5 when we have one of the following sets of painted layers:   1 green, 1 blue  1 red, 2 blues  2 reds, 1 green  3 reds, 1 blue   The total number of the ways to produce them is 40 .       Sample Input 2 2 5 6 0      Sample Output 2 1  The beauty of the tower is 0 only when all the layers are uncolored. Thus, the answer is 1 .       Sample Input 3 90081 33447 90629 6391049189      Sample Output 3 577742975

Solution (Python):
M=998244353;n,A,B,K=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==0)&(0<=r//B<=n):a+=F[i]*F[r//B]
print(a%M)

Solution (Python):
def modpow(a,n,p):
    if n==0:
        return 1
    x=modpow(a,n//2,p)
    x=(x*x)%p
    if (n%2)==1:
        x=(x*a)%p
    return x%p
def modinv(a,p):
    return modpow(a,p-2,p)
N,A,B,K=map(int,input().split())
p=998244353
fracmod=[1]*(N+1)
fracmodinv=[1]*(N+1)
ans=0
i=0
j=0
for i in range(1,N+1):
    fracmod[i]=(fracmod[i-1]*i)%p

for i in range(N+1):
    fracmodinv[i]=modinv(fracmod[i],p)
for i in range(N+1):
    forb=K-i*A      
    if (forb%B==0)and(0<=(forb//B)<=N):
        j=forb//B
        
        ans=(ans+((((fracmod[N]*fracmodinv[N-i]*fracmodinv[i])%p)*(fracmod[N]*fracmodinv[N-j]*fracmodinv[j]))%p))%p
print(ans)

Solution (Python):
from fractions import gcd
n,a,b,k=map(int,input().split())
def cmb(n, r):
    if (r < 0) or (n < r):
        return 0
    r = min(r, n - r)
    return fact[n] * factinv[r] * factinv[n-r] % mod

N = n 
mod = 998244353
fact = [1, 1]  
factinv = [1, 1]  
inv = [0, 1]  

for i in range(2, N + 1):
    fact.append((fact[-1] * i) % mod)
    inv.append((-inv[mod % i] * (mod // i)) % mod)
    factinv.append((factinv[-1] * inv[-1]) % mod)

ls=[]
res=0

for i in range(n+1):
    
    if (k-a*i)%b==0:
      j= (k-a*i)//b
      
      
      if j<=n:
        ls.append((i,j))

for tup in ls:
  res+=cmb(n,tup[0])*cmb(n,tup[1])
  
print(res%mod)

Solution (Python):
n, a, b, k = map(int, input().split())
mod = 998244353

MAX = 3 * 10 ** 5

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

inv = [1] * (MAX + 1)
for i in range(2, MAX + 1):
    inv[i] = inv[mod % i] * (mod - mod // i) % mod

fact_inv = [1] * (MAX + 1)
for i in range(1, MAX + 1):
    fact_inv[i] = fact_inv[i-1] * inv[i] % mod


def comb(n, k):
    if n < k:
        return 0
    return fact[n] * fact_inv[n-k] * fact_inv[k] % mod


ans = 0
for ea in range(n + 1):
    ebb = k - ea * a
    if ebb % b == 0:
        eb = ebb // b
        if 0 <= eb <= n:
            ans += comb(n, ea) * comb(n, eb)
            ans %= mod

print(ans)

Solution (Python):
from fractions import gcd

N, A, B, K = map(int, input().split())
if K == 0:
    print(1)
    exit()
c = gcd(A, B)
if K % c != 0:
    print(0)
    exit()
A //= c; B //= c; K //= c
mod = 998244353

fact = [1] * (N+1)
fact_inv = [1] * (N+1)
for i in range(1, N+1):
    fact[i] = i * fact[i-1] % mod
fact_inv[N] = pow(fact[N], mod-2, mod)
for i in range(1, N+1)[::-1]:
    fact_inv[i-1] = i * fact_inv[i] % mod
comb = lambda n, k: fact[n] * fact_inv[k] * fact_inv[n-k] % mod

ans = 0
for x in range(N+1):
    y, r = divmod(K - A * x, B)
    if r != 0:
        continue
    if 0 <= y <= N:
        ans = (ans + comb(N, x) * comb(N, y)) % mod
print(ans)

Solution (Python):
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 = 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)

Solution (Python):
n, a, b, k = map(int, input().split())

def cmb1(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 )

ans = 0
for i in range(n+1):
    x = a*i
    y = k-x
    if y < 0:
        continue
    if y%b == 0:
        j = y//b
        ans += cmb1(n, i, mod)*cmb1(n, j, mod)
print(ans%mod)

Solution (Python):
n, a, b, k = map(int, input().split())

def cmb1(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 )

ans = 0
for i in range(n+1):
    x = a*i
    y = k-x
    if y < 0:
        continue
    if y%b == 0:
        j = y//b
        ans += cmb1(n, i, mod)*cmb1(n, j, mod)
print(ans%mod)

Solution (Python):
N, A, B, K = map(int, input().split())
mod = 998244353





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


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)
ans = 0
for i in range(N+1):
    D = (K - A*i)
    if D % B == 0 and D//B <= N:
        ans += cmb(N, i, mod) * cmb(N, D//B, mod)
        ans %= mod

print(ans)

Solution (Python):
MOD = 998244353
n, a, b, k = map(int, input().split())
memo = [0]*(n//2+1)


def com(n, r):
    if n < r or r < 0:
        return 0
    if n-r < r:
        r = n-r
    if r == 0:
        memo[r] = 1
        return memo[r]
    if memo[r] != 0:
        return memo[r]
    X = Y = 1
    for i in range(1, r+1):
        Y = Y*i % MOD
        X = X*(n-i+1) % MOD
        memo[i] = X*pow(Y, MOD-2, MOD)
    return memo[r]


com(n, n//2)
ans = 0
for x in range(n+1):
    y, r = divmod(k-a*x, b)
    if y < 0:
        continue
    if r == 0:
        ans += com(n, x)*com(n, y)
        ans %= MOD
print(ans % MOD)

Solution (Python):
MOD = 998244353
n, a, b, k = map(int, input().split())
pair = []
for x in range(n+1):
    y, r = divmod(k-a*x, b)
    if y < 0:
        continue
    if r == 0:
        pair.append((x, y))


memo = [0]*n


def com(n, r):
    X = Y = 1
    if n < r:
        return 0
    if n-r < r:
        r = n-r
    if r == 0:
        memo[r] = 1
        return memo[r]
    if memo[r] != 0:
        return memo[r]
    for i in range(1, r+1):
        Y = Y*i % MOD
        X = X*(n-i+1) % MOD
        memo[i] = X*pow(Y, MOD-2, MOD)
    return memo[r]

com(n, n//2)

ans = 0
for x, y in pair:
    ans += com(n, x)*com(n, y) % MOD
print(ans % MOD)

Solution (Python):
n, a, b, k = map(int, input().split())

mod = 998244353
factorial = [1]
inverse = [1]
for i in range(1, n+1):
    factorial.append(factorial[-1] * i % mod)
    inverse.append(pow(factorial[-1], mod - 2, mod))

def comb(n, r, mod):
    if n < r or r < 0: return 0
    elif r == 0: return 1
    return factorial[n] * inverse[r] * inverse[n - r] % mod

ans = 0
for i in range(n+1):
    if a*i > k: break
    sum_b = k - a*i
    if sum_b % b != 0: continue
    j = sum_b // b
    if j > n: continue
    ans = (ans + comb(n, i, mod)*comb(n, j, mod)) % mod
print(ans)

Solution (Python):
n, a, b, k = map(int, input().split())

mod = 998244353
factorial = [1]
inverse = [1]
for i in range(1, n+1):
    factorial.append(factorial[-1] * i % mod)
    inverse.append(pow(factorial[-1], mod - 2, mod))

def comb(n, r, mod):
    if n < r or r < 0: return 0
    elif r == 0: return 1
    return factorial[n] * inverse[r] * inverse[n - r] % mod

ans = 0
for i in range(n+1):
    if a*i > k: break
    sum_b = k - a*i
    if sum_b % b != 0: continue
    j = sum_b // b
    if j > n: continue
    ans = (ans + comb(n, i, mod)*comb(n, j, mod)) % mod
print(ans)

Solution (Python):
import sys
input = sys.stdin.readline
from collections import *

MAX = 3*10**5+100
MOD = 998244353
fact = [0]*MAX 
inv = [0]*MAX 
finv = [0]*MAX 
fact[0] = 1
fact[1] = 1
finv[0] = 1
finv[1] = 1
inv[1] = 1
    
for i in range(2, MAX):
    fact[i] = fact[i-1]*i%MOD
    inv[i] = MOD-inv[MOD%i]*(MOD//i)%MOD
    finv[i] = finv[i-1]*inv[i]%MOD

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

N, A, B, K = map(int, input().split())
ans = 0

for i in range(N+1):
    if (K-i*A)%B!=0:
        continue
    
    j = (K-i*A)//B
    ans += C(N, i)*C(N, j)
    ans %= MOD

print(ans)

Solution (Python):
N,A,B,K = 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)

Solution (Python):
class Combination:
    """階乗とその逆元のテーブルをO(N)で事前作成し、組み合わせの計算をO(1)で行う"""
    def __init__(self, n, MOD):
        self.fact = [1]
        for i in range(1, n + 1):
            self.fact.append(self.fact[-1] * i % MOD)
        self.inv_fact = [0] * (n + 1)
        self.inv_fact[n] = pow(self.fact[n], MOD - 2, MOD)
        for i in reversed(range(n)):
            self.inv_fact[i] = self.inv_fact[i + 1] * (i + 1) % MOD
        self.MOD = MOD

    def factorial(self, k):
        """k!を求める O(1)"""
        return self.fact[k]

    def inverse_factorial(self, k):
        """k!の逆元を求める O(1)"""
        return self.inv_fact[k]

    def permutation(self, k, r):
        """kPrを求める O(1)"""
        if k < r:
            return 0
        return (self.fact[k] * self.inv_fact[k - r]) % self.MOD

    def combination(self, k, r):
        """kCrを求める O(1)"""
        if k < r:
            return 0
        return (self.fact[k] * self.inv_fact[k - r] * self.inv_fact[r]) % self.MOD

    def combination2(self, k, r):
        """kCrを求める O(r) kが大きいが、r <= nを満たしているときに使用"""
        if k < r:
            return 0
        res = 1
        for l in range(r):
            res *= (k - l)
            res %= self.MOD
        return (res * self.inv_fact[r]) % self.MOD


n, a, b, k = map(int, input().split())
MOD = 998244353
comb = Combination(10 ** 6, MOD)

ans = 0
for cnt_a in range(0, n + 1):
    point = a * cnt_a
    nokori = k - point
    if nokori >= 0 and nokori % b == 0:
        cnt_b = nokori // b
        if cnt_b > n:
            continue
        ans += comb.combination(n, cnt_a) * comb.combination(n, cnt_b)
        ans %= MOD
print(ans % MOD)

Solution (Python):
n,a,b,k=map(int,input().split())
mod=998244353
f=[1]
for i in range(n):f+=[f[-1]*(i+1)%mod]
def comb(a,b):return f[a]*pow(f[b],mod-2,mod)*pow(f[a-b],mod-2,mod)%mod
ans=0
for ac in range(k//a+1):
  asum=ac*a
  if (k-asum)%b:continue
  bc=(k-asum)//b
  if max(ac,bc)>n:continue
  ans=(ans+comb(n,ac)*comb(n,bc))%mod
print(ans)

Solution (Python):
n,a,b,k = map(int, input().split())

mod = 998244353
NUM = 10**6
g1 = [1, 1]
g2 = [1, 1]
inv = [0, 1]

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

def comb(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

ans=0
for i in range(n+1):
    if (k-i*a)%b == 0 and (k-i*a)>=0:
        j=(k-i*a)//b
        if j>n:
            continue
        ans+=(comb(n,i,mod)*comb(n,j,mod))
        ans%=mod
print(ans)

Solution (Python):
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 = 3*10**5
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=map(int,input().split())
ans=0
for i in range(0,N+1):
    if (K-A*i)%B==0 and N>=(K-A*i)//B>=0:
        j=(K-A*i)//B
        ans+=cmb(N,i,mod)*cmb(N,j,mod)
        ans%=mod

print(ans)

Solution (Python):
mod = 998244353


class Combination:
    """
    SIZEが10^6程度以下の二項係数を何回も呼び出したいときに使う
    使い方:
    comb = Combination(SIZE, MOD)
    comb(10, 3) => 120
    """

    def __init__(self, N, MOD=10 ** 9 + 7):
        self.MOD = MOD
        self.fact, self.inv = self._make_factorial_list(N)

    def __call__(self, n, k):
        if k < 0 or k > n:
            return 0
        res = self.fact[n] * self.inv[k] % self.MOD
        res = res * self.inv[n - k] % self.MOD
        return res

    def _make_factorial_list(self, N):
        fact = [1] * (N + 1)
        inv = [1] * (N + 1)
        MOD = self.MOD
        for i in range(1, N + 1):
            fact[i] = (fact[i - 1] * i) % MOD
        inv[N] = pow(fact[N], MOD - 2, MOD)
        for i in range(N, 0, -1):
            inv[i - 1] = (inv[i] * i) % MOD
        return fact, inv


if __name__ == "__main__":
    N, A, B, K = map(int, input().split())
    U = 3 * 10 ** 5 + 10
    comb = Combination(U, mod)

    ans = 0
    if A > B:
        A, B = B, A
    for i in range(N + 1):
        X = A * i
        if X > K:
            continue
        if (K - X) % B != 0:
            continue
        if (K - X) // B > N:
            continue
        j = (K - X) // B

        red = comb(N, i)
        blue = comb(N, j)
        ans = (ans + red * blue % mod) % mod
    print(ans)

Solution (Python):
N, A, B, K = map(int, input().split())
MOD = 998244353

g1 = [1, 1]
g2 = [1, 1]
inverse = [0, 1]
for i in range(2, 3*10**5+1):
    g1.append((g1[-1] * i) % MOD)
    inverse.append((-inverse[MOD % i] * (MOD // i)) % MOD)
    g2.append((g2[-1] * inverse[-1]) % MOD)

def comb(n, r, mod=MOD):
    if r < 0 or r > n:
        return 0
    r = min(r, n - r)
    return g1[n] * g2[r] * g2[n-r] % mod

ans = 0
for i in range(N+1):
  if (K - i*A) % B == 0:
    j = (K - i*A) // B
    ans += comb(N, i) * comb(N, j) % MOD
    ans %= MOD

print(ans)

Solution (Python):
import sys

readline = sys.stdin.readline
MOD = 10 ** 9 + 7
INF = float('INF')
sys.setrecursionlimit(10 ** 5)


def main():
    n, a, b, k = list(map(int, readline().split()))
    M = 998244353

    h = (n + 1) // 2
    l = h + 1
    com = [0] * l
    com[0] = 1

    for i in range(1, l):
        t = com[i - 1] * (n + 1 - i)
        t *= pow(i, M - 2, M)
        t %= M
        com[i] = t

    ans = 0
    rng = min(k // a + 1, n + 1)

    for x in range(rng):
        res = k - a * x
        if res % b == 0:
            y = res // b
            if y > n:
                continue
            if x > h:
                x = n - x
            if y > h:
                y = n - y
            add = com[x] * com[y]
            ans += add
            ans %= M

    print(ans)

if __name__ == '__main__':
    main()

Solution (Python):
N,A,B,K=map(int,input().split())


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 
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 )

ans=0
for i in range(N+1):
    if (K-(A*i))%B==0:
        ans+=cmb(N,i,mod)*cmb(N,(K-(A*i))//B,mod)
print(ans%mod)

Solution (Python):

import sys
import numpy as np

sr = lambda: sys.stdin.readline().rstrip()
ir = lambda: int(sr())
lr = lambda: list(map(int, sr().split()))

MOD = 998244353
N, A, B, K = lr()
answer = 0



def F(x, y):
    return comb[x] * comb[y] % MOD

def cumprod(arr, MOD):
    L = len(arr); Lsq = int(L**.5+1)
    arr = np.resize(arr, Lsq**2).reshape(Lsq, Lsq)
    for n in range(1, Lsq):
        arr[:, n] *= arr[:, n-1]; arr[:, n] %= MOD
    for n in range(1, Lsq):
        arr[n] *= arr[n-1, -1]; arr[n] %= MOD
    return arr.ravel()[:L]

def make_fact(U, MOD):
    x = np.arange(U, dtype=np.int64); x[0] = 1
    fact = cumprod(x, MOD)
    x = np.arange(U, 0, -1, dtype=np.int64); x[0] = pow(int(fact[-1]), MOD-2, MOD)
    fact_inv = cumprod(x, MOD)[::-1]
    return fact, fact_inv

U = N + 1 
fact, fact_inv = make_fact(U, MOD)
fact_N = fact[N]
comb = fact_N * fact_inv % MOD * fact_inv[::-1] % MOD

for x in range(N+1):
    temp = K-A*x
    q, r = divmod(temp, B)
    if r != 0:
        continue
    if q > N:
        continue
    if q < 0:
        break
    answer += F(x, q)

print(answer % MOD)

Solution (Python):

import sys
import numpy as np

sr = lambda: sys.stdin.readline().rstrip()
ir = lambda: int(sr())
lr = lambda: list(map(int, sr().split()))

MOD = 998244353
N, A, B, K = lr()
answer = 0



def F(x, y):
    return comb[x] * comb[y]

def cumprod(arr, MOD):
    L = len(arr); Lsq = int(L**.5+1)
    arr = np.resize(arr, Lsq**2).reshape(Lsq, Lsq)
    for n in range(1, Lsq):
        arr[:, n] *= arr[:, n-1]; arr[:, n] %= MOD
    for n in range(1, Lsq):
        arr[n] *= arr[n-1, -1]; arr[n] %= MOD
    return arr.ravel()[:L]

def make_fact(U, MOD):
    x = np.arange(U, dtype=np.int64); x[0] = 1
    fact = cumprod(x, MOD)
    x = np.arange(U, 0, -1, dtype=np.int64); x[0] = pow(int(fact[-1]), MOD-2, MOD)
    fact_inv = cumprod(x, MOD)[::-1]
    return fact, fact_inv

U = N + 1 
fact, fact_inv = make_fact(U, MOD)
fact_N = fact[N]
comb = fact_N * fact_inv % MOD * fact_inv[::-1] % MOD

for x in range(N+1):
    temp = K-A*x
    q, r = divmod(temp, B)
    if r != 0:
        continue
    if q > N:
        continue
    if q < 0:
        break
    answer += F(x, q)
    answer %= MOD

print(answer % MOD)

Solution (Python):

import sys
import numpy as np

sr = lambda: sys.stdin.readline().rstrip()
ir = lambda: int(sr())
lr = lambda: list(map(int, sr().split()))

MOD = 998244353
N, A, B, K = lr()
answer = 0



def F(x, y):
    return comb[x] * comb[y] % MOD

def cumprod(arr, MOD):
    L = len(arr); Lsq = int(L**.5+1)
    arr = np.resize(arr, Lsq**2).reshape(Lsq, Lsq)
    for n in range(1, Lsq):
        arr[:, n] *= arr[:, n-1]; arr[:, n] %= MOD
    for n in range(1, Lsq):
        arr[n] *= arr[n-1, -1]; arr[n] %= MOD
    return arr.ravel()[:L]

def make_fact(U, MOD):
    x = np.arange(U, dtype=np.int64); x[0] = 1
    fact = cumprod(x, MOD)
    x = np.arange(U, 0, -1, dtype=np.int64); x[0] = pow(int(fact[-1]), MOD-2, MOD)
    fact_inv = cumprod(x, MOD)[::-1]
    return fact, fact_inv

U = N + 1 
fact, fact_inv = make_fact(U, MOD)
fact_N = fact[N]
comb = fact_N * fact_inv % MOD * fact_inv[::-1] % MOD

for x in range(N+1):
    temp = K-A*x
    q, r = divmod(temp, B)
    if r != 0:
        continue
    if q > N:
        continue
    if q < 0:
        break
    answer += F(x, q)
    answer %= MOD

print(answer % MOD)

Solution (Python):

import sys
import numpy as np

sr = lambda: sys.stdin.readline().rstrip()
ir = lambda: int(sr())
lr = lambda: list(map(int, sr().split()))

MOD = 998244353
N, A, B, K = lr()
answer = 0



def cmb(n, k):
    if k < 0 or k > n: return 0 
    return fact_N * fact_inv[k] % MOD * fact_inv[n-k] % MOD

def cumprod(arr, MOD):
    L = len(arr); Lsq = int(L**.5+1)
    arr = np.resize(arr, Lsq**2).reshape(Lsq, Lsq)
    for n in range(1, Lsq):
        arr[:, n] *= arr[:, n-1]; arr[:, n] %= MOD
    for n in range(1, Lsq):
        arr[n] *= arr[n-1, -1]; arr[n] %= MOD
    return arr.ravel()[:L]

def make_fact(U, MOD):
    x = np.arange(U, dtype=np.int64); x[0] = 1
    fact = cumprod(x, MOD)
    x = np.arange(U, 0, -1, dtype=np.int64); x[0] = pow(int(fact[-1]), MOD-2, MOD)
    fact_inv = cumprod(x, MOD)[::-1]
    return fact, fact_inv

U = N + 1 
fact, fact_inv = make_fact(U, MOD)
fact_N = fact[N]
comb = fact_N * fact_inv % MOD * fact_inv[::-1] % MOD

for x in range(N+1):
    temp = K-A*x
    if temp%B != 0:
        continue
    y = temp // B
    if y > N:
        continue
    if y < 0:
        break
    answer += comb[x] * comb[y]
    answer %= MOD

print(answer % MOD)

Solution (Python):

import sys
import numpy as np

sr = lambda: sys.stdin.readline().rstrip()
ir = lambda: int(sr())
lr = lambda: list(map(int, sr().split()))

MOD = 998244353
N, A, B, K = lr()
answer = 0



def cmb(n, k):
    if k < 0 or k > n: return 0 
    return fact_N * fact_inv[k] % MOD * fact_inv[n-k] % MOD

def cumprod(arr, MOD):
    L = len(arr); Lsq = int(L**.5+1)
    arr = np.resize(arr, Lsq**2).reshape(Lsq, Lsq)
    for n in range(1, Lsq):
        arr[:, n] *= arr[:, n-1]; arr[:, n] %= MOD
    for n in range(1, Lsq):
        arr[n] *= arr[n-1, -1]; arr[n] %= MOD
    return arr.ravel()[:L]

def make_fact(U, MOD):
    x = np.arange(U, dtype=np.int64); x[0] = 1
    fact = cumprod(x, MOD)
    x = np.arange(U, 0, -1, dtype=np.int64); x[0] = pow(int(fact[-1]), MOD-2, MOD)
    fact_inv = cumprod(x, MOD)[::-1]
    return fact, fact_inv

U = N + 100 
fact, fact_inv = make_fact(U, MOD)
fact_N = fact[N]

for x in range(N+1):
    temp = K-A*x
    if temp%B != 0:
        continue
    y = temp // B
    if y > N:
        continue
    if y < 0:
        break
    answer += cmb(N, x) * cmb(N, y)
    answer %= MOD

print(answer % MOD)

Solution (Python):
mod = 998244353

def comb(n, r):
    if 2 * r > n: return comb(n, n-r)
    return fac[n] * inv[r] * inv[n-r] % mod

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

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

ans = 0
for i in range(n+1):
    Num = (k-a*i)
    if (0 <= Num//b) and (Num//b <= n) and (Num%b == 0):
      j = Num//b 
      ans += comb(n, i) * comb(n, j) 
print(ans%mod)

Solution (Python):
mod = 998244353

def comb(n, r):
    if 2 * r > n: return comb(n, n-r)
    return fac[n] * inv[r] * inv[n-r] % mod

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

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

ans = 0
for i in range(n+1):
    Num = (k-a*i)
    if (0 <= Num//b) and (Num//b <= n) and (Num%b == 0):
      j = Num//b 
      ans += comb(n, i) * comb(n, j) 
print(ans%mod)

Solution (Python):
import sys
import math
from collections import defaultdict
from bisect import bisect_left, bisect_right

sys.setrecursionlimit(10**7)
def input():
    return sys.stdin.readline()[:-1]

mod = 10**9 + 7

def I(): return int(input())
def LI(): return list(map(int, input().split()))
def LIR(row,col):
    if row <= 0:
        return [[] for _ in range(col)]
    elif col == 1:
        return [I() for _ in range(row)]
    else:
        read_all = [LI() for _ in range(row)]
        return map(list, zip(*read_all))



mod = 998244353



def fact_all(n,M=mod):
    f = [1]*(n+1)
    ans = 1
    for i in range(1,n+1):
        ans = (ans*i)%M
        f[i] = ans
    return f


def fact_inv_all(fact_all,M=mod):
    N = len(fact_all)
    finv = [0]*N
    finv[-1] = pow(fact_all[-1],M-2,M)
    for i in range(N-1)[::-1]:
        finv[i] = finv[i+1]*(i+1)%M
    return finv


def nCk(n,k,fact_list=[],inv_list=[],M=mod):
    return (((fact_list[n]*inv_list[k])%M)*inv_list[n-k])%M

N,A,B,K = LI()



fl = fact_all(N)
il = fact_inv_all(fl)

ans = 0
for i in range(N+1):
    if (K-A*i)%B == 0:
        num = (K-A*i)//B
        if 0 <= num <= N:
            ans += nCk(N,i,fl,il)*nCk(N,num,fl,il)
            ans %= mod

print(ans)

Solution (Python):
import sys
input = sys.stdin.buffer.readline
import copy
 
def main():
    N,A,B,K = 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()

Solution (Python):
import sys
input = sys.stdin.buffer.readline
import copy
 
def main():
    N,A,B,K = 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()

Solution (Python):

class Combination():
    def __init__(self, N:int, P:int):
        self.N = N
        self.P = P
        
        self.fact = [1, 1]
        
        self.factinv = [1, 1]
        
        self.inv = [0, 1]

        for i in range(2, N+1):
            self.fact.append((self.fact[-1] * i) % P)
            self.inv.append((-self.inv[P % i] * (P // i)) % P)
            self.factinv.append((self.factinv[-1] * self.inv[-1]) % P)
    
    def getComb(self, n:int, k:int):
        if (k < 0) or (n < k): return 0
        k = min(k, n - k)
        return self.fact[n] * self.factinv[k] * self.factinv[n-k] % self.P


N,A,B,K = map(int, input().split())

MOD = 998244353
C = Combination(N, MOD)

ans = 0
for a in range(N+1):
    
    if (K - a * A) % B == 0 and 0 <= (K - a * A) // B <= N:
        b = (K - a * A) // B
        ans = (ans + C.getComb(N, a) * C.getComb(N, b)) % MOD

print(ans)

Solution (Python):
n,a,b,k=map(int,input().split())
ans=0
comb=[1]
c=1
mod=998244353
for i in range(1,n+1):
  c*=n-i+1
  c*=pow(i,mod-2,mod)
  c%=mod
  comb.append(c)
for s in range(n+1):
  if k<a*s:
    break
  bt=k-a*s
  if bt%b==0 and 0<=bt//b<=n:
    t=bt//b
    ans=(ans+comb[s]*comb[t])%mod
print(ans)

Solution (Python):
N, A, B, K = map(int, input().split())
MOD = 998244353

fact_N  = 1
for i in range(1, N+1):
    fact_N = (fact_N*i)%MOD
fact_inv = [0 for i in range(N+1)]
fact_inv[N] = pow(fact_N, MOD-2, MOD)
for i in range(N-1, -1, -1):
    fact_inv[i] = (fact_inv[i+1] * (i+1))%MOD

comb = [ (((fact_N*fact_inv[i])%MOD)*fact_inv[N-i])%MOD for i in range(N+1)]

ans = 0
for a in range(N+1):
    tmp = K-a*A
    if (tmp%B == 0) and (0 <= tmp//B) and (tmp//B <= N):
        b = tmp//B
        ans = (ans+comb[a]*comb[b])%MOD
print(ans)

Solution (Python):
N,A,B,K = list(map(int,input().split()))


mod = 998244353

ans = 0

comb = [1]
c = 1

for i in range(1,N+1):
    c *= N - i + 1
    c *= pow(i,mod-2,mod)
    c %= mod
    comb.append(c)

for a in range(N+1):
    if K<a*A:
        break
    if (K-a*A)%B==0:
        b = (K-a*A)//B
        if b <= N:
            ans += comb[a]*comb[b]
            ans %= mod

print(ans)

Solution (Python):
N,A,B,K=map(int,input().split())
r=0
m=998244353
max_n=3*10**5
fac=[1]*(max_n+1)
inv=[1]*(max_n+1)
ifac=[1]*(max_n+1)
for n in range(2,max_n+1):
    fac[n]=(fac[n-1]*n)%m
    inv[n]=m-inv[m%n]*(m//n)%m
    ifac[n]=(ifac[n-1]*inv[n])%m
def comb(n,k):
    if n<k or n<0 or k<0:
        return 0
    return (fac[n]*ifac[k]*ifac[n-k])%m
for a in range(N+1):
    Bb=K-A*a
    if Bb>=0 and Bb%B==0:
        b=Bb//B
        if b<=N:
            r=(r+comb(N,a)*comb(N,b))%m
print(r)

Solution (Python):
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 = 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)

Solution (Python):
U = 3*10**5
MOD = 998244353
 
fact = [1]*(U+1)
fact_inv = [1]*(U+1)
 
for i in range(1,U+1):
    fact[i] = (fact[i-1]*i)%MOD
fact_inv[U] = pow(fact[U], MOD-2, MOD)
 
for i in range(U,0,-1):
    fact_inv[i-1] = (fact_inv[i]*i)%MOD

def comb(n, k):
    if k < 0 or k > n:
        return 0
    z = fact[n]
    z *= fact_inv[k]
    z %= MOD
    z *= fact_inv[n-k]
    z %= MOD
    return z
  
n, a, b, k = map(int, input().split())
ans = 0
for i in range(n+1):
  r = k - i*a
  if r%b==0 and 0 <= r//b <= n:
    ans += comb(n, i) * comb(n, r//b)
    ans %= MOD
print(ans)

Solution (Python):
mod = 998244353
def extgcd(a, b, c):
  x, y, u, v, k, l = 1, 0, 0, 1, a, b
  while l != 0:
    x, y, u, v = u, v, x - u * (k // l), y - (k // l)
    k, l = l, k % l
  if c == 1:
    return x
  elif c == 2:
    return y
  elif c == 3:
    return k  
def inved(x):
  return extgcd(x, mod, 1)

fact = [1 for i in range(300001)]
invf = [1 for i in range(300001)]
for i in range(300000):
  fact[i+1] = (fact[i] * (i + 1)) % mod
invf[300000] = inved(fact[300000])
for i in range(300000, 0, -1):
  invf[i-1] = (invf[i] * i) % mod



N, A, B, K = map(int, input().split())
g = extgcd(A, B, 3)
uppern = (fact[N] * fact[N]) % mod
S = 0
if K % g != 0:
  S = 0
else:
  left = max(-((B*N-K)//A), 0)
  right = min(K // A, N)
  for i in range(left, right + 1):
    if (K - A * i) % B == 0:
      l = (K - A * i) // B
      prod = (invf[i] * invf[N-i]) % mod
      prod *= (invf[l] * invf[N-l]) % mod
      prod %= mod
      prod *= uppern
      prod %= mod
      S += prod
      S %= mod
  print(S)