E - Four Square Tiles
Source: atcoder
Difficulty: 800

Problem Statement
You are given positive integers
N
,
H
, and
W
, with
H,W \le 3N-1
.
Find the number, modulo
998244353
, of ways to place four
N\times N
square tiles on an
H\times W
grid that satisfy all of the following conditions.
Each tile exactly covers
N^2
cells of the grid.
No cell is covered by more than one tile.
Here, the tiles are indistinguishable.
There are
T
test cases; solve each one.
Constraints
1\le T\le 2\times 10^5
1\le N,H,W\le 10^9
H,W\le 3N-1
All input values are integers.
Input
The input is given from Standard Input in the following format:
T
\text{case}_1
\vdots
\text{case}_T
Each case is given in the following format:
N
H
W
Output
Print
T
lines. The
i
-th line should contain the number, modulo
998244353
, of valid ways to place the tiles for the
i
-th test case.
Sample Input 1
4
2 4 5
2 5 5
1000 1000 1000
1000 2222 2025
Sample Output 1
9
79
0
262210557
For the first test case, there are
9
ways as illustrated in the following figure:

C++ solution:
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
typedef long long ll;
using namespace std;
namespace fast_IO {
#define IOSIZE 100000
    int precision = 3, POW[10] = {1, 10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000, 1000000000};    // 向下取整到小数点后第 precision 位
    char ibuf[IOSIZE], obuf[IOSIZE], *p1 = ibuf, *p2 = ibuf, *p3 = obuf; 
#ifdef ONLINE_JUDGE
#define getchar() ((p1==p2)and(p2=(p1=ibuf)+fread(ibuf,1,IOSIZE,stdin),p1==p2)?(EOF):(*p1++))
#endif 
#define putchar(x) ((p3==obuf+IOSIZE)&&(fwrite(obuf,p3-obuf,1,stdout),p3=obuf),*p3++=x)
#define isdigit(ch) (ch>47&&ch<58)
#define isspace(ch) (ch<33)
    template<typename T> inline T read() { T s = 0; int w = 1; char ch; while (ch = getchar(), !isdigit(ch) and (ch != EOF)) if (ch == '-') w = -1; if (ch == EOF) return false; while (isdigit(ch)) s = s * 10 + ch - 48, ch = getchar(); return s * w; }
    template<typename T> inline bool read(T &s) { s = 0; int w = 1; char ch; while (ch = getchar(), !isdigit(ch) and (ch != EOF)) if (ch == '-') w = -1; if (ch == EOF) return false; while (isdigit(ch)) s = s * 10 + ch - 48, ch = getchar(); return s *= w, true; }
    template<typename T> inline void print(T x) { if (x < 0) putchar('-'), x = -x; if (x > 9) print(x / 10); putchar(x % 10 + 48); }
    inline bool read(char &s) { while (s = getchar(), isspace(s) and s != EOF); return s != EOF; }
    inline void print(char x) { putchar(x); }
    inline bool read(char *s) { char ch; while (ch = getchar(), isspace(ch) && ch != EOF); if (ch == EOF) return false; while (!isspace(ch) and (ch != EOF)) *s++ = ch, ch = getchar(); *s = '\0'; return true; }
    inline void print(char *x) { while (*x) putchar(*x++); }
    inline void print(const char *x) { for (int i = 0; x[i]; i++) putchar(x[i]); }
    inline bool read(std::string& s) { s = ""; char ch; while (ch = getchar(), isspace(ch) and ch != EOF); if (ch == EOF) return false; while (!isspace(ch)) s += ch, ch = getchar(); return true; }
    inline void print(std::string x) { for (int i = 0, n = x.size(); i < n; i++) putchar(x[i]); }
    inline bool read(bool &b) { char ch; while(ch=getchar(), isspace(ch) and ch != EOF); return ch == EOF ?false :(b=ch^48, true); }
    inline void print(bool b) { putchar(b+48); }
    inline bool read(double &x) { char ch = getchar(); int f = 1; for(; (ch<48 or 57<ch) and (ch != EOF); ch=getchar()) if(ch == '-')    f = -1; if(ch == EOF)    return false; for(x=0; 47<ch and ch<58; ch=getchar())    x = x*10 + (ch^48); if(ch != '.')    return x *= f, true; double y = 0.1; for(ch=getchar(); 47<ch and ch<58; ch=getchar())    x += y*(ch^48), y /= 10; return x *= f, true; }
    inline void print(double x) { if(x < 0)    putchar('-'), x = -x; print((unsigned long long)(x)); if(precision) putchar('.'), print((unsigned long long)((x-(unsigned long long)(x))*POW[precision])); }
    // double 暂时只支持向下取整
    template<typename T, typename... T1> inline int read(T& a, T1&... other) { return read(a) + read(other...); }
    template<typename T, typename... T1> inline void print(T a, T1... other) { print(a), print(other...); }
    struct Fast_IO { 
        bool flag = 1;
        inline ~Fast_IO() { fwrite(obuf, p3 - obuf, 1, stdout); } 
        inline void SetPrecision(int x) { if(x > 9)    throw runtime_error("Precision too high!"); else if(x < 0) throw runtime_error("Precision too low!"); else precision = x; } // 浮点数精度设为 x，即精确到小数点后 x 位。
        inline operator bool() { return flag; }
    } io;
    template<typename T> Fast_IO& operator >> (Fast_IO &io, T &b) { return io.flag &= read(b), io; }
    template<typename T> Fast_IO& operator << (Fast_IO &io, T b) { return print(b), io; }
#define cout io
#define cin io
#define endl '\n'
} using namespace fast_IO;
namespace furina{
const int N = 2e5, mod = 998244353, inf = 0x7fffffff >> 1;
int inv2, inv6;
inline int add(int a, int b){
    return a + b >= mod ? a + b - mod : a + b;
}
inline int sub(int a, int b){
    return a - b < 0 ? a - b + mod : a - b;
}
inline int pow(int a, int b){
    int ans = 1;
    for(; b; b >>= 1, a = 1ll * a * a % mod)
        if(b & 1)
            ans = 1ll * ans * a % mod;
    return ans;
}
int n, h, w;
inline int sum1(int n){
    return ((1ll * (n + 1) * n) >> 1) % mod;
}
inline int sum2(int n){
    return 1ll * n * (n + 1) % mod * ((n << 1) + 1) % mod * inv6 % mod;
}
inline int sum3(int n){
    return 1ll * sum1(n) * sum1(n) % mod;
}
inline int get(int m){
    m -= n << 1;
    int ans = 0;
    ans = add(ans, 1ll * inv2 * add(sum2(m), sum1(m)) % mod);
    ans = add(ans, 1ll * inv6 * add(add(2ll * sum3(m) % mod, 3ll * sum2(m) % mod), sum1(m)) % mod);
    return 1ll * ans * inv2 % mod;
}
inline void solve(){
    cin >> n >> h >> w;
    if(h < (n << 1) || w < (n << 1)){
        cout << "0\n";
        return;
    }
    int x = sub(1ll * (h - (n << 1) + 1) * (h - (n << 1) + 2) % mod, sum1(h - (n << 1) + 1));
    int y = sub(1ll * (w - (n << 1) + 1) * (w - (n << 1) + 2) % mod, sum1(w - (n << 1) + 1));
    int ans = 1ll * x * x % mod * y % mod * y % mod;
    x = get(h);
    y = get(w);
    cout << sub(ans, 2ll * x * y % mod) << '\n';
}
int main(){
    inv2 = pow(2, mod - 2), inv6 = pow(6, mod - 2);
    int T;
    cin >> T;
    while(T--)
        solve();
    return 0;
}
}
signed main(){
    return furina::main();
}

Go solution:
package main

import (
	"bufio"
	"fmt"
	"os"
	"strconv"
	"strings"
)

const MOD = 998244353

func main() {
	scanner := bufio.NewScanner(os.Stdin)
	scanner.Split(bufio.ScanLines)
	scanner.Scan()
	T, _ := strconv.Atoi(scanner.Text())

	output := make([]string, 0, T)

	// Precompute the modular inverse of 18 using Fermat's Little Theorem
	inv18 := pow(18, MOD-2, MOD)

	for i := 0; i < T; i++ {
		scanner.Scan()
		line := scanner.Text()
		parts := strings.Fields(line)
		N, _ := strconv.Atoi(parts[0])
		H, _ := strconv.Atoi(parts[1])
		W, _ := strconv.Atoi(parts[2])

		h := H - 2*N
		w := W - 2*N

		if h < 0 || w < 0 {
			output = append(output, "0")
			continue
		}

		f := func(x int) int {
			return x * (x + 1) / 2 % MOD
		}

		a := f(h+1) * f(w+1) % MOD
		term := f(h) * f(w) % MOD
		term = term * f(h+2) % MOD
		term = term * f(w+2) % MOD
		term = term * inv18 % MOD
		res := (a*a%MOD - term + MOD) % MOD

		output = append(output, strconv.Itoa(res))
	}

	fmt.Println(strings.Join(output, "\n"))
}

// Function to compute (base^exp) % mod efficiently using Fermat's Little Theorem
func pow(base, exp, mod int) int {
	result := 1
	base = base % mod
	for exp > 0 {
		if exp%2 == 1 {
			result = (result * base) % mod
		}
		base = (base * base) % mod
		exp = exp / 2
	}
	return result
}

Java solution:
import java.io.*;
import java.util.*;

public class Main {
    static final int MOD = 998244353;

    public static void main(String[] args) throws IOException {
        BufferedReader br = new BufferedReader(new InputStreamReader(System.in));
        PrintWriter pw = new PrintWriter(System.out);
        int T = Integer.parseInt(br.readLine());
        while (T-- > 0) {
            StringTokenizer st = new StringTokenizer(br.readLine());
            long N = Long.parseLong(st.nextToken());
            long H = Long.parseLong(st.nextToken());
            long W = Long.parseLong(st.nextToken());
            long h = H - 2 * N;
            long w = W - 2 * N;
            if (h < 0 || w < 0) {
                pw.println(0);
                continue;
            }
            long combH2 = comb(h + 2, 2);
            long combW2 = comb(w + 2, 2);
            long term1 = (combH2 * combH2) % MOD;
            long term2 = (combW2 * combW2) % MOD;
            long total = (term1 * term2) % MOD;
            long combH4 = comb(h + 3, 4);
            long combW4 = comb(w + 3, 4);
            long subtract = (2 * combH4 % MOD) * combW4 % MOD;
            total = (total - subtract + MOD) % MOD;
            pw.println(total);
        }
        pw.flush();
    }

    static long comb(long n, long k) {
        if (k == 2) {
            return n * (n - 1) / 2 % MOD;
        } else if (k == 4) {
            return n * (n - 1) % MOD * (n - 2) % MOD * (n - 3) % MOD * inv(24) % MOD;
        }
        return 0;
    }

    static long inv(long a) {
        return modPow(a, MOD - 2);
    }

    static long modPow(long a, long b) {
        long res = 1;
        a = a % MOD;
        while (b > 0) {
            if ((b & 1) != 0) {
                res = res * a % MOD;
            }
            a = a * a % MOD;
            b >>= 1;
        }
        return res;
    }
}

JavaScript solution:
const MOD = 998244353;

function solve() {
    const input = require('fs').readFileSync('/dev/stdin', 'utf8').trim().split('\n');
    const T = parseInt(input[0], 10);
    const cases = input.slice(1, T + 1);
    
    const results = [];
    
    for (let i = 0; i < T; i++) {
        const [n, h, w] = cases[i].split(' ').map(Number);
        const h_adj = h - 2 * n;
        const w_adj = w - 2 * n;
        
        if (h_adj < 0 || w_adj < 0) {
            results.push(0);
            continue;
        }
        
        // Calculate comb(h_adj + 2, 2) mod MOD
        const comb_h2 = ((h_adj + 2) * (h_adj + 1) / 2) % MOD;
        const comb_w2 = ((w_adj + 2) * (w_adj + 1) / 2) % MOD;
        const term1 = (comb_h2 * comb_w2) % MOD;
        const term1_sq = (term1 * term1) % MOD;
        
        // Calculate comb(h_adj + 3, 4) mod MOD
        let comb_h4 = 0;
        if (h_adj + 3 >= 4) {
            const numerator_h = (h_adj + 3) * (h_adj + 2) * (h_adj + 1) * h_adj;
            const denom_h = 24;
            comb_h4 = (numerator_h / denom_h) % MOD;
        }
        
        // Calculate comb(w_adj + 3, 4) mod MOD
        let comb_w4 = 0;
        if (w_adj + 3 >= 4) {
            const numerator_w = (w_adj + 3) * (w_adj + 2) * (w_adj + 1) * w_adj;
            const denom_w = 24;
            comb_w4 = (numerator_w / denom_w) % MOD;
        }
        
        const term2 = (comb_h4 * comb_w4) % MOD;
        const term2_final = (2 * term2) % MOD;
        
        let result = (term1_sq - term2_final) % MOD;
        if (result < 0) result += MOD;
        
        results.push(result);
    }
    
    console.log(results.join('\n'));
}

solve();

Python solution:
f=lambda x:x*x+x>>1
for i in[*open(0)][1:]:n,h,w=map(int,i.split());h-=n+n;w-=n+n;a=f(h+1)*f(w+1);print((h>-1<w)*(a*a-f(h)*f(w)*f(h+2)*f(w+2)//18)%998244353)

Ruby solution:
MOD = 998244353

def solve
  t = gets.to_i
  t.times do
    n, h, w = gets.split.map(&:to_i)
    h_prime = h - 2 * n
    w_prime = w - 2 * n
    if h_prime < 0 || w_prime < 0
      puts 0
      next
    end
    
    # Calculate comb(h_prime + 2, 2)
    c_h2 = (h_prime + 2) * (h_prime + 1) / 2 % MOD
    # Calculate comb(w_prime + 2, 2)
    c_w2 = (w_prime + 2) * (w_prime + 1) / 2 % MOD
    term1 = (c_h2 * c_w2) % MOD
    term1 = (term1 * term1) % MOD
    
    # Calculate comb(h_prime + 3, 4)
    if h_prime + 3 >= 4
      c_h4 = (h_prime + 3) * (h_prime + 2) * (h_prime + 1) * h_prime / 24 % MOD
    else
      c_h4 = 0
    end
    # Calculate comb(w_prime + 3, 4)
    if w_prime + 3 >= 4
      c_w4 = (w_prime + 3) * (w_prime + 2) * (w_prime + 1) * w_prime / 24 % MOD
    else
      c_w4 = 0
    end
    term2 = (2 * c_h4 * c_w4) % MOD
    
    result = (term1 - term2) % MOD
    result += MOD if result < 0
    puts result
  end
end

solve