From 685950d3f2ac8e893a23443f9e087294cce2c6fa Mon Sep 17 00:00:00 2001
From: Marc Sunet <msunet@shellblade.net>
Date: Wed, 11 May 2022 09:56:33 -0700
Subject: Add random variable.

---
 random/CMakeLists.txt          |  3 ++-
 random/include/random/normal.h |  9 +++++++++
 random/include/random/random.h |  2 ++
 random/src/normal.c            | 17 +++++++++++++++++
 4 files changed, 30 insertions(+), 1 deletion(-)
 create mode 100644 random/include/random/normal.h
 create mode 100644 random/src/normal.c

diff --git a/random/CMakeLists.txt b/random/CMakeLists.txt
index ec80b1d..188d16b 100644
--- a/random/CMakeLists.txt
+++ b/random/CMakeLists.txt
@@ -3,7 +3,8 @@ cmake_minimum_required(VERSION 3.0)
 project(random)
 
 add_library(random
-  src/mt19937-64.c)
+  src/mt19937-64.c
+  src/normal.c)
 
 target_include_directories(random PUBLIC include)
 
diff --git a/random/include/random/normal.h b/random/include/random/normal.h
new file mode 100644
index 0000000..bee32a9
--- /dev/null
+++ b/random/include/random/normal.h
@@ -0,0 +1,9 @@
+#pragma once
+
+/// Generate two samples from the standard normal distribution.
+///
+/// |u| and |v| must be uniformly distributed in (0,1).
+void normal2(double u, double v, double* z0, double* z1);
+
+/// Map a sample from a standard normal distribution to an arbitrary normal.
+double normal_transform(double z, double mu, double sigma);
diff --git a/random/include/random/random.h b/random/include/random/random.h
index 5499f62..1f4a48d 100644
--- a/random/include/random/random.h
+++ b/random/include/random/random.h
@@ -1,3 +1,5 @@
 #pragma once
 
 #include <random/mt19937-64.h>
+#include <random/normal.h>
+
diff --git a/random/src/normal.c b/random/src/normal.c
new file mode 100644
index 0000000..d4bcac1
--- /dev/null
+++ b/random/src/normal.c
@@ -0,0 +1,17 @@
+#include <random/normal.h>
+
+#include <math.h>
+
+// Generate two samples in the standard normal distribution using the
+// Box-Muller transform.
+// https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform
+void normal2(double u, double v, double* z0, double* z1) {
+  const double r = sqrt(-2 * log(u));
+  const double x = 2 * M_PI * v;
+  *z0 = r * cos(x);
+  *z1 = r * sin(x);
+}
+
+double normal_transform(double z, double mu, double sigma) {
+  return z*sigma + mu;
+}
-- 
cgit v1.2.3