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):
import numpy as np

mod = 998244353


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


def ii():
    return int(input())


def lii():
    return list(map(int, input().split(' ')))


def lvi(N):
    l = []
    for _ in range(N):
        l.append(ii())
    return l


def lv(N):
    l = []
    for _ in range(N):
        l.append(input())
    return l


def lcm(a, b):
    """
    最小公倍数
    """
    return a * b / gcd(a,b)

















def solve_eucl(A, B, K, N):
    for _x in range(N):
        if (K - A * _x) % B == 0:
            x = _x
            y = (K - A * _x) // B
            return x, y


def modinv(a, m):
    orgm = m
    x, lastx, y, lasty = 0, 1, 1, 0
    while m:
        a, (quotient, m) = m, divmod(a, m)
        x, lastx = lastx - quotient * x, x
        y, lasty = lasty - quotient * y, y
    if a != 1:
        raise ValueError
    return lastx % orgm


def fastcombmod(m, n):
    nn = min(n, m - n)
    num = 1
    div = 1
    for i in range(nn):
        num *= m - i
        if num > mod:
            num %= mod
        div *= i + 1
        if div > mod:
            div %= mod
    return (num * modinv(div, mod)) % mod


def combmod_table(m):
    num = 1
    div = 1
    t = [1]
    for i in range(m // 2 ):
        num *= m - i
        if num > mod:
            num %= mod
        div *= i + 1
        if div > mod:
            div %= mod

        t.append((num * modinv(div, mod)) % mod)

    return t


def main():
    N, A, B, K = lii()
    if K == 0: return 1

    d = gcd(gcd(A, B), K)
    if d > 1:
        A = A // d
        B = B // d
        K = K // d

    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    ans = 0
    ct = combmod_table(N)
    for x in range(N+1):
        y, m = divmod(K - A * x, B)
        if m == 0 and 0 <= y <= N:
            ans = ans + ct[min(x, N-x)] * ct[min(y, N-y)]
            ans = ans % mod
    return ans

if __name__ == '__main__':
    print(main())

Solution (Python):
N,A,B,K = map(int,input().split())
sprime = 998244353 
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
 

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)

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

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

num=1
for i in range(1,N+1):
    num = (num*i) %mod


number =0
for i in range(N+1):
    if (K-A*i)%B==0 and K-A*i>=0:
        j=(K-A*i)//B
        if 0<=j<=N:
            t = (num*num)%mod
            t = (t*table[i])%mod
            t = (t*table[j])%mod
            t = (t*table[N-i])%mod
            t = (t*table[N-j])%mod
            number = (number+t)%mod
print(number)

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

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

num=1
for i in range(1,N+1):
    num = (num*i) %mod


number =0
for i in range(N+1):
    if (K-A*i)%B==0 and K-A*i>=0:
        j=(K-A*i)//B
        if 0<=j<=N:
            t = (num*num)%mod
            t = (t*table[i])%mod
            t = (t*table[j])%mod
            t = (t*table[N-i])%mod
            t = (t*table[N-j])%mod
            number = (number+t)%mod
print(number)

Solution (Python):
n,a,b,k=[int(i) for i in input().split()]
c=a+b
mod=998244353
ans=0

g1 = [1, 1]  
g2 = [1, 1]  
inverse = [0, 1]  

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


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

for i in range(n+1):
    if (k-b*i)%a==0 and k-b*i>=0:
        n_a=(k-b*i)//a
        n_b=i
        if  max(n_a, n_b) <=n:
            ans+=com(n, n_a, mod) * com(n, n_b, mod)  %mod

print(ans%mod)

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

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

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

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

Solution (Python):
def main():
    n, a, b, k = map(int, raw_input().split())
    mod = 998244353
    f = [1] * 300010
    invf = [1] * 300010
    for i in xrange(300000):
        f[i+1] = f[i] * (i + 1) % mod
    invf[300000] = pow(f[300000], mod - 2, mod)
    for i in xrange(300000, 1, -1):
        invf[i-1] = invf[i] * i % mod
    ans = 0
    for i in xrange(n + 1):
        if k < a * i:
            break
        if (k - a * i) % b == 0:
            t = (k - a * i) / b
            if t > n:
                continue
            ans += f[n] * invf[i] * invf[n-i] % mod * f[n] * invf[t] * invf[n-t] % mod
            if ans >= mod:
                ans -= mod
    print ans
main()

Solution (Python):
def gcd(a,b):
	if a > b:
		n,m = a,b
	else:
		n,m = b,a
	while m > 0:
		t = n % m
		n = m
		m = t
	return n

def gcd_3(a,b,c):
	return gcd(gcd(a,b),c)

g1 = [1, 1] 
g2 = [1, 1] 

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

N,A,B,K = map(int,input().split())
g = gcd_3(A,B,K)
A = A // g
B = B // g
K = K // g
mod = 998244353 
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
n0,m0=-1,-1
for i in range(K//A + 1):
	if (K - A*i) % B == 0:
		n0,m0 = i, (K - A*i)//B
		break
k = 0
while k*B + n0 <= N and m0 - k*A >= 0:
	n,m = k*B + n0, m0 - k*A
	if m <= N:
		ans = (ans + cmb(N,n,mod) * cmb(N,m,mod)) % mod
	k += 1
print(ans)

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 = 998244353

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 inv(x):
    return pow(x, mod - 2, mod)

nm = {}
nm[0] = 1
def nCr_(n):
    b = n // 2 + 1
    r = 1
    d = 1
    for k in range(1,b+1):
        r = r * (n-k+1) % mod
        d = d * k % mod
        nm[k] = r * inv(d) % mod

def nCr(n, b):
    if b > n - b:
        b = n - b
    return nm[b]

def main():
    n,a,b,k = LI()
    r = 0
    nCr_(n)
    for i in range(n+1):
        ai = a * i
        bi = k - ai
        if bi < 0 or bi % b > 0 or bi // b > n:
            continue
        bb = bi // b
        r += nCr(n,i) * nCr(n,bb)
        r %= mod

    return r


print(main())

Solution (Python):


import sys,collections,itertools,re,math,fractions,decimal,random,array,bisect,heapq



MOD = 998244353

fact = [1] * 300001
ifact = [1] * 300001
for i in range(1, 300001):
    fact[i] = (fact[i-1] * i) % MOD
    ifact[i] = pow(fact[i], MOD-2, MOD)


def nCk(n, k):
    if k < 0 or k > n: return 0
    return (fact[n] * ifact[n-k] * ifact[k]) % MOD

def solve(f):
    n, a, b, k = f.read_int_list()
    ans = 0

    for ai in xrange(n+1):
        res = k - a*ai
        if res < 0 or res%b != 0: continue
        bi = res/b
        ans = (ans + nCk(n, ai)*nCk(n, bi))%MOD

    return ans


class Reader(object):
    def __init__(self, filename=None):
        self.file = open(filename) if filename is not None else None
        self.case = 1

    def __readline(self):
        return self.file.next().strip() if self.file else raw_input()

    def next_case(self):
        self.file.next()
        self.case += 1

    def read_int(self): return int(self.__readline())
    def read_float(self): return float(self.__readline())
    def read_long(self): return long(self.__readline())
    def read_decimal(self): return decimal.Decimal(self.__readline())
    def read_str(self): return self.__readline()

    def read_int_list(self): return map(int, self.__readline().split())
    def read_float_list(self): return map(float, self.__readline().split())
    def read_long_list(self): return map(long, self.__readline().split())
    def read_decimal_list(self): return map(decimal.Decimal, self.__readline().split())
    def read_str_list(self): return self.__readline().split()

if __name__ == '__main__':
    filename = sys.argv[1] if len(sys.argv) > 1 else None
    f = Reader(filename)
    if f.file:
        while True:
            print "Case 
            try:
                f.next_case()
            except StopIteration:
                break
    else:
        print solve(f)

Solution (Python):
MOD = 998244353

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

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 k in range(1,N//2+1):
  v = memo[k-1]*(N-k+1)*modinv(k,MOD)
  v %= MOD
  memo[k] = memo[N-k] = v

cnt = 0

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

print(cnt)

Solution (Python):
MOD = 998244353

f_list = [1] * 300000
f_r_list = [1] * 300000

for i in range(299999):
	f_list[i + 1] = int((f_list[i] * (i + 2)) % MOD)

def power(n, x):
	if x == 1:
		return n
	elif x % 2 == 0:
		return power(int((n * n) % MOD), int(x / 2))
	else:
		return int((n * power(n, x - 1)) % MOD)

f_r_list[-1] = power(f_list[-1], MOD - 2)

for i in range(2, 300001):
	f_r_list[-i] = int((f_r_list[-i + 1] * (300002 - i)) % MOD)

def comb(n, r):
	if n == 0 or r == 0 or n == r:
		return 1
	else:
		return (((f_list[n - 1] * f_r_list[n - r - 1]) % MOD) * f_r_list[r - 1]) % MOD 


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

cmb = []

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

if cmb == []:
	print(0)
else:

	ans = 0

	for c in cmb:
		ans = int((ans + (comb(n, c[0]) * comb(n, c[1])) % MOD) % MOD)

	print(ans)

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

Solution (Python):

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

LARGE = 998244353
INVS = [0]*(N+1)
for i in range(1, N+1):
    INVS[i] = pow(i, LARGE-2, LARGE)
NCKL = [1]*(N+1)
for i in range(1, N+1):
    NCKL[i] = (NCKL[i-1] * (N+1-i) * INVS[i]) % LARGE

res = 0
for i in range(N+1):
    if (K - A*i) % B == 0:
        j = (K - A*i) // B
        if j >= 0 and j <= N:
            res += (NCKL[i] * NCKL[j]) % LARGE
print(res % LARGE)

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

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

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