CodeEff/ECCO huggingface dataset:
# Original solution
MOD = 998244353
HALF = (MOD + 1) // 2


def solve(a, b, c, d):
    dp1 = [[0] * (d + 1) for _ in range(c + 1)]
    dp2 = [[0] * (d + 1) for _ in range(c + 1)]
    dp3 = [[0] * (d + 1) for _ in range(c + 1)]
    dp1[a][b] = 1
    dp2[a][b] = 0
    dp3[a][b] = 0
    for i in range(a + 1, c + 1):
        dp1[i][b] = (pow(b - 1, i - a, MOD) + pow(b - 1, i - a - 1, MOD) * 1 * (i - a)) % MOD
        dp3[i][b] = (pow(b, i - a, MOD) - dp1[i][b]) % MOD
    for j in range(b + 1, d + 1):
        dp1[a][j] = (pow(a - 1, j - b, MOD) + pow(a - 1, j - b - 1, MOD) * 1 * (j - b)) % MOD
        dp2[a][j] = (pow(a, j - b, MOD) - dp1[a][j]) % MOD
    for i in range(a + 1, c + 1):
        for j in range(b + 1, d + 1):
            dp1[i][j] = (((dp1[i - 1][j] + dp2[i - 1][j]) * (j - 1) + (dp1[i][j - 1] + dp3[i][j - 1]) * (i - 1)) * HALF) % MOD
            dp2[i][j] = (dp2[i][j - 1] * i + dp1[i][j - 1] + dp3[i][j - 1]) % MOD
            dp3[i][j] = (dp3[i - 1][j] * j + dp1[i - 1][j] + dp2[i - 1][j]) % MOD

    # print(dp1, dp2, dp3)
    # print((dp1[c][d] + dp2[c][d] + dp3[c][d]) % MOD)
    return (dp1[c][d] + dp2[c][d] + dp3[c][d]) % MOD


def main():
    a, b, c, d = list(map(int, input().split()))
    res = solve(a, b, c, d)
    print(res)


def test():
    assert solve(1, 1, 2, 2) == 3
    assert solve(1, 1, 3, 3) == 30
    assert solve(2, 1, 3, 4) == 65
    # assert solve(31, 41, 59, 265) == 387222020


if __name__ == "__main__":
    test()
    main()


# Optimized solution
MOD = 998244353


def solve(a, b, c, d):
    dp1 = [[0] * (d + 1) for _ in range(c + 1)]
    dp2 = [[0] * (d + 1) for _ in range(c + 1)]
    dp3 = [[0] * (d + 1) for _ in range(c + 1)]
    dp1[a][b] = 1
    dp2[a][b] = 0
    dp3[a][b] = 0
    for i in range(a + 1, c + 1):
        dp1[i][b] = (pow(b - 1, i - a, MOD) + pow(b - 1, i - a - 1, MOD) * 1 * (i - a)) % MOD
        dp3[i][b] = (pow(b, i - a, MOD) - dp1[i][b]) % MOD
    for j in range(b + 1, d + 1):
        dp1[a][j] = (pow(a - 1, j - b, MOD) + pow(a - 1, j - b - 1, MOD) * 1 * (j - b)) % MOD
        dp2[a][j] = (pow(a, j - b, MOD) - dp1[a][j]) % MOD
    for i in range(a + 1, c + 1):
        for j in range(b + 1, d + 1):
            dp1[i][j] = ((dp1[i - 1][j] + dp2[i - 1][j]) * (j - 1)) % MOD
            dp2[i][j] = (dp2[i][j - 1] * i + dp1[i][j - 1] + dp3[i][j - 1]) % MOD
            dp3[i][j] = (dp3[i - 1][j] * j + dp1[i - 1][j] + dp2[i - 1][j]) % MOD

    # print(dp1, dp2, dp3)
    # print((dp1[c][d] + dp2[c][d] + dp3[c][d]) % MOD)
    return (dp1[c][d] + dp2[c][d] + dp3[c][d]) % MOD


def main():
    a, b, c, d = list(map(int, input().split()))
    res = solve(a, b, c, d)
    print(res)


def test():
    assert solve(1, 1, 2, 2) == 3
    assert solve(1, 1, 3, 3) == 30
    assert solve(2, 1, 3, 4) == 65
    # assert solve(31, 41, 59, 265) == 387222020


if __name__ == "__main__":
    test()
    main()
# Original solution
a,b,c,d = input().strip().split(' ')
a = int(a)
b = int(b)
c = int(c)
d = int(d)
ans = []
for i in range(a, c + 1):
  ans.append([])
  for j in range(b, d + 1):
    if a == i and j == b:
      ans[i - a].append(1)
    elif a == i:
      ans[i - a].append(i * ans[i - a][j - b - 1])
    elif b == j:
      ans[i - a].append(j * ans[i - a - 1][j - b])
    else:
      ans[i - a].append(i * ans[i - a][j - b - 1] + j * ans[i - a - 1][j - b] - (i - 1) * (j - 1) * ans[i - a - 1][j - b - 1])
print((ans[c - a][d - b] % 998244353))

# Optimized solution
a,b,c,d = input().strip().split(' ')
a = int(a)
b = int(b)
c = int(c)
d = int(d)
ans = [[], []]
for i in range(a, c + 1):
  i1 = 0 if i % 2 == 0 else 1
  i2 = 0 if i1 == 1 else 1
  ans[i1] = []
  for j in range(b, d + 1):
    if a == i and j == b:
      ans[i1].append(1)
    elif a == i:
      ans[i1].append(i * ans[i1][j - b - 1])
    elif b == j:
      ans[i1].append(j * ans[i2][j - b])
    else:
      ans[i1].append(i * ans[i1][j - b - 1] + j * ans[i2][j - b] - (i - 1) * (j - 1) * ans[i2][j - b - 1])
t = 0 if c % 2 == 0 else 1
print((ans[t][d - b] % 998244353))
# Original solution
import sys
input = sys.stdin.readline

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

dp = [[0] * (D+1) for _ in range(C+1)]
dp[A][B] = 1
for i in range(A+1, C+1):
    dp[i][B] = dp[i-1][B] * B % mod
for j in range(B, D):
    for i in range(A, C+1):
        x = dp[i][j] * i
        dp[i][j+1] = (dp[i][j+1] + x) % mod
        for k in range(i + 1, C+1):
            dp[k][j+1] = (dp[k][j+1] + x) % mod
            x = x * (j + 1) % mod
print((dp[C][D]))



# Optimized solution
import sys
input = sys.stdin.readline

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

dp = [[0] * (D+1) for _ in range(C+1)]
dp[A][B] = 1
for i in range(A, C+1):
    if i > A:
        dp[i][B] = dp[i-1][B] * B % mod
    for j in range(B+1, 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]))
# Original solution
import sys
sys.setrecursionlimit(10**7)

MOD = 998244353
A,B,C,D= list(map(int,input().split()))
#A,B,C,D = 1,1,100,100

def dpinit(ps, val=0):
    res = [val for i in [0]*ps[-1]]
    for i in ps[:-1][::-1]:
        res = [res[:] for k in [0]*i]
    return res

memo = dpinit((C+1,D+1))
memo[A][B] = 1
for Y in range(A,C+1):
    memo[Y][1] = 1
for X in range(B,D+1):
    memo[1][X] = 1

def cnt(C,D):
    if memo[C][D] != 0:
        return memo[C][D]
    res = 0
    if C > A:
        res += (D * cnt(C-1,D) % MOD)
    if D > B:
        res += (C * cnt(C,D-1) % MOD)
    if C > A and D > B:
        res -= (cnt(C-1,D-1)*(C-1)*(D-1) % MOD)
    res %= MOD
    memo[C][D] = res 
    return res

print((cnt(C,D)))






# Optimized solution
import sys
sys.setrecursionlimit(10**7)

MOD = 998244353
A,B,C,D= list(map(int,input().split()))
#A,B,C,D = 1,1,3000,3000

def dpinit(ps, val=0):
    res = [val for i in [0]*ps[-1]]
    for i in ps[:-1][::-1]:
        res = [res[:] for k in [0]*i]
    return res

memo = dpinit((C+1,D+1))

memo[A][B] = 1
for Y in range(A,C+1):
    memo[Y][B] = pow(B,Y-A,MOD)
for X in range(B,D+1):
    memo[A][X] = pow(A,X-B,MOD)

for y in range(A+1,C+1):
    for x in range(B+1,D+1):
        memo[y][x] = (x*memo[y-1][x] + y*memo[y][x-1] - (x-1)*(y-1)*memo[y-1][x-1]) % MOD

print((memo[C][D]))
# Original solution
import sys

sys.setrecursionlimit(10 ** 4)

MOD = 998244353


def solve(A, B, C, D):
    cache = {}

    def numWays(rows, cols):
        if rows == A and cols == B:
            return 1

        key = rows * 4096 + cols
        if key not in cache:

            assert rows >= A and cols >= B
            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
            cache[key] = ways % MOD
        return cache[key]

    return numWays(C, D) % MOD


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


# Optimized solution
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)))
# Original solution
A,B,C,D = list(map(int,input().split()))
MOD = 998244353

dp = [[0 for _ in range(D+1)] for _ in range(C+1)]
dp[A][B] = 1
#dp[i][j] ヨコi番目でそこまでの最大の高さがj

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

ans = dp[C][D]
print(ans)
      

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

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 or j > B:
      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

ans = dp[C][D]
print(ans)
# Original solution
#import numpy as np
import math
import collections
import bisect


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


dp = [[[1 for k in range(2)] for i in range(d + 1)] for j in range(c + 1)]
for k in range(2):
    for i in range(c + 1):
        for j in range(d + 1):
            dp[i][j][1] = 0
            if j < b or i < a:
                dp[i][j][k] = 0

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

print(((dp[c][d][0] + dp[c][d][1]) % 998244353))


# Optimized solution
#import numpy as np
import math
import collections
import bisect


def main():
    a, b, c, d = list(map(int, input().split()))

    dp0 = [[1 for i in range(d + 1)] for j in range(c + 1)]
    dp1 = [[0 for i in range(d + 1)] for j in range(c + 1)]
    for i in range(c + 1):
        for j in range(d + 1):
            if j < b or i < a:
                dp0[i][j] = 0
                dp1[i][j] = 0
    for i in range(a, c + 1):
        for j in range(b, d + 1):
            # print(i, j, dp[i][j])
            if a == i and b == j:
                continue
            dp0[i][j] = ((i * dp0[i][j - 1]) % 998244353 + dp1[i][j - 1]) % 998244353
            dp1[i][j] = (j * (dp0[i - 1][j] + dp1[i - 1][j]) % 998244353) % 998244353

    print(((dp0[c][d] + dp1[c][d]) % 998244353))


if __name__ == '__main__':
    main()
# Original solution
a,b,c,d=list(map(int,input().split()))
m=998244353
D=[[0]*(d+1)for _ in range(c+1)]
D[a][b]=1
for i in range(a,c+1):
 for j in range(b,d+1):
  D[i][j]+=(D[i][j-1]*i+D[i-1][j]*j-D[i-1][j-1]*(i-1)*(j-1))%m
print((D[c][d]))

# Optimized solution
a,b,c,d=list(map(int,input().split()));m=998244353;D=[[0]*(d+1)for _ in range(c+1)];D[a][b]=1
for i in range(a,c+1):
 for j in range(b,d+1):
  x=D[i-1]
  D[i][j]+=(D[i][j-1]*i+x[j]*j-x[j-1]*(i-1)*(j-1))%m
print((D[c][d]))
# Original solution
A,B,C,D=list(map(int,input().split()))
mod=998244353
dp=[[[0 for k in range(2)] for i in range(D+1)] for i in range(C+1)]
dp[C][D][0]=1
dp[C][D][1]=1
for i in range(D-1,B-1,-1):
    dp[C][i][0]=(dp[C][i+1][1]*C)%mod
    dp[C][i][1]=(dp[C][i+1][1]*C)%mod
for j in range(C-1,A-1,-1):
    dp[j][D][0]=(dp[j+1][D][0]*D)%mod
    dp[j][D][1]=(dp[j+1][D][0])%mod

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

print((dp[A][B][0]))


# Optimized solution
A,B,C,D=list(map(int,input().split()))
mod=998244353
dpx=[[0 for i in range(D+1)] for j in range(C+1)]
dpy=[[0 for i in range(D+1)] for j in range(C+1)]
dpx[C][D]=1
dpy[C][D]=1
for i in range(D-1,B-1,-1):
    dpx[C][i]=(dpy[C][i+1]*C)%mod
    dpy[C][i]=(dpy[C][i+1]*C)%mod
for j in range(C-1,A-1,-1):
    dpx[j][D]=(dpx[j+1][D]*D)%mod
    dpy[j][D]=(dpx[j+1][D])%mod

for i in range(C-1,A-1,-1):
    for j in range(D-1,B-1,-1):
        dpx[i][j]=((dpx[i+1][j]*j)%mod+(dpy[i][j+1]*i)%mod)%mod
        dpy[i][j]=(dpx[i+1][j]+(dpy[i][j+1]*i)%mod)%mod

print((dpx[A][B]))
# Original solution
A,B,C,D=list(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]))


# Optimized solution
A,B,C,D=list(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+1)+dp[x][y+1]*(B+y+1)-dp[x][y]*(A+x)*(B+y))%mod
print((dp[C-A][D-B]))
# Original solution
# coding: utf-8
import sys

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

"""
dp、無理な組み合わせを引く、追加された黒マスの数は決まっている
各マスのdp
下から来たものと左から来たものをdpで保持
下から来たものは次だけ制約あり
dp = [下から、左から]
"""
MOD = 998244353
A, B, C, D = lr()  # Cは縦
dp = [[[0, 0] for _ in range(D+1)] for _ in range(C+1)]   # 1-indexed
dp[A][B] = [0, 1]
for i in range(A, C+1):
    for j in range(B, D+1):
        # iが縦
        if i == A and j == B:
            continue
        if i > A:
            dp[i][j][0] = sum(dp[i-1][j]) * j
        if j > B:
            dp[i][j][1] = dp[i][j-1][0] + dp[i][j-1][1] * i
        dp[i][j][0] %= MOD
        dp[i][j][1] %= MOD

answer = sum(dp[C][D])
print((answer%MOD))



# Optimized solution
# coding: utf-8
import sys

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

"""
各マスのdp
下から来たものと左から来たものを足す
その後、重なっているところを引く
"""
MOD = 998244353
A, B, C, D = lr()  # Cは縦
dp = [[0] * (D+1) for _ in range(C+1)]   # 1-indexed
dp[A][B] = 1
for i in range(A, C+1):
    for j in range(B, D+1):
        # iが縦
        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

answer = dp[C][D]
print((answer%MOD))
# Original solution
# coding: utf-8
import sys

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

"""
各マスのdp
下から来たものと左から来たものを足す
その後、重なっているところを引く
"""
MOD = 998244353
A, B, C, D = lr()  # Cは縦
dp = [[0] * (D+1) for _ in range(C+1)]   # 1-indexed
dp[A][B] = 1
for i in range(A, C+1):
    for j in range(B, D+1):
        # iが縦
        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

answer = dp[C][D]
print((answer%MOD))


# Optimized solution
# coding: utf-8
import sys

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

"""
各マスのdp
下から来たものと左から来たものを足す
その後、重なっているところを引く
"""
MOD = 998244353
A, B, C, D = lr()  # Cは縦
dp = [[0] * (D+1) for _ in range(C+1)]   # 1-indexed
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][j] -= dp[i-1][j-1] * (i-1) * (j-1)
        dp[i][j] %= MOD

answer = dp[C][D]
print((answer%MOD))
# Original solution
import sys
from collections import defaultdict

def solve():
    input = sys.stdin.readline
    A, B, C, D = list(map(int, input().split()))
    mod = 998244353
    DP = defaultdict(int)
    DP[(A, B)] = 1
    for cb in range(B + 1, D + 1):
        DP[(A, cb)] = (DP[(A, cb - 1)] * A) % mod

    for ca in range(A + 1, C + 1):
        DP[(ca, B)] = (DP[(ca - 1, B)] * B) % mod
        for nb in range(B + 1, D + 1):
            DP[(ca, nb)] = (DP[(ca, nb - 1)] * ca + DP[(ca - 1, nb)] * nb - DP[(ca - 1, nb - 1)] * (ca - 1) * (nb - 1)) % mod

    print((DP[(C, D)]))
    #print(DP)

    return 0

if __name__ == "__main__":
    solve()

# Optimized solution
import sys
from collections import defaultdict

def solve():
    input = sys.stdin.readline
    A, B, C, D = list(map(int, input().split()))
    mod = 998244353

    DP = [[0 for j in range(D + 1)] for i in range(C + 1)]
    DP[A][B] = 1
    for b in range(B + 1, D + 1):
        DP[A][b] = (DP[A][b-1] * A) % mod
    
    for a in range(A + 1, C + 1):
        DP[a][B] = (DP[a-1][B] * B) % mod
        for b in range(B + 1, D + 1):
            DP[a][b] = ((DP[a][b-1] * a) % mod + (DP[a-1][b] * b) % mod - (DP[a-1][b-1] * (a-1) * (b-1)) % mod) % mod
    print((DP[C][D]))
    #print(DP)

    return 0

if __name__ == "__main__":
    solve()
# Original solution
import math,string,itertools,fractions,heapq,collections,re,array,bisect,sys,copy,functools
import time,random,resource

sys.setrecursionlimit(10**7)
inf = 10**20
eps = 1.0 / 10**10
mod = 10**9+7
mod2 = 998244353
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 list(map(int, sys.stdin.readline().split()))
def LLI(): return [list(map(int, l.split())) for l in sys.stdin.readlines()]
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 pe(s): return print(str(s), file=sys.stderr)
def JA(a, sep): return sep.join(map(str, a))
def JAA(a, s, t): return s.join(t.join(map(str, b)) for b in a)


def main():
    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 a < i:
                dp[i][j] += dp[i-1][j] * j % mod2
            if b < j:
                dp[i][j] += dp[i][j-1] * i % mod2
            if a < i and b < j:
                dp[i][j] -= dp[i-1][j-1] * (i-1) * (j-1)

            dp[i][j] %= mod2

    return dp[c][d]


print(main())




# Optimized solution
import math,string,itertools,fractions,heapq,collections,re,array,bisect,sys,copy,functools
import time,random,resource

sys.setrecursionlimit(10**7)
inf = 10**20
eps = 1.0 / 10**10
mod = 10**9+7
mod2 = 998244353
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 list(map(int, sys.stdin.readline().split()))
def LLI(): return [list(map(int, l.split())) for l in sys.stdin.readlines()]
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 pe(s): return print(str(s), file=sys.stderr)
def JA(a, sep): return sep.join(map(str, a))
def JAA(a, s, t): return s.join(t.join(map(str, b)) for b in a)


def main():
    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 a < i:
                dp[i][j] += dp[i-1][j] * j % mod2
            if b < j:
                dp[i][j] += dp[i][j-1] * i % mod2
            if a < i and b < j:
                dp[i][j] -= dp[i-1][j-1] * (i-1) * (j-1) % mod2

    return dp[c][d] % mod2


print(main())
# Original solution
a,b,c,d=list(map(int,input().split()))
mod=998244353
dct=[[[-1]*3 for _ in range(d+1)] for _ in range(c+1)]

import sys
sys.setrecursionlimit(10**7)
def dp(i,j,bl):
  #print(i,j,bl)
  if dct[i][j][bl]>0:
    return dct[i][j][bl]
  if (i,j)==(a,b):
    return 1 if bl==0 else 0
  if i<a or j<b:
    return 0
  if i==a:
    if bl==2:
      ret=0 if j-b<2 else pow(i,j-b,mod)-pow(i-1,j-b-1,mod)*(j-b)-pow(i-1,j-b,mod)
    elif bl==1:
      ret=0 if j-b<1 else pow(i-1,j-b-1,mod)*(j-b)
    elif bl==0:
      ret=pow(i-1,j-b,mod)
  elif bl==2:
    ret=dp(i,j-1,1)+dp(i,j-1,2)*i
  elif bl==1:
    ret=(dp(i-1,j,0)+dp(i-1,j,1)+dp(i-1,j,2))*j
  else:
    ret=dp(i,j-1,0)*(i-1)
  ret%=mod
  dct[i][j][bl]=ret
  return ret

print(((dp(c,d,0)+dp(c,d,1)+dp(c,d,2))%mod))


# Optimized solution
a,b,c,d=list(map(int,input().split()))
mod=998244353

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

for i in range(b+1,d+1):
  if a==1:
    dp1[a][i]=1 if i==b+1 else 0
    dp2[a][i]=0 if i==b+1 else 1
  else:
    dp1[a][i]=(pow(a-1,i-b-1,mod)*(i-b))%mod
    dp2[a][i]=(pow(a,i-b,mod)-dp1[a][i])%mod
for i in range(a+1,c+1):
  dp1[i][b]=pow(b,i-a,mod)
  dp2[i][b]=0
dp2[a][b]=1
for i in range(a+1,c+1):
  for j in range(b+1,d+1):
    dp1[i][j]=(j*(dp1[i-1][j]+dp2[i-1][j]))%mod
    dp2[i][j]=(dp1[i][j-1]+dp2[i][j-1]*i)%mod
print(((dp1[c][d]+dp2[c][d])%mod))

"""
print('dp1')
for x in dp1:
  print(x)
print('dp2')
for x in dp2:
  print(x)
"""
# Original solution
def main():
    A, B, C, D = list(map(int, input().split()))

    c = C - A
    d = D - B
    mod = 998244353
    dp = [[0]*(d+1) for i in range(c+1)]
    dp[0][0] = 1

    for i in range(c+1):
        prev = dp[i-1]
        crr = dp[i]
        if i < c:
            next = dp[i+1]
        for j in range(d+1):
            if i > 0 and j > 0:
                if i == j == 1:
                    crr[j] -= (A+i-1) * (B+j-1)
                else:
                    tmp = prev[j-1]
                    #print(tmp, i, j)
                    tmp *= (A+i-1)*(B+j-1)

                    # print(crr[j], tmp)
                    crr[j] -= tmp

            if i < c:
                next[j] += crr[j]*(B+j)
            if j < d:
                crr[j+1] += crr[j]*(A+i)
            # print(dp)
            crr[j] %= mod

    print((dp[c][d]))


main()


# Optimized solution
def main():
    A, B, C, D = list(map(int, input().split()))

    c = C - A
    d = D - B
    mod = 998244353
    dp = [[0]*(d+1) for i in range(c+1)]
    dp[0][0] = 1

    for i in range(c+1):

        for j in range(d+1):
            ret = 0
            if i == 0 and j == 0:
                continue
            if i > 0 and j > 0:
                if i == j == 1:
                    ret -= (A+i-1) * (B+j-1)
                else:
                    ret -= dp[i-1][j-1]
                    ret *= (A+i-1)*(B+j-1)
            ret += dp[i-1][j]*(B+j) + dp[i][j-1]*(A+i)

            ret %= mod
            dp[i][j] = ret
    print((dp[c][d]))


main()
# Original solution
def main():
    A, B, C, D = list(map(int, input().split()))
    mod = 998244353
    dp = [[0] * (D + 1) for _ in range(C + 1)]
    dp[A][B] = 1
    for b in range(B + 1, D + 1):
        dp[A][b] = dp[A][b - 1] * A
    for a in range(A + 1, C + 1):
        dp[a][B] = dp[a - 1][B] * B % mod
        for b in range(B + 1, D + 1):
            dp[a][b] = (dp[a - 1][b] * b + dp[a][b - 1] * a - dp[a - 1][b - 1] * (a - 1) * (b - 1)) % mod

    return dp[C][D]

print((main()))


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