Use float constants directly; cosmetic changes; initialize largest
correctly; allow test(N) to set number of calls in the tests.
diff --git a/Lib/random.py b/Lib/random.py
index 1fa1377..0d2f06b 100644
--- a/Lib/random.py
+++ b/Lib/random.py
@@ -30,7 +30,7 @@
# -------------------- normal distribution --------------------
-NV_MAGICCONST = 4*exp(-0.5)/sqrt(2)
+NV_MAGICCONST = 4*exp(-0.5)/sqrt(2.0)
verify('NV_MAGICCONST', 1.71552776992141)
def normalvariate(mu, sigma):
# mu = mean, sigma = standard deviation
@@ -44,7 +44,7 @@
u1 = random()
u2 = random()
z = NV_MAGICCONST*(u1-0.5)/u2
- zz = z*z/4
+ zz = z*z/4.0
if zz <= -log(u2):
break
return mu+z*sigma
@@ -75,7 +75,7 @@
# -------------------- von Mises distribution --------------------
-TWOPI = 2*pi
+TWOPI = 2.0*pi
verify('TWOPI', 6.28318530718)
def vonmisesvariate(mu, kappa):
@@ -86,15 +86,15 @@
if kappa <= 1e-6:
return TWOPI * random()
- a = 1.0 + sqrt(1 + 4 * kappa * kappa)
- b = (a - sqrt(2 * a))/(2 * kappa)
- r = (1 + b * b)/(2 * b)
+ a = 1.0 + sqrt(1.0 + 4.0 * kappa * kappa)
+ b = (a - sqrt(2.0 * a))/(2.0 * kappa)
+ r = (1.0 + b * b)/(2.0 * b)
while 1:
u1 = random()
z = cos(pi * u1)
- f = (1 + r * z)/(r + z)
+ f = (1.0 + r * z)/(r + z)
c = kappa * (r - f)
u2 = random()
@@ -112,15 +112,15 @@
# -------------------- gamma distribution --------------------
-LOG4 = log(4)
+LOG4 = log(4.0)
verify('LOG4', 1.38629436111989)
def gammavariate(alpha, beta):
# beta times standard gamma
- ainv = sqrt(2 * alpha - 1)
+ ainv = sqrt(2.0 * alpha - 1.0)
return beta * stdgamma(alpha, ainv, alpha - LOG4, alpha + ainv)
-SG_MAGICCONST = 1+log(4.5)
+SG_MAGICCONST = 1.0 + log(4.5)
verify('SG_MAGICCONST', 2.50407739677627)
def stdgamma(alpha, ainv, bbb, ccc):
@@ -140,11 +140,11 @@
while 1:
u1 = random()
u2 = random()
- v = log(u1/(1-u1))/ainv
+ v = log(u1/(1.0-u1))/ainv
x = alpha*exp(v)
z = u1*u1*u2
r = bbb+ccc*v-x
- if r + SG_MAGICCONST - 4.5*z >= 0 or r >= log(z):
+ if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= log(z):
return x
elif alpha == 1.0:
@@ -176,17 +176,21 @@
# -------------------- Gauss (faster alternative) --------------------
-# When x and y are two variables from [0, 1), uniformly distributed, then
-#
-# cos(2*pi*x)*log(1-y)
-# sin(2*pi*x)*log(1-y)
-#
-# are two *independent* variables with normal distribution (mu = 0, sigma = 1).
-# (Lambert Meertens)
-
gauss_next = None
def gauss(mu, sigma):
+
+ # When x and y are two variables from [0, 1), uniformly
+ # distributed, then
+ #
+ # cos(2*pi*x)*log(1-y)
+ # sin(2*pi*x)*log(1-y)
+ #
+ # are two *independent* variables with normal distribution
+ # (mu = 0, sigma = 1).
+ # (Lambert Meertens)
+
global gauss_next
+
if gauss_next != None:
z = gauss_next
gauss_next = None
@@ -195,23 +199,30 @@
log1_y = log(1.0 - random())
z = cos(x2pi) * log1_y
gauss_next = sin(x2pi) * log1_y
+
return mu + z*sigma
# -------------------- beta --------------------
def betavariate(alpha, beta):
+
+ # Discrete Event Simulation in C, pp 87-88.
+
y = expovariate(alpha)
z = expovariate(1.0/beta)
return z/(y+z)
# -------------------- test program --------------------
-def test():
+def test(*args):
print 'TWOPI =', TWOPI
print 'LOG4 =', LOG4
print 'NV_MAGICCONST =', NV_MAGICCONST
print 'SG_MAGICCONST =', SG_MAGICCONST
N = 200
+ if args:
+ if args[1:]: print 'Excess test() arguments ignored'
+ N = args[0]
test_generator(N, 'random()')
test_generator(N, 'normalvariate(0.0, 1.0)')
test_generator(N, 'lognormvariate(0.0, 1.0)')
@@ -234,7 +245,7 @@
sum = 0.0
sqsum = 0.0
smallest = 1e10
- largest = 1e-10
+ largest = -1e10
t0 = time.time()
for i in range(n):
x = eval(code)