```go
package main

import (
	"fmt"
	"io"
	"os"
)

const MOD = 998244353
const MAX = 3000005

var phi [MAX]int
var isPrime [MAX]bool
var primes []int
var fact [MAX]int64
var invFact [MAX]int64

func power(base, exp int64) int64 {
	res := int64(1)
	base %= MOD
	for exp > 0 {
		if exp&1 == 1 {
			res = res * base % MOD
		}
		base = base * base % MOD
		exp >>= 1
	}
	return res
}

func precompute() {
	fact[0] = 1
	invFact[0] = 1
	for i := 1; i < MAX; i++ {
		fact[i] = fact[i-1] * int64(i) % MOD
	}
	invFact[MAX-1] = power(fact[MAX-1], MOD-2)
	for i := MAX - 2; i >= 1; i-- {
		invFact[i] = invFact[i+1] * int64(i+1) % MOD
	}

	for i := 2; i < MAX; i++ {
		isPrime[i] = true
	}
	phi[1] = 1
	for i := 2; i < MAX; i++ {
		if isPrime[i] {
			primes = append(primes, i)
			phi[i] = i - 1
		}
		for _, p := range primes {
			if i*p >= MAX {
				break
			}
			isPrime[i*p] = false
			if i%p == 0 {
				phi[i*p] = phi[i] * p
				break
			} else {
				phi[i*p] = phi[i] * (p - 1)
			}
		}
	}
}

func gcd(a, b int) int {
	for b != 0 {
		a, b = b, a%b
	}
	return a
}

func lcm(a, b int) int {
	if a == 0 || b == 0 {
		return 0
	}
	return (a / gcd(a, b)) * b
}

func getDivisors(n int) []int {
	var res []int
	for i := 1; i*i <= n; i++ {
		if n%i == 0 {
			res = append(res, i)
			if i*i != n {
				res = append(res, n/i)
			}
		}
	}
	return res
}

func main() {
	precompute()

	in, err := io.ReadAll(os.Stdin)
	if err != nil {
		return
	}
	pos := 0

	nextInt := func() int {
		for pos < len(in) && in[pos] <= 32 {
			pos++
		}
		if pos >= len(in) {
			return 0
		}
		res := 0
		for pos < len(in) && in[pos] > 32 {
			res = res*10 + int(in[pos]-'0')
			pos++
		}
		return res
	}

	t := nextInt()
	if t == 0 {
		return
	}

	for tc := 0; tc < t; tc++ {
		a := nextInt()
		b := nextInt()
		c := nextInt()
		k := nextInt()

		d := make([]int, k)
		g := 0
		for i := 0; i < k; i++ {
			d[i] = nextInt()
			if i == 0 {
				g = d[i]
			} else {
				g = gcd(g, d[i])
			}
		}

		divA := getDivisors(a)
		divB := getDivisors(b)
		divC := getDivisors(c)

		weights := make(map[int]int64)
		for _, A := range divA {
			for _, B := range divB {
				lcmAB := lcm(A, B)
				for _, C := range divC {
					m := lcm(lcmAB, C)
					if g%m == 0 {
						phim := int64(phi[A]) * int64(phi[B]) % MOD
						phim = phim * int64(phi[C]) % MOD
						weights[m] = (weights[m] + phim) % MOD
					}
				}
			}
		}

		N := a * b * c
		ans := int64(0)

		for m, w := range weights {
			ways := fact[N/m]
			for j := 0; j < k; j++ {
				ways = ways * invFact[d[j]/m] % MOD
			}
			ans = (ans + w*ways) % MOD
		}

		ans = ans * power(int64(N), MOD-2) % MOD
		fmt.Println(ans)
	}
}
```