AtCoder Grand Contest 046 - Extension

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.

Examples

Input

1 1 2 2


Output

3


Input

2 1 3 4


Output

65


Input

31 41 59 265


Output

387222020

Source: atcoder | Difficulty: UNKNOWN_DIFFICULTY


Solutions:
(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])


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


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


(A, B, C, D) = map(int, input().split())
MOD = 998244353
dp = [[0] * (D + 2) for _ in range(C + 2)]
dp[A][B] = 1
for h in range(A, C + 1):
	for w in range(B, D + 1):
		d = dp[h][w] % MOD
		dp[h + 1][w] += d * w
		dp[h][w + 1] += d * h
		dp[h + 1][w + 1] -= d * w * h % MOD
print(dp[C][D] % MOD)


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


(A, B, C, D) = 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])


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


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


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)


(A, B, C, D) = map(int, input().split())
mod = 998244353
dp = [[0] * (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:
			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])


import sys, re
from collections import deque, defaultdict, Counter
from math import ceil, sqrt, hypot, factorial, pi, sin, cos, radians
from itertools import accumulate, permutations, combinations, product, groupby, combinations_with_replacement
from operator import itemgetter, mul
from copy import deepcopy
from string import ascii_lowercase, ascii_uppercase, digits
from bisect import bisect, bisect_left
from fractions import gcd
from heapq import heappush, heappop
from functools import reduce

def input():
	return sys.stdin.readline().strip()

def INT():
	return int(input())

def MAP():
	return map(int, input().split())

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

def ZIP(n):
	return zip(*(MAP() for _ in range(n)))
sys.setrecursionlimit(10 ** 9)
INF = float('inf')
mod = 998244353
(A, B, C, D) = MAP()
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 j > B:
			dp[i][j] += dp[i][j - 1] * i
		if i > A:
			dp[i][j] += dp[i - 1][j] * j
		if A < i and B < j:
			dp[i][j] -= (i - 1) * (j - 1) * dp[i - 1][j - 1]
		dp[i][j] %= mod
print(dp[C][D])


(a, b, c, d) = map(int, input().split())
mod = 998244353
if a == c and b == d:
	print(1)
	exit()
else:
	dp = [[0] * (d + 1) for _ 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
			else:
				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])


import sys
input = sys.stdin.buffer.readline
sys.setrecursionlimit(10 ** 7)
from collections import defaultdict
con = 998244353
INF = float('inf')

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

def main():
	(A, B, C, D) = getlist()
	DP = [[0] * (D + 1) for i in range(C + 1)]
	DP[A][B] = 1
	for i in range(B + 1, D + 1):
		DP[A][i] = DP[A][i - 1] * A % con
	for i in range(A + 1, C + 1):
		DP[i][B] = DP[i - 1][B] * B % con
	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 - DP[i - 1][j - 1] * (i - 1) * (j - 1)) % con
	ans = DP[C][D]
	print(ans)
main()


import bisect
import copy
import heapq
import sys
import itertools
import math
import queue
from functools import lru_cache
input = sys.stdin.readline
sys.setrecursionlimit(1000000)
mod = 998244353

def read_values():
	return map(int, input().split())

def read_index():
	return map(lambda x: int(x) - 1, input().split())

def read_list():
	return list(read_values())

def read_lists(N):
	return [read_list() for n in range(N)]

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


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


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)


import math
import queue
import bisect
from collections import deque, defaultdict
import heapq as hpq
from sys import stdin, setrecursionlimit
ipt = stdin.readline
setrecursionlimit(10 ** 7)
mod = 998244353
dir = [(-1, 0), (1, 0), (0, -1), (0, 1)]
alp = 'abcdefghijklmnopqrstuvwxyz'

def main(a, b, c, d):
	dp = [[0] * (d + 1) for i 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
			dpi = dp[i + 1][j] * (i + 1) + dp[i][j + 1] * (j + 1) - dp[i][j] * (i * j)
			dp[i + 1][j + 1] = dpi % mod
	return dp[c][d]
	return None
(a, b, c, d) = map(int, ipt().split())
print(main(a, b, c, d))


def main():
	import sys
	input = sys.stdin.readline
	mod = 998244353
	(A, B, C, D) = map(int, input().split())
	dp = [[1] * (D + 1) for _ in range(C + 1)]
	a = A
	for (i, w) in enumerate(range(B + 1, D + 1)):
		dp[A][w] = a
		a *= A
		a %= mod
	b = B
	for (i, h) in enumerate(range(A + 1, C + 1)):
		dp[h][B] = b
		b *= B
		b %= mod
	for h in range(A + 1, C + 1):
		for w in range(B + 1, D + 1):
			dp[h][w] = (dp[h - 1][w] * w % mod + dp[h][w - 1] * h % mod - dp[h - 1][w - 1] * (h - 1) * (w - 1) % mod) % mod
	print(dp[C][D])
main()


import sys
input = sys.stdin.readline
(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(c - a + 1):
	for j in range(d - b + 1):
		if i > 0:
			dp[i][j] += dp[i - 1][j] * (b + j)
		if j > 0:
			dp[i][j] += dp[i][j - 1] * (a + i)
		if i > 0 and j > 0:
			dp[i][j] -= dp[i - 1][j - 1] * (b + j - 1) * (a + i - 1)
		dp[i][j] %= 998244353
print(dp[-1][-1])


(a, b, c, d) = map(int, input().split())
dp = [[0] * d for _ in range(c)]
dp[a - 1][b - 1] = 1
mod = 998244353
if a == c and b == d:
	print(1)
	exit()
if a == c:
	print(a ** (d - b) % mod)
	exit()
if b == d:
	print(b ** (c - a) % mod)
	exit()
for y in range(a - 1, c):
	for x in range(b - 1, d):
		dp[y][x] += dp[y][x - 1] * (y + 1) + dp[y - 1][x] * (x + 1) - dp[y - 1][x - 1] * x * y
		dp[y][x] = dp[y][x] % mod
print(dp[-1][-1] % mod)


import sys
input = sys.stdin.readline
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):
	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])


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

def main():
	MOD = 998244353
	(a, b, c, d) = map(int, readline().split())
	(h, w) = (c - a, d - b)
	dp = [[1] * (d + 1) for _ in range(c + 1)]
	for row in range(a + 1, c + 1):
		dp[row][b] = dp[row - 1][b] * b % MOD
	for col in range(b + 1, d + 1):
		dp[a][col] = dp[a][col - 1] * a % MOD
	for row in range(a + 1, c + 1):
		for col in range(b + 1, d + 1):
			x = dp[row - 1][col] * col + dp[row][col - 1] * row
			y = dp[row - 1][col - 1] * (row - 1) * (col - 1)
			dp[row][col] = (x - y) % MOD
	print(dp[c][d])
main()


import sys
input = sys.stdin.readline

def I():
	return int(input())

def MI():
	return map(int, input().split())

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

def main():
	mod = 998244353
	(A, B, C, D) = MI()
	if A == C and B == D:
		print(1)
		exit()
	F = [[0] * (D + 1) for _ in range(C + 1)]

	def calc(a, b):
		temp = F[a][b - 1] * a + F[a - 1][b] * b - F[a - 1][b - 1] * (a - 1) * (b - 1)
		return temp % mod
	F[A][B] = 1
	for i in range(A + 1, C + 1):
		F[i][B] = F[i - 1][B] * B % mod
	for j in range(B + 1, D + 1):
		F[A][j] = F[A][j - 1] * A % mod
	for i in range(A + 1, C + 1):
		for j in range(B + 1, D + 1):
			F[i][j] = calc(i, j)
	print(F[C][D])
main()


(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):
		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])


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

def amari(a):
	return a % 998244353
X0 = [[0] * d for i in range(c)]
X1 = [[0] * d for i in range(c)]
X2 = [[0] * d for i in range(c)]
X0[a - 1][b - 1] = 1
for j in range(b, d):
	X0[a - 1][j] = amari(X0[a - 1][j - 1] * (a - 1))
	X1[a - 1][j] = amari(X0[a - 1][j - 1] + X1[a - 1][j - 1] * (a - 1))
	X2[a - 1][j] = amari(X1[a - 1][j - 1] + X2[a - 1][j - 1] * a)
for i in range(a, c):
	for j in range(b, d + 1):
		X1[i][j - 1] = amari((X0[i - 1][j - 1] + X1[i - 1][j - 1] + X2[i - 1][j - 1]) * j)
		if j >= 2:
			X2[i][j - 1] = amari(X1[i][j - 2] + X2[i][j - 2] * (i + 1))
print(amari(X0[c - 1][d - 1] + X1[c - 1][d - 1] + X2[c - 1][d - 1]))


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


mod = 998244353
(A, B, C, D) = map(int, input().split())
dp = [0] * (D + 1)
dp[B:D + 1] = [pow(A, i, mod) for i in range(D + 1 - B)]
for a in range(A + 1, C + 1):
	T = [0] * (D + 1)
	for b in range(B, D + 1):
		T[b] = (a * T[b - 1] + b * dp[b] - (a - 1) * (b - 1) * dp[b - 1]) % mod
	dp = T
print(dp[D])


(a, b, c, d) = map(int, input().split())
mod = 998244353
DP = [[0] * (d + 1) for _ in range(c + 1)]
DP[a][b] = 1
for y in range(a + 1, c + 1):
	DP[y][b] = DP[y - 1][b] * b % mod
for x in range(b + 1, d + 1):
	DP[a][x] = DP[a][x - 1] * a % mod
for y in range(a + 1, c + 1):
	for x in range(b + 1, d + 1):
		DP[y][x] = (DP[y - 1][x] * x + (DP[y][x - 1] - DP[y - 1][x - 1] * (x - 1)) * y + DP[y - 1][x - 1] * (x - 1)) % mod
print(DP[c][d])


(A, B, C, D) = (int(x) for x in input().split())
MOD = 998244353
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
		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])


(a, b, c, d) = [int(x) for x in input().split()]
MOD = 998244353
dp = [[0 for i in range(d + 10)] for j in range(c + 10)]
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)
		dp[i][j] %= MOD
print(dp[c][d])


(a, b, c, d) = 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])


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


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


(a, b, c, d) = map(int, input().split())
p = [[0] * (d + 2) for i in range(c + 1)]
p[a][b] = 1
m = 998244353
for i in range(a, c + 1):
	for j in range(b, d + 1):
		p[i][j] += (j * p[i - 1][j] + i * p[i][j - 1] - (i - 1) * (j - 1) * p[i - 1][j - 1]) % m
		p[i][j] %= m
print(p[-1][-2] % m)


import sys
import math
import fractions
from collections import deque
from collections import defaultdict
sys.setrecursionlimit(10 ** 7)
(A, B, C, D) = map(int, input().split())
mod = 998244353
dp = [[0] * (D + 1) for _ in range(C + 1)]
dp[A][B] = 1
for i in range(A, C):
	dp[i + 1][B] = dp[i][B] * B % mod
for i in range(B, D):
	dp[A][i + 1] = dp[A][i] * A % mod
for i in range(A, C):
	for j in range(B, D):
		dp[i + 1][j + 1] += dp[i + 1][j] * (i + 1)
		dp[i + 1][j + 1] %= mod
		dp[i + 1][j + 1] += dp[i][j + 1] * (j + 1)
		dp[i + 1][j + 1] %= mod
		dp[i + 1][j + 1] -= dp[i][j] * i % mod * j % mod
		if dp[i + 1][j + 1] < 0:
			dp[i + 1][j + 1] += mod
print(dp[C][D])


import sys

def input():
	return sys.stdin.readline().strip()

def mapint():
	return map(int, input().split())
sys.setrecursionlimit(10 ** 9)
(A, B, C, D) = mapint()
mod = 998244353
dp = [[0] * D for _ in range(C)]
dp[A - 1][B - 1] = 1
for i in range(A - 1, C):
	for j in range(B - 1, D):
		if i == A - 1 and j == B - 1:
			continue
		val = 0
		if i > A - 1:
			val += dp[i - 1][j] * (j + 1)
		if j > B - 1:
			val += dp[i][j - 1] * (i + 1)
		if i > A - 1 and j > B - 1:
			val -= dp[i - 1][j - 1] * i * j
		dp[i][j] = val % mod
print(dp[C - 1][D - 1])


import sys
sr = lambda : sys.stdin.readline().rstrip()
ir = lambda : int(sr())
lr = lambda : list(map(int, sr().split()))
MOD = 998244353
(A, B, C, D) = lr()
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):
		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)


(a, b, h, w) = map(int, input().split())
dp = [[0 for _ in range(w + 1)] for _ in range(h + 1)]
mod = 998244353
for i in range(1, h + 1):
	if i < a:
		continue
	for j in range(1, w + 1):
		if j < b:
			continue
		if i == a and j == b:
			dp[i][j] = 1
			continue
		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[-1][-1])


MOD = 998244353
inpl = lambda : list(map(int, input().split()))
(A, B, C, D) = inpl()
dp = [[0] * (D + 1) for _ in range(C + 1)]
dp[A][B] = 1
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
for a in range(A + 1, C + 1):
	for b in range(B + 1, D + 1):
		dp[a][b] = (b * dp[a - 1][b] + a * dp[a][b - 1] - (a - 1) * (b - 1) * dp[a - 1][b - 1]) % MOD
print(dp[C][D])


def main():
	(A, B, C, D) = map(int, input().split())
	mod = 998244353
	dp = [[0] * (D + 2) for i in range(C + 2)]
	dp[A][B] = 1
	for i in range(A, C + 1):
		for j in range(B, D + 1):
			ret = dp[i][j]
			if i - A == j - B == 1:
				ret -= (i - 1) * (j - 1)
			else:
				ret -= dp[i - 1][j - 1] * (j - 1) * (i - 1)
			ret %= mod
			dp[i][j] = ret
			dp[i + 1][j] += ret * j
			dp[i][j + 1] += ret * i
	print(dp[-2][-2])
main()


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])
main()


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


(A, B, C, D) = map(int, input().split())
MOD = 998244353
dp = [[1] * C for _ in [0] * D]
if B < D:
	dp[B][A - 1] = A
for i in range(B + 1, D):
	j = A - 1
	dp[i][j] = dp[i - 1][j] * A % MOD
if A < C:
	dp[B - 1][A] = B
for j in range(A + 1, C):
	i = B - 1
	dp[i][j] = dp[i][j - 1] * B % MOD
for i in range(B, D):
	for j in range(A, C):
		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])


(a, b, c, d) = map(int, input().split())
mod = 998244353
dp = []
for i in range(c + 1):
	temp = [0] * (d + 1)
	dp.append(temp)
dp[a][b] = 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
		else:
			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])


(A, B, C, D) = 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])


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


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


import sys
from collections import defaultdict

def solve():
	input = sys.stdin.readline
	(A, B, C, D) = 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])
	return 0
solve()


import math
import collections
import bisect
(a, b, c, d) = 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):
		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)


p = 998244353
(a, b, c, d) = map(int, input().split())
dp = [[0 for _ in range(d + 1)] for k in range(c + 1)]
for i in range(a, c + 1):
	for j in range(b, d + 1):
		if i == a:
			if j == b:
				dp[i][j] = 1
			else:
				dp[i][j] = dp[i][j - 1] * i % p
		elif j == b:
			dp[i][j] = dp[i - 1][j] * j % p
		else:
			dp[i][j] = (dp[i - 1][j] * j + dp[i][j - 1] * i - dp[i - 1][j - 1] * (i - 1) * (j - 1)) % p
print(dp[c][d])


(a, b, c, d) = map(int, input().split())
p = 998244353
n = c - a + 2
m = d - b + 2
data = [[0] * m for i in range(n)]
data[1][1] = 1
for i in range(1, n):
	gyou = a + i - 1
	for j in range(1, m):
		retu = b + j - 1
		if i * j == 1:
			continue
		data[i][j] = data[i - 1][j] * retu % p + data[i][j - 1] * gyou % p - data[i - 1][j - 1] * (gyou - 1) % p * (retu - 1) % p
		data[i][j] = data[i][j] % p
print(data[n - 1][m - 1] % p)


(A, B, C, D) = map(int, input().split())
MOD = 998244353
dp0 = [[0] * (D + 1) for _ in range(C + 1)]
dp1 = [[0] * (D + 1) for _ in range(C + 1)]
dp0[A][B] = 1
for x in range(A, C + 1):
	for y in range(B, D + 1):
		if x == A and y == B:
			continue
		dp0[x][y] = (dp0[x][y - 1] * x + dp1[x][y - 1] * 1) % MOD
		dp1[x][y] = (dp0[x - 1][y] * y + dp1[x - 1][y] * y) % MOD
ans = (dp0[C][D] + dp1[C][D]) % MOD
print(ans)


(a, b, c, d) = map(int, input().split())
dp = [[1 for i in range(d + 1)] for j in range(c + 1)]
mod = 998244353
for y in range(c + 1):
	if y < a:
		continue
	for x in range(d + 1):
		if x < b:
			continue
		if y == a and x == b:
			continue
		py = 0
		px = 0
		m = 0
		if y > a:
			py = dp[y - 1][x] * x
		if x > b:
			px = dp[y][x - 1] * y
		if y > a and x > b:
			m = dp[y - 1][x - 1] * (y - 1) * (x - 1)
		dp[y][x] = py + px - m
		dp[y][x] %= mod
print(dp[c][d])


import sys
sys.setrecursionlimit(1 << 25)
read = sys.stdin.readline
ra = range
enu = enumerate

def mina(*argv, sub=1):
	return list(map(lambda x: x - sub, argv))

def read_a_int():
	return int(read())

def read_ints():
	return list(map(int, read().split()))

def read_col(H):
	ret = []
	for _ in range(H):
		ret.append(list(map(int, read().split())))
	return tuple(map(list, zip(*ret)))

def read_tuple(H):
	ret = []
	for _ in range(H):
		ret.append(tuple(map(int, read().split())))
	return ret

def read_matrix(H):
	ret = []
	for _ in range(H):
		ret.append(list(map(int, read().split())))
	return ret

def combination_mod(n, r, mod=998244353):
	if r > n:
		return 0
	r = min(r, n - r)
	nf = rf = 1
	for i in range(r):
		nf = nf * (n - i) % mod
		rf = rf * (i + 1) % mod
	return nf * pow(rf, mod - 2, mod) % mod
MOD = 998244353
INF = 2 ** 31
from collections import defaultdict, Counter, deque
from operator import itemgetter
from itertools import product, permutations, combinations
from bisect import bisect_left, bisect_right
from math import gcd
(A, B, C, D) = read_ints()
i_max = C - A
j_max = D - B
dp = [[0] * (j_max + 1) for _ in ra(i_max + 1)]
for i in ra(i_max + 1):
	dp[i][0] = pow(B, i, MOD)
for j in ra(j_max + 1):
	dp[0][j] = pow(A, j, MOD)
for (i, j) in product(range(1, i_max + 1), range(1, j_max + 1)):
	dp[i][j] = dp[i][j - 1] * (A + i) + dp[i - 1][j] * (B + j) - (A + i - 1) * (B + j - 1)
	dp[i][j] = dp[i][j - 1] * (A + i) + dp[i - 1][j] * (B + j) - (A + i - 1) * (B + j - 1) * dp[i - 1][j - 1]
	dp[i][j] %= MOD
print(dp[-1][-1])


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


import sys

def input():
	return sys.stdin.readline().strip()

def list2d(a, b, c):
	return [[c] * b for i in range(a)]

def list3d(a, b, c, d):
	return [[[d] * c for j in range(b)] for i in range(a)]

def list4d(a, b, c, d, e):
	return [[[[e] * d for j in range(c)] for j in range(b)] for i in range(a)]

def ceil(x, y=1):
	return int(-(-x // y))

def INT():
	return int(input())

def MAP():
	return map(int, input().split())

def LIST(N=None):
	return list(MAP()) if N is None else [INT() for i in range(N)]

def Yes():
	print('Yes')

def No():
	print('No')

def YES():
	print('YES')

def NO():
	print('NO')
sys.setrecursionlimit(10 ** 9)
INF = 10 ** 19
MOD = 998244353
(a, b, c, d) = MAP()
dp = list2d(c + 1, d + 1, 0)
dp[a][b] = 1
for i in range(a, c + 1):
	for j in range(b, d + 1):
		dp[i][j] += dp[i - 1][j] * j
		dp[i][j] += dp[i][j - 1] * i
		dp[i][j] -= dp[i - 1][j - 1] * (i - 1) * (j - 1)
		dp[i][j] %= MOD
ans = dp[c][d]
print(ans)


(a, b, c, d) = map(int, input().split())
dp = [[1] * (d - b + 1) for _ in range(c - a + 1)]
mod = 998244353
for h in range(1, c - a + 1):
	dp[h][0] = pow(b, h, mod)
for w in range(1, d - b + 1):
	dp[0][w] = pow(a, w, mod)
for w in range(1, d - b + 1):
	for h in range(1, c - a + 1):
		dp[h][w] = dp[h][w - 1] * (h + a) + dp[h - 1][w] * (w + b) - dp[h - 1][w - 1] * (w + b - 1) * (h + a - 1)
		dp[h][w] %= mod
print(dp[c - a][d - b])


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)


mod = 998244353
(A, B, C, D) = map(int, input().split())
X = D - B
Y = C - A
if X == 0:
	print(pow(B, Y, mod))
	exit()
if Y == 0:
	print(pow(A, X, mod))
	exit()
dp = [B, A]
for x in range(1, X):
	dp.append(dp[-1] * A % mod)
for y in range(1, Y):
	dp_new = [None] * (X + 1)
	dp_new[0] = dp[0] * B % mod
	a = 0
	for x in range(1, X + 1):
		a = (a + dp[x - 1]) * (A + y) % mod
		b = dp[x] * (B + x) % mod
		dp_new[x] = (a + b) % mod
	dp = dp_new
k = 1
ans = 0
for d in dp[::-1]:
	ans += d * k
	ans %= mod
	k *= C
	k %= mod
print(ans)


import sys
sys.setrecursionlimit(10 ** 7)
MOD = 998244353
(A, B, C, D) = map(int, input().split())

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


def main():
	(A, B, C, D) = 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()


import sys
input = sys.stdin.readline
import math
import bisect
import heapq
from collections import defaultdict
from collections import deque
from collections import Counter
from functools import lru_cache
MOD = 10 ** 9 + 7
INF = float('inf')
AZ = 'abcdefghijklmnopqrstuvwxyz'

def I():
	return int(input().strip())

def S():
	return input().strip()

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

def SL():
	return list(map(str, input().split()))

def ILs(n):
	return list((int(input()) for _ in range(n)))

def SLs(n):
	return list((input().strip() for _ in range(n)))

def ILL(n):
	return [list(map(int, input().split())) for _ in range(n)]

def SLL(n):
	return [list(map(str, input().split())) for _ in range(n)]

def P(arg):
	print(arg)
	return

def Y():
	print('Yes')
	return

def N():
	print('No')
	return

def E():
	exit()

def PE(arg):
	print(arg)
	exit()

def YE():
	print('Yes')
	exit()

def NE():
	print('No')
	exit()

def DD(arg):
	return defaultdict(arg)

def inv(n):
	return pow(n, MOD - 2, MOD)
kaijo_memo = []

def kaijo(n):
	if len(kaijo_memo) > n:
		return kaijo_memo[n]
	if len(kaijo_memo) == 0:
		kaijo_memo.append(1)
	while len(kaijo_memo) <= n:
		kaijo_memo.append(kaijo_memo[-1] * len(kaijo_memo) % MOD)
	return kaijo_memo[n]
gyaku_kaijo_memo = []

def gyaku_kaijo(n):
	if len(gyaku_kaijo_memo) > n:
		return gyaku_kaijo_memo[n]
	if len(gyaku_kaijo_memo) == 0:
		gyaku_kaijo_memo.append(1)
	while len(gyaku_kaijo_memo) <= n:
		gyaku_kaijo_memo.append(gyaku_kaijo_memo[-1] * pow(len(gyaku_kaijo_memo), MOD - 2, MOD) % MOD)
	return gyaku_kaijo_memo[n]

def nCr(n, r):
	if n == r:
		return 1
	if n < r or r < 0:
		return 0
	ret = 1
	ret = ret * kaijo(n) % MOD
	ret = ret * gyaku_kaijo(r) % MOD
	ret = ret * gyaku_kaijo(n - r) % MOD
	return ret

def factorization(n):
	arr = []
	temp = n
	for i in range(2, int(-(-n ** 0.5 // 1)) + 1):
		if temp % i == 0:
			cnt = 0
			while temp % i == 0:
				cnt += 1
				temp //= i
			arr.append([i, cnt])
	if temp != 1:
		arr.append([temp, 1])
	if arr == []:
		arr.append([n, 1])
	return arr

def make_divisors(n):
	divisors = []
	for i in range(1, int(n ** 0.5) + 1):
		if n % i == 0:
			divisors.append(i)
			if i != n // i:
				divisors.append(n // i)
	return divisors

def make_primes(N):
	max = int(math.sqrt(N))
	seachList = [i for i in range(2, N + 1)]
	primeNum = []
	while seachList[0] <= max:
		primeNum.append(seachList[0])
		tmp = seachList[0]
		seachList = [i for i in seachList if i % tmp != 0]
	primeNum.extend(seachList)
	return primeNum

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

def lcm(a, b):
	return a * b // gcd(a, b)

def count_bit(n):
	count = 0
	while n:
		n &= n - 1
		count += 1
	return count

def base_10_to_n(X, n):
	if X // n:
		return base_10_to_n(X // n, n) + [X % n]
	return [X % n]

def base_n_to_10(X, n):
	return sum((int(str(X)[-i - 1]) * n ** i for i in range(len(str(X)))))

def int_log(n, a):
	count = 0
	while n >= a:
		n //= a
		count += 1
	return count
(A, B, C, D) = IL()
(H, W) = (C - A + 1, D - B + 1)
dp = [[0 for w in range(W)] for h in range(H)]
for w in range(W):
	dp[0][w] = pow(A, w, 998244353)
for h in range(H):
	dp[h][0] = pow(B, h, 998244353)
for w in range(1, W):
	for h in range(1, H):
		dp[h][w] = (dp[h - 1][w] + dp[h][w - 1] + dp[h - 1][w] * (B + w - 1) + dp[h][w - 1] * (A + h - 1) - dp[h - 1][w - 1] * (B + w - 1) * (A + h - 1)) % 998244353
print(dp[-1][-1])


import sys
sys.setrecursionlimit(10 ** 7)
input = sys.stdin.readline
f_inf = float('inf')
mod = 998244353

def resolve():
	(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][j - 1] * i + dp[i - 1][j] * j - dp[i - 1][j - 1] * (i - 1) * (j - 1)
			dp[i][j] %= mod
	print(dp[-1][-1])
resolve()


import sys, bisect, string, math, time, functools, random, fractions
from heapq import heappush, heappop, heapify
from collections import deque, defaultdict, Counter
from itertools import permutations, combinations, groupby
rep = range

def Golf():
	(n, *t) = map(int, open(0).read().split())

def I():
	return int(input())

def S_():
	return input()

def IS():
	return input().split()

def LS():
	return [i for i in input().split()]

def MI():
	return map(int, input().split())

def LI():
	return [int(i) for i in input().split()]

def LI_():
	return [int(i) - 1 for i in input().split()]

def NI(n):
	return [int(input()) for i in range(n)]

def NI_(n):
	return [int(input()) - 1 for i in range(n)]

def StoLI():
	return [ord(i) - 97 for i in input()]

def ItoS(n):
	return chr(n + 97)

def LtoS(ls):
	return ''.join([chr(i + 97) for i in ls])

def RA():
	return map(int, open(0).read().split())

def GI(V, E, ls=None, Directed=False, index=1):
	org_inp = []
	g = [[] for i in range(V)]
	FromStdin = True if ls == None else False
	for i in range(E):
		if FromStdin:
			inp = LI()
			org_inp.append(inp)
		else:
			inp = ls[i]
		if len(inp) == 2:
			(a, b) = inp
			c = 1
		else:
			(a, b, c) = inp
		if index == 1:
			a -= 1
			b -= 1
		aa = (a, c)
		bb = (b, c)
		g[a].append(bb)
		if not Directed:
			g[b].append(aa)
	return (g, org_inp)

def GGI(h, w, search=None, replacement_of_found='.', mp_def={'#': 1, '.': 0}, boundary=1):
	mp = [boundary] * (w + 2)
	found = {}
	for i in range(h):
		s = input()
		for char in search:
			if char in s:
				found[char] = (i + 1) * (w + 2) + s.index(char) + 1
				mp_def[char] = mp_def[replacement_of_found]
		mp += [boundary] + [mp_def[j] for j in s] + [boundary]
	mp += [boundary] * (w + 2)
	return (h + 2, w + 2, mp, found)

def TI(n):
	return GI(n, n - 1)

def accum(ls):
	rt = [0]
	for i in ls:
		rt += [rt[-1] + i]
	return rt

def bit_combination(n, base=2):
	rt = []
	for tb in range(base ** n):
		s = [tb // base ** bt % base for bt in range(n)]
		rt += [s]
	return rt

def gcd(x, y):
	if y == 0:
		return x
	if x % y == 0:
		return y
	while x % y != 0:
		(x, y) = (y, x % y)
	return y

def YN(x):
	print(['NO', 'YES'][x])

def Yn(x):
	print(['No', 'Yes'][x])

def show(*inp, end='\n'):
	if show_flg:
		print(*inp, end=end)
mo = 10 ** 9 + 7
inf = float('inf')
FourNb = [(-1, 0), (1, 0), (0, 1), (0, -1)]
EightNb = [(-1, 0), (1, 0), (0, 1), (0, -1), (1, 1), (-1, -1), (1, -1), (-1, 1)]
compas = dict(zip('WENS', FourNb))
cursol = dict(zip('LRUD', FourNb))
l_alp = string.ascii_lowercase
read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
input = lambda : sys.stdin.readline().rstrip()

class Tree:

	def __init__(self, inp_size=None, ls=None, init=True, index=1):
		self.LCA_init_stat = False
		self.ETtable = []
		if init:
			if ls == None:
				self.stdin(inp_size, index=index)
			else:
				self.node_size = len(ls) + 1
				(self.edges, _) = GI(self.node_size, self.node_size - 1, ls, index=index)
		return

	def stdin(self, inp_size=None, index=1):
		if inp_size == None:
			self.node_size = int(input())
		else:
			self.node_size = inp_size
		(self.edges, _) = GI(self.node_size, self.node_size - 1, index=index)
		return

	def listin(self, ls, index=0):
		self.node_size = len(ls) + 1
		(self.edges, _) = GI(self.node_size, self.node_size - 1, ls, index=index)
		return

	def dfs(self, x, func=lambda pr, prv, nx, dist: prv + dist, root_v=0):
		q = deque([x])
		v = [None] * self.node_size
		v[x] = root_v
		while q:
			c = q.pop()
			for (nb, d) in self.edges[c]:
				if v[nb] == None:
					q.append(nb)
					v[nb] = func(c, v[c], nb, d)
		return v

	def bfs(self, x, func=lambda pr, prv, nx, dist: prv + dist, root_v=0):
		q = deque([x])
		v = [None] * self.node_size
		v[x] = root_v
		while q:
			c = q.popleft()
			for (nb, d) in self.edges[c]:
				if v[nb] == None:
					q.append(nb)
					v[nb] = func(c, v[c], nb, d)
		return v

	def parent(self, x):
		return self.dfs(0, func=lambda pr, prv, nx, dist: pr, root_v=-1)

	def topological_sort(self, x):
		tps = []
		q = deque([x])
		v = [None] * self.node_size
		v[x] = 0
		while q:
			c = q.popleft()
			tps.append(c)
			for (nb, d) in self.edges[c]:
				if v[nb] == None:
					q.append(nb)
					v[nb] = 0
		return tps

	def EulerTour(self, x):
		q = deque()
		q.append(x)
		self.depth = [None] * self.node_size
		self.depth[x] = 0
		self.ETtable = []
		self.ETdepth = []
		self.ETin = [-1] * self.node_size
		self.ETout = [-1] * self.node_size
		cnt = 0
		while q:
			c = q.pop()
			if c < 0:
				ce = ~c
			else:
				ce = c
				for (nb, d) in self.edges[ce]:
					if self.depth[nb] == None:
						q.append(~ce)
						q.append(nb)
						self.depth[nb] = self.depth[ce] + 1
			self.ETtable.append(ce)
			self.ETdepth.append(self.depth[ce])
			if self.ETin[ce] == -1:
				self.ETin[ce] = cnt
			else:
				self.ETout[ce] = cnt
			cnt += 1
		return

	def LCA_init(self, root):
		self.EulerTour(root)
		self.st = SparseTable(self.ETdepth, init_func=min, init_idl=inf)
		self.LCA_init_stat = True
		return

	def LCA(self, root, x, y):
		if self.LCA_init_stat == False:
			self.LCA_init(root)
		(xin, xout) = (self.ETin[x], self.ETout[x])
		(yin, yout) = (self.ETin[y], self.ETout[y])
		a = min(xin, yin)
		b = max(xout, yout, xin, yin)
		id_of_min_dep_in_et = self.st.query_id(a, b + 1)
		return self.ETtable[id_of_min_dep_in_et]

	def __str__(self):
		return str(self.edges)

	def show(self):
		if all([all([d == 1 for (nd, d) in e]) for e in self.edges]):
			print([[nd for (nd, d) in e] for e in self.edges])
		else:
			print(self)

class SparseTable:

	def __init__(self, ls, init_func=min, init_idl=float('inf')):
		self.func = init_func
		self.idl = init_idl
		self.size = len(ls)
		self.N0 = self.size.bit_length()
		self.table = [ls[:]]
		self.index = [list(range(self.size))]
		self.lg = [0] * (self.size + 1)
		for i in range(2, self.size + 1):
			self.lg[i] = self.lg[i >> 1] + 1
		for i in range(self.N0):
			tmp = [self.func(self.table[i][j], self.table[i][min(j + (1 << i), self.size - 1)]) for j in range(self.size)]
			tmp_id = [self.index[i][j] if self.table[i][j] == self.func(self.table[i][j], self.table[i][min(j + (1 << i), self.size - 1)]) else self.index[i][min(j + (1 << i), self.size - 1)] for j in range(self.size)]
			self.table += [tmp]
			self.index += [tmp_id]

	def query(self, l, r):
		if r > self.size:
			r = self.size
		N = self.lg[r - l]
		return self.func(self.table[N][l], self.table[N][max(0, r - (1 << N))])

	def query_id(self, l, r):
		if r > self.size:
			r = self.size
		N = self.lg[r - l]
		(a, b) = (self.index[N][l], self.index[N][max(0, r - (1 << N))])
		if self.table[0][a] == self.func(self.table[N][l], self.table[N][max(0, r - (1 << N))]):
			b = a
		return b

	def right_bound(self, pos, x):
		k = (self.size - pos).bit_length()
		for j in range(k)[::-1]:
			nx = pos + (1 << j)
			if nx < n and self.query(pos, nx + 1) <= x:
				pos += 1 << j
		return pos

	def left_bound(self, pos, x):
		k = pos.bit_length()
		for j in range(k)[::-1]:
			nx = pos - (1 << j)
			if 0 <= nx and self.query(nx, pos + 1) <= x:
				pos -= 1 << j
		return pos

	def __str__(self):
		return str(self.table[0])

	def print(self):
		for i in self.table:
			print(*i)
show_flg = False
show_flg = True
ans = 0

def LI():
	return [int(i) for i in input().split()]
mo = 998244353
(a, b, c, d) = LI()
h = c - a + 1
w = d - b + 1
dp = [[0] * (d + 1) for i in range(c + 1)]
dp[a][b] = 1
for i in range(a, c):
	dp[i + 1][b] = dp[i][b] * b
	dp[i + 1][b] %= mo
for i in range(b, d):
	dp[a][i + 1] = dp[a][i] * a
	dp[a][i + 1] %= mo
for i in range(a, c):
	for j in range(b, d):
		dp[i + 1][j + 1] = dp[i][j + 1] + dp[i + 1][j] + dp[i][j + 1] * j + dp[i + 1][j] * i - dp[i][j] * i * j
		dp[i + 1][j + 1] %= mo
ans = dp[c][d]
print(ans)


from sys import stdin
input = stdin.readline
(a, b, c, d) = map(int, input().split())
MAX = 998244353
rowcol = [[0] * (d + 1) for _ in range(c + 1)]
rowcol[a][b] = 1
row2 = [[0] * (d + 1) for _ in range(c + 1)]
col2 = [[0] * (d + 1) for _ in range(c + 1)]
for i in range(b + 1, d + 1):
	rowcol[a][i] = rowcol[a][i - 1] * a % MAX
for j in range(a + 1, c + 1):
	for i in range(b, d + 1):
		if i == b:
			rowcol[j][b] = rowcol[j - 1][b] * b % MAX
		else:
			rowcol[j][i] = (rowcol[j - 1][i] + row2[j - 1][i]) * (i - 1) % MAX
			row2[j][i] = (rowcol[j][i - 1] + col2[j][i - 1] + row2[j][i - 1] * j) % MAX
			col2[j][i] = (rowcol[j - 1][i] + col2[j - 1][i] * i + row2[j - 1][i]) % MAX
print((rowcol[c][d] + row2[c][d] + col2[c][d]) % MAX)


MOD = 998244353
(A, B, C, D) = map(int, input().split())
dp = [[1] * (D + 1) for _ in range(C + 1)]
powA = 1
for i in range(A + 1, C + 1):
	powA *= B
	powA %= MOD
	dp[i][B] = powA
powB = 1
for i in range(B + 1, D + 1):
	powB *= A
	powB %= MOD
	dp[A][i] = powB
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][j] -= (i - 1) * (j - 1) * dp[i - 1][j - 1]
		dp[i][j] %= MOD
print(dp[C][D])


(a, b, c, d) = map(int, input().split())
m = c - a
n = d - b
t = [[1]]
for i in range(m):
	t[0].append(b * t[0][i] % 998244353)
for i in range(n):
	t.append([a * t[i][0] % 998244353] + [0] * m)
for i in range(1, n + 1):
	for j in range(1, m + 1):
		t[i][j] = t[i][j - 1] * (b + i) + t[i - 1][j] * (a + j) - t[i - 1][j - 1] * ((a + j - 1) * (b + i - 1))
		t[i][j] %= 998244353
print(t[n][m])


MOD = 998244353
(A, B, C, D) = map(int, input().split())
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
		elif i == A:
			dp[i][j] = dp[i][j - 1] * i % MOD
		elif j == B:
			dp[i][j] = dp[i - 1][j] * j % MOD
		else:
			dp[i][j] = (dp[i][j - 1] * i + dp[i - 1][j] * j - dp[i - 1][j - 1] * (i - 1) * (j - 1)) % MOD
print(dp[C][D])


M = 998244353
(A, B, C, D) = map(int, input().split())
table_a = [None] * (D - B + 1)
table_h = [None] * (D - B + 1)
table_v = [None] * (D - B + 1)
tmp_table_a = [None] * (D - B + 1)
tmp_table_h = [None] * (D - B + 1)
tmp_table_v = [None] * (D - B + 1)
table_a[0] = 1
c = 1
for x in range(1, D - B + 1):
	c = c * A % M
	table_a[x] = c
	table_v[x] = c
c = 1
for y in range(1, C - A + 1):
	c = c * B % M
	tmp_table_a[0] = c
	tmp_table_h[0] = c
	for x in range(1, D - B + 1):
		tmp_table_v[x] = tmp_table_a[x - 1] * (A + y) % M
		tmp_table_h[x] = table_a[x] * (B + x) % M
		tmp_table_a[x] = (tmp_table_h[x] + (tmp_table_a[x - 1] - tmp_table_h[x - 1]) * (A + y) + tmp_table_h[x - 1]) % M
	(table_a, tmp_table_a) = (tmp_table_a, table_a)
	(table_v, tmp_table_v) = (tmp_table_v, table_v)
	(table_h, tmp_table_h) = (tmp_table_h, table_h)
print(table_a[D - B])


import sys
input = sys.stdin.readline
mod = 998244353
NN = 10 ** 4
inv = [0]
for i in range(1, NN):
	inv.append(pow(i, mod - 2, mod))

def main():
	(A, B, C, D) = map(int, input().split())
	dp = [[0] * (D + 1) for _ in range(C + 1)]
	dp2 = [[0] * (D + 1) for _ in range(C + 1)]
	for w in range(B - 1, D + 1):
		dp[A][w] = 0
		dp2[A][w] = pow(A, D - B, mod)
	for c in range(A + 1, C + 1):
		if D - B - 1 >= 0:
			t = pow(c * inv[c - 1] % mod, D - B - 1, mod)
		kakeru = (c - 1) * inv[c] % mod
		dp[c][B] = pow(B, c - A, mod) * pow(c, D - B, mod) % mod
		dp2[c][B] = dp[c][B]
		for d in range(B + 1, D + 1):
			dp[c][d] = (dp2[c - 1][d - 1] + dp[c - 1][d] * d % mod) * t % mod
			t = t * kakeru % mod
			dp2[c][d] = (dp2[c][d - 1] + dp[c][d]) % mod
	print(dp2[C][D])
main()


(A, B, C, D) = map(int, input().split())
MOD = 998244353
a = 1
arr = []
for j in range(D - B + 1):
	arr.append(a)
	a *= A
	a %= MOD
dp = [arr]
b = B
for i in range(1, C - A + 1):
	arr = [b]
	for j in range(1, D - B + 1):
		x = arr[-1]
		y = dp[-1][j]
		z = dp[-1][j - 1]
		h = A + i
		w = B + j
		arr.append((x * h + y * w - z * ((h - 1) * (w - 1))) % MOD)
	b *= B
	b %= MOD
	dp.append(arr)
print(dp[C - A][D - B])


MOD = 998244353

def print_dp(dp, c):
	for i in range(1, c + 1):
		print(dp[i][1:])

def main():
	(a, b, c, d) = map(int, input().split())
	dp = [[0] * (d + 1) for _ in range(c + 1)]
	t = [[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):
			dp[i][j] += dp[i - 1][j] * j
			dp[i][j] += dp[i][j - 1] * i
			dp[i][j] -= t[i][j - 1] * (i - 1)
			dp[i][j] %= MOD
			t[i][j] += dp[i - 1][j] * j
			t[i][j] %= MOD
	print(dp[c][d])
main()


import sys
readline = sys.stdin.buffer.readline

def even(n):
	return 1 if n % 2 == 0 else 0
(A, B, C, D) = map(int, readline().split())
mod = 998244353
dp0 = [[0] * (D + 1) for _ in range(C + 1)]
dp1 = [[0] * (D + 1) for _ in range(C + 1)]
dp0[A][B] = 1
for x in range(A, C + 1):
	for y in range(B, D + 1):
		if x == A and y == B:
			continue
		dp0[x][y] = (dp0[x][y - 1] * x + dp1[x][y - 1] * 1) % mod
		dp1[x][y] = (dp0[x - 1][y] * y + dp1[x - 1][y] * y) % mod
ans = (dp0[C][D] + dp1[C][D]) % mod
print(ans)


(a, b, c, d) = 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)


(A, B, C, D) = map(int, input().split())
MOD = 998244353
MAX_SIZE = 3010
dp = [[0] * MAX_SIZE for _ in range(MAX_SIZE)]
dp[A][B] = 1
for h in range(A, C + 1):
	for w in range(B, D + 1):
		if h == A and w == B:
			continue
		dp[h][w] = dp[h - 1][w] * w + dp[h][w - 1] * h - dp[h - 1][w - 1] * (h - 1) * (w - 1)
		dp[h][w] %= MOD
print(dp[C][D])


import sys
sys.setrecursionlimit(10 ** 8)
from functools import lru_cache
INF = float('inf')
MOD = 998244353

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


(A, B, C, D) = list(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 + 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 - 1][j] * j + DP[i][j - 1] * i - DP[i - 1][j - 1] * (i - 1) * (j - 1)) % mod
print(DP[C][D] % mod)


(A, B, C, D) = map(int, input().split())
p = 998244353
dp = [[0] * (D + 1) for _ in range(C + 1)]
for h in range(A, C + 1):
	for w in range(B, D + 1):
		if h == A and w == B:
			dp[A][B] = 1
			continue
		if h > A and w > B:
			x = dp[h][w - 1] * h
			y = dp[h - 1][w] * w
			z = dp[h - 1][w - 1] * (h - 1) * (w - 1)
			sub = x + y - z
			dp[h][w] = sub % p
		elif h == A:
			dp[h][w] = dp[h][w - 1] * h % p
		else:
			dp[h][w] = dp[h - 1][w] * w % p
print(dp[C][D])


from sys import stdin
input = stdin.readline
MOD = 998244353

def main():
	(A, B, C, D) = map(int, input().rstrip().split())
	DP = [[0 for _ in range(D + 1)] for _ 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)) % MOD
	print(DP[C][D] % MOD)
main()


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


import sys
stdin = sys.stdin

def ni():
	return int(ns())

def na():
	return list(map(int, stdin.readline().split()))

def naa(N):
	return [na() for _ in range(N)]

def ns():
	return stdin.readline().rstrip()
(A, B, C, D) = na()
mod = 998244353
dp = [[0] * (D + 1) for _ in range(C + 1)]
dp[A][B] = 1
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
for a in range(A + 1, C + 1):
	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
print(dp[C][D])


(a, b, c, d) = [int(_) for _ in input().split()]
mod = 998244353
dp = [0] * 3005
dp[b] = 1
for y in range(b + 1, d + 1):
	dp[y] = dp[y - 1] * a % mod
for x in range(a + 1, c + 1):
	dp2 = [0] * 3005
	dp2[b] = dp[b] * b % mod
	for y in range(b + 1, d + 1):
		dp2[y] = (dp[y] * y + dp2[y - 1] * x - dp[y - 1] * (x - 1) * (y - 1)) % mod
	dp = dp2
print(dp[d])


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


from collections import defaultdict, deque, Counter
from heapq import heappush, heappop, heapify
import math
import bisect
import random
from itertools import permutations, accumulate, combinations, product
import sys
import string
from bisect import bisect_left, bisect_right
from math import factorial, ceil, floor, gcd
from operator import mul
from functools import reduce
from operator import mul
import pprint
sys.setrecursionlimit(2147483647)
INF = 10 ** 20

def LI():
	return list(map(int, sys.stdin.buffer.readline().split()))

def I():
	return int(sys.stdin.buffer.readline())

def LS():
	return sys.stdin.buffer.readline().rstrip().decode('utf-8').split()

def S():
	return sys.stdin.buffer.readline().rstrip().decode('utf-8')

def IR(n):
	return [I() for i in range(n)]

def LIR(n):
	return [LI() for i in range(n)]

def SR(n):
	return [S() for i in range(n)]

def LSR(n):
	return [LS() for i in range(n)]

def SRL(n):
	return [list(S()) for i in range(n)]

def MSRL(n):
	return [[int(j) for j in list(S())] for i in range(n)]
mod = 998244353
(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):
		dp[i][j] += dp[i - 1][j] * j
		dp[i][j] += dp[i][j - 1] * i
		dp[i][j] -= dp[i - 1][j - 1] * (i - 1) * (j - 1)
		dp[i][j] %= mod
print(dp[c][d])


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
	return (dp1[c][d] + dp2[c][d] + dp3[c][d]) % MOD

def main():
	(a, b, c, d) = 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
test()
main()


from sys import stdin
input = stdin.readline
printn = lambda x: print(x, end='')
inn = lambda : int(input())
inl = lambda : list(map(int, input().split()))
inm = lambda : map(int, input().split())
ins = lambda : input().strip()
DBG = True
BIG = 10 ** 18
R = 998244353

def ddprint(x):
	if DBG:
		print(x)
(a, b, c, d) = inm()
dp = [[0] * (d + 1) for i in range(c + 1)]
dp[a][b] = 1
for i in range(a + 1, c + 1):
	dp[i][b] = dp[i - 1][b] * b % R
for j in range(b + 1, d + 1):
	dp[a][j] = dp[a][j - 1] * a % R
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)) % R
print(dp[c][d])


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


mod = 998244353
(a, b, c, d) = map(int, input().split())
dp0 = [0] * (c - a + 1)
dp1 = [1]
for i in range(c - a):
	dp1.append(dp1[-1] * b % mod)
for j in range(1, d - b + 1):
	new_dp0 = [(dp0[0] * a + dp1[0]) % mod]
	new_dp1 = [dp1[0] * (a - 1) % mod]
	for i in range(1, c - a + 1):
		new_dp0.append((dp0[i] * (a + i) + dp1[i]) % mod)
		new_dp1.append((new_dp0[-2] + new_dp1[-1]) * (b + j) % mod)
	dp0 = new_dp0
	dp1 = new_dp1
print((dp0[-1] + dp1[-1]) % mod)


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]
print(main())


(a, b, c, d) = map(int, input().split())
mod = 998244353
dp = [[0 for _ in range(c - a + 1)] for _ in range(d - b + 1)]
dp[0][0] = 1
for j in range(c - a + 1):
	for i in range(d - b + 1):
		if i == j == 0:
			dp[i][j] = 1
		elif i == 0:
			dp[i][j] = dp[i][j - 1] * b % mod
		elif j == 0:
			dp[i][j] = dp[i - 1][j] * a % mod
		else:
			tmp1 = (j + a - 1) * dp[i - 1][j] % mod
			tmp2 = (i + b - 1) * dp[i][j - 1] % mod
			tmp3 = -dp[i - 1][j - 1] * (i + b - 1) * (j + a - 1) % mod
			tmp4 = (dp[i - 1][j] + dp[i][j - 1]) % mod
			dp[i][j] = (tmp1 + tmp2 + tmp3 + tmp4) % mod
print(dp[-1][-1])


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


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


Q = 998244353

def main():
	(A, B, C, D) = map(int, input().split())
	dp = [[0] * 3002 for _ in range(3001)]
	dp[C][D] = 1
	for a in range(C - 1, A - 1, -1):
		dp[a][D] = dp[a + 1][D] * D % Q
	for b in range(D - 1, B - 1, -1):
		dp[C][b] = dp[C][b + 1] * C % Q
	for a in range(C - 1, A - 1, -1):
		z = dp[a + 1][D] * a % Q
		for b in range(D - 1, B - 1, -1):
			dp[a][b] = (b * dp[a + 1][b] + z) % Q
			z = (z + dp[a + 1][b]) * a % Q
	print(dp[A][B])
main()


import sys
sys.setrecursionlimit(10 ** 8)
input = sys.stdin.readline

def main():
	(A, B, C, D) = [int(x) for x in input().split()]
	MOD = 998244353
	dp = [[0] * (D - B + 1) for j in range(C - A + 1)]
	for j in range(C - A + 1):
		if j == 0:
			dp[j][0] = 1
		else:
			dp[j][0] = pow(B, j, MOD)
		for i in range(1, D - B + 1):
			if j == 0:
				dp[j][i] = dp[j][i - 1] * (A + j) % MOD
			else:
				dp[j][i] = (dp[j][i - 1] * (A + j) + (B + i) * dp[j - 1][i] - dp[j - 1][i - 1] * ((A + j - 1) * (B + i - 1))) % MOD
	print(dp[-1][-1] % MOD)
main()


import sys, collections as cl, bisect as bs
sys.setrecursionlimit(100000)
input = sys.stdin.readline
mod = 998244353
Max = sys.maxsize

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

def m():
	return map(int, input().split())

def onem():
	return int(input())

def s(x):
	a = []
	if len(x) == 0:
		return []
	aa = x[0]
	su = 1
	for i in range(len(x) - 1):
		if aa != x[i + 1]:
			a.append([aa, su])
			aa = x[i + 1]
			su = 1
		else:
			su += 1
	a.append([aa, su])
	return a

def jo(x):
	return ' '.join(map(str, x))

def max2(x):
	return max(map(max, x))

def In(x, a):
	k = bs.bisect_left(a, x)
	if k != len(a) and a[k] == x:
		return True
	else:
		return False

def pow_k(x, n):
	ans = 1
	while n:
		if n % 2:
			ans *= x
		x *= x
		n >>= 1
	return ans
(a, b, c, d) = m()
ans = 0
po = c - a
poo = d - b
dp = [[0 for i in range(poo + 1)] for j in range(po + 1)]
dp[0][0] = 1
for i in range(po + 1):
	for j in range(poo + 1):
		if i == 0:
			if j != 0:
				dp[i][j] = dp[i][j - 1] * a
				dp[i][j] %= mod
		elif j == 0:
			dp[i][j] = dp[i - 1][j] * b
			dp[i][j] %= mod
		else:
			dp[i][j] = dp[i - 1][j] * (b + j) + dp[i][j - 1] * (a + i) - (a + i - 1) * (b + j - 1) * dp[i - 1][j - 1]
			dp[i][j] %= mod
print(dp[-1][-1])


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

def doubling(n, m, modulo=0):
	y = 1
	base = n
	tmp = m
	while tmp:
		if tmp % 2 == 1:
			y *= base
			if modulo > 0:
				y %= modulo
		base *= base
		if modulo > 0:
			base %= modulo
		tmp //= 2
	return y
for i in range(b, d):
	dp[a][i + 1] = dp[a][i] * a
	dp[a][i + 1] %= mod
for i in range(a, c):
	dp[i + 1][b] = dp[i][b] * b
	dp[i + 1][b] %= mod
for i in range(a + 1, c + 1):
	for j in range(b + 1, d + 1):
		dp[i][j] = j * dp[i - 1][j] + i * dp[i][j - 1] - (i - 1) * (j - 1) * dp[i - 1][j - 1]
		dp[i][j] %= mod
print(dp[c][d])


import sys
MOD = 998244353

def main(A, B, C, D):
	(I, J) = (C - A, D - B)
	dp = [[0] * (I + 1) for _ in range(I + J + 1)]
	dp[0][0] = 1
	for k in range(1, I + J + 1):
		for i in range(min(I, k) + 1):
			j = k - i
			v = dp[k - 1][i - 1] * (j + B) % MOD if i >= 1 else 0
			h = dp[k - 1][i] * (i + A) % MOD
			d = dp[k - 2][i - 1] * (i + A - 1) * (j + B - 1) if k >= 2 and i >= 1 else 0
			dp[k][i] = (v + h - d) % MOD
	print(dp[I + J][I])
input = sys.stdin.readline
(A, B, C, D) = map(int, input().split())
main(A, B, C, D)


import sys
import math
input = sys.stdin.readline
getS = lambda : input().strip()
getN = lambda : int(input())
getList = lambda : list(map(int, input().split()))
getZList = lambda : [int(x) - 1 for x in input().split()]
INF = 10 ** 20
MOD = 998244353

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

def main():
	n = getN()
	for _ in range(n):
		solve()
solve()


(a, b, c, d) = map(int, input().split())
dp = [[1 for _ in range(d - b + 1)] for _ in range(c - a + 1)]
mod = 998244353
for x in range(c - a + 1):
	for y in range(d - b + 1):
		if x == 0:
			if y > 0:
				dp[x][y] = dp[x][y - 1] * (a + x) % mod
		elif y == 0:
			dp[x][y] = dp[x - 1][y] * (b + y) % mod
		else:
			dp[x][y] = (dp[x - 1][y] * (b + y) + dp[x][y - 1] * (a + x) - dp[x - 1][y - 1] * (a + x - 1) * (b + y - 1)) % mod
print(dp[-1][-1])


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])
main()


MOD = 998244353

def solve(A, B, C, D):
	(H, W) = (C - A, D - B)
	dp = [[0] * (W + 1) for _ in range(H + 1)]
	t = 1
	for i in range(W + 1):
		dp[0][i] = t
		t *= A
		t %= MOD
	for j in range(1, H + 1):
		dp[j][0] = dp[j - 1][0] * B % MOD
		for i in range(1, W + 1):
			dp[j][i] = dp[j - 1][i] * (B + i) + dp[j][i - 1] * (A + j) - dp[j - 1][i - 1] * (B + i - 1) * (A + j - 1)
			dp[j][i] %= MOD
	return dp[-1][-1]
print(solve(*map(int, input().split())))


(a, b, c, d) = map(int, input().split())
mod = 998244353
dp = [[0] * (d - b + 1) for i in range(c - a + 1)]
if a == c and b == d:
	print(1)
elif a == c:
	print(pow(a, d - b, mod))
elif b == d:
	print(pow(b, c - a, mod))
else:
	for i in range(c - a):
		p = 0
		for j in range(d - b + 1):
			if i == j == 0:
				dp[i][j] = 1
			dp[i + 1][j] = (p + dp[i][j] * (b + j)) % mod
			p += dp[i][j]
			p %= mod
			p *= a + i
			p %= mod
	ans = 0
	for d in dp[-1]:
		ans *= c
		ans %= mod
		ans += d
		ans %= mod
	print(ans)


def main():
	MOD = 998244353
	(A, B, H, W) = map(int, input().split())
	dp = [[0] * (W + 1) for _ in range(H + 1)]
	dp[A][B] = 1
	for h in range(A, H + 1):
		for w in range(B, W + 1):
			if h == A and w == B:
				continue
			t = 0
			if h > A:
				t += dp[h - 1][w] * w
			if w > B:
				t += dp[h][w - 1] * h
			if h > A and w > B:
				t -= dp[h - 1][w - 1] * (h - 1) * (w - 1)
			dp[h][w] = t % MOD
	print(dp[H][W])
main()


import sys
input = lambda : sys.stdin.readline().rstrip()
sys.setrecursionlimit(max(1000, 10 ** 9))
write = lambda x: sys.stdout.write(x + '\n')
(a, b, c, d) = map(int, input().split())
dp = [[0] * (d + 1) for _ in range(c + 1)]
dp[a][b] = 1
M = 998244353
for x in range(a + 1, c + 1):
	dp[x][b] = b * dp[x - 1][b]
	dp[x][b] %= M
for y in range(b + 1, d + 1):
	dp[a][y] = a * dp[a][y - 1]
	dp[a][y] %= M
for x in range(a + 1, c + 1):
	for y in range(b + 1, d + 1):
		dp[x][y] = dp[x - 1][y] + dp[x][y - 1] + (x - 1) * dp[x][y - 1] + (y - 1) * dp[x - 1][y] - (x - 1) * (y - 1) * dp[x - 1][y - 1]
		dp[x][y] %= M
print(dp[c][d])


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])
main()


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


import sys
sys.setrecursionlimit(10 ** 7)
input = sys.stdin.readline
f_inf = float('inf')
mod = 998244353

def resolve():
	(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)
			dp[i][j] %= mod
	print(dp[-1][-1] % mod)
resolve()


(a, s, d, f) = map(int, input().split())
dp = [[0] * (f + 1) for i in range(d + 1)]
dp[a][s] = 1
mod = 998244353
for x in range(a, d + 1):
	for y in range(s, f + 1):
		if x == a and y == s:
			continue
		dp[x][y] = (dp[x - 1][y] * y + dp[x][y - 1] * x - dp[x - 1][y - 1] * (x - 1) * (y - 1)) % mod
print(dp[d][f])


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
		if i == a:
			dp[i][j] = dp[i][j - 1] * i
		elif j == b:
			dp[i][j] = dp[i - 1][j] * j
		else:
			dp[i][j] = dp[i - 1][j] + dp[i][j - 1]
			dp[i][j] += dp[i - 1][j - 1] * (i - 1) * (j - 1)
			dp[i][j] += (dp[i - 1][j] - dp[i - 1][j - 1] * (i - 1)) * (j - 1)
			dp[i][j] += (dp[i][j - 1] - dp[i - 1][j - 1] * (j - 1)) * (i - 1)
		dp[i][j] %= MOD
ans = dp[i][j]
print(ans)


(A, B, C, D) = map(int, input().split())
MOD = 998244353
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 i in range(1, min(C - A, D - B) + 1):
	for j in range(i, C - A + 1):
		DP[j][i] = DP[j - 1][i] * (B + i) + DP[j][i - 1] * (A + j) - DP[j - 1][i - 1] * (B + i - 1) * (A + j - 1)
		DP[j][i] %= MOD
	for j in range(i + 1, D - B + 1):
		DP[i][j] = DP[i - 1][j] * (B + j) + DP[i][j - 1] * (A + i) - DP[i - 1][j - 1] * (B + j - 1) * (A + i - 1)
		DP[i][j] %= MOD
print(DP[-1][-1])