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

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

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

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

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

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

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

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

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

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

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

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

dp[a][b] = 1
if d > b:
    dp[a][b+1] = a
if c > a:
    dp[a+1][b] = b
for i in range(b+2,d+1):
    dp[a][i] = dp[a][i-1]*a%mod

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

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


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

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

print(dp[c][d])

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

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

if __name__ == "__main__":
    main()

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 % MOD

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

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

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

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

Solution (Python):
def solve():
    MOD = 998244353

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

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

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

    print(dp[C][D])


solve()

Solution (Python):






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

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

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

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 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,C):
    for b in range(B,D):
        dp[a+1][b+1] = (dp[a+1][b]+dp[a][b+1]+a*dp[a+1][b]+b*dp[a][b+1]-a*b*dp[a][b])%mod
print(dp[C][D])

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

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

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

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

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

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

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

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

Solution (Python):

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


if __name__ == '__main__':
    main()

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

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

"""
def nibu(x,n,r):
	ll = 0
	rr = r
	while True:
		mid = (ll+rr)//2

	if rr == mid:
		return ll
	if (ここに評価入れる):
		rr = mid
	else:
		ll = mid+1
"""

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
		else:
			if 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])

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

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

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

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

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

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

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

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

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

if __name__ == '__main__':
	main()

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


if __name__ == '__main__':
    print(solve(*map(int,input().split())))

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

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

if __name__ == "__main__":
    solve()