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
(limited to 'random')
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