this is my attempt of the Wichmann-Hill Pseudo Random Number Generator, if you can find any improvement please help:
Public Static Function RndX() As Double
' Wichmann, B.A. and I.D. Hill, Algorithm AS 183: An Efficient and Portable
' Pseudo-Random Number Generator, Applied Statistics, 31, 188-190, 1982.
' C IX, IY, IZ SHOULD BE SET TO INTEGER VALUES BETWEEN 1 AND 30000 BEFORE FIRST ENTRY
' IX = MOD(171 * IX, 30269)
' IY = MOD(172 * IY, 30307)
' IZ = MOD(170 * IZ, 30323)
' RANDOM = AMOD(FLOAT(IX) / 30269# + FLOAT(IY) / 30307# + FLOAT(IZ) / 30323#, 1#)

' Note: the FORTRAN function AMOD(A,P) = A - (INT(A / P) * P)
' VBA Mod function rounds non-integer values of A,P FIRST, so it won't work in the last line!

Dim ix As Long, iy As Long, iz As Long
Dim F_ix As Single, F_iy As Single, F_iz As Single
Dim sum As Single
Dim seeded As Boolean ' initialized to 'False' during first call to function

If seeded = False Then
ix = 123 ' use a fixed initial seed
iy = 234
iz = 345
seeded = True ' static variable set true so we don't 're-seed' every time
End If

' the main part of the linear congruential algorithm is only three lines long
ix = (171 * ix) Mod 30269
iy = (172 * iy) Mod 30307
iz = (170 * iz) Mod 30323

'CSng convert the integer values to single precision reals for arithmetic
F_ix = CSng(ix) / 30269! ' use ! to ensure denominator is also single precision real
F_iy = CSng(iy) / 30307!
F_iz = CSng(iz) / 30323!
sum = F_ix + F_iy + F_iz

' take the fractional part of the sum as the random number
RndX = sum - Int(sum)
End Function

Then i use the random number to use in the box muller fuction i attempted to make:

Public Function RndBM()

' User function to transform uniform random numbers from randBM (x,y)
' into standard normal variates (Z_1, Z_2) using the Box-Muller transformation
' Box-Muller Ann. Math. Stat. 29, 610-611 (1958)

Dim x As Double, y As Double, s As Double, Z_1 As Double
Const TWO_PI As Double = 2# * 3.14159265358979

x = RndX() ' generate two uniform random numbers, x and y
y = RndX()

s = Sqr(-2# * Log(x)) ' transform to normal random numbers, r1 and r2
Z_1 = s * Cos(TWO_PI * y)

RndBM = Z_1 ' return the cosine term
End Function