diff --git a/.gitmodules b/.gitmodules
new file mode 100644
index 00000000..12b39b25
--- /dev/null
+++ b/.gitmodules
@@ -0,0 +1,6 @@
+[submodule "docs/py/Benchmarks/discrete-gaussian-differential-privacy"]
+	path = docs/py/Benchmarks/discrete-gaussian-differential-privacy
+	url = https://github.com/IBM/discrete-gaussian-differential-privacy/
+[submodule "docs/py/Benchmarks/differential-privacy-library"]
+	path = docs/py/Benchmarks/differential-privacy-library
+	url = https://github.com/IBM/differential-privacy-library/
diff --git a/build/py/run_gaussian_benchmarks.sh b/build/py/run_gaussian_benchmarks.sh
new file mode 100755
index 00000000..f5ad6af7
--- /dev/null
+++ b/build/py/run_gaussian_benchmarks.sh
@@ -0,0 +1,3 @@
+#1/bin/bash
+
+PYTHONPATH=.:build/py/DafnyVMC-py:docs/py/benchmarks/differential-privacy-library:docs/py/benchmarks/discrete-gaussian-differential-privacy python3 docs/py/benchmarks/gaussian_benchmarks.py
\ No newline at end of file
diff --git a/build/py/run_gaussian_diagrams.sh b/build/py/run_gaussian_diagrams.sh
new file mode 100755
index 00000000..0a86feee
--- /dev/null
+++ b/build/py/run_gaussian_diagrams.sh
@@ -0,0 +1,3 @@
+#1/bin/bash
+
+PYTHONPATH=.:build/py/DafnyVMC-py:docs/py/benchmarks/differential-privacy-library:docs/py/benchmarks/discrete-gaussian-differential-privacy python3 docs/py/benchmarks/gaussian_diagrams.py
\ No newline at end of file
diff --git a/docs/py/Benchmarks/differential-privacy-library b/docs/py/Benchmarks/differential-privacy-library
new file mode 160000
index 00000000..2ec58659
--- /dev/null
+++ b/docs/py/Benchmarks/differential-privacy-library
@@ -0,0 +1 @@
+Subproject commit 2ec58659f8a1795a02f73d1fb0506477b3345a55
diff --git a/docs/py/Benchmarks/discrete-gaussian-differential-privacy b/docs/py/Benchmarks/discrete-gaussian-differential-privacy
new file mode 160000
index 00000000..cb190d2a
--- /dev/null
+++ b/docs/py/Benchmarks/discrete-gaussian-differential-privacy
@@ -0,0 +1 @@
+Subproject commit cb190d2a990a78eff6e21159203bc888e095f01b
diff --git a/docs/py/Benchmarks/gaussian_benchmarks.py b/docs/py/Benchmarks/gaussian_benchmarks.py
new file mode 100644
index 00000000..40bc8b7f
--- /dev/null
+++ b/docs/py/Benchmarks/gaussian_benchmarks.py
@@ -0,0 +1,90 @@
+import timeit
+import secrets
+import numpy
+import matplotlib.pyplot as plt
+from decimal import Decimal
+import DafnyVMC
+from diffprivlib.mechanisms import GaussianDiscrete
+import discretegauss
+from datetime import datetime
+import tqdm
+
+vmc_mean = []
+vmc_std = []
+ibm_dgdp_mean = []
+ibm_dgdp_std = []
+ibm_dpl_mean = []
+ibm_dpl_std = []
+
+fig,ax1 = plt.subplots()
+
+rng = secrets.SystemRandom()
+r = DafnyVMC.Random()
+
+sigmas = []
+for epsilon_times_100 in tqdm.tqdm(range(1, 500, 2)):
+    vmc = []
+    ibm_dgdp = []
+    ibm_dpl = []
+
+    # The GaussianDiscrete class does not expose the sampler directly, and needs to be instantiated with `(epsilon, delta)`.
+    # We access its `_scale` member to get the values `sigma`'s needed by `DafnyVMC` and `discretegauss`.
+    g = GaussianDiscrete(epsilon=0.01 * epsilon_times_100, delta=0.00001)
+    sigma = g._scale
+    sigmas += [sigma]
+
+    sigma_num, sigma_denom = Decimal(sigma).as_integer_ratio()
+    sigma_squared = sigma ** 2
+
+    for i in range(1100):
+        start_time = timeit.default_timer()
+        r.DiscreteGaussianSample(sigma_num, sigma_denom)
+        elapsed = timeit.default_timer() - start_time
+        vmc.append(elapsed)
+
+    for i in range(1100):
+        start_time = timeit.default_timer()
+        discretegauss.sample_dgauss(sigma_squared, rng)
+        elapsed = timeit.default_timer() - start_time
+        ibm_dgdp.append(elapsed)
+
+    for i in range(1100):
+        start_time = timeit.default_timer()
+        # The sampler is not directly accessible, so we call `.randomise(0)` instead, as it adds a noise drawn according to a discrete Gaussian to `0`.
+        g.randomise(0)
+        elapsed = timeit.default_timer() - start_time
+        ibm_dpl.append(elapsed)
+
+    vmc = numpy.array(vmc[-1000:])
+    ibm_dgdp = numpy.array(ibm_dgdp[-1000:])
+    ibm_dpl = numpy.array(ibm_dpl[-1000:])
+
+    vmc_mean.append(vmc.mean()*1000.0)
+    vmc_std.append(vmc.std()*1000.0)
+    ibm_dgdp_mean.append(ibm_dgdp.mean()*1000.0)
+    ibm_dgdp_std.append(ibm_dgdp.std()*1000.0)
+    ibm_dpl_mean.append(ibm_dpl.mean()*1000.0)
+    ibm_dpl_std.append(ibm_dpl.std()*1000.0)
+
+
+ax1.plot(sigmas, vmc_mean, color='green', linewidth=1.0, label='VMC')
+ax1.fill_between(sigmas, numpy.array(vmc_mean)-0.5*numpy.array(vmc_std), numpy.array(vmc_mean)+0.5*numpy.array(vmc_std),
+    alpha=0.2, facecolor='k',
+    linewidth=2, linestyle='dashdot', antialiased=True)
+
+ax1.plot(sigmas, ibm_dgdp_mean, color='red', linewidth=1.0, label='IBM-DGDP')
+ax1.fill_between(sigmas, numpy.array(ibm_dgdp_mean)-0.5*numpy.array(ibm_dgdp_std), numpy.array(ibm_dgdp_mean)+0.5*numpy.array(ibm_dgdp_std),
+    alpha=0.2,  facecolor='y',
+    linewidth=2, linestyle='dashdot', antialiased=True)
+
+ax1.plot(sigmas, ibm_dpl_mean, color='purple', linewidth=1.0, label='IBM-DPL')
+ax1.fill_between(sigmas, numpy.array(ibm_dpl_mean)-0.5*numpy.array(ibm_dpl_std), numpy.array(ibm_dpl_mean)+0.5*numpy.array(ibm_dpl_std),
+    alpha=0.2,  facecolor='y',
+    linewidth=2, linestyle='dashdot', antialiased=True)
+
+ax1.set_xlabel("Sigma")
+ax1.set_ylabel("Sampling Time (ms)")
+plt.legend(loc = 'best')
+now = datetime.now()
+filename = 'GaussianBenchmarks' + now.strftime("%H%M%S") + '.pdf'
+plt.savefig(filename)
\ No newline at end of file
diff --git a/docs/py/Benchmarks/gaussian_diagrams.py b/docs/py/Benchmarks/gaussian_diagrams.py
new file mode 100644
index 00000000..2e460e72
--- /dev/null
+++ b/docs/py/Benchmarks/gaussian_diagrams.py
@@ -0,0 +1,45 @@
+import matplotlib.pyplot as plt
+import secrets
+from decimal import Decimal
+from datetime import datetime
+import DafnyVMC
+import discretegauss
+from diffprivlib.mechanisms import GaussianDiscrete
+
+fig, axs = plt.subplots(8, 3, figsize=(20, 20))
+
+rng = secrets.SystemRandom()
+r = DafnyVMC.Random()
+
+for i in range(8):
+    vmc_data = []
+    ibm_dgdp_data = []
+    ibm_dpl_data = []
+
+    epsilon_times_100 = 1 + (i**2)*2.5
+    g = GaussianDiscrete(epsilon=0.01 * epsilon_times_100, delta=0.00001)
+
+    sigma = g._scale
+    sigma_squared = sigma ** 2
+    sigma_num, sigma_denom = Decimal(sigma).as_integer_ratio()
+
+    title_vmc = 'VMC, Sigma = ' + str(sigma)
+    title_ibm_dgdp = 'IBM-DGDP, Sigma = ' + str(sigma)
+    title_ibm_dpl = 'IBM-DPL, Sigma = ' + str(sigma)
+
+    for _ in range(100000):
+        vmc_data.append(r.DiscreteGaussianSample(sigma_num, sigma_denom))
+        ibm_dgdp_data.append(discretegauss.sample_dgauss(sigma_squared, rng))
+        ibm_dpl_data.append(g.randomise(0))
+
+    axs[i, 0].hist(vmc_data, color='lightgreen', ec='black', bins=50)
+    axs[i, 0].set_title(title_vmc)
+    axs[i, 1].hist(ibm_dgdp_data, color='lightgreen', ec='black', bins=50)
+    axs[i, 1].set_title(title_ibm_dgdp)
+    axs[i, 2].hist(ibm_dpl_data, color='lightgreen', ec='black', bins=50)
+    axs[i, 2].set_title(title_ibm_dpl)
+
+now = datetime.now()
+filename = 'GaussianDiagrams' + now.strftime("%H%M%S") + '.pdf'
+plt.subplots_adjust(wspace=0.4, hspace=0.4)
+plt.savefig(filename)
\ No newline at end of file