I don't totally understand the math, but davmre is suggesting calling randn(), not rand() -- so sampling from a normal distribution instead of a uniform distribution.
But that pdf is also the joint pdf of N i.i.d. gaussians, evident by decomposing f=a * exp(-x1^2) * ... * exp(-xN^2) [2], which is the joint pdf x1,...xN s.t. fx1=c * exp(-x1^2), ..., fxn=c * exp(-xN^2).
[1] Since exp(-R^2) does not depend on direction but only on distance from the origin
[2] The fact that f(x1,...xN)=f1 * ... * fN if x1,...xN are independent follows directly from the fact that P(A & B) = P(A)*P(B) if A and B are independent events.