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

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

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

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

dp = [[0]*(d+2) 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][j-1]*i+dp[i-1][j]*j-dp[i-1][j-1]*(i-1)*(j-1))%mod
        dp[i][j]%=mod
print(dp[-1][-2]%mod)

Solution (Python):

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

Solution (Python):

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)

Solution (Python):

import math
import collections
import bisect


def main():
    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)


if __name__ == '__main__':
    main()

Solution (Python):
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()
 
 
if __name__ == "__main__":
    
    solve()

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

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

Solution (Python):
a,b,c,d=map(int,input().split());m=998244353;x=range(c+d);D=[[0]*6011for _ in[0]*6011];D[a][b]=1
for i in x:
 for j in x:D[i+1][j+1]+=(D[i+1][j]*(i+1)+D[i][j+1]*(j+1)-D[i][j]*i*j)%m
print(D[c][d])

Solution (Python):
a,b,c,d=map(int,input().split());m=998244353;n=6011;x=range(c+d);D=[[0]*6011for _ in[0]*6011];D[a][b]=1
for i in x:
 for j in x:D[i+1][j+1]+=(D[i+1][j]*(i+1)+D[i][j+1]*(j+1)-D[i][j]*i*j)%m
print(D[c][d])

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

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

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

Solution (Python):
a,b,c,d=map(int,input().split());m=998244353;x=range(3010);D=[[0]*3011for _ in[0]*3011];D[a][b]=1
for i in x:
 for j in x:D[i+1][j+1]+=(D[i+1][j]*(i+1)+D[i][j+1]*(j+1)-D[i][j]*i*j)%m
print(D[c][d])

Solution (Python):
a,b,c,d=map(int,input().split());m=998244353;x=range(3010);D=[[0]*3011for _ in [0]*3011];D[a][b]=1
for i in x:
 for j in x:D[i+1][j+1]+=(D[i+1][j]*(i+1)+D[i][j+1]*(j+1)-D[i][j]*i*j)%m
print(D[c][d])

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

Solution (Python):
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(c):
 for j in range(d):D[i+1][j+1]+=(D[i+1][j]*(i+1)+D[i][j+1]*(j+1)-D[i][j]*i*j)%m
print(D[c][d])

Solution (Python):
A,B,C,D=map(int,input().split())
dp=[[0 for _ in range(D+1)] for _ in range(C+1)]
dp[A][B]=1
p=998244353

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

print(dp[C][D])

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

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

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

Solution (Python):

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()  
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][j] -= dp[i-1][j-1] * (i-1) * (j-1)
        dp[i][j] %= MOD

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

Solution (Python):
import sys
import numpy as np

 
read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines
 
MOD = 998_244_353
 
def main(A, B, C, D):
	d = np.zeros((C + 1, D + 1), np.int64)
	d[A,B]=1
	for i in range(A+1,C+1):
  		d[i,B] = B*d[i-1,B] % MOD
	for j in range(B+1,D+1):
		d[A,j] = A*d[A,j-1] % MOD
 
	for i in range(A+1,C+1):
		for j in range(B+1,D+1):
			d[i,j]=((((j*d[i-1,j]) % MOD + i*d[i,j-1] % MOD) % MOD) -((i-1)*((j-1)*d[i-1,j-1] % MOD)% MOD)% MOD)
    
	return d[C,D] % MOD
          
"""    
    print(d[C][D] % MOD)
    dp = np.zeros((C + 1, D + 1), np.int64)
    dp[C, D] = 1
    for d in range(D, 0, -1):
        dp[C, d - 1] = dp[C, d] * C % MOD
    for c in range(C - 1, -1, -1):
        x = 0
        for d in range(D, -1, -1):
            dp[c, d] = (d * dp[c + 1, d] + x) % MOD
            x += dp[c + 1, d]
            x *= c
            x %= MOD
    return dp[A, B]
"""

if sys.argv[-1] == 'ONLINE_JUDGE':
    import numba
    from numba.pycc import CC
    i8 = numba.from_dtype(np.int64)
    signature = (i8, i8, i8, i8)
 
    cc = CC('my_module')
    cc.export('main', signature)(main)
    cc.compile()
 
from my_module import main
 
A, B, C, D = map(int, read().split())
print(main(A, B, C, D))

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

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

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

Solution (Python):
import sys

def main():
    input = sys.stdin.readline
    mod = 998244353
    
    a, b, c, d = map(int, input().split())
    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
                continue
            if j == b:
                dp[i][j] = dp[i-1][j]*j%mod
                continue
            if i == a:
                dp[i][j] = dp[i][j-1]*i%mod
                continue
            dp[i][j] = (dp[i-1][j]*j+dp[i][j-1]*i+mod-dp[i-1][j-1]*(i-1)*(j-1)%mod)%mod
    print(dp[c][d])        
            

if __name__ == "__main__":
    main()

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


if __name__ == '__main__':
    main()

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

Solution (Python):


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):
    '''H is number of rows
    A列、B列が与えられるようなとき
    ex1)A,B=read_col(H)    ex2) A,=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):
    '''H is number of rows'''
    ret = []
    for _ in range(H):
        ret.append(tuple(map(int, read().split())))
    return ret


def read_matrix(H):
    '''H is number of rows'''
    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

'''
A+i,B+jまで広げたときを考える
このときの通りの総数は

更新則



これちがう↓
(1+(A+i-1,B+j)まで広げたの通りの総数)*(1+(A+i,B+j-1)の通りの総数)-1 となる
'''


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

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

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

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

if __name__ == '__main__':
    main()

Solution (Python):

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

Solution (Python):
a,b,c,d=map(int,input().split())
MOD=998244353
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][j]+=dp[i-1][j]*j
      dp[i][j]-=dp[i-1][j-1]*(i-1)*(j-1)
      dp[i][j]%=MOD
print(dp[c][d]%MOD)

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

print(X[C][D])

Solution (Python):
import bisect
import copy
import heapq
import math
import sys
from collections import *
from functools import lru_cache
from itertools import accumulate, combinations, permutations, product
def input():
    return sys.stdin.readline()[:-1]
def ruiseki(lst):
    return [0]+list(accumulate(lst))
sys.setrecursionlimit(500000)
mod=998244353
al=[chr(ord('a') + i) for i in range(26)]
direction=[[1,0],[0,1],[-1,0],[0,-1]]

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

Solution (Python):
a,b,c,d=map(int,input().split())
mod=998244353
dp=[(d+1)*[0]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
    if i==a:
      dp[i][j]=dp[i][j-1]*i%mod
      continue
    if j==b:
      dp[i][j]=dp[i-1][j]*j%mod
      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])

Solution (Python):
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+1):
    for j in range(b,d+1):
        dp[i][j] += dp[i][j-1]*i+dp[i-1][j]*j-dp[i-1][j-1]*(i-1)*(j-1)
        dp[i][j] %= mod
print(dp[c][d])

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

Solution (Python):

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

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

Solution (Python):
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+1,c+1):
    dp[i][b] = dp[i-1][b]*b
    dp[i][b] %= mod
for i in range(a,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)
        dp[i][j] %= mod
print(dp[c][d])

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

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

print(dp[C][D])

Solution (Python):
if __name__ == '__main__':
    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)

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

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

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

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

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

Solution (Python):

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

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


if __name__ == "__main__":
    main()

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


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

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

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

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