Problem:
 Score : 600 points    Problem Statement We have a grid with A horizontal rows and B vertical columns, with the squares painted white. On this grid, we will repeatedly apply the following operation:   Assume that the grid currently has a horizontal rows and b vertical columns. Choose "vertical" or "horizontal".  If we choose "vertical", insert one row at the top of the grid, resulting in an (a+1) \times b grid.  If we choose "horizontal", insert one column at the right end of the grid, resulting in an a \times (b+1) grid.    Then, paint one of the added squares black, and the other squares white.   Assume the grid eventually has C horizontal rows and D vertical columns. Find the number of ways in which the squares can be painted in the end, modulo 998244353 .      Constraints  1 \leq A \leq C \leq 3000  1 \leq B \leq D \leq 3000  A , B , C , and D are integers.         Input Input is given from Standard Input in the following format:  A  B  C  D       Output Print the number of ways in which the squares can be painted in the end, modulo 998244353 .        Sample Input 1 1 1 2 2      Sample Output 1 3  Any two of the three squares other than the bottom-left square can be painted black.       Sample Input 2 2 1 3 4      Sample Output 2 65       Sample Input 3 31 41 59 265      Sample Output 3 387222020

Solution (Python):
import sys
readline = sys.stdin.readline

MOD = 998244353
A, B, C, D = map(int, readline().split())

dp = [[0]*D for _ in range(C)]
dp[A-1][B-1] = 1
for h in range(A, C):
    dp[h][B-1] = dp[h-1][B-1]*B%MOD

for w in range(B, D):
    dp[A-1][w] = dp[A-1][w-1]*A%MOD

for s in range(A+B, C+D-1):
    for h in range(A, C):
        w = s - h
        if not B-1 < w < D:
            continue
        dp[h][w] = (dp[h][w-1] + dp[h-1][w] + dp[h-1][w]*w + dp[h][w-1]*h - dp[h-1][w-1]*h*w)%MOD

print(dp[C-1][D-1])

Solution (Python):
MOD = 998244353

A, B, C, D = map(int, input().split())

dp = [[0] * (D-B+1) for _ in range(C-A+1)]
dp[0][0] = 1
for i in range(1, C-A+1):
    dp[i][0] = dp[i-1][0] * B
    dp[i][0] %= MOD
for i in range(1, D-B+1):
    dp[0][i] = dp[0][i-1] * A
    dp[0][i] %= MOD

for a in range(A, C):
    for b in range(B, D):
        dp[a-A+1][b-B+1] = dp[a-A][b-B+1] * (b+1) + dp[a-A+1][b-B] * (a+1) - dp[a-A][b-B] * a * b
        dp[a-A+1][b-B+1] %= MOD

print(dp[C-A][D-B])

Solution (Python):
import sys
def input():return sys.stdin.readline()[:-1]
def N(): return int(input())
def NM():return map(int,input().split())
def L():return list(NM())
def LN(n):return [N() for i in range(n)]
def LL(n):return [L() for i in range(n)]

a,b,c,d=NM()
dp=[[0]*(d+1) for i in range(c+1)]
dp[a][b]=1
mod=998244353
for i in range(a+1,c+1):
    dp[i][b]=dp[i-1][b]*b%mod
for j in range(b+1,d+1):
    dp[a][j]=dp[a][j-1]*a%mod
for i in range(a+1,c+1):
    for j in range(b+1,d+1):
        dp[i][j]=(dp[i][j-1]*i+dp[i-1][j]*j-(i-1)*(j-1)*dp[i-1][j-1])%mod
print(dp[c][d])

Solution (Python):
A,B,C,D=map(int,input().split())
mod=998244353
dp=[[0 for y in range(D-B+1)]for x in range(C-A+1)]
dp[0][0]=1
for x in range(C-A):
    dp[x+1][0]=(dp[x][0]*B)%mod
for y in range(D-B):
    dp[0][y+1]=(dp[0][y]*A)%mod
for x in range(C-A):
    for y in range(D-B):
        dp[x+1][y+1]=(dp[x+1][y]*(A+x))%mod
        dp[x+1][y+1]+=(dp[x][y+1]*(B+y))%mod
        dp[x+1][y+1]%=mod
        dp[x+1][y+1]-=(dp[x][y]*(A+x)*(B+y))%mod
        dp[x+1][y+1]%=mod
        dp[x+1][y+1]+=(dp[x+1][y]+dp[x][y+1])%mod
        dp[x+1][y+1]%=mod

print(dp[C-A][D-B])

Solution (Python):
def main():
    A, B, C, D = list(map(int, input().split()))
    mod = 998244353
    u = [0] * (D + 1)
    v = [0] * (D + 1)
    u[B] = 1
    for b in range(B + 1, D + 1):
        u[b] = u[b - 1] * A
    for a in range(A + 1, C + 1):
        v[B] = u[B] * B % mod
        for b in range(B + 1, D + 1):
            v[b] = (u[b] * b + v[b - 1] * a - u[b - 1] * (a - 1) * (b - 1)) % mod
        u, v = v, u

    return u[D]

print(main())

Solution (Python):
def main():
    mod=998244353
    a,b,c,d=map(int,input().split())
    dx,dy=c-a,d-b
    if dx==0:
        return pow(a,dy,mod)
    if dy==0:
        return pow(b,dx,mod)
    DP=[[0]*(d+1) for _ in range(c+1)]
    DP[a][b]=1
    for x in range(a+1,c+1):
        DP[x][b]=pow(b,x-a,mod)
    for y in range(b+1,d+1):
        DP[a][y]=pow(a,y-b,mod)
    
    for x in range(a+1,c+1):
        for y in range(b+1,d+1):
            DP[x][y]=DP[x-1][y]*y+DP[x][y-1]*x
            DP[x][y]-=DP[x-1][y-1]*(x-1)*(y-1)
            DP[x][y]%=mod
    
    return DP[-1][-1]

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

Solution (Python):
ma = 3010
dp = [[0] * ma for _ in range(ma)]
mod = 998244353


a, b, c, d = map(int, input().split())
dp[a][b] = 1
for i in range(b+1, ma):
    dp[a][i] = dp[a][i - 1] * a % mod
for i in range(a+1, ma):
    dp[i][b] = dp[i - 1][b] * b % mod


for i in range(a+1, ma):
    for j in range(b+1, ma):
        dp[i][j] = (-dp[i - 1][j - 1] * (i-1)*(j-1) +
                    dp[i - 1][j] * j + dp[i][j - 1] * i) % mod
print(dp[c][d])

Solution (Python):

import sys
sys.setrecursionlimit(10**9)
INF=10**18
MOD=998244353
input=lambda: sys.stdin.readline().rstrip()
YesNo=lambda b: bool([print('Yes')] if b else print('No'))
YESNO=lambda b: bool([print('YES')] if b else print('NO'))
int1=lambda x:int(x)-1

a,b,c,d=map(int,input().split())
MAX=3050
dp=[[0]*MAX for _ in range(MAX)]
dp[a][b]=1
for i in range(a,c+1):
    for j in range(b,d+1):
        dp[i][j]+=dp[i][j-1]*i+dp[i-1][j]*j-dp[i-1][j-1]*(i-1)*(j-1)
        dp[i][j]%=MOD
print(dp[c][d])

Solution (Python):
MOD = 998244353
A,B,C,D = map(int,input().split())
dp = [[1 for i in range(D+1)] for j in range(C+1)]
for i in range(B+1, D+1):
	dp[A][i] = dp[A][i-1] * A % MOD
for j in range(A+1, C+1):
	dp[j][B] = dp[j-1][B] * B % MOD

for i in range(A+1, C+1):
	for j in range(B+1, D+1):
		dp[i][j] = (dp[i-1][j] * j + dp[i][j-1] * i - dp[i-1][j-1] * (i-1)*(j-1)) % MOD
print(dp[C][D])

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

dp = [[0] * (D+1) for _ in range(C+1)]

x, y = A, B
dp[x][y] = 1

for i in range(x, C+1):
    for j in range(y, D+1):
        if i == x and j == y:
            continue

        if i > x:
            dp[i][j] += dp[i-1][j] * j
        if j > y:
            dp[i][j] += dp[i][j-1] * i

        if i > x and j > y:
            dp[i][j] -= (i-1) * (j-1) * dp[i-1][j-1]

        dp[i][j] %= mod

print(dp[C][D])

Solution (Python):


import sys
read = sys.stdin.read
readline = sys.stdin.readline

a,b,c,d = map(int,readline().split())


dp1 = [[0]*(d+1) for _ in range(c+1)] 
dp2 = [[0]*(d+1) for _ in range(c+1)] 

dp1[a][b] = 1



MOD = 998244353
for i in range(a,c+1):
    for j in range(b,d+1):
        if i != c:
            dp1[i+1][j] += j*dp1[i][j]
            dp1[i+1][j] += dp2[i][j]
            dp1[i+1][j] %= MOD
        if j != d:
            dp2[i][j+1] += i*dp2[i][j]
            dp2[i][j+1] += i*dp1[i][j]
            dp2[i][j+1] %= MOD
            



print((dp1[c][d]+dp2[c][d])%MOD)

Solution (Python):
MOD = 998244353
a, b, c, d = map(int, input().split())
dp = [[0] * (d + 1) for _ in range(c + 1)]
dp[a][b] = 1
for i in range(a, c + 1):
    for j in range(b, d + 1):
        if i > a:
            dp[i][j] += dp[i - 1][j] * j
        if j > b:
            dp[i][j] += dp[i][j - 1] * i
        if i > a and j > b:
            dp[i][j] -= dp[i - 1][j - 1] * \
                (i - 1) * (j - 1)
        dp[i][j] %= MOD
print(dp[c][d])

Solution (Python):
MOD = 998244353

def fastExponentiation(a,b,N): 
    ans = 1
    bDiv2 = b
    aPow2 = a
    while bDiv2 > 0:
        if bDiv2 % 2 != 0:
            ans *= aPow2
            ans %= N
        bDiv2 //= 2
        aPow2 = aPow2 * aPow2
        aPow2 %= N
    return ans

a,b,c,d = map(int,input().split())
DP = []
for i in range(c-a+1):
    DP.append([])
    for j in range(d-b+1):
        DP[i].append(0)

for i in range(c-a+1):
    partialSum = 0
    if i == 0:
        DP[i][0] = 1
        continue
    for j in range(d-b+1):
        if j == 0:
            DP[i][j] = (DP[i-1][j] * b) % MOD
        else:
            partialSum += DP[i-1][j-1]
            partialSum *= (a + i - 1)
            partialSum %= MOD
            DP[i][j] = (DP[i-1][j] * (j + b) + partialSum) % MOD

ans = 0
for j in range(d-b+1):
    ans = ans + DP[c-a][j] * fastExponentiation(c,d-b-j,MOD)
    ans %= MOD
print(ans)

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

A,B,C,D=map(int,input().split())
mod=998244353

DP=[[0]*(D+1) for i in range(C+1)]
DP[A][B]=1

for i in range(A,C+1):
    for j in range(B,D+1):
        if i==A and j==B:
            continue
        DP[i][j]=(DP[i-1][j]*j+DP[i][j-1]*i-DP[i-1][j-1]*(i-1)*(j-1))%mod

print(DP[C][D]%mod)

Solution (Python):
def main():
    mod = 998244353
    a, b, c, d = map(int, input().split())
    dp = [[0]*(d+1) for _ in [0]*(c+1)]
    dp[a][b] = 1
    for i in range(a+b+1, c+d+1):
        for p in range(c, a-1, -1):
            q = i-p
            if q > d:
                break
            if q < 0:
                continue
            dp[p][q] = (q*dp[p-1][q]+p*dp[p][q-1] -
                        (p-1)*(q-1)*dp[p-1][q-1]) % mod
    print(dp[c][d])


main()

Solution (Python):
MOD = 998244353


def solve(A, B, C, D):
    numWays = [[None for c in range(D + 1)] for r in range(C + 1)]

    numWays[A][B] = 1

    for rows in range(A, C + 1):
        for cols in range(B, D + 1):
            if rows == A and cols == B:
                continue
            ways = 0
            if rows > A:
                ways += cols * numWays[rows - 1][cols] % MOD
            if cols > B:
                ways += rows * numWays[rows][cols - 1] % MOD
            if rows > A and cols > B:
                ways -= (cols - 1) * (rows - 1) * numWays[rows - 1][cols - 1] % MOD
            numWays[rows][cols] = ways % MOD

    return numWays[C][D] % MOD


A, B, C, D = [int(x) for x in input().split()]
print(solve(A, B, C, D))

Solution (Python):
import sys
sys.setrecursionlimit(10**7) 
import math
from copy import copy, deepcopy
from copy import deepcopy as dcp
from operator import itemgetter
from bisect import bisect_left, bisect, bisect_right

from collections import deque


from collections import Counter

from itertools import accumulate,combinations,permutations

from heapq import heapify,heappop,heappush


from functools import lru_cache

from decimal import Decimal

def input(): 
    x=sys.stdin.readline()
    return x[:-1] if x[-1]=="\n" else x
def printl(li): _=print(*li, sep="\n") if li else None
def argsort(s, return_sorted=False): 
    inds=sorted(range(len(s)), key=lambda k: s[k])
    if return_sorted: return inds, [s[i] for i in inds]
    return inds
def alp2num(c,cap=False): return ord(c)-97 if not cap else ord(c)-65
def num2alp(i,cap=False): return chr(i+97) if not cap else chr(i+65)
def matmat(A,B):
    K,N,M=len(B),len(A),len(B[0])
    return [[sum([(A[i][k]*B[k][j]) for k in range(K)]) for j in range(M)] for i in range(N)]
def matvec(M,v):
    N,size=len(v),len(M)
    return [sum([M[i][j]*v[j] for j in range(N)]) for i in range(size)]
def T(M):
    n,m=len(M),len(M[0])
    return [[M[j][i] for j in range(n)] for i in range(m)]
def main():
    mod = 998244353
    

    
    A,B,C,D = map(int, input().split())
    
    
    
    dp=[[0]*(D+1) for _ in range(C+1)]
    dp[A][B]=1
    for i in range(A,C+1):
        for j in range(B,D+1):
            if i==A and j==B:continue
            
            dp[i][j]=((dp[i-1][j]*j+dp[i][j-1]*i)-dp[i-1][j-1]*(i-1)*(j-1))%mod
    print(dp[-1][-1])





if __name__ == "__main__":
    main()

Solution (Python):
A, B, C, D = map(int, input().split())
MOD = 998244353
dp = [[0]*(D+1) for i in range(C+1)]
for i in range(A, C+1):
    for j in range(B, D+1):
        if (i, j) == (A, B):
            dp[i][j] = 1
        else:
            dp[i][j] = (dp[i-1][j]*j + dp[i][j-1]*i - dp[i-1][j-1]*(i-1)*(j-1)) % MOD

print(dp[-1][-1])

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

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



A,B,C,D = LI()

dp = [[0]*(D+1) for _ in range(C+1)]
dp[A][B] = 1

for i in range(A,C+1):
    for j in range(B,D+1):
        if i == A and j == B:
            continue
        dp[i][j] = (dp[i-1][j]*j)%mod + (dp[i][j-1]*i)%mod\
                     - (dp[i-1][j-1]*(i-1)*(j-1))%mod
        dp[i][j] %= mod

print(dp[C][D])

Solution (Python):
a,b,c,d = map(int,input().split())

mod = 998244353

dp = [[0] * d for _ in range(c)]
dp[a-1][b-1] = 1

for i in range(c):
    for j in range(d):
        if i < a and j < b or i < a-1 or j < b-1:
            pass
        else:
            if i == 0:
                dp[i][j] = 1
            elif j == 0:
                dp[i][j] = 1
            else:
                dp[i][j] = dp[i-1][j] * (j+1) + dp[i][j-1] * (i+1) - dp[i-1][j-1] * i*j
        dp[i][j] %= mod

print(dp[-1][-1])

Solution (Python):
import sys
import numpy as np
from functools import lru_cache

read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines

MOD = 998_244_353

def main(A, B, C, D):
    dp = np.zeros((C + 1, D + 1), np.int64)
    dp[C, D] = 1
    for d in range(D, 0, -1):
        dp[C, d - 1] = dp[C, d] * C % MOD
    for c in range(C - 1, -1, -1):
        x = 0
        for d in range(D, -1, -1):
            dp[c, d] = (d * dp[c + 1, d] + x) % MOD
            x += dp[c + 1, d]
            x *= c
            x %= MOD
    return dp[A, B]

if sys.argv[-1] == 'ONLINE_JUDGE':
    import numba
    from numba.pycc import CC
    i8 = numba.from_dtype(np.int64)
    signature = (i8, i8, i8, i8)

    cc = CC('my_module')
    cc.export('main', signature)(main)
    cc.compile()

from my_module import main

A, B, C, D = map(int, read().split())
print(main(A, B, C, D))

Solution (Python):
"""

ありうるものは？
差の分だけ塗れる

ｄｐ？
dp[i][j] = dp[i-1][j] * j + dp[i][j-1] * i
i,jの中に 丁度i-A,j-B個ある置き方は全部できる

dp[i][j] = dp[i-1][j] * j + dp[i][j-1] * i - dp[i-1][j-1] * (i-1) * (j-1)

i,j置いた時に
i+j個の黒が入っているはず
→ i-1,jにi+j-1こはいってる * j + i-1,jにi+j-1こはいってる* i - どっちも満たしてる

"""

A,B,C,D = map(int,input().split())

dp = [[0] * (D+1) for i in range(C+1)]

dp[A][B] = 1
mod = 998244353

for i in range(A,C+1):

    for j in range(B,D+1):

        if i == A and j == B:
            continue

        dp[i][j] = dp[i-1][j] * j + dp[i][j-1] * i - dp[i-1][j-1] * (i-1) * (j-1)
        dp[i][j] %= mod

print (dp[C][D]%mod)

Solution (Python):
import sys
readline = sys.stdin.readline

P = 998244353
def mod_prod(a,b):
    return (a % P) * (b % P) % P

def mod_fact(a):
    ret = 1
    for i in range(a):
        ret = ret * (i + 1) % P
    return ret

def mod_pow(a, n):
    ret = 1
    while n > 0:
        if n % 2 == 1:
            ret = ret * a % P
        a = a * a % P
        n = n >> 1
    return ret

def mod_div(a, b):
    b_inv = mod_pow(b, P - 2)
    return a * b_inv % P

def mod_comb(a,b):
    ret = 1
    b = min(b, a - b)
    for i in range(b):
        ret *= a - i
        ret %= P
    ret = mod_div(ret, mod_fact(b))
    return ret

def solve():
    A, B, C, D = map(int, readline().split())

    dp = [[0] * D for i in range(C)]

    dp[A - 1][B - 1] = 1

    for i in range(C):
        for j in range(D):
            if i >= A or j >= B:
                dp[i][j] = (
                    (dp[i - 1][j] * (j + 1) if i > 0 else 0)
                    + (dp[i][j - 1] * (i + 1) if j > 0 else 0)
                    - (dp[i - 1][j - 1] * i * j if i > 0 and j > 0 else 0)
                   ) % P
    print(dp[-1][-1])

solve()

Solution (Python):
import sys
input = sys.stdin.readline
input = sys.stdin.buffer.readline





def RD(): return sys.stdin.read()
def II(): return int(input())
def MI(): return map(int,input().split())
def MF(): return map(float,input().split())
def LI(): return list(map(int,input().split()))
def LF(): return list(map(float,input().split()))
def TI(): return tuple(map(int,input().split()))





def main():
	a,b,c,d=MI()
	dp=[[0]*(d+1) for _ in range(c+1)]
	dp[a][b]=1

	for i in range(a-1,c):
		for j in range(b-1,d):
			if i+1==a and j+1==b:
				continue
			
			dp[i+1][j+1]=dp[i][j+1]*(j+1)+dp[i+1][j]*(i+1)-i*j*dp[i][j]
			dp[i+1][j+1]%=998244353

	print(dp[c][d])
















if __name__ == "__main__":
	main()

Solution (Python):
A,B,C,D = map(int,input().split())
MOD = 998244353
dp = [[0]*D for i in range(C)]
dp[A-1][B-1] = 1

for i in range(C):
    for j in range(D):
        if i!=A-1 or j!=B-1:
            val1 = dp[i-1][j]*(j+1)
            if i==0:
                val1 = 0
            val2 = dp[i][j-1]*(i+1)
            if j==0:
                val2 = 0
            val3 = 0
            if val1*val2!=0:
                val3 = i*j
            dp[i][j] = (val1+val2 - val3*dp[i-1][j-1])%MOD
            

print(dp[C-1][D-1])

Solution (Python):
import sys

readline = sys.stdin.readline
readall = sys.stdin.read
ns = lambda: readline().rstrip()
ni = lambda: int(readline().rstrip())
nm = lambda: map(int, readline().split())
nl = lambda: list(map(int, readline().split()))
prn = lambda x: print(*x, sep='\n')

def solve():
    mod = 998244353
    a, b, c, d = nm()
    n, m = c-a, d-b
    dp = [[0]*(m+1) for _ in range(n+1)]
    dp[0][0] = 1
    for i in range(n):
        dp[i+1][0] = dp[i][0] * b % mod
    for i in range(m):
        dp[0][i+1] = dp[0][i] * a % mod
    for i in range(1, n+1):
        for j in range(1, m+1):
            dp[i][j] = (dp[i][j-1] * (a + i) + dp[i-1][j] * (b + j) - dp[i-1][j-1] * (a + i - 1) * (b + j - 1)) % mod
    
    
    print(dp[n][m])

    return

solve()

Solution (Python):
import sys
sys.setrecursionlimit(10000000)
MOD = 998244353
INF = 10 ** 9

def main():
    A,B,C,D = map(int,input().split())
    dp = [[0]*(1+D) for _ in range(C + 1)]

    dp[A][B] = 1
    for i in range(A - 1,C):
        for j in range(B - 1,D):
            if i == A - 1 and j == B - 1:
                continue
            dp[i + 1][j + 1] = (i + 1)*dp[i + 1][j] + (j + 1)*dp[i][j + 1] - dp[i][j] * i*j
            dp[i + 1][j + 1] %= MOD
    print(dp[C][D])
if __name__ == '__main__':
    main()

Solution (Python):
a, b, c, d = map(int, input().split())
MOD = 998244353


dp = [[0] * (d + 1) for i in range(c + 1)]
dp[a][b] = 1

for i in range(a, c):
    dp[i + 1][b] = b * dp[i][b]
    dp[i + 1][b] %= MOD
for j in range(b, d):
    dp[a][j + 1] = a * dp[a][j]
    dp[a][j + 1] %= MOD
for i in range(a, c):
    for j in range(b, d):
        dp[i + 1][j + 1] = dp[i][j + 1] * (j + 1) + dp[i + 1][j] * (i + 1) - dp[i][j] * (i * j)
        dp[i + 1][j + 1] %= MOD
print(dp[-1][-1])

Solution (Python):
a, b, c, d = map(int,input().split())
mod = 998244353
A = [[0] * (d+1) for i in range(c+1)]

A[a][b] = 1

for i in range(a, c+1):
    for j in range(b, d+1):
        A[i][j] = (A[i][j] + A[i][j-1] * i + A[i-1][j] * j ) % mod
        if i!=a and j!=b:
            A[i][j] = (A[i][j] - A[i-1][j-1] * (i-1)*(j-1)) % mod
print(A[c][d])

Solution (Python):
P = 998244353
A, B, C, D = map(int, input().split())
X = [[0] * 3001 for _ in range(3001)]
for a in range(3001):
    for b in range(3001):
        if a < A or b < B: continue
        if a == A and b == B:
            X[a][b] = 1
        elif a == A:
            X[a][b] = X[a][b-1] * a % P
        elif b == B:
            X[a][b] = X[a-1][b] * b % P
        else:
            X[a][b] = (X[a-1][b] * b + X[a][b-1] * a - X[a-1][b-1] * (a-1) * (b-1)) % P

print(X[C][D])

Solution (Python):
a, b, c, d = map(int,input().split())
dp = [[0 for j in range(d+1)] for i in range(c+1)]
for i in range(a, c+1):
  for j in range(b, d+1):
    if i == a and j == b:
      dp[i][j] = 1
      continue
    dp[i][j] = (dp[i-1][j] * j + dp[i][j-1] * i - dp[i-1][j-1] * (i-1) * (j-1) ) % 998244353
    dp[i][j] += 998244353
    dp[i][j] %= 998244353
print(dp[c][d])

Solution (Python):
A,B,C,D=map(int,input().split())
dp=[[0]*(D+1) for i in range(C+1)]
m=998244353
dp[A][B]=1
for i in range(C-A):
    dp[A+i+1][B]=dp[A+i][B]*B%m
for j in range(D-B):
    dp[A][B+j+1]=dp[A][B+j]*A%m
for i in range(C-A):
    for j in range(D-B):
        dp[A+i+1][B+j+1]=(dp[A+i][B+j+1]*(B+j+1)+dp[A+i+1][B+j]*(A+i+1)-dp[A+i][B+j]*(A+i)*(B+j))%m
print(dp[-1][-1])

Solution (Python):
MOD = 998244353
A, B, C, D = map(int, input().split())
dp = [[0 for _ in range(D+1)] for _ in range(C+1)]
dp[A][B] = 1
for i in range(A, C+1):
	for j in range(B, D+1):
		if i > A:
			dp[i][j] += j * dp[i-1][j]
		if j > B:
			dp[i][j] += i * dp[i][j-1]
		if i>A and j>B:
			dp[i][j] -= (i-1) * (j-1) * dp[i-1][j-1]
		dp[i][j] %= MOD
print(dp[C][D])

Solution (Python):
a,b,c,d = map(int,input().split())
mod = 998244353
dp = [[0 for i in range(d+1)] for j in range(c+1)]
dp[a][b] = 1
for i in range(a,c+1):
  for j in range(b,d+1):
    if i == a and j == b:
      continue
    if i > a:
      dp[i][j] += dp[i-1][j]*j
    if j > b:
      dp[i][j] += dp[i][j-1]*i
    if i != a and j != b:
      dp[i][j] -= dp[i-1][j-1]*(i-1)*(j-1)
    dp[i][j] %= mod
print(dp[c][d])