Problem:
 Score : 1800 points    Problem Statement We have a square grid with N rows and M columns. Takahashi will write an integer in each of the squares, as follows:   First, write 0 in every square.  For each i=1,2,...,N , choose an integer k_i  (0\leq k_i\leq M) , and add 1 to each of the leftmost k_i squares in the i -th row.  For each j=1,2,...,M , choose an integer l_j  (0\leq l_j\leq N) , and add 1 to each of the topmost l_j squares in the j -th column.   Now we have a grid where each square contains 0 , 1 , or 2 . Find the number of different grids that can be made this way, modulo 998244353 .
We consider two grids different when there exists a square with different integers.      Constraints  1 \leq N,M \leq 5\times 10^5  N and M are integers.         Input Input is given from Standard Input in the following format:  N  M       Output Print the number of different grids that can be made, modulo 998244353 .        Sample Input 1 1 2      Sample Output 1 8  Let (a,b) denote the grid where the square to the left contains a and the square to the right contains b . Eight grids can be made: (0,0),(0,1),(1,0),(1,1),(1,2),(2,0),(2,1), and (2,2) .       Sample Input 2 2 3      Sample Output 2 234       Sample Input 3 10 7      Sample Output 3 995651918       Sample Input 4 314159 265358      Sample Output 4 70273732

Solution (Python):
N, M = map(int, input().split())
mod = 998244353

if N > M:
  N, M = M, N

fact = [0] * (M + 1)
powm = [0] * (M + 1)
pown = [0] * (M + 1)
fact[0] = 1
powm[0] = 1
pown[0] = 1
for i in range(1, M + 1):
  fact[i] = fact[i - 1] * i % mod
  powm[i] = powm[i - 1] * (M + 1) % mod
  pown[i] = pown[i - 1] * (N + 1) % mod

def pow(n, m):
  if m == 0:
    return 1
  elif m == 1:
    return n
  elif m % 2 == 0:
    return pow(n, m // 2)**2 % mod
  else:
    return pow(n, m // 2)**2 % mod * n % mod

inv_fact = [0] * (M + 1)
inv_fact[M] = pow(fact[M], mod-2)
for i in reversed(range(0, M)):
  inv_fact[i] = inv_fact[i + 1] * (i + 1) % mod
    
def C(n, r):
  return fact[n] * inv_fact[r] % mod * inv_fact[n - r] % mod
  
ans = 0
for i in range(N+1):
  ans += (-1)**i * C(N, i) * C(M, i) * fact[i] * powm[N - i] * pown[M - i]
  ans = ans % mod

print(ans)