How to generate huge dimension uniform distribution over small range of int in C++ fastly?

My problem is generating N dimensional array, N is huge, more than 10^8, each dimension lies at i.i.d uniform distribution over 1~q(or 0~q-1), where q is very small, q<=2^4. And array must be good statistically. So the basic soltuion is

constexpr auto N = 1000000000UZ;
constexpr auto q = 12;
std::array<std::uint8_t, N> Array{};
std::random_device Device{};
std::mt19937_64 Eng{Device()};
std::uniform_int_distribution<std::uint8_t> Dis(0, q);
std::ranges::generate(Array, [&]{return Dis(Eng);});

But the problem lies in performance, I have several plans to improve it:

1. because s=q^8<=2^32 so, we use

std::uniform_int_distribution<std::uint8_t> Dis(0, s);

add decompose the result t<q^8 into 8 different t_i<q, but this decomposition is not straightforward, and may have performance flaw in decomposition.

2. use boost::random::uniform_smallint, but I don't know how much improvement will be? this cann't be used together with method 1.

3. use multi-threading like openmp or <thread>, but C++ PRNG may not be thread-safe, so it's hard to write to my knowledge.

4. use other generator such as pcg32 or anything else, but these are not thread-safe as well.

Does anyone offer any suggestions?