Define the new Rand and Source types to allow creating

isolated sources of random values.

Add normal and exponential distributions.

Add some tests for the normal and exponential distributions.

R=rsc
APPROVED=rsc
DELTA=1005  (904 added, 80 deleted, 21 changed)
OCL=35501
CL=35779

src/pkg/rand/Makefile

@@ -6,6 +6,9 @@ include $(GOROOT)/src/Make.$(GOARCH)
 
 TARG=rand
 GOFILES=\
+	exp.go\
+	normal.go\
+	rand.go\
 	rng.go\
 
 include $(GOROOT)/src/Make.pkg

src/pkg/rand/normal.go

+// Copyright 2009 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package rand
+
+import (
+	"math";
+)
+
+/*
+ * Normal distribution
+ *
+ * See "The Ziggurat Method for Generating Random Variables"
+ * (Marsaglia & Tsang, 2000)
+ * http://www.jstatsoft.org/v05/i08/paper [pdf]
+ */
+
+const (
+	rn = 3.442619855899;
+)
+
+func absInt32(i int32) uint32 {
+	if i < 0 {
+		return uint32(-i);
+	}
+	return uint32(i);
+}
+
+// NormFloat64 returns a normally distributed float64 in the range
+// [-math.MaxFloat64, +math.MaxFloat64] with
+// standard normal distribution (mean = 0, stddev = 1).
+// To produce a different normal distribution, callers can
+// adjust the output using:
+//
+//  sample = NormFloat64() * desiredStdDev + desiredMean
+//
+func (r *Rand) NormFloat64() float64 {
+	for {
+		j := int32(r.Uint32());	// Possibly negative
+		i := j&0x7F;
+		x := float64(j)*float64(wn[i]);
+		if absInt32(j) < kn[i] {
+			// This case should be hit better than 99% of the time.
+			return x;
+		}
+
+		if i == 0 {
+			// This extra work is only required for the base strip.
+			for {
+				x = -math.Log(r.Float64()) * (1.0/rn);
+				y := -math.Log(r.Float64());
+				if y+y >= x*x {
+					break;
+				}
+			}
+			if j > 0 {
+				return rn+x;
+			}
+			return -rn - x;
+		}
+		if fn[i] + float32(r.Float64())*(fn[i-1]-fn[i]) < float32(math.Exp(-.5 * x * x)) {
+			return x;
+		}
+	}
+	panic("unreachable");
+}
+
+var kn = [128]uint32{
+	0x76ad2212, 0x0, 0x600f1b53, 0x6ce447a6, 0x725b46a2,
+	0x7560051d, 0x774921eb, 0x789a25bd, 0x799045c3, 0x7a4bce5d,
+	0x7adf629f, 0x7b5682a6, 0x7bb8a8c6, 0x7c0ae722, 0x7c50cce7,
+	0x7c8cec5b, 0x7cc12cd6, 0x7ceefed2, 0x7d177e0b, 0x7d3b8883,
+	0x7d5bce6c, 0x7d78dd64, 0x7d932886, 0x7dab0e57, 0x7dc0dd30,
+	0x7dd4d688, 0x7de73185, 0x7df81cea, 0x7e07c0a3, 0x7e163efa,
+	0x7e23b587, 0x7e303dfd, 0x7e3beec2, 0x7e46db77, 0x7e51155d,
+	0x7e5aabb3, 0x7e63abf7, 0x7e6c222c, 0x7e741906, 0x7e7b9a18,
+	0x7e82adfa, 0x7e895c63, 0x7e8fac4b, 0x7e95a3fb, 0x7e9b4924,
+	0x7ea0a0ef, 0x7ea5b00d, 0x7eaa7ac3, 0x7eaf04f3, 0x7eb3522a,
+	0x7eb765a5, 0x7ebb4259, 0x7ebeeafd, 0x7ec2620a, 0x7ec5a9c4,
+	0x7ec8c441, 0x7ecbb365, 0x7ece78ed, 0x7ed11671, 0x7ed38d62,
+	0x7ed5df12, 0x7ed80cb4, 0x7eda175c, 0x7edc0005, 0x7eddc78e,
+	0x7edf6ebf, 0x7ee0f647, 0x7ee25ebe, 0x7ee3a8a9, 0x7ee4d473,
+	0x7ee5e276, 0x7ee6d2f5, 0x7ee7a620, 0x7ee85c10, 0x7ee8f4cd,
+	0x7ee97047, 0x7ee9ce59, 0x7eea0eca, 0x7eea3147, 0x7eea3568,
+	0x7eea1aab, 0x7ee9e071, 0x7ee98602, 0x7ee90a88, 0x7ee86d08,
+	0x7ee7ac6a, 0x7ee6c769, 0x7ee5bc9c, 0x7ee48a67, 0x7ee32efc,
+	0x7ee1a857, 0x7edff42f, 0x7ede0ffa, 0x7edbf8d9, 0x7ed9ab94,
+	0x7ed7248d, 0x7ed45fae, 0x7ed1585c, 0x7ece095f, 0x7eca6ccb,
+	0x7ec67be2, 0x7ec22eee, 0x7ebd7d1a, 0x7eb85c35, 0x7eb2c075,
+	0x7eac9c20, 0x7ea5df27, 0x7e9e769f, 0x7e964c16, 0x7e8d44ba,
+	0x7e834033, 0x7e781728, 0x7e6b9933, 0x7e5d8a1a, 0x7e4d9ded,
+	0x7e3b737a, 0x7e268c2f, 0x7e0e3ff5, 0x7df1aa5d, 0x7dcf8c72,
+	0x7da61a1e, 0x7d72a0fb, 0x7d30e097, 0x7cd9b4ab, 0x7c600f1a,
+	0x7ba90bdc, 0x7a722176, 0x77d664e5,
+}
+var wn = [128]float32{
+	1.7290405e-09, 1.2680929e-10, 1.6897518e-10, 1.9862688e-10,
+	2.2232431e-10, 2.4244937e-10, 2.601613e-10, 2.7611988e-10,
+	2.9073963e-10, 3.042997e-10, 3.1699796e-10, 3.289802e-10,
+	3.4035738e-10, 3.5121603e-10, 3.616251e-10, 3.7164058e-10,
+	3.8130857e-10, 3.9066758e-10, 3.9975012e-10, 4.08584e-10,
+	4.1719309e-10, 4.2559822e-10, 4.338176e-10, 4.418672e-10,
+	4.497613e-10, 4.5751258e-10, 4.651324e-10, 4.7263105e-10,
+	4.8001775e-10, 4.87301e-10, 4.944885e-10, 5.015873e-10,
+	5.0860405e-10, 5.155446e-10, 5.2241467e-10, 5.2921934e-10,
+	5.359635e-10, 5.426517e-10, 5.4928817e-10, 5.5587696e-10,
+	5.624219e-10, 5.6892646e-10, 5.753941e-10, 5.818282e-10,
+	5.882317e-10, 5.946077e-10, 6.00959e-10, 6.072884e-10,
+	6.135985e-10, 6.19892e-10, 6.2617134e-10, 6.3243905e-10,
+	6.386974e-10, 6.449488e-10, 6.511956e-10, 6.5744005e-10,
+	6.6368433e-10, 6.699307e-10, 6.7618144e-10, 6.824387e-10,
+	6.8870465e-10, 6.949815e-10, 7.012715e-10, 7.075768e-10,
+	7.1389966e-10, 7.202424e-10, 7.266073e-10, 7.329966e-10,
+	7.394128e-10, 7.4585826e-10, 7.5233547e-10, 7.58847e-10,
+	7.653954e-10, 7.719835e-10, 7.7861395e-10, 7.852897e-10,
+	7.920138e-10, 7.987892e-10, 8.0561924e-10, 8.125073e-10,
+	8.194569e-10, 8.2647167e-10, 8.3355556e-10, 8.407127e-10,
+	8.479473e-10, 8.55264e-10, 8.6266755e-10, 8.7016316e-10,
+	8.777562e-10, 8.8545243e-10, 8.932582e-10, 9.0117996e-10,
+	9.09225e-10, 9.174008e-10, 9.2571584e-10, 9.341788e-10,
+	9.427997e-10, 9.515889e-10, 9.605579e-10, 9.697193e-10,
+	9.790869e-10, 9.88676e-10, 9.985036e-10, 1.0085882e-09,
+	1.0189509e-09, 1.0296151e-09, 1.0406069e-09, 1.0519566e-09,
+	1.063698e-09, 1.0758702e-09, 1.0885183e-09, 1.1016947e-09,
+	1.1154611e-09, 1.1298902e-09, 1.1450696e-09, 1.1611052e-09,
+	1.1781276e-09, 1.1962995e-09, 1.2158287e-09, 1.2369856e-09,
+	1.2601323e-09, 1.2857697e-09, 1.3146202e-09, 1.347784e-09,
+	1.3870636e-09, 1.4357403e-09, 1.5008659e-09, 1.6030948e-09,
+}
+var fn = [128]float32{
+	1, 0.9635997, 0.9362827, 0.9130436, 0.89228165, 0.87324303,
+	0.8555006, 0.8387836, 0.8229072, 0.8077383, 0.793177,
+	0.7791461, 0.7655842, 0.7524416, 0.73967725, 0.7272569,
+	0.7151515, 0.7033361, 0.69178915, 0.68049186, 0.6694277,
+	0.658582, 0.6479418, 0.63749546, 0.6272325, 0.6171434,
+	0.6072195, 0.5974532, 0.58783704, 0.5783647, 0.56903,
+	0.5598274, 0.5507518, 0.54179835, 0.5329627, 0.52424055,
+	0.5156282, 0.50712204, 0.49871865, 0.49041483, 0.48220766,
+	0.4740943, 0.46607214, 0.4581387, 0.45029163, 0.44252872,
+	0.43484783, 0.427247, 0.41972435, 0.41227803, 0.40490642,
+	0.39760786, 0.3903808, 0.3832238, 0.37613547, 0.36911446,
+	0.3621595, 0.35526937, 0.34844297, 0.34167916, 0.33497685,
+	0.3283351, 0.3217529, 0.3152294, 0.30876362, 0.30235484,
+	0.29600215, 0.28970486, 0.2834622, 0.2772735, 0.27113807,
+	0.2650553, 0.25902456, 0.2530453, 0.24711695, 0.241239,
+	0.23541094, 0.22963232, 0.2239027, 0.21822165, 0.21258877,
+	0.20700371, 0.20146611, 0.19597565, 0.19053204, 0.18513499,
+	0.17978427, 0.17447963, 0.1692209, 0.16400786, 0.15884037,
+	0.15371831, 0.14864157, 0.14361008, 0.13862377, 0.13368265,
+	0.12878671, 0.12393598, 0.119130544, 0.11437051, 0.10965602,
+	0.104987256, 0.10036444, 0.095787846, 0.0912578, 0.08677467,
+	0.0823389, 0.077950984, 0.073611505, 0.06932112, 0.06508058,
+	0.06089077, 0.056752663, 0.0526674, 0.048636295, 0.044660863,
+	0.040742867, 0.03688439, 0.033087887, 0.029356318,
+	0.025693292, 0.022103304, 0.018592102, 0.015167298,
+	0.011839478, 0.008624485, 0.005548995, 0.0026696292,
+}

src/pkg/rand/rand.go

+// Copyright 2009 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// Package rand implements pseudo-random number generators.
+package rand
+
+// A Source represents a source of uniformly-distributed
+// pseudo-random int64 values in the range [0, 1<<63).
+type Source interface {
+	Int63() int64;
+	Seed(seed int64);
+}
+
+// NewSource returns a new pseudo-random Source seeded with the given value.
+func NewSource(seed int64) Source {
+	var rng rngSource;
+	rng.Seed(seed);
+	return &rng;
+}
+
+// A Rand is a source of random numbers.
+type Rand struct {
+	src Source;
+}
+
+// New returns a new Rand that uses random values from src
+// to generate other random values.
+func New(src Source) *Rand {
+	return &Rand{src};
+}
+
+// Seed uses the provided seed value to initialize the generator to a deterministic state.
+func (r *Rand) Seed(seed int64) {
+	r.src.Seed(seed);
+}
+
+// Int63 returns a non-negative pseudo-random 63-bit integer as an int64.
+func (r *Rand) Int63() int64 {
+	return r.src.Int63();
+}
+
+// Uint32 returns a pseudo-random 32-bit value as a uint32.
+func (r *Rand) Uint32() uint32 {
+	return uint32(r.Int63() >> 31);
+}
+
+// Int31 returns a non-negative pseudo-random 31-bit integer as an int32.
+func (r *Rand) Int31() int32 {
+	return int32(r.Int63() >> 32);
+}
+
+// Int returns a non-negative pseudo-random int.
+func (r *Rand) Int() int {
+	u := uint(r.Int63());
+	return int(u<<1>>1);	// clear sign bit if int == int32
+}
+
+// Int63n returns, as an int64, a non-negative pseudo-random number in [0,n).
+func (r *Rand) Int63n(n int64) int64 {
+	if n <= 0 {
+		return 0;
+	}
+	max := int64((1<<63) - 1 - (1<<63)%uint64(n));
+	v := r.Int63();
+	for v > max {
+		v = r.Int63();
+	}
+	return v%n;
+}
+
+// Int31n returns, as an int32, a non-negative pseudo-random number in [0,n).
+func (r *Rand) Int31n(n int32) int32 {
+	return int32(r.Int63n(int64(n)));
+}
+
+// Intn returns, as an int, a non-negative pseudo-random number in [0,n).
+func (r *Rand) Intn(n int) int {
+	return int(r.Int63n(int64(n)));
+}
+
+// Float64 returns, as a float64, a pseudo-random number in [0.0,1.0).
+func (r *Rand) Float64() float64 {
+	return float64(r.Int63())/(1<<63);
+}
+
+// Float32 returns, as a float32, a pseudo-random number in [0.0,1.0).
+func (r *Rand) Float32() float32 {
+	return float32(r.Float64());
+}
+
+// Float returns, as a float, a pseudo-random number in [0.0,1.0).
+func (r *Rand) Float() float {
+	return float(r.Float64());
+}
+
+// Perm returns, as a slice of n ints, a pseudo-random permutation of the integers [0,n).
+func (r *Rand) Perm(n int) []int {
+	m := make([]int, n);
+	for i := 0; i < n; i++ {
+		m[i] = i;
+	}
+	for i := 0; i < n; i++ {
+		j := r.Intn(i+1);
+		m[i], m[j] = m[j], m[i];
+	}
+	return m;
+}
+
+/*
+ * Top-level convenience functions
+ */
+
+var globalRand = New(NewSource(1))
+
+// Seed uses the provided seed value to initialize the generator to a deterministic state.
+func Seed(seed int64) {
+	globalRand.Seed(seed);
+}
+
+// Int63 returns a non-negative pseudo-random 63-bit integer as an int64.
+func Int63() int64 {
+	return globalRand.Int63();
+}
+
+// Uint32 returns a pseudo-random 32-bit value as a uint32.
+func Uint32() uint32 {
+	return globalRand.Uint32();
+}
+
+// Int31 returns a non-negative pseudo-random 31-bit integer as an int32.
+func Int31() int32 {
+	return globalRand.Int31();
+}
+
+// Int returns a non-negative pseudo-random int.
+func Int() int {
+	return globalRand.Int();
+}
+
+// Int63n returns, as an int64, a non-negative pseudo-random number in [0,n).
+func Int63n(n int64) int64 {
+	return globalRand.Int63n(n);
+}
+
+// Int31n returns, as an int32, a non-negative pseudo-random number in [0,n).
+func Int31n(n int32) int32 {
+	return globalRand.Int31n(n);
+}
+
+// Intn returns, as an int, a non-negative pseudo-random number in [0,n).
+func Intn(n int) int {
+	return globalRand.Intn(n);
+}
+
+// Float64 returns, as a float64, a pseudo-random number in [0.0,1.0).
+func Float64() float64 {
+	return globalRand.Float64();
+}
+
+// Float32 returns, as a float32, a pseudo-random number in [0.0,1.0).
+func Float32() float32 {
+	return globalRand.Float32();
+}
+
+// Float returns, as a float, a pseudo-random number in [0.0,1.0).
+func Float() float {
+	return globalRand.Float();
+}
+
+// Perm returns, as a slice of n ints, a pseudo-random permutation of the integers [0,n).
+func Perm(n int) []int {
+	return globalRand.Perm(n);
+}
+
+// NormFloat64 returns a normally distributed float64 in the range
+// [-math.MaxFloat64, +math.MaxFloat64] with
+// standard normal distribution (mean = 0, stddev = 1).
+// To produce a different normal distribution, callers can
+// adjust the output using:
+//
+//  sample = NormFloat64() * desiredStdDev + desiredMean
+//
+func NormFloat64() float64 {
+	return globalRand.NormFloat64();
+}
+
+// ExpFloat64 returns an exponentially distributed float64 in the range
+// (0, +math.MaxFloat64] with an exponential distribution whose rate parameter
+// (lambda) is 1 and whose mean is 1/lambda (1).
+// To produce a distribution with a different rate parameter,
+// callers can adjust the output using:
+//
+//  sample = ExpFloat64() / desiredRateParameter
+//
+func ExpFloat64() float64 {
+	return globalRand.ExpFloat64();
+}

src/pkg/rand/rand_test.go

+// Copyright 2009 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package rand
+
+import (
+	"math";
+	"fmt";
+	"os";
+	"testing";
+)
+
+const (
+	numTestSamples = 10000;
+)
+
+type statsResults struct {
+	mean		float64;
+	stddev		float64;
+	closeEnough	float64;
+	maxError	float64;
+}
+
+func max(a, b float64) float64 {
+	if a > b {
+		return a;
+	}
+	return b;
+}
+
+func nearEqual(a, b, closeEnough, maxError float64) bool {
+	absDiff := math.Fabs(a-b);
+	if absDiff < closeEnough {	// Necessary when one value is zero and one value is close to zero.
+		return true;
+	}
+	return absDiff / max(math.Fabs(a), math.Fabs(b)) < maxError;
+}
+
+var testSeeds = []int64{1, 1754801282, 1698661970, 1550503961}
+
+// checkSimilarDistribution returns success if the mean and stddev of the
+// two statsResults are similar.
+func (this *statsResults) checkSimilarDistribution(expected *statsResults) os.Error {
+	if !nearEqual(this.mean, expected.mean, expected.closeEnough, expected.maxError) {
+		s := fmt.Sprintf("mean %v != %v (allowed error %v, %v)", this.mean, expected.mean, expected.closeEnough, expected.maxError);
+		fmt.Println(s);
+		return os.ErrorString(s);
+	}
+	if !nearEqual(this.stddev, expected.stddev, 0, expected.maxError) {
+		s := fmt.Sprintf("stddev %v != %v (allowed error %v, %v)", this.stddev, expected.stddev, expected.closeEnough, expected.maxError);
+		fmt.Println(s);
+		return os.ErrorString(s);
+	}
+	return nil;
+}
+
+func getStatsResults(samples []float64) *statsResults {
+	res := new(statsResults);
+	var sum float64;
+	for i := range samples {
+		sum += samples[i];
+	}
+	res.mean = sum/float64(len(samples));
+	var devsum float64;
+	for i := range samples {
+		devsum += math.Pow(samples[i] - res.mean, 2);
+	}
+	res.stddev = math.Sqrt(devsum/float64(len(samples)));
+	return res;
+}
+
+func checkSampleDistribution(t *testing.T, samples []float64, expected *statsResults) {
+	actual := getStatsResults(samples);
+	err := actual.checkSimilarDistribution(expected);
+	if err != nil {
+		t.Errorf(err.String());
+	}
+}
+
+func checkSampleSliceDistributions(t *testing.T, samples []float64, nslices int, expected *statsResults) {
+	chunk := len(samples)/nslices;
+	for i := 0; i < nslices; i++ {
+		low := i*chunk;
+		var high int;
+		if i == nslices-1 {
+			high = len(samples)-1;
+		} else {
+			high = (i+1)*chunk;
+		}
+		checkSampleDistribution(t, samples[low:high], expected);
+	}
+}
+
+//
+// Normal distribution tests
+//
+
+func generateNormalSamples(nsamples int, mean, stddev float64, seed int64) []float64 {
+	r := New(NewSource(seed));
+	samples := make([]float64, nsamples);
+	for i := range samples {
+		samples[i] = r.NormFloat64() * stddev + mean;
+	}
+	return samples;
+}
+
+func testNormalDistribution(t *testing.T, nsamples int, mean, stddev float64, seed int64) {
+	//fmt.Printf("testing nsamples=%v mean=%v stddev=%v seed=%v\n", nsamples, mean, stddev, seed);
+
+	samples := generateNormalSamples(nsamples, mean, stddev, seed);
+	errorScale := max(1.0, stddev);	// Error scales with stddev
+	expected := &statsResults{mean, stddev, 0.10 * errorScale, 0.08 * errorScale};
+
+	// Make sure that the entire set matches the expected distribution.
+	checkSampleDistribution(t, samples, expected);
+
+	// Make sure that each half of the set matches the expected distribution.
+	checkSampleSliceDistributions(t, samples, 2, expected);
+
+	// Make sure that each 7th of the set matches the expected distribution.
+	checkSampleSliceDistributions(t, samples, 7, expected);
+}
+
+// Actual tests
+
+func TestStandardNormalValues(t *testing.T) {
+	for _, seed := range testSeeds {
+		testNormalDistribution(t, numTestSamples, 0, 1, seed);
+	}
+}
+
+func TestNonStandardNormalValues(t *testing.T) {
+	for sd := float64(0.5); sd < 1000; sd *= 2 {
+		for m := float64(0.5); m < 1000; m *= 2 {
+			for _, seed := range testSeeds {
+				testNormalDistribution(t, numTestSamples, m, sd, seed);
+			}
+		}
+	}
+}
+
+//
+// Exponential distribution tests
+//
+
+func generateExponentialSamples(nsamples int, rate float64, seed int64) []float64 {
+	r := New(NewSource(seed));
+	samples := make([]float64, nsamples);
+	for i := range samples {
+		samples[i] = r.ExpFloat64() / rate;
+	}
+	return samples;
+}
+
+func testExponentialDistribution(t *testing.T, nsamples int, rate float64, seed int64) {
+	//fmt.Printf("testing nsamples=%v rate=%v seed=%v\n", nsamples, rate, seed);
+
+	mean := 1/rate;
+	stddev := mean;
+
+	samples := generateExponentialSamples(nsamples, rate, seed);
+	errorScale := max(1.0, 1/rate);	// Error scales with the inverse of the rate
+	expected := &statsResults{mean, stddev, 0.10 * errorScale, 0.20 * errorScale};
+
+	// Make sure that the entire set matches the expected distribution.
+	checkSampleDistribution(t, samples, expected);
+
+	// Make sure that each half of the set matches the expected distribution.
+	checkSampleSliceDistributions(t, samples, 2, expected);
+
+	// Make sure that each 7th of the set matches the expected distribution.
+	checkSampleSliceDistributions(t, samples, 7, expected);
+}
+
+// Actual tests
+
+func TestStandardExponentialValues(t *testing.T) {
+	for _, seed := range testSeeds {
+		testExponentialDistribution(t, numTestSamples, 1, seed);
+	}
+}
+
+func TestNonStandardExponentialValues(t *testing.T) {
+	for rate := float64(0.05); rate < 10; rate *= 2 {
+		for _, seed := range testSeeds {
+			testExponentialDistribution(t, numTestSamples, rate, seed);
+		}
+	}
+}
+
+//
+// Table generation tests
+//
+
+func initNorm() (testKn []uint32, testWn, testFn []float32) {
+	const m1 = 1<<31;
+	var (
+		dn	float64	= rn;
+		tn		= dn;
+		vn	float64	= 9.91256303526217e-3;
+	)
+
+	testKn = make([]uint32, 128);
+	testWn = make([]float32, 128);
+	testFn = make([]float32, 128);
+
+	q := vn / math.Exp(-0.5 * dn * dn);
+	testKn[0] = uint32((dn/q)*m1);
+	testKn[1] = 0;
+	testWn[0] = float32(q/m1);
+	testWn[127] = float32(dn/m1);
+	testFn[0] = 1.0;
+	testFn[127] = float32(math.Exp(-0.5 * dn * dn));
+	for i := 126; i >= 1; i-- {
+		dn = math.Sqrt(-2.0 * math.Log(vn/dn + math.Exp(-0.5 * dn * dn)));
+		testKn[i+1] = uint32((dn/tn)*m1);
+		tn = dn;
+		testFn[i] = float32(math.Exp(-0.5 * dn * dn));
+		testWn[i] = float32(dn/m1);
+	}
+	return;
+}
+
+func initExp() (testKe []uint32, testWe, testFe []float32) {
+	const m2 = 1<<32;
+	var (
+		de	float64	= re;
+		te		= de;
+		ve	float64	= 3.9496598225815571993e-3;
+	)
+
+	testKe = make([]uint32, 256);
+	testWe = make([]float32, 256);
+	testFe = make([]float32, 256);
+
+	q := ve / math.Exp(-de);
+	testKe[0] = uint32((de/q)*m2);
+	testKe[1] = 0;
+	testWe[0] = float32(q/m2);
+	testWe[255] = float32(de/m2);
+	testFe[0] = 1.0;
+	testFe[255] = float32(math.Exp(-de));
+	for i := 254; i >= 1; i-- {
+		de = -math.Log(ve/de + math.Exp(-de));
+		testKe[i+1] = uint32((de/te)*m2);
+		te = de;
+		testFe[i] = float32(math.Exp(-de));
+		testWe[i] = float32(de/m2);
+	}
+	return;
+}
+
+// compareUint32Slices returns the first index where the two slices
+// disagree, or <0 if the lengths are the same and all elements
+// are identical.
+func compareUint32Slices(s1, s2 []uint32) int {
+	if len(s1) != len(s2) {
+		if len(s1) > len(s2) {
+			return len(s2)+1;
+		}
+		return len(s1)+1;
+	}
+	for i := range s1 {
+		if s1[i] != s2[i] {
+			return i;
+		}
+	}
+	return -1;
+}
+
+// compareFloat32Slices returns the first index where the two slices
+// disagree, or <0 if the lengths are the same and all elements
+// are identical.
+func compareFloat32Slices(s1, s2 []float32) int {
+	if len(s1) != len(s2) {
+		if len(s1) > len(s2) {
+			return len(s2)+1;
+		}
+		return len(s1)+1;
+	}
+	for i := range s1 {
+		if !nearEqual(float64(s1[i]), float64(s2[i]), 0, 1e-7) {
+			return i;
+		}
+	}
+	return -1;
+}
+
+func TestNormTables(t *testing.T) {
+	testKn, testWn, testFn := initNorm();
+	if i := compareUint32Slices(kn[0:len(kn)], testKn); i >= 0 {
+		t.Errorf("kn disagrees at index %v; %v != %v\n", i, kn[i], testKn[i]);
+	}
+	if i := compareFloat32Slices(wn[0:len(wn)], testWn); i >= 0 {
+		t.Errorf("wn disagrees at index %v; %v != %v\n", i, wn[i], testWn[i]);
+	}
+	if i := compareFloat32Slices(fn[0:len(fn)], testFn); i >= 0 {
+		t.Errorf("fn disagrees at index %v; %v != %v\n", i, fn[i], testFn[i]);
+	}
+}
+
+func TestExpTables(t *testing.T) {
+	testKe, testWe, testFe := initExp();
+	if i := compareUint32Slices(ke[0:len(ke)], testKe); i >= 0 {
+		t.Errorf("ke disagrees at index %v; %v != %v\n", i, ke[i], testKe[i]);
+	}
+	if i := compareFloat32Slices(we[0:len(we)], testWe); i >= 0 {
+		t.Errorf("we disagrees at index %v; %v != %v\n", i, we[i], testWe[i]);
+	}
+	if i := compareFloat32Slices(fe[0:len(fe)], testFe); i >= 0 {
+		t.Errorf("fe disagrees at index %v; %v != %v\n", i, fe[i], testFe[i]);
+	}
+}

src/pkg/rand/rng.go

@@ -2,10 +2,11 @@
 // Use of this source code is governed by a BSD-style
 // license that can be found in the LICENSE file.
 
-// Uniformly distributed pseudo-random numbers.
 package rand
 
 /*
+ * Uniform distribution
+ *
  * algorithm by
  * DP Mitchell and JA Reeds
  */
@@ -22,10 +23,6 @@ const (
 )
 
 var (
-	rng_tap		int;		// index into vector
-	rng_feed	int;		// index into vector
-	rng_vec		[_LEN]int64;	// current feedback register
-
 	// cooked random numbers
 	// the state of the rng
 	// after 780e10 iterations
@@ -185,6 +182,12 @@ var (
 	};
 )
 
+type rngSource struct {
+	tap	int;		// index into vec
+	feed	int;		// index into vec
+	vec	[_LEN]int64;	// current feedback register
+}
+
 // seed rng x[n+1] = 48271 * x[n] mod (2**31 - 1)
 func seedrand(x int32) int32 {
 	hi := x/_Q;
@@ -197,9 +200,9 @@ func seedrand(x int32) int32 {
 }
 
 // Seed uses the provided seed value to initialize the generator to a deterministic state.
-func Seed(seed int32) {
-	rng_tap = 0;
-	rng_feed = _LEN-_TAP;
+func (rng *rngSource) Seed(seed int64) {
+	rng.tap = 0;
+	rng.feed = _LEN-_TAP;
 
 	seed = seed%_M;
 	if seed < 0 {
@@ -209,7 +212,7 @@ func Seed(seed int32) {
 		seed = 89482311;
 	}
 
-	x := seed;
+	x := int32(seed);
 	for i := -20; i < _LEN; i++ {
 		x = seedrand(x);
 		if i >= 0 {
@@ -220,99 +223,24 @@ func Seed(seed int32) {
 			x = seedrand(x);
 			u ^= int64(x);
 			u ^= rng_cooked[i];
-			rng_vec[i] = u&_MASK;
+			rng.vec[i] = u&_MASK;
 		}
 	}
 }
 
 // Int63 returns a non-negative pseudo-random 63-bit integer as an int64.
-func Int63() int64 {
-	rng_tap--;
-	if rng_tap < 0 {
-		rng_tap += _LEN;
+func (rng *rngSource) Int63() int64 {
+	rng.tap--;
+	if rng.tap < 0 {
+		rng.tap += _LEN;
 	}
 
-	rng_feed--;
-	if rng_feed < 0 {
-		rng_feed += _LEN;
+	rng.feed--;
+	if rng.feed < 0 {
+		rng.feed += _LEN;
 	}
 
-	x := (rng_vec[rng_feed]+rng_vec[rng_tap])&_MASK;
-	rng_vec[rng_feed] = x;
+	x := (rng.vec[rng.feed] + rng.vec[rng.tap])&_MASK;
+	rng.vec[rng.feed] = x;
 	return x;
 }
-
-// Uint32 returns a pseudo-random 32-bit value as a uint32.
-func Uint32() uint32 {
-	return uint32(Int63()>>31);
-}
-
-// Int31 returns a non-negative pseudo-random 31-bit integer as an int32.
-func Int31() int32 {
-	return int32(Int63()>>32);
-}
-
-// Int returns a non-negative pseudo-random int.  All bits but the top bit are random.
-func Int() int {
-	u := uint(Int63());
-	return int(u<<1>>1);	// clear sign bit if int == int32
-}
-
-// Int63n returns, as an int64, a non-negative pseudo-random number in [0,n).
-func Int63n(n int64) int64 {
-	if n <= 0 {
-		return 0;
-	}
-	max := int64((1<<63) - 1 - (1<<63)%uint64(n));
-	v := Int63();
-	for v > max {
-		v = Int63();
-	}
-	return v%n;
-}
-
-// Int31n returns, as an int32, a non-negative pseudo-random number in [0,n).
-func Int31n(n int32) int32 {
-	return int32(Int63n(int64(n)));
-}
-
-// Intn returns, as an int, a non-negative pseudo-random number in [0,n).
-func Intn(n int) int {
-	return int(Int63n(int64(n)));
-}
-
-// Float64 returns, as a float64, a pseudo-random number in [0.0,1.0).
-func Float64() float64 {
-	x := float64(Int63())/float64(_MAX);
-	for x >= 1 {
-		x = float64(Int63())/float64(_MAX);
-	}
-	return x;
-}
-
-// Float32 returns, as a float32, a pseudo-random number in [0.0,1.0).
-func Float32() float32 {
-	return float32(Float64());
-}
-
-// Float returns, as a float, a pseudo-random number in [0.0,1.0).
-func Float() float {
-	return float(Float64());
-}
-
-// Perm returns, as a slice of n ints, a pseudo-random permutation of the integers [0,n).
-func Perm(n int) []int {
-	m := make([]int, n);
-	for i := 0; i < n; i++ {
-		m[i] = i;
-	}
-	for i := 0; i < n; i++ {
-		j := Intn(n);
-		m[i], m[j] = m[j], m[i];
-	}
-	return m;
-}
-
-func init() {
-	Seed(1);
-}