A simple way:
It is unlikely in this sort of game that you will need 3 Gaussian variates just once.
You need some store variable that can contains either a triplet of Gaussian variates or nothing (Null, Nothing, Empty, whatever that is in your programming language, you didn't tell us which one).
Initially, the store contains nothing (empty).
When asked for a triplet:
if the store contains a triplet, just return that triplet. And mark the store as empty.
if the store is empty, run Box-Muller 3 times. That gives you 2 triplets. Put the second triplet in the store. Return the first triplet.
An alternative way for the mathematically inclined programmer:
If one just tries to adapt Box-Muller to 3 dimensions, the sole tricky part is to get the norm of the random 3D vector. The rest is about the 2 spherical angles θ (theta) and φ (phi), which is easy stuff.
It turns out that in 3 dimensions, that norm involves the inverse of the incomplete gamma function.
And if you have Python and Numpy/Scipy, this is function scipy.special.gammaincinv.
We can thus write this code:
import math
import numpy.random as rd
import scipy.special as sp
# convert 3 uniform [0,1) variates into 3 unit Gaussian variates:
def boxMuller3d(u3):
u0,u1,u2 = u3 # 3 uniform random numbers in [0,1)
gamma = u0
norm2 = 2.0 * sp.gammaincinv(1.5, gamma) # "regularized" versions
norm = math.sqrt(norm2)
zr = (2.0 * u1) - 1.0 # sin(theta)
hr = math.sqrt(1.0 - zr*zr) # cos(theta)
phi = 2.0 * math.pi * u2
xr = hr * math.cos(phi)
yr = hr * math.sin(phi)
g3 = list(map(lambda c: c*norm, [xr, yr, zr]))
return g3
# generate 3 uniform variates and convert them into 3 unit Gaussian variates:
def gauss3(rng):
u3 = rng.uniform(0.0, 1.0, 3)
g3 = boxMuller3d(u3)
return g3
To (partly) check correctness, we can have this small main program, which displays the statistical moments of order 1 to 4 of the resulting random serie:
randomSeed = 42
rng = rd.default_rng(randomSeed)
count = 3000000 # (X,Y,Z) triplet count
variates = []
for i in range(count):
g3 = gauss3(rng)
variates += g3
ln = len(variates)
print("length=%d\n" % ln)
# Checking statistical moments of order 1 to 4:
m1 = sum(variates) / ln
m2 = sum( map(lambda x: x*x, variates) ) / ln
m3 = sum( map(lambda x: x**3, variates) ) / ln
m4 = sum( map(lambda x: x**4, variates) ) / ln
print("m1=%g m2=%g m3=%g m4=%g\n" % (m1,m2,m3,m4))
Test program output:
length=9000000
m1=-0.000455911 m2=1.00025 m3=-0.000563454 m4=3.00184
We thus can see that these moments are reasonably close to their mathematically expected values, respectively 0,1,0,3.