From 8d7b109462d9ff09075dd481bea3f9d3a1218514 Mon Sep 17 00:00:00 2001
From: didier amyot <didier@miaou>
Date: Tue, 6 Oct 2015 23:10:21 -0400
Subject: [PATCH] this is a bad commit message

---
 aes.go         | 44 +++++++++++++++++++++++++++++++++++
 buxMuller.go   | 22 ++++++++++++++++++
 marsaglia.go   | 32 +++++++++++++++++++++++++
 random.go      | 14 +++++++++++
 random_test.go | 63 ++++++++++++++++++++++++++++++++++++++++++++++++++
 readme.md      | 13 +++++++++++
 6 files changed, 188 insertions(+)
 create mode 100644 aes.go
 create mode 100644 buxMuller.go
 create mode 100644 marsaglia.go
 create mode 100644 random.go
 create mode 100644 random_test.go
 create mode 100644 readme.md

diff --git a/aes.go b/aes.go
new file mode 100644
index 0000000..77f8932
--- /dev/null
+++ b/aes.go
@@ -0,0 +1,44 @@
+package main
+
+import "crypto/aes"
+
+var seed = []byte{1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5, 6}
+var idx = 0
+
+var block, _ = aes.NewCipher([]byte{1, 4, 2, 5, 1, 4, 2, 5, 1, 4, 5, 2, 1, 4, 5, 2})
+
+func getBytesForInt() []byte {
+	if idx < len(seed) {
+		idx += 8
+		return seed[idx-8 : idx]
+	}
+	block.Encrypt(seed, seed)
+	idx = 8
+	return seed[0:8]
+}
+
+func aesInt63() int64 {
+	intAsBytes := getBytesForInt()
+	ret := int64(intAsBytes[0])
+	ret += 265 * int64(intAsBytes[1])
+	ret += 265 * 265 * int64(intAsBytes[2])
+	ret += 265 * 265 * 265 * int64(intAsBytes[3])
+	ret += 265 * 265 * 265 * 256 * int64(intAsBytes[4])
+	ret += 265 * 265 * 265 * 256 * 256 * int64(intAsBytes[5])
+	ret += 265 * 265 * 265 * 256 * 256 * 256 * int64(intAsBytes[6])
+	ret += 265 * 265 * 265 * 256 * 256 * 256 * 256 * int64(intAsBytes[7])
+
+	if ret < 0 {
+		return -ret
+	}
+	return ret
+
+}
+
+func aesRandom() float64 {
+	f := float64(aesInt63()) / (1 << 63)
+	if f == 1 {
+		f = 0
+	}
+	return f
+}
diff --git a/buxMuller.go b/buxMuller.go
new file mode 100644
index 0000000..4639a49
--- /dev/null
+++ b/buxMuller.go
@@ -0,0 +1,22 @@
+package main
+
+import (
+	"math"
+	"math/rand"
+)
+
+var normalBoxMullerAlreadyComputer bool = false
+var normalBoxMullerAlreadyComputedValue float64 = 0
+
+func NormalBoxMuller(r *rand.Rand) float64 {
+	if normalBoxMullerAlreadyComputer {
+		normalBoxMullerAlreadyComputer = false
+		return normalBoxMullerAlreadyComputedValue
+	}
+	radius := math.Sqrt(-2 * math.Log(r.Float64()))
+	theta := 2 * math.Pi * r.Float64()
+
+	normalBoxMullerAlreadyComputer = true
+	normalBoxMullerAlreadyComputedValue = radius * math.Cos(theta)
+	return radius * math.Cos(theta)
+}
diff --git a/marsaglia.go b/marsaglia.go
new file mode 100644
index 0000000..62b6ca2
--- /dev/null
+++ b/marsaglia.go
@@ -0,0 +1,32 @@
+package main
+
+import (
+	"math"
+	"math/rand"
+)
+
+var normalMarsagliaAlreadyComputer bool = false
+var normalMarsagliaAlreadyComputedValue float64 = 0
+
+func NormalMarsaglia(r *rand.Rand) float64 {
+	if normalMarsagliaAlreadyComputer {
+		normalMarsagliaAlreadyComputer = false
+		return normalMarsagliaAlreadyComputedValue
+	}
+	x, y := getPointInUnitDisc(r)
+	radius := math.Sqrt(x*x + y*y)
+
+	normalMarsagliaAlreadyComputer = true
+	coeff := math.Sqrt(-4 * math.Log(radius) / radius)
+	normalMarsagliaAlreadyComputedValue = x * coeff
+	return y * coeff
+}
+
+func getPointInUnitDisc(r *rand.Rand) (float64, float64) {
+	x := 2 * (r.Float64() - 0.5)
+	y := 2 * (r.Float64() - 0.5)
+	if x*x+y*y < 1 {
+		return x, y
+	}
+	return getPointInUnitDisc(r)
+}
diff --git a/random.go b/random.go
new file mode 100644
index 0000000..933cb42
--- /dev/null
+++ b/random.go
@@ -0,0 +1,14 @@
+package main
+
+import (
+	"math"
+	"math/rand"
+)
+
+func ExpFloat64(r *rand.Rand) float64 {
+	return -math.Log(r.Float64())
+}
+
+func main() {
+
+}
diff --git a/random_test.go b/random_test.go
new file mode 100644
index 0000000..b1c318f
--- /dev/null
+++ b/random_test.go
@@ -0,0 +1,63 @@
+package main
+
+import (
+	"math/rand"
+	"testing"
+
+	"github.com/bszcz/mt19937_64"
+)
+
+func BenchmarkMersenneTwister(b *testing.B) {
+	r := mt19937_64.New()
+	for i := 0; i < b.N; i++ {
+		r.Real1()
+	}
+}
+
+func BenchmarkStdUniforme01(b *testing.B) {
+	r := rand.New(rand.NewSource(99))
+	for i := 0; i < b.N; i++ {
+		r.Float64()
+	}
+}
+
+func BenchmarkAESRandom(b *testing.B) {
+	for i := 0; i < b.N; i++ {
+		aesRandom()
+	}
+}
+
+func BenchmarkStdExpoVariate(b *testing.B) {
+	r := rand.New(rand.NewSource(99))
+	for i := 0; i < b.N; i++ {
+		r.ExpFloat64()
+	}
+}
+
+func BenchmarkExpoVariate(b *testing.B) {
+	r := rand.New(rand.NewSource(99))
+	for i := 0; i < b.N; i++ {
+		ExpFloat64(r)
+	}
+}
+
+func BenchmarkNormalBoxMuller(b *testing.B) {
+	r := rand.New(rand.NewSource(99))
+	for i := 0; i < b.N; i++ {
+		NormalBoxMuller(r)
+	}
+}
+
+func BenchmarkStdNormal(b *testing.B) {
+	r := rand.New(rand.NewSource(99))
+	for i := 0; i < b.N; i++ {
+		r.NormFloat64()
+	}
+}
+
+func BenchmarkNormalMarsaglia(b *testing.B) {
+	r := rand.New(rand.NewSource(99))
+	for i := 0; i < b.N; i++ {
+		NormalMarsaglia(r)
+	}
+}
diff --git a/readme.md b/readme.md
new file mode 100644
index 0000000..635d2ae
--- /dev/null
+++ b/readme.md
@@ -0,0 +1,13 @@
+# uqam simulation
+
+Exploratory project to look into performance of some method to simulate
+most common random variable in comparison with the ones provided by Go's
+standard library.
+
+As one can expect the conclusion is to use the functions provided by Go
+and not to write your own.
+
+### running the benchmarks
+Just use Go's benchmarking capabilities with the bellow line of code.
+
+    go test -bench .