diff --git a/audit.log b/audit.log
index 015176fc..6408e696 100644
--- a/audit.log
+++ b/audit.log
@@ -1,30 +1,6 @@
 src/DafnyVMC.dfy(23,6): UniformPowerOfTwoSample: Definition has `assume {:axiom}` statement in body. 
-src/Distributions/Bernoulli/Correctness.dfy(59,17): BernoulliCorrectnessCaseFalse: Declaration has explicit `{:axiom}` attribute. 
-src/Distributions/BernoulliExpNeg/Correctness.dfy(229,6): CorrectnessLe1: Definition has `assume {:axiom}` statement in body. 
-src/Distributions/BernoulliExpNeg/Correctness.dfy(346,17): Le1SeriesConvergesToExpNeg: Declaration has explicit `{:axiom}` attribute. 
-src/Distributions/BernoulliExpNeg/Correctness.dfy(411,6): Le1LoopCutDecomposeProb: Definition has `assume {:axiom}` statement in body. 
-src/Distributions/BernoulliExpNeg/Correctness.dfy(412,6): Le1LoopCutDecomposeProb: Definition has `assume {:axiom}` statement in body. 
-src/Distributions/BernoulliExpNeg/Correctness.dfy(419,6): Le1LoopCutDecomposeProb: Definition has `assume {:axiom}` statement in body. 
-src/Distributions/BernoulliExpNeg/Correctness.dfy(420,6): Le1LoopCutDecomposeProb: Definition has `assume {:axiom}` statement in body. 
-src/Distributions/BernoulliExpNeg/Correctness.dfy(424,6): Le1LoopCutDecomposeProb: Definition has `assume {:axiom}` statement in body. 
-src/Distributions/BernoulliExpNeg/Correctness.dfy(425,6): Le1LoopCutDecomposeProb: Definition has `assume {:axiom}` statement in body. 
-src/Distributions/BernoulliExpNeg/Correctness.dfy(615,8): Le1LoopCorrectnessEq: Definition has `assume {:axiom}` statement in body. 
-src/Distributions/BernoulliExpNeg/Correctness.dfy(616,8): Le1LoopCorrectnessEq: Definition has `assume {:axiom}` statement in body. 
-src/Distributions/BernoulliExpNeg/Model.dfy(112,10): Le1LoopTerminatesAlmostSurely: Definition has `assume {:axiom}` statement in body. 
-src/Distributions/BernoulliExpNeg/Model.dfy(148,4): Le1LoopDivergenceProbabilityBound: Definition has `assume {:axiom}` statement in body. 
-src/Distributions/BernoulliExpNeg/Model.dfy(161,8): Le1LoopDivergenceProbabilityBound: Definition has `assume {:axiom}` statement in body. 
-src/Distributions/BernoulliExpNeg/Model.dfy(162,8): Le1LoopDivergenceProbabilityBound: Definition has `assume {:axiom}` statement in body. 
-src/Distributions/BernoulliExpNeg/Model.dfy(185,10): Le1LoopDivergenceProbabilityBound: Definition has `assume {:axiom}` statement in body. 
-src/Distributions/BernoulliExpNeg/Model.dfy(71,17): Le1LoopIterCorrectness: Declaration has explicit `{:axiom}` attribute. 
-src/Distributions/DiscreteLaplace/Correctness.dfy(37,17): Correctness: Declaration has explicit `{:axiom}` attribute. 
-src/Distributions/DiscreteLaplace/Model.dfy(89,17): SampleInnerLoopTerminatesAlmostSurely: Declaration has explicit `{:axiom}` attribute. 
-src/Distributions/DiscreteLaplace/Model.dfy(93,17): SampleLoopTerminatesAlmostSurely: Declaration has explicit `{:axiom}` attribute. 
 src/Math/Analysis/Limits.dfy(126,17): Sandwich: Declaration has explicit `{:axiom}` attribute. 
 src/Math/Analysis/Reals.dfy(35,17): LeastUpperBoundProperty: Declaration has explicit `{:axiom}` attribute. 
-src/Math/Exponential.dfy(55,17): ExpSeriesConverges: Declaration has explicit `{:axiom}` attribute. 
-src/Math/Exponential.dfy(58,17): FunctionalEquation: Declaration has explicit `{:axiom}` attribute. 
-src/Math/Exponential.dfy(61,17): Increasing: Declaration has explicit `{:axiom}` attribute. 
-src/Math/Exponential.dfy(65,17): EvalOne: Declaration has explicit `{:axiom}` attribute. 
 src/Math/Measures.dfy(125,17): ProbabilityLe1: Declaration has explicit `{:axiom}` attribute. 
 src/Math/Measures.dfy(130,17): IsMonotonic: Declaration has explicit `{:axiom}` attribute. 
 src/ProbabilisticProgramming/Independence.dfy(105,17): BindIsIndep: Declaration has explicit `{:axiom}` attribute. 
@@ -35,12 +11,10 @@ src/ProbabilisticProgramming/Independence.dfy(77,17): IsIndepImpliesMeasurableNa
 src/ProbabilisticProgramming/Independence.dfy(82,17): IsIndepImpliesIsIndepFunction: Declaration has explicit `{:axiom}` attribute. 
 src/ProbabilisticProgramming/Independence.dfy(87,17): CoinIsIndep: Declaration has explicit `{:axiom}` attribute. 
 src/ProbabilisticProgramming/Independence.dfy(91,17): ReturnIsIndep: Declaration has explicit `{:axiom}` attribute. 
-src/ProbabilisticProgramming/Loops.dfy(362,17): EnsureWhileTerminates: Declaration has explicit `{:axiom}` attribute. 
-src/ProbabilisticProgramming/Loops.dfy(367,17): EnsureWhileTerminatesAlmostSurelyViaLimit: Declaration has explicit `{:axiom}` attribute. 
-src/ProbabilisticProgramming/Loops.dfy(385,17): WhileProbabilityViaLimit: Declaration has explicit `{:axiom}` attribute. 
-src/ProbabilisticProgramming/Loops.dfy(402,17): UntilProbabilityFraction: Declaration has explicit `{:axiom}` attribute. 
-src/ProbabilisticProgramming/Loops.dfy(435,4): EnsureUntilTerminatesAndForAll: Definition has `assume {:axiom}` statement in body. 
-src/ProbabilisticProgramming/Loops.dfy(458,4): WhileIsIndep: Definition has `assume {:axiom}` statement in body. 
+src/ProbabilisticProgramming/Loops.dfy(361,17): EnsureWhileTerminates: Declaration has explicit `{:axiom}` attribute. 
+src/ProbabilisticProgramming/Loops.dfy(371,17): UntilProbabilityFraction: Declaration has explicit `{:axiom}` attribute. 
+src/ProbabilisticProgramming/Loops.dfy(404,4): EnsureUntilTerminatesAndForAll: Definition has `assume {:axiom}` statement in body. 
+src/ProbabilisticProgramming/Loops.dfy(427,4): WhileIsIndep: Definition has `assume {:axiom}` statement in body. 
 src/ProbabilisticProgramming/RandomSource.dfy(50,17): ProbIsProbabilityMeasure: Declaration has explicit `{:axiom}` attribute. 
 src/ProbabilisticProgramming/RandomSource.dfy(54,17): CoinHasProbOneHalf: Declaration has explicit `{:axiom}` attribute. 
 src/ProbabilisticProgramming/RandomSource.dfy(61,17): MeasureHeadDrop: Declaration has explicit `{:axiom}` attribute. 
diff --git a/docs/dafny/ExamplesRandom.dfy b/docs/dafny/ExamplesRandom.dfy
index b33d1f89..041f1ad0 100644
--- a/docs/dafny/ExamplesRandom.dfy
+++ b/docs/dafny/ExamplesRandom.dfy
@@ -4,7 +4,6 @@
  *******************************************************************************/
 
 module Examples {
-  import Rationals
   import DafnyVMC
   import Helper
 
@@ -15,13 +14,6 @@ module Examples {
     var r: DafnyVMC.Random := new DafnyVMC.Random();
 
     var t := 0;
-    for i := 0 to n {
-      var b := r.CoinSample();
-      if b {
-        t := t + 1;
-      }
-    }
-    print "Estimated parameter for CoinSample(): ", (t as real) / (n as real), " (should be around 0.5)\n";
 
     var a := 0;
     var b := 0;
@@ -55,7 +47,7 @@ module Examples {
 
     t := 0;
     for i := 0 to n {
-      var b := r.BernoulliSample(Rationals.Rational(1, 5));
+      var b := r.BernoulliSample(1, 5);
       if b {
         t := t + 1;
       }
@@ -65,7 +57,7 @@ module Examples {
 
     t := 0;
     for i := 0 to n {
-      var b := r.BernoulliSample(Rationals.Rational(0, 5));
+      var b := r.BernoulliSample(0, 5);
       if b {
         t := t + 1;
       }
@@ -75,7 +67,7 @@ module Examples {
 
     t := 0;
     for i := 0 to n {
-      var b := r.BernoulliSample(Rationals.Rational(5, 5));
+      var b := r.BernoulliSample(5, 5);
       if b {
         t := t + 1;
       }
@@ -85,7 +77,7 @@ module Examples {
 
     t := 0;
     for i := 0 to n {
-      var u := r.BernoulliExpNegSample(Rationals.Rational(12381, 5377)); // about -ln(0.1)
+      var u := r.BernoulliExpNegSample(12381, 5377); // about -ln(0.1)
       if u {
         t := t + 1;
       }
@@ -96,7 +88,7 @@ module Examples {
     var count1 := 0;
     var countneg1 := 0;
     for i := 0 to n {
-      var u := r.DiscreteLaplaceSample(Rationals.Rational(7, 5));
+      var u := r.DiscreteLaplaceSample(7, 5);
       match u {
         case -1 => countneg1 := countneg1 + 1;
         case 0 => count0 := count0 + 1;
@@ -113,7 +105,7 @@ module Examples {
     count1 := 0;
     countneg1 := 0;
     for i := 0 to n {
-      var u := r.DiscreteGaussianSample(Rationals.Rational(7, 5));
+      var u := r.DiscreteGaussianSample(7, 5);
       match u {
         case -1 => countneg1 := countneg1 + 1;
         case 0 => count0 := count0 + 1;
diff --git a/docs/java/BuildTest.java b/docs/java/BuildTest.java
index f048cdf1..2c64a04d 100644
--- a/docs/java/BuildTest.java
+++ b/docs/java/BuildTest.java
@@ -13,30 +13,28 @@ public static void main(String[] args) {
         boolean[] arr5 = {true, false, false, true};
         long[] arr6 = {111111L, 333333L, 999999L};
         short[] arr7 = {-3, 0, 3};
-        Rationals.Rational gamma = new Rationals.Rational(BigInteger.valueOf(3), BigInteger.valueOf(5));
+        BigInteger num = BigInteger.valueOf(3);
+        BigInteger den = BigInteger.valueOf(5);
         
         /* STANDARD RNG */
         System.out.println("\nSTANDARD RNG TESTS\n");
         
         DafnyVMC.Random r = new DafnyVMC.Random();
 
-        System.out.println("Example of Coin sampling");
-        System.out.println(r.CoinSample());
-
         System.out.println("Example of Uniform sampling");
         System.out.println(r.UniformSample(BigInteger.valueOf(4)));
 
         System.out.println("Example of Bernoulli sampling");
-        System.out.println(r.BernoulliSample(gamma));
+        System.out.println(r.BernoulliSample(num,den));
 
         System.out.println("Example of BernoulliExpNeg sampling");
-        System.out.println(r.BernoulliExpNegSample(gamma));
+        System.out.println(r.BernoulliExpNegSample(num,den));
 
         System.out.println("Example of DiscreteGaussian sampling");
-        System.out.println(r.DiscreteGaussianSample(gamma));
+        System.out.println(r.DiscreteGaussianSample(num,den));
 
         System.out.println("Example of DiscreteLaplace sampling");
-        System.out.println(r.DiscreteLaplaceSample(gamma));
+        System.out.println(r.DiscreteLaplaceSample(num,den));
 
         System.out.println("Example of Fisher-Yates: BigInteger");
         r.Shuffle(arr1);
@@ -72,23 +70,20 @@ public static void main(String[] args) {
         SecureRandom rng = new SecureRandom();
         DafnyVMC.Random t = new DafnyVMC.Random(() -> rng);
 
-        System.out.println("Example of Coin sampling");
-        System.out.println(t.CoinSample());
-
         System.out.println("Example of Uniform sampling");
         System.out.println(t.UniformSample(BigInteger.valueOf(4)));
 
         System.out.println("Example of Bernoulli sampling");
-        System.out.println(t.BernoulliSample(gamma));
+        System.out.println(t.BernoulliSample(num,den));
 
         System.out.println("Example of BernoulliExpNeg sampling");
-        System.out.println(r.BernoulliExpNegSample(gamma));
+        System.out.println(r.BernoulliExpNegSample(num,den));
 
         System.out.println("Example of DiscreteGaussian sampling");
-        System.out.println(t.DiscreteGaussianSample(gamma));
+        System.out.println(t.DiscreteGaussianSample(num,den));
 
         System.out.println("Example of DiscreteLaplace sampling");
-        System.out.println(t.DiscreteLaplaceSample(gamma));
+        System.out.println(t.DiscreteLaplaceSample(num,den));
 
         System.out.println("Example of Fisher-Yates: BigInteger");
         t.Shuffle(arr1);
diff --git a/scripts/audit.sh b/scripts/audit.sh
index 6815b666..12de8579 100755
--- a/scripts/audit.sh
+++ b/scripts/audit.sh
@@ -8,4 +8,4 @@ then
   exit 1
 fi
 
-$DAFNY audit dfyconfig.toml | grep -v '{:termination false}\|{:extern}\|decreases *\|Dafny auditor completed\|Dafny program verifier\|No terms found to trigger on\|Compiled declaration has no body' | sed 's/Warning: //' | sed 's/Possible mitigation: .*//' | sed '/^$/d' | sort > audit.log
+$DAFNY audit dfyconfig.toml | grep -v '{:termination false}\|{:extern}\|decreases *\|DafnyVMCTrait\|Dafny auditor completed\|Dafny program verifier\|No terms found to trigger on\|Compiled declaration has no body' | sed 's/Warning: //' | sed 's/Possible mitigation: .*//' | sed '/^$/d' | sort > audit.log
diff --git a/scripts/checkformat.sh b/scripts/checkformat.sh
index b9be73a3..1f255980 100755
--- a/scripts/checkformat.sh
+++ b/scripts/checkformat.sh
@@ -10,4 +10,5 @@ fi
 
 echo "Checking the formatting of all Dafny files..."
 echo "Run \`DAFNY=dafny scripts/format.sh\` to fix any errors."
+$DAFNY format src/DafnyVMCTrait.dfy
 $DAFNY format --check dfyconfig.toml
diff --git a/src/DafnyVMCTrait.dfy b/src/DafnyVMCTrait.dfy
index f5365a0c..cf60dfd0 100644
--- a/src/DafnyVMCTrait.dfy
+++ b/src/DafnyVMCTrait.dfy
@@ -1,19 +1,227 @@
-/*******************************************************************************
- *  Copyright by the contributors to the Dafny Project
- *  SPDX-License-Identifier: MIT
- *******************************************************************************/
 
 module DafnyVMCTrait {
-  import Coin
-  import Bernoulli
+
   import UniformPowerOfTwo
-  import Uniform
-  import BernoulliExpNeg
-  import DiscreteGaussian
-  import DiscreteLaplace
+
   import FisherYates
-  import Monad
 
-  trait {:termination false} RandomTrait extends Coin.Implementation.Trait, UniformPowerOfTwo.Implementation.Trait, Bernoulli.Implementation.Trait, Uniform.Implementation.Trait, BernoulliExpNeg.Implementation.Trait, DiscreteGaussian.Implementation.Trait, DiscreteLaplace.Implementation.Trait, FisherYates.Implementation.Trait  {
-  }
-}
\ No newline at end of file
+  import opened Pos
+
+  trait {:termination false} RandomTrait extends  UniformPowerOfTwo.Implementation.Trait, FisherYates.Implementation.Trait {
+
+    method {:verify false} UniformSample (n: pos)
+      returns (o: nat)
+      modifies this
+      decreases *
+    {
+      var x := UniformPowerOfTwoSample(2 * n);
+      while ! (x < n)
+        decreases *
+      {
+        x := UniformPowerOfTwoSample(2 * n);
+      }
+      var r := x;
+      o := r;
+    }
+
+    method {:verify false} BernoulliSample (num: nat, den: pos)
+      returns (o: bool)
+      requires num <= den 
+      modifies this
+      decreases *
+    {
+      var d := UniformSample(den);
+      o := d < num;
+    }
+
+    method {:verify false} BernoulliExpNegSampleUnitLoop (num: nat, den: pos, state: (bool,pos))
+      returns (o: (bool,pos))
+      requires num <= den 
+      modifies this
+      decreases *
+    {
+      var A := BernoulliSample(num, state.1 * den);
+      o := (A,state.1 + 1);
+    }
+
+    method {:verify false} BernoulliExpNegSampleUnitAux (num: nat, den: pos)
+      returns (o: nat)
+      requires num <= den 
+      modifies this
+      decreases *
+    {
+      var state := (true,1);
+      while state.0
+        decreases *
+      {
+        state := BernoulliExpNegSampleUnitLoop(num, den, state);
+      }
+      var r := state;
+      o := r.1;
+    }
+
+    method {:verify false} BernoulliExpNegSampleUnit (num: nat, den: pos)
+      returns (o: bool)
+      requires num <= den 
+      modifies this
+      decreases *
+    {
+      var K := BernoulliExpNegSampleUnitAux(num, den);
+      if K % 2 == 0 {
+        o := true;
+      } else {
+        o := false;
+      }
+    }
+
+    method {:verify false} BernoulliExpNegSampleGenLoop (iter: nat)
+      returns (o: bool)
+      modifies this
+      decreases *
+    {
+      if iter == 0 {
+        o := true;
+      } else {
+        var B := BernoulliExpNegSampleUnit(1, 1);
+        var R := BernoulliExpNegSampleGenLoop(iter - 1);
+        o := B == true && R == true;
+      }
+    }
+
+    method {:verify false} BernoulliExpNegSample (num: nat, den: pos)
+      returns (o: bool)
+      modifies this
+      decreases *
+    {
+      if num <= den {
+        var X := BernoulliExpNegSampleUnit(num, den);
+        o := X;
+      } else {
+        var gamf := num / den;
+        var B := BernoulliExpNegSampleGenLoop(gamf);
+        if B == true {
+          var num := num - gamf * den;
+          var X := BernoulliExpNegSampleUnit(num, den);
+          o := X;
+        } else {
+          o := false;
+        }
+      }
+    }
+
+    method {:verify false} DiscreteLaplaceSampleLoopIn1Aux (t: pos)
+      returns (o: (nat,bool))
+      modifies this
+      decreases *
+    {
+      var U := UniformSample(t);
+      var D := BernoulliExpNegSample(U, t);
+      o := (U,D);
+    }
+
+    method {:verify false} DiscreteLaplaceSampleLoopIn1 (t: pos)
+      returns (o: nat)
+      modifies this
+      decreases *
+    {
+      var x := DiscreteLaplaceSampleLoopIn1Aux(t);
+      while ! (x.1)
+        decreases *
+      {
+        x := DiscreteLaplaceSampleLoopIn1Aux(t);
+      }
+      var r1 := x;
+      o := r1.0;
+    }
+
+    method {:verify false} DiscreteLaplaceSampleLoopIn2Aux (num: nat, den: pos, K: (bool,pos))
+      returns (o: (bool,pos))
+      requires num <= den 
+      modifies this
+      decreases *
+    {
+      var A := BernoulliExpNegSampleUnit(num, den);
+      o := (A,K.1 + 1);
+    }
+
+    method {:verify false} DiscreteLaplaceSampleLoopIn2 (num: nat, den: pos)
+      returns (o: pos)
+      modifies this
+      decreases *
+    {
+      var K := (true,1);
+      while K.0
+        decreases *
+      {
+        K := DiscreteLaplaceSampleLoopIn2Aux(1, 1, K);
+      }
+      var r2 := K;
+      o := r2.1;
+    }
+
+    method {:verify false} DiscreteLaplaceSampleLoop (num: pos, den: pos)
+      returns (o: (bool,nat))
+      modifies this
+      decreases *
+    {
+      var U := DiscreteLaplaceSampleLoopIn1(num);
+      var v := DiscreteLaplaceSampleLoopIn2(1, 1);
+      var V := v - 2;
+      var X := U + num * V;
+      var Y := X / den;
+      var B := BernoulliSample(1, 2);
+      o := (B,Y);
+    }
+
+    method {:verify false} DiscreteLaplaceSample (num: pos, den: pos)
+      returns (o: int)
+      modifies this
+      decreases *
+    {
+      var x := DiscreteLaplaceSampleLoop(num, den);
+      while ! (! (x.0 == true && x.1 == 0))
+        decreases *
+      {
+        x := DiscreteLaplaceSampleLoop(num, den);
+      }
+      var r := x;
+      var Z := if r.0 == true then - (r.1) else r.1;
+      o := Z;
+    }
+
+    method {:verify false} DiscreteGaussianSampleLoop (num: pos, den: pos, t: pos)
+      returns (o: (int,bool))
+      modifies this
+      decreases *
+    {
+      var Y := DiscreteLaplaceSample(t, 1);
+      var y := (if Y < 0 then -Y else Y);
+      var n := (y * t * den - num) * (y * t * den - num);
+      var d := 2 * num * (t) * (t) * den;
+      var C := BernoulliExpNegSample(n, d);
+      o := (Y,C);
+    }
+
+    method {:verify false} DiscreteGaussianSample (num: pos, den: pos)
+      returns (o: int)
+      modifies this
+      decreases *
+    {
+      var ti := num / den;
+      var t := ti + 1;
+      var num := (num) * (num);
+      var den := (den) * (den);
+      var x := DiscreteGaussianSampleLoop(num, den, t);
+      while ! (x.1)
+        decreases *
+      {
+        x := DiscreteGaussianSampleLoop(num, den, t);
+      }
+      var r := x;
+      o := r.0;
+    }
+
+
+}
+
+}
diff --git a/src/Distributions/Bernoulli/Correctness.dfy b/src/Distributions/Bernoulli/Correctness.dfy
deleted file mode 100644
index 04662453..00000000
--- a/src/Distributions/Bernoulli/Correctness.dfy
+++ /dev/null
@@ -1,152 +0,0 @@
-/*******************************************************************************
- *  Copyright by the contributors to the Dafny Project
- *  SPDX-License-Identifier: MIT
- *******************************************************************************/
-
-module Bernoulli.Correctness {
-  import Measures
-  import Helper
-  import Uniform
-  import Rand
-  import Independence
-  import Monad
-  import Model
-
-  /*******
-   Lemmas
-  *******/
-
-  lemma SampleIsIndep(m: nat, n: nat)
-    requires n != 0
-    requires m <= n
-    ensures Independence.IsIndep(Model.Sample(m, n))
-  {
-    var f := Uniform.Model.Sample(n);
-    var g := (k: nat) => k < m;
-
-    assert Independence.IsIndep(f) by {
-      Uniform.Correctness.SampleIsIndep(n);
-    }
-
-    Independence.MapIsIndep(f, g);
-    reveal Model.Sample();
-  }
-
-  // The probability mass function (PMF) of the bernoulli distribution
-  opaque function BernoulliMass(numer: nat, denom: nat): bool -> real
-    requires denom > 0
-    requires numer <= denom
-  {
-    b => if b then numer as real / denom as real else 1.0 - numer as real / denom as real
-  }
-
-  lemma BernoulliCorrectness(m: nat, n: nat, b: bool)
-    requires n != 0
-    requires m <= n
-    ensures
-      var event := iset s | Model.Sample(m, n)(s).Equals(b);
-      && event in Rand.eventSpace
-      && Rand.prob(event) == BernoulliMass(m, n)(b)
-  {
-    reveal BernoulliMass();
-    if b {
-      BernoulliCorrectnessCaseTrue(m, n);
-    } else {
-      BernoulliCorrectnessCaseFalse(m, n);
-    }
-  }
-
-  lemma {:axiom} BernoulliCorrectnessCaseFalse(m: nat, n: nat)
-    requires n != 0
-    requires m <= n
-    ensures
-      var falseEvent := iset s | Model.Sample(m, n)(s).Equals(false);
-      && falseEvent in Rand.eventSpace
-      && Rand.prob(falseEvent) == 1.0 - m as real / n as real
-
-
-  lemma BernoulliCorrectnessCaseTrue(m: nat, n: nat)
-    requires n != 0
-    requires m <= n
-    ensures
-      var trueEvent := iset s | Model.Sample(m, n)(s).Equals(true);
-      && trueEvent in Rand.eventSpace
-      && Rand.prob(trueEvent) == m as real / n as real
-  {
-    reveal Model.Sample();
-    var e := iset s | Model.Sample(m, n)(s).Equals(true);
-    var ltM: nat -> bool := x => x < m;
-    var ltM1: nat -> bool := x => x < m - 1;
-
-    if m == 0 {
-      assert e == iset{} by {
-        forall s ensures !Model.Sample(m, n)(s).Equals(true) {
-          calc {
-            Model.Sample(m, n)(s).Equals(true);
-            Uniform.Model.Sample(n)(s).Satisfies(x => x < 0);
-            false;
-          }
-        }
-      }
-
-      assert e in Rand.eventSpace && Rand.prob(e) == 0.0 by {
-        Rand.ProbIsProbabilityMeasure();
-        assert Measures.IsSigmaAlgebra(Rand.eventSpace);
-        assert Measures.IsPositive(Rand.eventSpace, Rand.prob);
-      }
-    } else {
-      var e1 := iset s | Uniform.Model.Sample(n)(s).Equals(m - 1);
-      var e2 := (iset s {:trigger Uniform.Model.Sample(n)(s).value} | Model.Sample(m-1, n)(s).Equals(true));
-
-      assert (iset s | Uniform.Model.Sample(n)(s).Satisfies(ltM1)) == e2 by {
-        assert (iset s | Uniform.Model.Sample(n)(s).Satisfies(ltM1)) == (iset s {:trigger Uniform.Model.Sample(n)(s).value} | Model.Sample(m-1, n)(s).Equals(true)) by {
-          forall s ensures Uniform.Model.Sample(n)(s).Satisfies(ltM1) <==> Model.Sample(m-1, n)(s).Equals(true) {
-            assert Uniform.Model.Sample(n)(s).Satisfies(ltM1) <==> Model.Sample(m-1, n)(s).Equals(true);
-          }
-        }
-      }
-
-      assert e == e1 + e2 by {
-        calc {
-          e;
-          iset s | Uniform.Model.Sample(n)(s).Satisfies(ltM);
-          iset s | Uniform.Model.Sample(n)(s).Equals(m-1) || Uniform.Model.Sample(n)(s).Satisfies(ltM1);
-          (iset s | Uniform.Model.Sample(n)(s).Equals(m-1)) + (iset s | Uniform.Model.Sample(n)(s).Satisfies(ltM1));
-          e1 + e2;
-        }
-      }
-
-      assert A1: e1 in Rand.eventSpace && Rand.prob(e1) == 1.0 / (n as real) by {
-        Uniform.Correctness.UniformFullCorrectness(n, m-1);
-        assert e1 == Uniform.Correctness.SampleEquals(n, m-1);
-      }
-
-      assert A2: e2 in Rand.eventSpace && Rand.prob(e2) == (m-1) as real / n as real by {
-        BernoulliCorrectnessCaseTrue(m-1, n);
-      }
-
-      calc {
-        Rand.prob(e);
-        Rand.prob(e1 + e2);
-        { reveal A1; reveal A2; assert e1 * e2 == iset{}; Rand.ProbIsProbabilityMeasure(); Measures.PosCountAddImpliesAdd(Rand.eventSpace, Rand.prob); assert Measures.IsAdditive(Rand.eventSpace, Rand.prob); }
-        Rand.prob(e1) + Rand.prob(e2);
-        { reveal A1; reveal A2; }
-        (1.0 / n as real) + ((m - 1) as real / n as real);
-        { AdditionOfFractions(m, n); }
-        m as real / n as real;
-      }
-
-      assert e in Rand.eventSpace by {
-        Rand.ProbIsProbabilityMeasure();
-        reveal A1;
-        reveal A2;
-        Measures.BinaryUnionIsMeasurable(Rand.eventSpace, e1, e2);
-      }
-    }
-  }
-
-  lemma AdditionOfFractions(m: nat, n: nat)
-    requires n > 0
-    ensures (1.0 / n as real) + ((m - 1) as real / n as real) == (m as real) / n as real
-  {}
-}
diff --git a/src/Distributions/Bernoulli/Implementation.dfy b/src/Distributions/Bernoulli/Implementation.dfy
deleted file mode 100644
index 4d213504..00000000
--- a/src/Distributions/Bernoulli/Implementation.dfy
+++ /dev/null
@@ -1,27 +0,0 @@
-/*******************************************************************************
- *  Copyright by the contributors to the Dafny Project
- *  SPDX-License-Identifier: MIT
- *******************************************************************************/
-
-module Bernoulli.Implementation {
-  import Rationals
-  import Monad
-  import Model
-  import Interface
-
-  trait {:termination false} Trait extends Interface.Trait {
-
-    method BernoulliSample(p: Rationals.Rational) returns (c: bool)
-      modifies `s
-      decreases *
-      requires 0 <= p.numer <= p.denom
-      ensures Model.Sample(p.numer, p.denom)(old(s)) == Monad.Result(c, s)
-    {
-      reveal Model.Sample();
-      var k := UniformSample(p.denom);
-      c := k < p.numer;
-      reveal Model.Sample();
-    }
-
-  }
-}
diff --git a/src/Distributions/Bernoulli/Interface.dfy b/src/Distributions/Bernoulli/Interface.dfy
deleted file mode 100644
index c640b6e5..00000000
--- a/src/Distributions/Bernoulli/Interface.dfy
+++ /dev/null
@@ -1,21 +0,0 @@
-/*******************************************************************************
- *  Copyright by the contributors to the Dafny Project
- *  SPDX-License-Identifier: MIT
- *******************************************************************************/
-
-module Bernoulli.Interface {
-  import Rationals
-  import Monad
-  import Uniform
-  import Model
-
-  trait {:termination false} Trait extends Uniform.Interface.Trait {
-
-    method BernoulliSample(p: Rationals.Rational) returns (c: bool)
-      modifies `s
-      decreases *
-      requires 0 <= p.numer <= p.denom
-      ensures Model.Sample(p.numer, p.denom)(old(s)) == Monad.Result(c, s)
-
-  }
-}
diff --git a/src/Distributions/Bernoulli/Model.dfy b/src/Distributions/Bernoulli/Model.dfy
deleted file mode 100644
index 7782a987..00000000
--- a/src/Distributions/Bernoulli/Model.dfy
+++ /dev/null
@@ -1,18 +0,0 @@
-/*******************************************************************************
- *  Copyright by the contributors to the Dafny Project
- *  SPDX-License-Identifier: MIT
- *******************************************************************************/
-
-module Bernoulli.Model {
-  import Uniform
-  import Monad
-
-  // Footnote 5, p. 82
-  opaque ghost function Sample(numer: nat, denom: nat): (f: Monad.Hurd<bool>)
-    requires denom != 0
-    requires numer <= denom
-  {
-    Monad.Map(Uniform.Model.Sample(denom), (k: nat) => k < numer)
-  }
-
-}
diff --git a/src/Distributions/BernoulliExpNeg/Correctness.dfy b/src/Distributions/BernoulliExpNeg/Correctness.dfy
deleted file mode 100644
index 4723da43..00000000
--- a/src/Distributions/BernoulliExpNeg/Correctness.dfy
+++ /dev/null
@@ -1,628 +0,0 @@
-/*******************************************************************************
- *  Copyright by the contributors to the Dafny Project
- *  SPDX-License-Identifier: MIT
- *******************************************************************************/
-
-module BernoulliExpNeg.Correctness {
-  import Measures
-  import Rationals
-  import Sequences
-  import Limits
-  import Series
-  import Helper
-  import NatArith
-  import RealArith
-  import Exponential
-  import Rand
-  import Monad
-  import Independence
-  import Loops
-  import Bernoulli
-  import Model
-  import Uniform
-
-  lemma Correctness(gamma: Rationals.Rational)
-    requires 0 <= gamma.numer
-    requires 0 < gamma.denom
-    decreases gamma.numer
-    ensures Rand.prob(iset s | Model.Sample(gamma)(s).Equals(true)) == Exponential.Exp(-gamma.ToReal())
-  {
-    if gamma.numer <= gamma.denom {
-      assert (iset s | Model.Sample(gamma)(s).Equals(true)) == (iset s | Model.SampleLe1(gamma)(s).Equals(true)) by {
-        reveal Model.Sample();
-      }
-      calc {
-        Rand.prob(iset s | Model.SampleLe1(gamma)(s).Equals(true));
-        { CorrectnessLe1(gamma); }
-        Exponential.Exp(-gamma.ToReal());
-      }
-    } else {
-      var trueEvent := iset s | Model.Sample(gamma)(s).Equals(true);
-      var firstSampleTrueEvent := iset s | Model.SampleLe1(Rationals.Int(1))(s).Equals(true);
-      var gamma' := Rationals.Rational(gamma.numer - gamma.denom, gamma.denom);
-      var secondSampleTrueEvent := iset s | Model.Sample(gamma')(s).Equals(true);
-      assert decomposeTrueEvent: trueEvent == firstSampleTrueEvent * Monad.BitstreamsWithRestIn(Model.SampleLe1(Rationals.Int(1)), secondSampleTrueEvent) by {
-        forall s ensures s in trueEvent <==> s in firstSampleTrueEvent && s in Monad.BitstreamsWithRestIn(Model.SampleLe1(Rationals.Int(1)), secondSampleTrueEvent) {
-          SampleTrueIff(gamma, gamma', s, secondSampleTrueEvent);
-        }
-      }
-      assert independence: Rand.prob(trueEvent) == Rand.prob(firstSampleTrueEvent) * Rand.prob(secondSampleTrueEvent) by {
-        assert Independence.IsIndepFunction(Model.SampleLe1(Rationals.Int(1))) by {
-          SampleLe1IsIndep(Rationals.Int(1));
-          Independence.IsIndepImpliesIsIndepFunction(Model.SampleLe1(Rationals.Int(1)));
-        }
-        assert secondSampleTrueEvent in Rand.eventSpace by {
-          assert secondSampleTrueEvent == Monad.BitstreamsWithValueIn(Model.Sample(gamma'), iset{true});
-          SampleIsIndep(gamma');
-          Independence.IsIndepImpliesMeasurablePreImageBool(Model.Sample(gamma'), iset{true});
-        }
-        Independence.ResultsIndependent(
-          Model.SampleLe1(Rationals.Int(1)),
-          iset{true},
-          secondSampleTrueEvent);
-        assert Monad.BitstreamsWithValueIn(Model.SampleLe1(Rationals.Int(1)), iset{true}) == firstSampleTrueEvent;
-        reveal decomposeTrueEvent;
-      }
-      assert firstProb: Rand.prob(firstSampleTrueEvent) == Exponential.Exp(-1.0) by {
-        calc {
-          Rand.prob(firstSampleTrueEvent);
-          { CorrectnessLe1(Rationals.Int(1)); }
-          Exponential.Exp(-1.0);
-        }
-      }
-      assert gammaMinusOne: gamma'.ToReal() == gamma.ToReal() - 1.0 by {
-        RationalMinusOneFact(gamma, gamma');
-      }
-      assert secondProb: Rand.prob(secondSampleTrueEvent) == Exponential.Exp(-gamma.ToReal() + 1.0) by {
-        calc {
-          Rand.prob(secondSampleTrueEvent);
-          { Correctness(gamma'); }
-          Exponential.Exp(-gamma'.ToReal());
-          { reveal gammaMinusOne; }
-          Exponential.Exp(-(gamma.ToReal() - 1.0));
-          Exponential.Exp(-gamma.ToReal() + 1.0);
-        }
-      }
-      calc {
-        Rand.prob(iset s | Model.Sample(gamma)(s).Equals(true));
-        Rand.prob(trueEvent);
-        { reveal independence; }
-        Rand.prob(firstSampleTrueEvent) * Rand.prob(secondSampleTrueEvent);
-        { reveal firstProb; reveal secondProb; }
-        Exponential.Exp(-1.0) * Exponential.Exp(-gamma.ToReal() + 1.0);
-        { Exponential.FunctionalEquation(-1.0, -gamma.ToReal() + 1.0); }
-        Exponential.Exp(-gamma.ToReal());
-      }
-    }
-  }
-
-  lemma RationalMinusOneFact(gamma: Rationals.Rational, gamma': Rationals.Rational)
-    requires gamma.numer > gamma.denom
-    requires gamma.denom == gamma'.denom
-    requires gamma'.numer == gamma.numer - gamma.denom
-    ensures gamma'.ToReal() == gamma.ToReal() - 1.0
-  {
-    var numer := gamma.numer as real;
-    var denom := gamma.denom as real;
-    assert denom != 0.0;
-    calc {
-      gamma'.ToReal();
-      gamma'.numer as real / gamma'.denom as real;
-      (gamma.numer - gamma.denom) as real / denom;
-      (numer - denom) / denom;
-      { RealArith.DivideSubtraction(numer, denom, denom); }
-      numer / denom - denom / denom;
-      { RealArith.FractionAsOne(denom); }
-      numer / denom - 1.0;
-      gamma.ToReal() - 1.0;
-    }
-  }
-
-  lemma SampleTrueIff(gamma: Rationals.Rational, gamma': Rationals.Rational, s: Rand.Bitstream, secondSampleTrueEvent: iset<Rand.Bitstream>)
-    requires 0 < gamma.denom < gamma.numer
-    requires gamma'.denom == gamma.denom
-    requires gamma'.numer == gamma.numer - gamma.denom
-    requires forall s' :: s' in secondSampleTrueEvent <==> Model.Sample(gamma')(s').Equals(true)
-    ensures Model.Sample(gamma)(s).Equals(true) <==> Model.SampleLe1(Rationals.Int(1))(s).Equals(true) && Model.SampleLe1(Rationals.Int(1))(s).RestIn(secondSampleTrueEvent)
-  {
-    reveal Model.Sample();
-    if Model.SampleLe1(Rationals.Int(1))(s).Equals(true) {
-      var s' := Model.SampleLe1(Rationals.Int(1))(s).rest;
-      assert Model.Sample(gamma)(s) == Model.Sample(gamma')(s');
-      assert Model.Sample(gamma')(s').Equals(true) <==> s' in secondSampleTrueEvent;
-    } else {
-      assert !Model.Sample(gamma)(s).Equals(true);
-    }
-  }
-
-  lemma SampleIsIndep(gamma: Rationals.Rational)
-    requires 0 <= gamma.numer
-    decreases gamma.numer
-    ensures Independence.IsIndep(Model.Sample(gamma))
-  {
-    reveal Model.Sample();
-    if gamma.numer <= gamma.denom {
-      SampleLe1IsIndep(gamma);
-    } else {
-      SampleLe1IsIndep(Rationals.Int(1));
-      var gamma' := Rationals.Rational(gamma.numer - gamma.denom, gamma.denom);
-      var f: bool -> Monad.Hurd<bool> :=
-        (b: bool) =>
-          var res: Monad.Hurd<bool> :=
-            if b
-            then Model.Sample(Rationals.Rational(gamma.numer - gamma.denom, gamma.denom))
-            else Monad.Return(false);
-          res;
-      forall b: bool ensures Independence.IsIndep(f(b)) {
-        if b {
-          SampleIsIndep(gamma');
-        } else {
-          Independence.ReturnIsIndep(false);
-        }
-      }
-      assert Independence.IsIndep(Model.Sample(gamma)) by {
-        assert Model.Sample(gamma) == Monad.Bind(Model.SampleLe1(Rationals.Int(1)), f);
-        Independence.BindIsIndep(Model.SampleLe1(Rationals.Int(1)), f);
-      }
-    }
-  }
-
-  // Proves the correctness of Model.SampleLe1
-  lemma CorrectnessLe1(gamma: Rationals.Rational)
-    requires 0 <= gamma.numer <= gamma.denom
-    ensures Rand.prob(iset s | Model.SampleLe1(gamma)(s).Equals(true)) == Exponential.Exp(-gamma.ToReal())
-  {
-    var event := Monad.BitstreamsWithValueIn(Model.SampleLe1(gamma), iset{true});
-    var resultSeq := (k: nat) => if k % 2 == 0 then iset{} else iset ak: (bool, nat) | ak.1 == k;
-    var resultSet := Measures.CountableUnion(resultSeq);
-    var partEvent := (k: nat) => Monad.BitstreamsWithValueIn(Model.Le1Loop(gamma)((true, 0)), resultSeq(k));
-    var probPartEvent := (k: nat) => Rand.prob(partEvent(k));
-    // The strategy is to decompose event into parts depending on the `k` returned from `Model.Le1Loop`.
-    // Then we sum over all possible `k`s and obtain a series converging to the exponential function at -gamma.
-    assert decomposeEvent: event == Measures.CountableUnion(partEvent) by {
-      forall s ensures s in event <==> s in Measures.CountableUnion(partEvent) {
-        var oddNats := iset ak: (bool, nat) | ak.1 % 2 == 1;
-        var resultSeqUnion := iset k: nat, ak <- resultSeq(k) :: ak;
-        assert oddNatsEquiv: oddNats == resultSeqUnion by {
-          forall ak: (bool, nat) ensures ak in oddNats <==> ak in resultSeqUnion {
-            calc {
-              ak in oddNats;
-              ak.1 % 2 == 1;
-              ak in resultSeq(ak.1);
-            }
-          }
-        }
-        assert partEventCharacterization: forall k: nat :: s in partEvent(k) <==> Model.Le1Loop(gamma)((true, 0))(s).In(resultSeq(k)) by {
-          forall k: nat ensures s in partEvent(k) <==> Model.Le1Loop(gamma)((true, 0))(s).In(resultSeq(k)) {}
-        }
-        calc {
-          s in event;
-          Model.SampleLe1(gamma)(s).Equals(true);
-          Model.Le1Loop(gamma)((true, 0))(s).Satisfies((ak: (bool, nat)) => ak.1 % 2 == 1);
-          Model.Le1Loop(gamma)((true, 0))(s).In(oddNats);
-          { reveal oddNatsEquiv; }
-          Model.Le1Loop(gamma)((true, 0))(s).In(resultSeqUnion);
-          exists k: nat :: Model.Le1Loop(gamma)((true, 0))(s).In(resultSeq(k));
-          { reveal partEventCharacterization; }
-          exists k: nat :: s in partEvent(k);
-        }
-      }
-    }
-    assert probPart: forall k: nat :: probPartEvent(k) == Le1SeriesTerm(gamma.ToReal())(k) by {
-      forall k: nat ensures probPartEvent(k) == Le1SeriesTerm(gamma.ToReal())(k) {
-        if k % 2 == 0 {
-          assert partEvent(k) == iset{};
-          assert probPartEvent(k) == 0.0 by {
-            Rand.ProbIsProbabilityMeasure();
-          }
-        } else {
-          assert probPartEvent(k) == Le1SeriesTerm(gamma.ToReal())(k) by {
-            Le1LoopCorrectnessEq(gamma, k);
-          }
-        }
-      }
-    }
-    assert probPartsSum1: Series.SumsTo(probPartEvent, Rand.prob(event)) by {
-      reveal decomposeEvent;
-      Rand.ProbIsProbabilityMeasure();
-      // TODO: prove measurability
-      assume {:axiom} forall k: nat :: partEvent(k) in Rand.eventSpace;
-      assert Measures.PairwiseDisjoint(partEvent) by {
-        forall m: nat, n: nat | m != n ensures partEvent(m) * partEvent(n) == iset{} {}
-      }
-      Measures.MeasureOfCountableDisjointUnionIsSum(Rand.eventSpace, Rand.prob, partEvent, probPartEvent);
-    }
-    assert probPartsSum2: Series.SumsTo(probPartEvent, Exponential.Exp(-gamma.ToReal())) by {
-      assert Series.SumsTo(Le1SeriesTerm(gamma.ToReal()), Exponential.Exp(-gamma.ToReal())) by {
-        Le1SeriesConvergesToExpNeg(gamma.ToReal());
-      }
-      reveal probPart;
-      Series.SumOfEqualsIsEqual(Le1SeriesTerm(gamma.ToReal()), probPartEvent, Exponential.Exp(-gamma.ToReal()));
-    }
-    assert Rand.prob(iset s | Model.SampleLe1(gamma)(s).Equals(true)) == Exponential.Exp(-gamma.ToReal()) by {
-      assert Rand.prob(iset s | Model.SampleLe1(gamma)(s).Equals(true)) == Rand.prob(event) by {
-        assert event == iset s | Model.SampleLe1(gamma)(s).Equals(true);
-      }
-      reveal probPartsSum1;
-      reveal probPartsSum2;
-      Series.SumIsUnique(probPartEvent, Rand.prob(event), Exponential.Exp(-gamma.ToReal()));
-    }
-  }
-
-  lemma SampleLe1IsIndep(gamma: Rationals.Rational)
-    requires 0 <= gamma.numer <= gamma.denom
-    ensures Independence.IsIndep(Model.SampleLe1(gamma))
-  {
-    Le1LoopIsIndep(gamma);
-    Independence.MapIsIndep(Model.Le1Loop(gamma)((true, 0)), (ak: (bool, nat)) => ak.1 % 2 == 1);
-  }
-
-  lemma Le1LoopIsIndep(gamma: Rationals.Rational)
-    requires 0 <= gamma.numer <= gamma.denom
-    ensures Independence.IsIndep(Model.Le1Loop(gamma)((true, 0)))
-  {
-    Model.Le1LoopTerminatesAlmostSurely(gamma);
-    forall ak: (bool, nat) ensures Independence.IsIndep(Model.Le1LoopIter(gamma)(ak)) {
-      Le1LoopIterIsIndep(gamma, ak);
-    }
-    Loops.WhileIsIndep(
-      Model.Le1LoopCondition,
-      Model.Le1LoopIter(gamma),
-      (true, 0)
-    );
-    assert Independence.IsIndep(Model.Le1Loop(gamma)((true, 0))) by {
-      reveal Model.Le1Loop();
-    }
-  }
-
-  lemma Le1LoopIterIsIndep(gamma: Rationals.Rational, ak: (bool, nat))
-    requires 0 <= gamma.numer <= gamma.denom
-    ensures Independence.IsIndep(Model.Le1LoopIter(gamma)(ak))
-  {
-    var k': nat := ak.1 + 1;
-    Bernoulli.Correctness.SampleIsIndep(gamma.numer, k' * gamma.denom);
-    Independence.MapIsIndep(
-      Bernoulli.Model.Sample(gamma.numer, k' * gamma.denom),
-      a => (a, k')
-    );
-  }
-
-  // Proves the correctness of `Model.Le1LoopIter`.
-  // Correctness means that when run on an initial value of `(true, k)` it produces
-  // (true, k + 1) with probability gamma / (k + 1)
-  // (false, k + 1) with probability 1 - gamma / (k + 1)
-  // and (implicitly) all other probabilities are zero.
-  lemma Le1LoopIterCorrectness(gamma: Rationals.Rational, k: nat)
-    requires 0 <= gamma.numer <= gamma.denom
-    ensures Rand.prob(Monad.BitstreamsWithValueIn(Model.Le1LoopIter(gamma)((true, k)), iset{(true, k + 1)})) == gamma.ToReal() / (k + 1) as real
-    ensures Rand.prob(Monad.BitstreamsWithValueIn(Model.Le1LoopIter(gamma)((true, k)), iset{(false, k + 1)})) == 1.0 - gamma.ToReal() / (k + 1) as real
-  {
-    var k': nat := k + 1;
-    var denom: nat := k' * gamma.denom;
-    assert denom >= 1;
-    var eventTrue := Monad.BitstreamsWithValueIn(Model.Le1LoopIter(gamma)((true, k)), iset{(true, k')});
-    var eventTrue2 := iset s | Bernoulli.Model.Sample(gamma.numer, denom)(s).Equals(true);
-    assert eventTrue == eventTrue2 by {
-      forall s
-        ensures s in eventTrue <==> s in eventTrue2
-      {
-        reveal Bernoulli.Model.Sample();
-        if Uniform.Model.Sample(denom)(s).Result? {
-          Uniform.Model.SampleBound(denom, s);
-        }
-      }
-    }
-    assert Rand.prob(eventTrue2) == gamma.numer as real / denom as real by {
-      Bernoulli.Correctness.BernoulliCorrectness(gamma.numer, denom, true);
-      reveal Bernoulli.Correctness.BernoulliMass();
-    }
-    var eventFalse := Monad.BitstreamsWithValueIn(Model.Le1LoopIter(gamma)((true, k)), iset{(false, k')});
-    var eventFalse2 := iset s | Bernoulli.Model.Sample(gamma.numer, denom)(s).Equals(false);
-    assert eventFalse == eventFalse2;
-    assert Rand.prob(eventFalse2) == 1.0 - gamma.numer as real / denom as real by {
-      Bernoulli.Correctness.BernoulliCorrectness(gamma.numer, denom, false);
-      reveal Bernoulli.Correctness.BernoulliMass();
-    }
-    assert gamma.numer as real / denom as real == gamma.ToReal() / k' as real by {
-      calc {
-        gamma.ToReal() / k' as real;
-        (gamma.numer as real / gamma.denom as real) / k' as real;
-        { RealArith.DivMulEqualsDivDiv(gamma.numer as real, gamma.denom as real, k' as real); }
-        gamma.numer as real / (gamma.denom as real * k' as real);
-        gamma.numer as real / (denom as real);
-      }
-    }
-  }
-
-  // Sequence of terms of the infinite series used in the correctness proof.
-  // It's like the exponential series of exp(-gamma) except that pairs of terms are moved to the odd positions,
-  // so it's 0.0 + (ExpTerm(gamma, 0) - ExpTerm(gamma, 1)) + 0.0 + (ExpTerm(gamma, 2) - ExpTerm(gamma, 3))
-  function Le1SeriesTerm(gamma: real): nat -> real {
-    (n: nat) => if n % 2 == 0 then 0.0 else Exponential.ExpTerm(gamma, n - 1) - Exponential.ExpTerm(gamma, n)
-  }
-
-  // Proves that the Le1SeriesTerm sequence is summable and has the same sum as the normal exponential series.
-  // On paper the proof is trivial (just move the summands around), but it's not as easy to formalize.
-  lemma {:axiom} Le1SeriesConvergesToExpNeg(gamma: real)
-    ensures Series.SumsTo(Le1SeriesTerm(gamma), Exponential.Exp(-gamma))
-
-  // A bounded version of `Model.Le1Loop`.
-  opaque ghost function Le1LoopCut(gamma: Rationals.Rational, ak: (bool, nat)): nat -> Monad.Hurd<(bool, nat)>
-    requires 0 <= gamma.numer <= gamma.denom
-  {
-    (fuel: nat) => Loops.WhileCut(
-        Model.Le1LoopCondition,
-        Model.Le1LoopIter(gamma),
-        ak,
-        fuel
-      )
-  }
-
-  // Decomposes the probability of `Le1LoopCut(gamma, (true, k))(fuel)` producing a value (false, m) with m <= k + n,
-  // by branching on the result of the sample `a` from the first loop iteration.
-  lemma Le1LoopCutDecomposeProb(gamma: Rationals.Rational, k: nat, n: int, fuel: nat)
-    requires 0 <= gamma.numer <= gamma.denom
-    requires fuel >= 1
-    ensures
-      Rand.prob(Monad.BitstreamsWithValueIn(Le1LoopCut(gamma, (true, k))(fuel), iset m: nat | m <= k + n :: (false, m)))
-      ==
-      Rand.prob(Monad.BitstreamsWithValueIn(Model.Le1LoopIter(gamma)((true, k)), iset{(false, k + 1)}))
-      * Rand.prob(Monad.BitstreamsWithValueIn(Le1LoopCut(gamma, (false, k + 1))(fuel - 1), iset m: nat | m <= k + n :: (false, m)))
-      + Rand.prob(Monad.BitstreamsWithValueIn(Model.Le1LoopIter(gamma)((true, k)), iset{(true, k + 1)}))
-      * Rand.prob(Monad.BitstreamsWithValueIn(Le1LoopCut(gamma, (true, k + 1))(fuel - 1), iset m: nat | m <= k + n :: (false, m)))
-  {
-    var resultSet := iset m: nat | m <= k + n :: (false, m);
-    var init: (bool, nat) := (true, k);
-    var event := Monad.BitstreamsWithValueIn(Le1LoopCut(gamma, init)(fuel), resultSet);
-    var k': nat := k + 1;
-    // first loop iteration returns `a == true`
-    var firstIterTrue := Monad.BitstreamsWithValueIn(Model.Le1LoopIter(gamma)(init), iset m: nat :: (true, m));
-    var firstIterTrue2 := Monad.BitstreamsWithValueIn(Model.Le1LoopIter(gamma)(init), iset{(true, k + 1)});
-    assert firstIterTrue == firstIterTrue2 by {
-      forall s ensures s in Monad.BitstreamsWithValueIn(Model.Le1LoopIter(gamma)(init), iset m: nat :: (true, m)) <==> s in Monad.BitstreamsWithValueIn(Model.Le1LoopIter(gamma)(init), iset{(true, k + 1)}) {}
-    }
-    var resultAfterFirstTrue := Monad.BitstreamsWithValueIn(Le1LoopCut(gamma, (true, k + 1))(fuel - 1), resultSet);
-    var seedsWithResultAfterFirstTrue := Monad.BitstreamsWithRestIn(Model.Le1LoopIter(gamma)(init), resultAfterFirstTrue);
-    // first loop iteration returns `a == false`
-    var firstIterFalse := Monad.BitstreamsWithValueIn(Model.Le1LoopIter(gamma)(init), iset m: nat :: (false, m));
-    var firstIterFalse2 := Monad.BitstreamsWithValueIn(Model.Le1LoopIter(gamma)(init), iset{(false, k + 1)});
-    assert firstIterFalse == firstIterFalse2 by {
-      forall s ensures s in Monad.BitstreamsWithValueIn(Model.Le1LoopIter(gamma)(init), iset m: nat :: (false, m)) <==> s in Monad.BitstreamsWithValueIn(Model.Le1LoopIter(gamma)(init), iset{(false, k + 1)}) {}
-    }
-    var resultAfterFirstFalse := Monad.BitstreamsWithValueIn(Le1LoopCut(gamma, (false, k + 1))(fuel - 1), resultSet);
-    var seedsWithResultAfterFirstFalse := Monad.BitstreamsWithRestIn(Model.Le1LoopIter(gamma)(init), resultAfterFirstFalse);
-    assert decomposeEvent: event == firstIterFalse * seedsWithResultAfterFirstFalse + firstIterTrue * seedsWithResultAfterFirstTrue by {
-      forall s ensures s in event <==> (s in firstIterFalse2 && s in seedsWithResultAfterFirstFalse) || (s in firstIterTrue2 && s in seedsWithResultAfterFirstTrue) {
-        reveal Le1LoopCut();
-        match Model.Le1LoopIter(gamma)(init)(s)
-        case Diverging => assert Le1LoopCut(gamma, init)(fuel)(s).Diverging?;
-        case Result((a', k'), s') =>
-          assert Le1LoopCut(gamma, init)(fuel)(s) == Le1LoopCut(gamma, (a', k'))(fuel - 1)(s');
-          calc {
-            s in event;
-            Le1LoopCut(gamma, init)(fuel)(s).In(resultSet);
-            (s in firstIterFalse2 && Le1LoopCut(gamma, (a', k'))(fuel - 1)(s').In(resultSet)) || (s in firstIterTrue2 && Le1LoopCut(gamma, (a', k'))(fuel - 1)(s').In(resultSet));
-            (s in firstIterFalse2 && a' == false && k' == k + 1 && Le1LoopCut(gamma, (a', k'))(fuel - 1)(s').In(resultSet)) || (s in firstIterTrue2 && a' == true && k' == k + 1 && Le1LoopCut(gamma, (a', k'))(fuel - 1)(s').In(resultSet));
-            (s in firstIterFalse2 && s in seedsWithResultAfterFirstFalse) || (s in firstIterTrue2 && s in seedsWithResultAfterFirstTrue);
-          }
-      }
-    }
-    assert decomposeProb: Rand.prob(event) == Rand.prob(firstIterFalse * seedsWithResultAfterFirstFalse) + Rand.prob(firstIterTrue * seedsWithResultAfterFirstTrue) by {
-      assume {:axiom} firstIterFalse * seedsWithResultAfterFirstFalse in Rand.eventSpace;
-      assume {:axiom} firstIterTrue * seedsWithResultAfterFirstTrue in Rand.eventSpace;
-      Rand.ProbIsProbabilityMeasure();
-      assert (firstIterFalse * seedsWithResultAfterFirstFalse) * (firstIterTrue * seedsWithResultAfterFirstTrue) == iset{};
-      Measures.MeasureOfDisjointUnionIsSum(Rand.eventSpace, Rand.prob, firstIterFalse * seedsWithResultAfterFirstFalse, firstIterTrue * seedsWithResultAfterFirstTrue);
-      reveal decomposeEvent;
-    }
-    assert trueIndependent: Rand.prob(firstIterTrue * seedsWithResultAfterFirstTrue) == Rand.prob(firstIterTrue) * Rand.prob(resultAfterFirstTrue) by {
-      assume {:axiom} Independence.IsIndepFunction(Model.Le1LoopIter(gamma)(init));
-      assume {:axiom} resultAfterFirstTrue in Rand.eventSpace;
-      Independence.ResultsIndependent(Model.Le1LoopIter(gamma)(init), iset{(true, k + 1)}, resultAfterFirstTrue);
-    }
-    assert falseIndependent: Rand.prob(firstIterFalse * seedsWithResultAfterFirstFalse) == Rand.prob(firstIterFalse) * Rand.prob(resultAfterFirstFalse) by {
-      assume {:axiom} Independence.IsIndepFunction(Model.Le1LoopIter(gamma)(init));
-      assume {:axiom} resultAfterFirstFalse in Rand.eventSpace;
-      Independence.ResultsIndependent(Model.Le1LoopIter(gamma)(init), iset{(false, k + 1)}, resultAfterFirstFalse);
-    }
-    calc {
-      Rand.prob(event);
-      { reveal decomposeProb; }
-      Rand.prob(firstIterFalse * seedsWithResultAfterFirstFalse) + Rand.prob(firstIterTrue * seedsWithResultAfterFirstTrue);
-      { reveal falseIndependent; reveal trueIndependent; }
-      Rand.prob(firstIterFalse) * Rand.prob(resultAfterFirstFalse) + Rand.prob(firstIterTrue) * Rand.prob(resultAfterFirstTrue);
-    }
-  }
-
-  // Proves correctness of `Le1LoopCut`:
-  // i.e. that the probability of `Le1LoopCut(gamma, (true, k))(fuel)` producing a value (false, m) with m <= k + n is
-  // (a) 0 if n <= 0 (intuitvely because `k` only increases, so the bound of `k + n` is already passed)
-  // (b) 1.0 - gamma^l / (k + 1)^(l) where l = min(n, fuel) and (k + 1)^(l) is the rising factorial (k + 1) * (k + 2) * ... * (k + l)
-  lemma Le1LoopCutCorrectness(gamma: Rationals.Rational, k: nat, n: int, fuel: nat)
-    decreases fuel
-    requires 0 <= gamma.numer <= gamma.denom
-    ensures
-      Rand.prob(Monad.BitstreamsWithValueIn(Le1LoopCut(gamma, (true, k))(fuel), iset m: nat | m <= k + n :: (false, m)))
-      == if n <= 0 then 0.0 else 1.0 - Exponential.ExpTerm(gamma.ToReal(), NatArith.Min(n, fuel), k + 1)
-  {
-    var resultSet := iset m: nat | m <= k + n :: (false, m);
-    var event := Monad.BitstreamsWithValueIn(Le1LoopCut(gamma, (true, k))(fuel), resultSet);
-    if fuel == 0 {
-      assert event == iset{} by {
-        forall s ensures s !in event {
-          reveal Le1LoopCut();
-          assert !Le1LoopCut(gamma, (true, k))(fuel)(s).In(resultSet);
-        }
-      }
-      assert Rand.prob(event) == 0.0 by {
-        Rand.ProbIsProbabilityMeasure();
-      }
-      assert Rand.prob(event) == if n <= 0 then 0.0 else 1.0 - Exponential.ExpTerm(gamma.ToReal(), NatArith.Min(n, fuel), k + 1) by {
-        reveal Exponential.ExpTerm();
-        reveal RealArith.Pow();
-        reveal NatArith.Factorial();
-      }
-    } else {
-      var k': nat := k + 1;
-      var firstIterTrue := Monad.BitstreamsWithValueIn(Model.Le1LoopIter(gamma)((true, k)), iset{(true, k + 1)});
-      var firstIterFalse := Monad.BitstreamsWithValueIn(Model.Le1LoopIter(gamma)((true, k)), iset{(false, k + 1)});
-      var desiredResultAfterFirstIterTrue := Monad.BitstreamsWithValueIn(Le1LoopCut(gamma, (true, k + 1))(fuel - 1), resultSet);
-      var desiredResultAfterFirstIterFalse := Monad.BitstreamsWithValueIn(Le1LoopCut(gamma, (false, k + 1))(fuel - 1), resultSet);
-
-      assert probFirstFalse: Rand.prob(firstIterFalse) == 1.0 - gamma.ToReal() / k' as real by {
-        Le1LoopIterCorrectness(gamma, k);
-      }
-      assert probResultAfterFirstFalse: Rand.prob(desiredResultAfterFirstIterFalse) == if 0 < n then 1.0 else 0.0 by {
-        if 0 < n {
-          assert desiredResultAfterFirstIterFalse == Measures.SampleSpace() by {
-            forall s ensures s in desiredResultAfterFirstIterFalse {
-              reveal Le1LoopCut();
-              assert Le1LoopCut(gamma, (false, k + 1))(fuel - 1)(s) == Monad.Return((false, k + 1))(s);
-            }
-          }
-          assert Rand.prob(desiredResultAfterFirstIterFalse) == 1.0 by {
-            Rand.ProbIsProbabilityMeasure();
-          }
-        } else {
-          assert desiredResultAfterFirstIterFalse == iset{} by {
-            forall s ensures s !in desiredResultAfterFirstIterFalse {
-              reveal Le1LoopCut();
-              assert Le1LoopCut(gamma, (false, k + 1))(fuel - 1)(s) == Monad.Return((false, k + 1))(s);
-              assert (false, k + 1) !in resultSet;
-            }
-          }
-          assert Rand.prob(desiredResultAfterFirstIterFalse) == 0.0 by {
-            Rand.ProbIsProbabilityMeasure();
-          }
-        }
-      }
-      assert probFirstTrue: Rand.prob(firstIterTrue) == gamma.ToReal() / k' as real by {
-        Le1LoopIterCorrectness(gamma, k);
-      }
-      assert probResultAfterFirstTrue: Rand.prob(desiredResultAfterFirstIterTrue) == if n <= 1 then 0.0 else 1.0 - Exponential.ExpTerm(gamma.ToReal(), NatArith.Min(n, fuel) - 1, k' + 1) by {
-        Le1LoopCutCorrectness(gamma, k', n - 1, fuel - 1);
-        assert n >= 1 ==> NatArith.Min(n, fuel) - 1 == NatArith.Min(n - 1, fuel - 1);
-      }
-      if n <= 0 {
-        calc {
-          Rand.prob(event);
-          { Le1LoopCutDecomposeProb(gamma, k, n, fuel); }
-          Rand.prob(firstIterFalse) * Rand.prob(desiredResultAfterFirstIterFalse) + Rand.prob(firstIterTrue) * Rand.prob(desiredResultAfterFirstIterTrue);
-          { reveal probResultAfterFirstFalse; reveal probResultAfterFirstTrue; }
-          0.0;
-        }
-      } else {
-        if n == 1 {
-          NatArith.FactoralPositive(1, k');
-          calc {
-            Rand.prob(event);
-            { Le1LoopCutDecomposeProb(gamma, k, n, fuel); }
-            Rand.prob(firstIterFalse) * Rand.prob(desiredResultAfterFirstIterFalse) + Rand.prob(firstIterTrue) * Rand.prob(desiredResultAfterFirstIterTrue);
-            { reveal probFirstFalse; reveal probResultAfterFirstFalse; reveal probFirstTrue; reveal probResultAfterFirstTrue; }
-            1.0 - gamma.ToReal() / k' as real;
-            { reveal RealArith.Pow(); reveal NatArith.Factorial(); }
-            1.0 - RealArith.Pow(gamma.ToReal(), 1) / NatArith.Factorial(1, k') as real;
-            { reveal Exponential.ExpTerm(); }
-            1.0 - Exponential.ExpTerm(gamma.ToReal(), n, k + 1);
-          }
-          assert Rand.prob(event) == if n <= 0 then 0.0 else 1.0 - Exponential.ExpTerm(gamma.ToReal(), NatArith.Min(n, fuel), k + 1);
-        } else {
-          calc {
-            Rand.prob(event);
-            { Le1LoopCutDecomposeProb(gamma, k, n, fuel); }
-            Rand.prob(firstIterFalse) * Rand.prob(desiredResultAfterFirstIterFalse) + Rand.prob(firstIterTrue) * Rand.prob(desiredResultAfterFirstIterTrue);
-            { reveal probFirstFalse; reveal probResultAfterFirstFalse; reveal probFirstTrue; reveal probResultAfterFirstTrue; }
-            (1.0 - gamma.ToReal() / k' as real) * 1.0 + (gamma.ToReal() / k' as real) * (1.0 - Exponential.ExpTerm(gamma.ToReal(), NatArith.Min(n, fuel) - 1, k' + 1));
-            1.0 - gamma.ToReal() / k' as real * Exponential.ExpTerm(gamma.ToReal(), NatArith.Min(n, fuel) - 1, k' + 1);
-            { Exponential.ExpTermStep(gamma.ToReal(), NatArith.Min(n, fuel), k'); }
-            1.0 - Exponential.ExpTerm(gamma.ToReal(), NatArith.Min(n, fuel), k + 1);
-          }
-        }
-        assert Rand.prob(event) == if n <= 0 then 0.0 else 1.0 - Exponential.ExpTerm(gamma.ToReal(), NatArith.Min(n, fuel), k + 1);
-      }
-      assert Rand.prob(event) == if n <= 0 then 0.0 else 1.0 - Exponential.ExpTerm(gamma.ToReal(), NatArith.Min(n, fuel), k + 1);
-    }
-    assert Rand.prob(event) == if n <= 0 then 0.0 else 1.0 - Exponential.ExpTerm(gamma.ToReal(), NatArith.Min(n, fuel), k + 1);
-  }
-
-  // Proves the first version of correctness of `Model.Le1Loop`:
-  // i.e. that the probability of `Model.Le1Loop(gamma)((true, 0))` producing a value (_, m) with m <= n is
-  // 1.0 - gamma^n / n!
-  lemma Le1LoopCorrectnessLe(gamma: Rationals.Rational, n: nat)
-    requires 0 <= gamma.numer <= gamma.denom
-    ensures
-      Rand.prob(Monad.BitstreamsWithValueIn(Model.Le1Loop(gamma)(((true, 0))), iset ak: (bool, nat) | ak.1 <= n))
-      == 1.0 - Exponential.ExpTerm(gamma.ToReal(), n)
-  {
-    var resultSet := iset ak: (bool, nat) | ak.1 <= n;
-    var resultSetRestricted := iset m: nat | m <= n :: (false, m);
-    var limit := 1.0 - Exponential.ExpTerm(gamma.ToReal(), n);
-    assert resultSetRestricted == iset a <- resultSet | !Model.Le1LoopCondition(a);
-    var sequence: nat -> real := Loops.WhileCutProbability(Model.Le1LoopCondition, Model.Le1LoopIter(gamma), (true, 0), resultSetRestricted);
-    assert Limits.ConvergesTo(sequence, limit) by {
-      forall fuel: nat | fuel >= n ensures sequence(fuel) == limit {
-        calc {
-          sequence(fuel);
-          Loops.WhileCutProbability(Model.Le1LoopCondition, Model.Le1LoopIter(gamma), (true, 0), resultSetRestricted)(fuel);
-          { reveal Le1LoopCut(); }
-          Rand.prob(Monad.BitstreamsWithValueIn(Le1LoopCut(gamma, (true, 0))(fuel), resultSetRestricted));
-          { Le1LoopCutCorrectness(gamma, 0, n, fuel); reveal Exponential.ExpTerm(); reveal RealArith.Pow(); reveal NatArith.Factorial(); }
-          1.0 - Exponential.ExpTerm(gamma.ToReal(), NatArith.Min(n, fuel));
-          limit;
-        }
-      }
-      Limits.ConstantSequenceConverges(sequence, limit, n);
-    }
-    reveal Model.Le1Loop();
-    assert Rand.prob(Monad.BitstreamsWithValueIn(Model.Le1Loop(gamma)(((true, 0))), resultSet)) == limit by {
-      Loops.WhileProbabilityViaLimit(Model.Le1LoopCondition, Model.Le1LoopIter(gamma), (true, 0), resultSet, resultSetRestricted, limit);
-    }
-  }
-
-  // Proves the second version of correctness of `Model.Le1Loop`:
-  // i.e. that the probability of `Model.Le1Loop(gamma)((true, 0))` producing a value (_, n) is
-  // 0 if n == 0
-  // gamma^(n - 1) / (n - 1)! - gamma^n / n!
-  lemma Le1LoopCorrectnessEq(gamma: Rationals.Rational, n: nat := 0)
-    requires 0 <= gamma.numer <= gamma.denom
-    ensures
-      Rand.prob(Monad.BitstreamsWithValueIn(Model.Le1Loop(gamma)(((true, 0))), iset ak: (bool, nat) | ak.1 == n))
-      == if n == 0 then 0.0 else Exponential.ExpTerm(gamma.ToReal(), n - 1) - Exponential.ExpTerm(gamma.ToReal(), n)
-  {
-    var resultEqN := iset ak: (bool, nat) | ak.1 == n;
-    var eventEqN := Monad.BitstreamsWithValueIn(Model.Le1Loop(gamma)(((true, 0))), resultEqN);
-    var resultLeN := iset ak: (bool, nat) | ak.1 <= n;
-    var eventLeN := Monad.BitstreamsWithValueIn(Model.Le1Loop(gamma)(((true, 0))), resultLeN);
-    if n == 0 {
-      assert resultEqN == resultLeN;
-      calc {
-        Rand.prob(Monad.BitstreamsWithValueIn(Model.Le1Loop(gamma)(((true, 0))), resultEqN));
-        Rand.prob(Monad.BitstreamsWithValueIn(Model.Le1Loop(gamma)(((true, 0))), resultLeN));
-        { Le1LoopCorrectnessLe(gamma, n); }
-        1.0 - Exponential.ExpTerm(gamma.ToReal(), n);
-        { reveal Exponential.ExpTerm(); reveal RealArith.Pow(); reveal NatArith.Factorial(); }
-        0.0;
-      }
-    } else {
-      var resultLeN1 := iset ak: (bool, nat) | ak.1 <= n - 1;
-      assert resultLeN1 * resultEqN == iset{};
-      assert resultLeN1 + resultEqN == resultLeN;
-      var eventLeN1 := Monad.BitstreamsWithValueIn(Model.Le1Loop(gamma)(((true, 0))), resultLeN1);
-      assert eventLeN1 * eventEqN == iset{};
-      assert eventLeN1 + eventEqN == eventLeN;
-      assert Rand.prob(eventLeN) == Rand.prob(eventLeN1) + Rand.prob(eventEqN) by {
-        Rand.ProbIsProbabilityMeasure();
-        assume {:axiom} eventLeN1 in Rand.eventSpace;
-        assume {:axiom} eventEqN in Rand.eventSpace;
-        Measures.MeasureOfDisjointUnionIsSum(Rand.eventSpace, Rand.prob, eventLeN1, eventEqN);
-      }
-      calc {
-        Rand.prob(eventEqN);
-        Rand.prob(eventLeN) - Rand.prob(eventLeN1);
-        { Le1LoopCorrectnessLe(gamma, n); Le1LoopCorrectnessLe(gamma, n - 1); }
-        (1.0 - Exponential.ExpTerm(gamma.ToReal(), n)) - (1.0 - Exponential.ExpTerm(gamma.ToReal(), n - 1));
-        Exponential.ExpTerm(gamma.ToReal(), n - 1) - Exponential.ExpTerm(gamma.ToReal(), n);
-      }
-    }
-  }
-}
diff --git a/src/Distributions/BernoulliExpNeg/Equivalence.dfy b/src/Distributions/BernoulliExpNeg/Equivalence.dfy
deleted file mode 100644
index 579ce284..00000000
--- a/src/Distributions/BernoulliExpNeg/Equivalence.dfy
+++ /dev/null
@@ -1,95 +0,0 @@
-/*******************************************************************************
- *  Copyright by the contributors to the Dafny Project
- *  SPDX-License-Identifier: MIT
- *******************************************************************************/
-
-module BernoulliExpNeg.Equivalence {
-  import Rationals
-  import Rand
-  import Monad
-  import Model
-  import Loops
-  import Bernoulli
-
-  /************
-   Definitions
-  ************/
-
-  opaque ghost predicate CaseLe1LoopInvariant(gamma: Rationals.Rational, oldS: Rand.Bitstream, a: bool, k: nat, s: Rand.Bitstream)
-    requires 0 <= gamma.numer <= gamma.denom
-  {
-    Model.Le1Loop(gamma)((true, 0))(oldS) == Model.Le1Loop(gamma)((a, k))(s)
-  }
-
-  /*******
-   Lemmas
-  *******/
-
-  lemma SampleUnfold(gamma: Rationals.Rational, s: Rand.Bitstream, prevGamma: Rationals.Rational, prevS: Rand.Bitstream, b: bool)
-    requires gamma.numer > 0
-    requires prevGamma.denom == gamma.denom
-    requires prevGamma.numer == gamma.numer + gamma.denom
-    requires Model.SampleLe1(Rationals.Int(1))(prevS) == Monad.Result(b, s)
-    ensures Model.Sample(prevGamma)(prevS) == if b then Model.Sample(gamma)(s) else Monad.Result(false, s)
-  {
-    reveal Model.Sample();
-  }
-
-  lemma Le1LoopUnroll(gamma: Rationals.Rational, ak: (bool, nat), s: Rand.Bitstream)
-    requires 0 <= gamma.numer <= gamma.denom
-    requires ak.0
-    ensures Model.Le1Loop(gamma)(ak)(s) == Monad.Bind(Model.Le1LoopIter(gamma)(ak), Model.Le1Loop(gamma))(s)
-  {
-    Model.Le1LoopTerminatesAlmostSurely(gamma);
-    calc {
-      Model.Le1Loop(gamma)(ak)(s);
-      { reveal Model.Le1Loop(); }
-      Loops.While(Model.Le1LoopCondition, Model.Le1LoopIter(gamma))(ak)(s);
-      { reveal Model.Le1Loop();
-        Loops.WhileUnroll(Model.Le1LoopCondition, Model.Le1LoopIter(gamma), ak, s); }
-      Monad.Bind(Model.Le1LoopIter(gamma)(ak), Model.Le1Loop(gamma))(s);
-    }
-  }
-
-  lemma EnsureCaseLe1LoopInvariantOnEntry(gamma: Rationals.Rational, s: Rand.Bitstream)
-    requires 0 <= gamma.numer <= gamma.denom
-    ensures CaseLe1LoopInvariant(gamma, s, true, 0, s)
-  {
-    reveal CaseLe1LoopInvariant();
-  }
-
-  lemma EnsureCaseLe1LoopInvariantMaintained(gamma: Rationals.Rational, oldS: Rand.Bitstream, k: nat, s: Rand.Bitstream, a': bool, k': nat, s': Rand.Bitstream)
-    requires 0 <= gamma.numer <= gamma.denom
-    requires k' == k + 1
-    requires inv: CaseLe1LoopInvariant(gamma, oldS, true, k, s)
-    requires bernoulli: Monad.Result(a' , s') == Bernoulli.Model.Sample(gamma.numer, k' * gamma.denom)(s)
-    ensures CaseLe1LoopInvariant(gamma, oldS, a', k', s')
-  {
-    assert iter: Monad.Result((a', k'), s') == Model.Le1LoopIter(gamma)((true, k))(s) by {
-      reveal bernoulli;
-    }
-    calc {
-      Model.Le1Loop(gamma)((true, 0))(oldS);
-      { reveal CaseLe1LoopInvariant(); reveal inv; }
-      Model.Le1Loop(gamma)((true, k))(s);
-      { reveal iter; Le1LoopUnroll(gamma, (true, k), s); }
-      Model.Le1Loop(gamma)((a', k'))(s');
-    }
-    reveal CaseLe1LoopInvariant();
-  }
-
-  lemma EnsureCaseLe1PostCondition(gamma: Rationals.Rational, oldS: Rand.Bitstream, k: nat, s: Rand.Bitstream, c: bool)
-    requires 0 <= gamma.numer <= gamma.denom
-    requires CaseLe1LoopInvariant(gamma, oldS, false, k, s)
-    requires c <==> (k % 2 == 1)
-    ensures Monad.Result(c, s) == Model.SampleLe1(gamma)(oldS)
-  {
-    calc {
-      Model.Le1Loop(gamma)((true, 0))(oldS);
-      { reveal CaseLe1LoopInvariant(); }
-      Model.Le1Loop(gamma)((false, k))(s);
-      { reveal Model.Le1Loop(); Loops.WhileInitialConditionNegated(Model.Le1LoopCondition, Model.Le1LoopIter(gamma), (false, k), s); }
-      Monad.Result((false, k), s);
-    }
-  }
-}
diff --git a/src/Distributions/BernoulliExpNeg/Implementation.dfy b/src/Distributions/BernoulliExpNeg/Implementation.dfy
deleted file mode 100644
index ea90cf75..00000000
--- a/src/Distributions/BernoulliExpNeg/Implementation.dfy
+++ /dev/null
@@ -1,69 +0,0 @@
-/*******************************************************************************
- *  Copyright by the contributors to the Dafny Project
- *  SPDX-License-Identifier: MIT
- *******************************************************************************/
-
-module BernoulliExpNeg.Implementation {
-  import NatArith
-  import Rationals
-  import Rand
-  import Monad
-  import Interface
-  import Model
-  import Bernoulli
-  import Equivalence
-
-  trait {:termination false} Trait extends Interface.Trait {
-
-    // Based on Algorithm 1 in https://arxiv.org/pdf/2004.00010.pdf; unverified
-    method BernoulliExpNegSample(gamma: Rationals.Rational) returns (c: bool)
-      modifies `s
-      requires gamma.numer >= 0
-      decreases *
-      ensures Monad.Result(c, s) == Model.Sample(gamma)(old(s))
-    {
-      var gamma' := gamma;
-      while gamma'.numer > gamma'.denom
-        decreases gamma'.numer
-        invariant gamma'.numer >= 0
-        invariant Model.Sample(gamma)(old(s)) == Model.Sample(gamma')(s)
-      {
-        ghost var prevGamma := gamma';
-        ghost var prevS := s;
-        var b := BernoulliExpNegSampleCaseLe1(Rationals.Int(1));
-        gamma' := Rationals.Rational(gamma'.numer - gamma'.denom, gamma'.denom);
-        Equivalence.SampleUnfold(gamma', s, prevGamma, prevS, b);
-        if !b {
-          return false;
-        }
-      }
-      c:= BernoulliExpNegSampleCaseLe1(gamma');
-      reveal Model.Sample();
-    }
-
-    method BernoulliExpNegSampleCaseLe1(gamma: Rationals.Rational) returns (c: bool)
-      modifies `s
-      requires 0 <= gamma.numer <= gamma.denom
-      decreases *
-      ensures Monad.Result(c, s) == Model.SampleLe1(gamma)(old(s))
-    {
-      var k: nat := 0;
-      var a := true;
-      Equivalence.EnsureCaseLe1LoopInvariantOnEntry(gamma, s);
-      while a
-        decreases *
-        invariant Equivalence.CaseLe1LoopInvariant(gamma, old(s), a, k, s)
-      {
-        ghost var prevK: nat := k;
-        ghost var prevS := s;
-        k := k + 1;
-        NatArith.MulMonotonic(1, gamma.denom, k, gamma.denom);
-        a := BernoulliSample(Rationals.Rational(gamma.numer, k * gamma.denom));
-        Equivalence.EnsureCaseLe1LoopInvariantMaintained(gamma, old(s), prevK, prevS, a, k, s);
-      }
-      c := k % 2 == 1;
-      Equivalence.EnsureCaseLe1PostCondition(gamma, old(s), k, s, c);
-    }
-  }
-
-}
diff --git a/src/Distributions/BernoulliExpNeg/Interface.dfy b/src/Distributions/BernoulliExpNeg/Interface.dfy
deleted file mode 100644
index 9d5ab259..00000000
--- a/src/Distributions/BernoulliExpNeg/Interface.dfy
+++ /dev/null
@@ -1,21 +0,0 @@
-/*******************************************************************************
- *  Copyright by the contributors to the Dafny Project
- *  SPDX-License-Identifier: MIT
- *******************************************************************************/
-
-module BernoulliExpNeg.Interface {
-  import Rationals
-  import Monad
-  import Bernoulli
-  import Model
-
-  trait {:termination false} Trait extends Bernoulli.Interface.Trait {
-
-    method BernoulliExpNegSample(gamma: Rationals.Rational) returns (c: bool)
-      modifies `s
-      decreases *
-      requires gamma.numer >= 0
-      ensures Monad.Result(c, s) == Model.Sample(gamma)(old(s))
-
-  }
-}
diff --git a/src/Distributions/BernoulliExpNeg/Model.dfy b/src/Distributions/BernoulliExpNeg/Model.dfy
deleted file mode 100644
index 0235358e..00000000
--- a/src/Distributions/BernoulliExpNeg/Model.dfy
+++ /dev/null
@@ -1,259 +0,0 @@
-/*******************************************************************************
- *  Copyright by the contributors to the Dafny Project
- *  SPDX-License-Identifier: MIT
- *******************************************************************************/
-
-module BernoulliExpNeg.Model {
-  import Rationals
-  import RealArith
-  import Sequences
-  import Limits
-  import Measures
-  import Rand
-  import Uniform
-  import Bernoulli
-  import Monad
-  import Independence
-  import Loops
-  import BernoulliModel = Bernoulli.Model
-
-  opaque ghost function Sample(gamma: Rationals.Rational): Monad.Hurd<bool>
-    requires gamma.denom != 0
-    requires gamma.numer >= 0
-    decreases gamma.numer
-  {
-    if gamma.numer <= gamma.denom
-    then SampleLe1(gamma)
-    else Monad.Bind(
-        SampleLe1(Rationals.Int(1)),
-        b =>
-          var res: Monad.Hurd<bool> :=
-            if b
-            then Sample(Rationals.Rational(gamma.numer - gamma.denom, gamma.denom))
-            else Monad.Return(false);
-          res
-      )
-  }
-
-  ghost function SampleLe1(gamma: Rationals.Rational): Monad.Hurd<bool>
-    requires 0 <= gamma.numer <= gamma.denom
-  {
-    Monad.Map(
-      Le1Loop(gamma)((true, 0)),
-      (ak: (bool, nat)) => ak.1 % 2 == 1
-    )
-  }
-
-  opaque ghost function Le1Loop(gamma: Rationals.Rational): ((bool, nat)) -> Monad.Hurd<(bool, nat)>
-    requires 0 <= gamma.numer <= gamma.denom
-  {
-    Le1LoopTerminatesAlmostSurely(gamma);
-    Loops.While(
-      Le1LoopCondition,
-      Le1LoopIter(gamma)
-    )
-  }
-
-  ghost function Le1LoopCondition(ak: (bool, nat)): bool {
-    ak.0
-  }
-
-  ghost function Le1LoopIter(gamma: Rationals.Rational): ((bool, nat)) -> Monad.Hurd<(bool, nat)>
-    requires 0 <= gamma.numer <= gamma.denom
-  {
-    (ak: (bool, nat)) =>
-      var k': nat := ak.1 + 1;
-      Monad.Map(
-        BernoulliModel.Sample(gamma.numer, k' * gamma.denom),
-        a => (a, k'))
-  }
-
-  lemma {:axiom} Le1LoopIterCorrectness(gamma: Rationals.Rational, ak: (bool, nat), k': nat) // TODO
-    requires k' == ak.1 + 1
-    requires 0 <= gamma.numer <= gamma.denom
-    ensures Rand.prob(Monad.BitstreamsWithValueIn(Le1LoopIter(gamma)(ak), iset{(true, k')})) == gamma.numer as real / (k' * gamma.denom) as real
-    ensures Rand.prob(Monad.BitstreamsWithValueIn(Le1LoopIter(gamma)(ak), iset{(false, k')})) == 1.0 - gamma.numer as real / (k' * gamma.denom) as real
-
-  ghost function Le1LoopDivergenceProbBound(k: nat): nat -> real {
-    (fuel: nat) => if fuel == 0 then 1.0 else 1.0 / (k + fuel) as real
-  }
-
-  /* Proves that Le1Loop terminates almost surely.
-   *
-   * It does so by bounding the probability of nontermination (divergence) after `fuel` iterations between 0.0 and `Le1LoopDivergenceProbBound(k)(fuel)`.
-   * Since both lower and upper bounds converge to zero as `fuel` tends to infinity, the divergence probability is zero.
-   */
-  lemma Le1LoopTerminatesAlmostSurely(gamma: Rationals.Rational)
-    requires 0 <= gamma.numer <= gamma.denom
-    ensures Loops.WhileTerminatesAlmostSurely(Le1LoopCondition, Le1LoopIter(gamma))
-  {
-    var divergenceLower := Sequences.Constant(0.0);
-    forall init: (bool, nat) ensures Loops.WhileTerminatesAlmostSurelyInit(Le1LoopCondition, Le1LoopIter(gamma), init) {
-      var divergenceProb := Loops.WhileCutDivergenceProbability(Le1LoopCondition, Le1LoopIter(gamma), init);
-      assert Sequences.IsLeq(divergenceProb, Le1LoopDivergenceProbBound(0)) by {
-        forall fuel: nat ensures divergenceProb(fuel) <= Le1LoopDivergenceProbBound(0)(fuel) {
-          Le1LoopDivergenceProbabilityBound(gamma, init, fuel);
-          assert divergenceProb(fuel) <= Le1LoopDivergenceProbBound(init.1)(fuel) <= Le1LoopDivergenceProbBound(0)(fuel) by {
-            if fuel > 0 {
-              RealArith.DivisionIsAntiMonotonic(1.0, (fuel + init.1) as real, fuel as real);
-            }
-          }
-        }
-      }
-      assert Limits.ConvergesTo(Le1LoopDivergenceProbBound(0), 0.0) by {
-        assert Sequences.IsSuffixOf(Sequences.OneOverNPlus1, Le1LoopDivergenceProbBound(0), 1);
-        Limits.OneOverNPlus1ConvergesToZero();
-        Limits.SuffixConvergesToSame(Le1LoopDivergenceProbBound(0), Sequences.OneOverNPlus1, 1, 0.0);
-      }
-      assert Sequences.IsLeq(divergenceLower, divergenceProb) by {
-        forall n: nat ensures divergenceProb(n) >= 0.0 {
-          Rand.ProbIsProbabilityMeasure();
-          var event := iset s | !Loops.WhileCutTerminatesWithFuel(Le1LoopCondition, Le1LoopIter(gamma), init, s)(n);
-          assume {:axiom} event in Rand.eventSpace; // TODO
-          assert Rand.prob(event) >= 0.0;
-        }
-      }
-      assert Limits.ConvergesTo(divergenceLower, 0.0) by {
-        Limits.ConstantSequenceConverges(divergenceLower, 0.0);
-      }
-      assert Limits.ConvergesTo(divergenceProb, 0.0) by {
-        Limits.Sandwich(divergenceLower, divergenceProb, Le1LoopDivergenceProbBound(0), 0.0);
-      }
-      Loops.EnsureWhileTerminatesAlmostSurelyViaLimit(Le1LoopCondition, Le1LoopIter(gamma), init);
-    }
-  }
-
-  /* Proves a bound on the nontermination (divergence) probability of Le1Loop(gamma)(ak) after `n` iterations.
-   *
-   * We could prove an even stronger bound of (gamma / (k + 1)) * (gamma / (k + 2)) * ... * (gamma / (k + n)),
-   * where k = ak.1.
-   * But such a bound is more complicated to express, so we instead prove the upper bound 1 if n = 0 and 1 / (k + n) else.
-   *
-   * We do this inductively.
-   * For n = 0 the bound of 1 is clear.
-   * For n >= 1, we reason as follows:
-   * If ak.0 is false, the loop stops and the divergence probability is zero.
-   * So the interesting case is ak.0 == true.
-   * If n == 1, the probability of sampling true from the Bernoulli(gamma / (k + 1)) distribution is <= 1 / (k + 1) == 1 / (k + n).
-   * Since the loop exits if we sample false, this is an upper bound on the divergence probability.
-   * If n > 1, the probability of divergence is at most that of the remaining n - 1 loop iterations diverging.
-   * By induction hypothesis, this is <= 1 / (k' + n - 1) == 1 / (k + n) as desired.
-   */
-  lemma Le1LoopDivergenceProbabilityBound(gamma: Rationals.Rational, ak: (bool, nat), n: nat)
-    requires 0 <= gamma.numer <= gamma.denom
-    decreases n
-    ensures Loops.WhileCutDivergenceProbability(Le1LoopCondition, Le1LoopIter(gamma), ak)(n) <= Le1LoopDivergenceProbBound(ak.1)(n)
-  {
-    var divergenceEvent := iset s | !Loops.WhileCutTerminatesWithFuel(Le1LoopCondition, Le1LoopIter(gamma), ak, s)(n);
-    assume {:axiom} divergenceEvent in Rand.eventSpace;
-    if n == 0 {
-      assert Loops.WhileCutDivergenceProbability(Le1LoopCondition, Le1LoopIter(gamma), ak)(n) == Rand.prob(divergenceEvent);
-      Rand.ProbIsProbabilityMeasure();
-      Measures.ProbabilityLe1(divergenceEvent, Rand.eventSpace, Rand.prob);
-      assert Loops.WhileCutDivergenceProbability(Le1LoopCondition, Le1LoopIter(gamma), ak)(n) <= Le1LoopDivergenceProbBound(ak.1)(n);
-    } else {
-      if Le1LoopCondition(ak) {
-        var (a, k) := ak;
-        var k': nat := k + 1;
-        var continueValues := iset ak': (bool, nat) | Le1LoopCondition(ak');
-        var firstIterationContinuesEvent := Monad.BitstreamsWithValueIn(Le1LoopIter(gamma)(ak), continueValues);
-        var restOfLoopDiverges := iset s | !Loops.WhileCutTerminatesWithFuel(Le1LoopCondition, Le1LoopIter(gamma), (true, k'), s)(n - 1);
-        assume {:axiom} firstIterationContinuesEvent in Rand.eventSpace; // TODO
-        assume {:axiom} restOfLoopDiverges in Rand.eventSpace; // TODO
-        assert divergenceEvent == firstIterationContinuesEvent * Monad.BitstreamsWithRestIn(Le1LoopIter(gamma)(ak), restOfLoopDiverges) by {
-          forall s ensures s in divergenceEvent <==> s in firstIterationContinuesEvent && s in Monad.BitstreamsWithRestIn(Le1LoopIter(gamma)(ak), restOfLoopDiverges) {
-            match Le1LoopIter(gamma)(ak)(s)
-            case Diverging => {}
-            case Result(ak', s') =>
-              assert s in firstIterationContinuesEvent <==> ak'.0;
-              calc {
-                s in divergenceEvent;
-                !Loops.WhileCutTerminatesWithFuel(Le1LoopCondition, Le1LoopIter(gamma), ak, s)(n);
-                { Loops.WhileCutTerminatesWithFuelUnroll(Le1LoopCondition, Le1LoopIter(gamma), ak, s, ak', s', n); }
-                !Loops.WhileCutTerminatesWithFuel(Le1LoopCondition, Le1LoopIter(gamma), ak', s')(n - 1);
-              }
-              if ak'.0 {
-                assert ak' == (true, k');
-                assert s' in restOfLoopDiverges <==> !Loops.WhileCutTerminatesWithFuel(Le1LoopCondition, Le1LoopIter(gamma), ak', s')(n - 1);
-              } else {
-                assert s !in firstIterationContinuesEvent;
-                assert s !in divergenceEvent;
-              }
-          }
-        }
-        assert Rand.prob(divergenceEvent) == Rand.prob(firstIterationContinuesEvent) * Rand.prob(restOfLoopDiverges) by {
-          assume {:axiom} Independence.IsIndepFunction(Le1LoopIter(gamma)(ak)); // TODO
-          Independence.ResultsIndependent(Le1LoopIter(gamma)(ak), continueValues, restOfLoopDiverges);
-        }
-        calc {
-          Loops.WhileCutDivergenceProbability(Le1LoopCondition, Le1LoopIter(gamma), ak)(n);
-          Rand.prob(divergenceEvent);
-          Rand.prob(firstIterationContinuesEvent) * Rand.prob(restOfLoopDiverges);
-        <= { Le1LoopDivergenceProbabilityBoundHelper(gamma, ak, k', n, firstIterationContinuesEvent, restOfLoopDiverges); }
-          Le1LoopDivergenceProbBound(ak.1)(n);
-        }
-      } else { // i.e. if !LeLoopCondition(ak)
-        forall s ensures Loops.WhileCutTerminatesWithFuel(Le1LoopCondition, Le1LoopIter(gamma), ak, s)(n) {}
-        assert Loops.WhileCutDivergenceProbability(Le1LoopCondition, Le1LoopIter(gamma), ak)(n) == 0.0 by {
-          assert iset{} == iset s | !Loops.WhileCutTerminatesWithFuel(Le1LoopCondition, Le1LoopIter(gamma), ak, s)(n);
-          Rand.ProbIsProbabilityMeasure();
-        }
-      }
-      assert Loops.WhileCutDivergenceProbability(Le1LoopCondition, Le1LoopIter(gamma), ak)(n) <= Le1LoopDivergenceProbBound(ak.1)(n);
-    }
-  }
-
-  lemma Le1LoopDivergenceProbabilityBoundHelper(gamma: Rationals.Rational, ak: (bool, nat), k': nat, n: nat, firstIterationContinuesEvent: iset<Rand.Bitstream>, restOfLoopDiverges: iset<Rand.Bitstream>)
-    requires 0 <= gamma.numer <= gamma.denom
-    requires n >= 1
-    requires k' == ak.1 + 1
-    requires firstIterationContinuesEvent == iset s | Le1LoopIter(gamma)(ak)(s).Satisfies(Le1LoopCondition)
-    requires restOfLoopDiverges == iset s | !Loops.WhileCutTerminatesWithFuel(Le1LoopCondition, Le1LoopIter(gamma), (true, k'), s)(n - 1)
-    requires firstIterationContinuesEvent in Rand.eventSpace
-    requires restOfLoopDiverges in Rand.eventSpace
-    decreases n, 1
-    ensures Rand.prob(firstIterationContinuesEvent) * Rand.prob(restOfLoopDiverges) <= Le1LoopDivergenceProbBound(ak.1)(n)
-  {
-    var k := ak.1;
-    var k' := k + 1;
-    assert 0.0 <= Rand.prob(restOfLoopDiverges) <= 1.0 by {
-      Rand.ProbIsProbabilityMeasure();
-      Measures.ProbabilityLe1(restOfLoopDiverges, Rand.eventSpace, Rand.prob);
-    }
-    assert 0.0 <= Rand.prob(firstIterationContinuesEvent) <= 1.0 by {
-      Rand.ProbIsProbabilityMeasure(); Measures.ProbabilityLe1(firstIterationContinuesEvent, Rand.eventSpace, Rand.prob);
-    }
-    if n == 1 {
-      assert Rand.prob(firstIterationContinuesEvent) <= 1.0 / k' as real by {
-        assert firstIterationContinuesEvent == Monad.BitstreamsWithValueIn(Le1LoopIter(gamma)(ak), iset{(true, k')}) by {
-          forall s ensures s in firstIterationContinuesEvent <==> s in Monad.BitstreamsWithValueIn(Le1LoopIter(gamma)(ak), iset{(true, k')}) {}
-        }
-        calc {
-          Rand.prob(firstIterationContinuesEvent);
-          { Le1LoopIterCorrectness(gamma, ak, k'); }
-          gamma.numer as real / (k' * gamma.denom) as real;
-        <=
-          1.0 / k' as real;
-        }
-      }
-      calc {
-        Rand.prob(firstIterationContinuesEvent) * Rand.prob(restOfLoopDiverges);
-      <= { RealArith.MulMonotonic(Rand.prob(firstIterationContinuesEvent), Rand.prob(restOfLoopDiverges), 1.0); }
-        Rand.prob(firstIterationContinuesEvent);
-      <=
-        1.0 / k' as real;
-        1.0 / (k + n) as real;
-      }
-    } else {
-      calc {
-        Rand.prob(firstIterationContinuesEvent) * Rand.prob(restOfLoopDiverges);
-      <= { RealArith.MulMonotonic(Rand.prob(restOfLoopDiverges), Rand.prob(firstIterationContinuesEvent), 1.0); }
-        Rand.prob(restOfLoopDiverges);
-        Loops.WhileCutDivergenceProbability(Le1LoopCondition, Le1LoopIter(gamma), (true, k'))(n - 1);
-      <= { Le1LoopDivergenceProbabilityBound(gamma, (true, k'), n - 1); }
-        1.0 / (k + n) as real;
-      }
-    }
-  }
-
-}
diff --git a/src/Distributions/DiscreteGaussian/Implementation.dfy b/src/Distributions/DiscreteGaussian/Implementation.dfy
deleted file mode 100644
index 503df92b..00000000
--- a/src/Distributions/DiscreteGaussian/Implementation.dfy
+++ /dev/null
@@ -1,37 +0,0 @@
-/*******************************************************************************
- *  Copyright by the contributors to the Dafny Project
- *  SPDX-License-Identifier: MIT
- *******************************************************************************/
-
-module DiscreteGaussian.Implementation {
-  import Rationals
-  import Interface
-
-  trait {:termination false} Trait extends Interface.Trait {
-
-    // Based on Algorithm 3 in https://arxiv.org/pdf/2004.00010.pdf; unverified
-    // Note that we take sigma as a parameter, not sigma^2, to avoid square roots.
-    method DiscreteGaussianSample(sigma: Rationals.Rational) returns (y: int)
-      modifies `s
-      requires sigma.numer >= 1
-      decreases *
-    {
-      var sigmaSquared := sigma.Mul(sigma);
-      var t := sigma.Floor() + 1;
-      while true
-        decreases *
-      {
-        y := DiscreteLaplaceSample(Rationals.Int(t));
-        var yAbs := if y < 0 then -y else y;
-        var numeratorTerm := Rationals.Int(yAbs).Sub(sigmaSquared.Div(Rationals.Int(t)));
-        var gamma := numeratorTerm.Mul(numeratorTerm)
-        .Div(Rationals.Int(2).Mul(sigmaSquared));
-        var c := BernoulliExpNegSample(gamma);
-        if c {
-          return;
-        }
-      }
-    }
-
-  }
-}
diff --git a/src/Distributions/DiscreteGaussian/Interface.dfy b/src/Distributions/DiscreteGaussian/Interface.dfy
deleted file mode 100644
index 1f937adf..00000000
--- a/src/Distributions/DiscreteGaussian/Interface.dfy
+++ /dev/null
@@ -1,18 +0,0 @@
-/*******************************************************************************
- *  Copyright by the contributors to the Dafny Project
- *  SPDX-License-Identifier: MIT
- *******************************************************************************/
-
-module DiscreteGaussian.Interface {
-  import Rationals
-  import BernoulliExpNeg
-  import DiscreteLaplace
-
-  trait {:termination false} Trait extends DiscreteLaplace.Interface.Trait, BernoulliExpNeg.Interface.Trait {
-    // Takes the sigma (not sigma^2!) as a fraction
-    method DiscreteGaussianSample(sigma: Rationals.Rational) returns (y: int)
-      modifies `s
-      requires sigma.numer >= 1
-      decreases *
-  }
-}
diff --git a/src/Distributions/DiscreteLaplace/Correctness.dfy b/src/Distributions/DiscreteLaplace/Correctness.dfy
deleted file mode 100644
index a36a92d2..00000000
--- a/src/Distributions/DiscreteLaplace/Correctness.dfy
+++ /dev/null
@@ -1,41 +0,0 @@
-/*******************************************************************************
- *  Copyright by the contributors to the Dafny Project
- *  SPDX-License-Identifier: MIT
- *******************************************************************************/
-
-module DiscreteLaplace.Correctness {
-
-  import Rationals
-  import Exponential
-  import Rand
-  import Model
-  import Monad
-  import RealArith
-
-  ghost function Numerator1(scale: Rationals.Rational): real
-    requires scale.numer >= 1
-  {
-    Exponential.Exp((1.0 / scale.ToReal()) - 1.0)
-  }
-
-  ghost function Numerator2(scale: Rationals.Rational, x: int): real
-    requires scale.numer >= 1
-  {
-    Exponential.Exp(-((RealArith.Abs(x as real)) / scale.ToReal()))
-  }
-
-  ghost function Denominator(scale: Rationals.Rational): (r: real)
-    requires scale.numer >= 1
-    ensures r > 0.0
-  {
-    var x := (1.0 / scale.ToReal()) + 1.0;
-    Exponential.Positive(x);
-    Exponential.Exp(x)
-  }
-
-  // Based on Definition 28 in https://arxiv.org/pdf/2004.00010.pdf
-  lemma {:axiom} Correctness(scale: Rationals.Rational, x: int)
-    requires scale.numer >= 1
-    ensures Rand.prob(iset s | Model.Sample(scale)(s).Equals(x)) == (Numerator1(scale) * Numerator2(scale, x)) / Denominator(scale)
-
-}
diff --git a/src/Distributions/DiscreteLaplace/Equivalence.dfy b/src/Distributions/DiscreteLaplace/Equivalence.dfy
deleted file mode 100644
index fddfb93e..00000000
--- a/src/Distributions/DiscreteLaplace/Equivalence.dfy
+++ /dev/null
@@ -1,397 +0,0 @@
-/*******************************************************************************
- *  Copyright by the contributors to the Dafny Project
- *  SPDX-License-Identifier: MIT
- *******************************************************************************/
-
-module DiscreteLaplace.Equivalence {
-  import Rand
-  import Monad
-  import Rationals
-  import BernoulliExpNeg
-  import Uniform
-  import Coin
-  import Model
-  import Loops
-
-  /*******
-   Lemmas
-  *******/
-
-  lemma SampleLiftToEnsures(scale: Rationals.Rational, s: Rand.Bitstream, t: Rand.Bitstream, x: (bool, int))
-    requires R1: scale.numer >= 1
-    requires R2: Monad.Result(x, s) == Model.SampleLoop(scale)(t)
-    ensures Model.Sample(scale)(t) == Monad.Result(if x.0 then -x.1 else x.1, s)
-  {
-    var f := (x: (bool, int)) => if x.0 then -x.1 else x.1;
-
-    calc {
-      Monad.Result(if x.0 then -x.1 else x.1, s);
-      Monad.Result(f(x), s);
-      Monad.Result(x, s).Map(f);
-      { reveal R1; reveal R2; }
-      Model.SampleLoop(scale)(t).Map(f);
-      Monad.Map(Model.SampleLoop(scale), f)(t);
-      { reveal Model.Sample(); }
-      Model.Sample(scale)(t);
-    }
-  }
-
-  lemma SampleLoopLiftToEnsures(scale: Rationals.Rational, s: Rand.Bitstream, t: Rand.Bitstream, bY: (bool, int))
-    requires R1: scale.numer >= 1
-    requires R2: Model.SampleLoop(scale)(s) == Model.SampleLoop(scale, bY)(t)
-    requires R3: !(bY.0 && (bY.1 == 0))
-    ensures Model.SampleLoop(scale)(s) == Monad.Result(bY, t)
-  {
-    calc {
-      Monad.Result(bY, t);
-      Monad.Return(bY)(t);
-      { reveal R1; SampleLoopTailRecursiveEquivalence(scale, t, bY); reveal R3; }
-      Model.SampleLoop(scale, bY)(t);
-      { reveal R1; reveal R2; }
-      Model.SampleLoop(scale)(s);
-    }
-  }
-
-  lemma SampleInnerLoopLiftToEnsures(s: Rand.Bitstream, t: Rand.Bitstream, a: bool, v: int)
-    requires R1: Model.SampleInnerLoop()(s) == Model.SampleInnerLoop((a, v))(t)
-    requires R2: !a
-    ensures Model.SampleInnerLoopFull()(s) == Monad.Result(v, t)
-  {
-    var f := (x: (bool, int)) => x.1;
-
-    assert A: Model.SampleInnerLoop()(s) == Monad.Result((a, v), t) by {
-      calc {
-        Model.SampleInnerLoop()(s);
-        { reveal R1; }
-        Model.SampleInnerLoop((a, v))(t);
-        { SampleInnerLoopTailRecursiveEquivalence(t, (a, v)); reveal R2; }
-        Monad.Return((a, v))(t);
-        Monad.Result((a, v), t);
-      }
-    }
-
-    calc {
-      Model.SampleInnerLoopFull()(s);
-      { reveal Model.SampleInnerLoopFull(); }
-      Model.SampleInnerLoop()(s).Map(f);
-      { reveal A; }
-      Monad.Result((a, v), t).Map(f);
-      Monad.Result(v, t);
-    }
-  }
-
-  lemma {:rlimit 100000} SampleLoopTailRecursiveEquivalence(scale: Rationals.Rational, s: Rand.Bitstream, bY: (bool, int) := (true, 0))
-    requires scale.numer >= 1
-    ensures
-      Model.SampleLoop(scale, bY)(s) ==
-      if bY.0 && (bY.1 == 0) then
-        Monad.Bind(
-          Uniform.Model.Sample(scale.numer),
-          (u: nat) =>
-            Monad.Bind(
-              BernoulliExpNeg.Model.Sample(Rationals.Rational(u, scale.numer)),
-              (d: bool) =>
-                if d then
-                  Monad.Bind(
-                    Model.SampleInnerLoopFull(),
-                    (v: int) =>
-                      Monad.Bind(
-                        Coin.Model.Sample,
-                        (b: bool) =>
-                          var x := u + scale.numer * v;
-                          var y := x / scale.denom;
-                          Model.SampleLoop(scale, (b, y))
-                      )
-                  )
-                else
-                  Model.SampleLoop(scale, bY)
-            )
-        )(s)
-      else
-        Monad.Return(bY)(s)
-  {
-    calc {
-      Model.SampleLoop(scale, bY)(s);
-    == { reveal Model.SampleLoop();
-         reveal Loops.While();
-         reveal Model.SampleLoopCondition();
-         reveal Model.SampleLoopBody(); }
-      Loops.While(Model.SampleLoopCondition, Model.SampleLoopBody(scale))(bY)(s);
-    == { Model.SampleLoopTerminatesAlmostSurely(scale);
-         Loops.WhileUnroll(Model.SampleLoopCondition, Model.SampleLoopBody(scale), bY, s);
-         reveal Model.SampleLoop();
-         reveal Model.SampleLoopCondition();
-         reveal Model.SampleLoopBody(); }
-      if Model.SampleLoopCondition(bY) then
-        Monad.Bind(Model.SampleLoopBody(scale)(bY), (bY: (bool, int)) => Model.SampleLoop(scale, bY))(s)
-      else
-        Monad.Return(bY)(s);
-    == { reveal Model.SampleLoopCondition();
-         reveal Model.SampleLoopBody(); }
-      if bY.0 && (bY.1 == 0) then
-        Monad.Bind(
-          Monad.Bind(
-            Uniform.Model.Sample(scale.numer),
-            (u: nat) =>
-              Monad.Bind(
-                BernoulliExpNeg.Model.Sample(Rationals.Rational(u, scale.numer)),
-                (d: bool) =>
-                  if d then
-                    Monad.Bind(
-                      Model.SampleInnerLoopFull(),
-                      (v: int) =>
-                        Monad.Bind(
-                          Coin.Model.Sample,
-                          (b: bool) =>
-                            var x := u + scale.numer * v;
-                            var y := x / scale.denom;
-                            Monad.Return((b, y))
-                        )
-                    )
-                  else
-                    Monad.Return(bY)
-              )
-          ),
-          (bY: (bool, int)) => Model.SampleLoop(scale, bY)
-        )(s)
-      else
-        Monad.Return(bY)(s);
-    ==
-      if bY.0 && (bY.1 == 0) then
-        Monad.Bind(
-          Uniform.Model.Sample(scale.numer),
-          (u: nat) =>
-            Monad.Bind(
-              Monad.Bind(
-                BernoulliExpNeg.Model.Sample(Rationals.Rational(u, scale.numer)),
-                (d: bool) =>
-                  if d then
-                    Monad.Bind(
-                      Model.SampleInnerLoopFull(),
-                      (v: int) =>
-                        Monad.Bind(
-                          Coin.Model.Sample,
-                          (b: bool) =>
-                            var x := u + scale.numer * v;
-                            var y := x / scale.denom;
-                            Monad.Return((b, y))
-                        )
-                    )
-                  else
-                    Monad.Return(bY)
-              ),
-              (bY: (bool, int)) => Model.SampleLoop(scale, bY)
-            )
-        )(s)
-      else
-        Monad.Return(bY)(s);
-    ==
-      if bY.0 && (bY.1 == 0) then
-        Monad.Bind(
-          Uniform.Model.Sample(scale.numer),
-          (u: nat) =>
-            Monad.Bind(
-              BernoulliExpNeg.Model.Sample(Rationals.Rational(u, scale.numer)),
-              (d: bool) =>
-                Monad.Bind(
-                  if d then
-                    Monad.Bind(
-                      Model.SampleInnerLoopFull(),
-                      (v: int) =>
-                        Monad.Bind(
-                          Coin.Model.Sample,
-                          (b: bool) =>
-                            var x := u + scale.numer * v;
-                            var y := x / scale.denom;
-                            Monad.Return((b, y))
-                        )
-                    )
-                  else
-                    Monad.Return(bY),
-                  (bY: (bool, int)) => Model.SampleLoop(scale, bY)
-                )
-            )
-        )(s)
-      else
-        Monad.Return(bY)(s);
-    ==
-      if bY.0 && (bY.1 == 0) then
-        Monad.Bind(
-          Uniform.Model.Sample(scale.numer),
-          (u: nat) =>
-            Monad.Bind(
-              BernoulliExpNeg.Model.Sample(Rationals.Rational(u, scale.numer)),
-              (d: bool) =>
-                if d then
-                  Monad.Bind(
-                    Monad.Bind(
-                      Model.SampleInnerLoopFull(),
-                      (v: int) =>
-                        Monad.Bind(
-                          Coin.Model.Sample,
-                          (b: bool) =>
-                            var x := u + scale.numer * v;
-                            var y := x / scale.denom;
-                            Monad.Return((b, y))
-                        )
-                    ),
-                    (bY: (bool, int)) => Model.SampleLoop(scale, bY)
-                  )
-                else
-                  Monad.Bind(
-                    Monad.Return(bY),
-                    (bY: (bool, int)) => Model.SampleLoop(scale, bY)
-                  )
-            )
-        )(s)
-      else
-        Monad.Return(bY)(s);
-    ==
-      if bY.0 && (bY.1 == 0) then
-        Monad.Bind(
-          Uniform.Model.Sample(scale.numer),
-          (u: nat) =>
-            Monad.Bind(
-              BernoulliExpNeg.Model.Sample(Rationals.Rational(u, scale.numer)),
-              (d: bool) =>
-                if d then
-                  Monad.Bind(
-                    Model.SampleInnerLoopFull(),
-                    (v: int) =>
-                      Monad.Bind(
-                        Monad.Bind(
-                          Coin.Model.Sample,
-                          (b: bool) =>
-                            var x := u + scale.numer * v;
-                            var y := x / scale.denom;
-                            Monad.Return((b, y))
-                        ),
-                        (bY: (bool, int)) => Model.SampleLoop(scale, bY)
-                      )
-                  )
-                else
-                  Model.SampleLoop(scale, bY)
-            )
-        )(s)
-      else
-        Monad.Return(bY)(s);
-    ==
-      if bY.0 && (bY.1 == 0) then
-        Monad.Bind(
-          Uniform.Model.Sample(scale.numer),
-          (u: nat) =>
-            Monad.Bind(
-              BernoulliExpNeg.Model.Sample(Rationals.Rational(u, scale.numer)),
-              (d: bool) =>
-                if d then
-                  Monad.Bind(
-                    Model.SampleInnerLoopFull(),
-                    (v: int) =>
-                      Monad.Bind(
-                        Coin.Model.Sample,
-                        (b: bool) =>
-                          var x := u + scale.numer * v;
-                          var y := x / scale.denom;
-                          Monad.Bind(
-                            Monad.Return((b, y)),
-                            (bY: (bool, int)) => Model.SampleLoop(scale, bY)
-                          )
-                      )
-                  )
-                else
-                  Model.SampleLoop(scale, bY)
-            )
-        )(s)
-      else
-        Monad.Return(bY)(s);
-    ==
-      if bY.0 && (bY.1 == 0) then
-        Monad.Bind(
-          Uniform.Model.Sample(scale.numer),
-          (u: nat) =>
-            Monad.Bind(
-              BernoulliExpNeg.Model.Sample(Rationals.Rational(u, scale.numer)),
-              (d: bool) =>
-                if d then
-                  Monad.Bind(
-                    Model.SampleInnerLoopFull(),
-                    (v: int) =>
-                      Monad.Bind(
-                        Coin.Model.Sample,
-                        (b: bool) =>
-                          var x := u + scale.numer * v;
-                          var y := x / scale.denom;
-                          Model.SampleLoop(scale, (b, y))
-                      )
-                  )
-                else
-                  Model.SampleLoop(scale, bY)
-            )
-        )(s)
-      else
-        Monad.Return(bY)(s);
-    }
-  }
-
-  lemma SampleInnerLoopTailRecursiveEquivalence(s: Rand.Bitstream, x: (bool, int) := (true, 0))
-    ensures
-      Model.SampleInnerLoop(x)(s) ==
-      if x.0 then
-        Monad.Bind(
-          BernoulliExpNeg.Model.Sample(Rationals.Int(1)),
-          (a: bool) => Model.SampleInnerLoop((a, if a then x.1 + 1 else x.1))
-        )(s)
-      else
-        Monad.Return(x)(s)
-  {
-    calc {
-      Model.SampleInnerLoop(x)(s);
-    == { reveal Loops.While();
-         reveal Model.SampleInnerLoop();
-         reveal Model.SampleInnerLoopBody();
-         reveal Model.SampleInnerLoopCondition(); }
-      Loops.While(Model.SampleInnerLoopCondition, Model.SampleInnerLoopBody)(x)(s);
-    == { Model.SampleInnerLoopTerminatesAlmostSurely();
-         Loops.WhileUnroll(Model.SampleInnerLoopCondition, Model.SampleInnerLoopBody, x, s);
-         reveal Model.SampleInnerLoop();
-         reveal Model.SampleInnerLoopBody();
-         reveal Model.SampleInnerLoopCondition(); }
-      if Model.SampleInnerLoopCondition(x) then
-        Monad.Bind(Model.SampleInnerLoopBody(x), (y: (bool, int)) => Model.SampleInnerLoop(y))(s)
-      else
-        Monad.Return(x)(s);
-    == { reveal Model.SampleInnerLoopBody();
-         reveal Model.SampleInnerLoopCondition(); }
-      if x.0 then
-        Monad.Bind(
-          Monad.Bind(
-            BernoulliExpNeg.Model.Sample(Rationals.Int(1)),
-            (a: bool) => Monad.Return((a, if a then x.1 + 1 else x.1))
-          ),
-          (y: (bool, int)) => Model.SampleInnerLoop(y)
-        )(s)
-      else
-        Monad.Return(x)(s);
-    ==
-      if x.0 then
-        Monad.Bind(
-          BernoulliExpNeg.Model.Sample(Rationals.Int(1)),
-          (a: bool) =>
-            Monad.Bind(
-              Monad.Return((a, if a then x.1 + 1 else x.1)),
-              (y: (bool, int)) => Model.SampleInnerLoop(y)
-            )
-        )(s)
-      else
-        Monad.Return(x)(s);
-    ==
-      if x.0 then
-        Monad.Bind(
-          BernoulliExpNeg.Model.Sample(Rationals.Int(1)),
-          (a: bool) => Model.SampleInnerLoop((a, if a then x.1 + 1 else x.1))
-        )(s)
-      else
-        Monad.Return(x)(s);
-    }
-  }
-}
\ No newline at end of file
diff --git a/src/Distributions/DiscreteLaplace/Implementation.dfy b/src/Distributions/DiscreteLaplace/Implementation.dfy
deleted file mode 100644
index a44cad09..00000000
--- a/src/Distributions/DiscreteLaplace/Implementation.dfy
+++ /dev/null
@@ -1,86 +0,0 @@
-/*******************************************************************************
- *  Copyright by the contributors to the Dafny Project
- *  SPDX-License-Identifier: MIT
- *******************************************************************************/
-
-module DiscreteLaplace.Implementation {
-  import Rationals
-  import Interface
-  import Monad
-  import Model
-  import Uniform
-  import BernoulliExpNeg
-  import Bernoulli
-  import Coin
-  import Equivalence
-  import Loops
-
-  trait {:termination false} Trait extends Interface.Trait {
-
-    // Based on Algorithm 2 in https://arxiv.org/pdf/2004.00010.pdf; unverified
-    method DiscreteLaplaceSample(scale: Rationals.Rational) returns (z: int)
-      modifies `s
-      requires scale.numer >= 1
-      decreases *
-      ensures Model.Sample(scale)(old(s)) == Monad.Result(z, s)
-    {
-      var x := DiscreteLaplaceSampleLoop(scale);
-      Equivalence.SampleLiftToEnsures(scale, s, old(s), x);
-      z := if x.0 then -x.1 else x.1;
-    }
-
-    method {:rlimit 100000} DiscreteLaplaceSampleLoop(scale: Rationals.Rational) returns (bY: (bool, int))
-      modifies `s
-      requires scale.numer >= 1
-      decreases *
-      ensures Model.SampleLoop(scale)(old(s)) == Monad.Result(bY, s)
-    {
-      var b := true;
-      var y := 0;
-
-
-      while b && y == 0
-        decreases *
-        invariant Model.SampleLoop(scale)(old(s)) == Model.SampleLoop(scale, (b, y))(s)
-      {
-        Equivalence.SampleLoopTailRecursiveEquivalence(scale, s, (b, y));
-
-        var u := UniformSample(scale.numer);
-        var d := BernoulliExpNegSample(Rationals.Rational(u, scale.numer));
-        if d {
-          var v := DisceteLaplaceSampleInnerLoop();
-          b := CoinSample();
-          var x := u + scale.numer * v;
-          y := x / scale.denom;
-        }
-      }
-
-      Equivalence.SampleLoopLiftToEnsures(scale, old(s), s, (b, y));
-
-      return (b, y);
-    }
-
-    method DisceteLaplaceSampleInnerLoop() returns (v: int)
-      modifies `s
-      decreases *
-      ensures Model.SampleInnerLoopFull()(old(s)) == Monad.Result(v, s)
-    {
-      var a := true;
-      v := 0;
-
-      while a
-        decreases *
-        invariant Model.SampleInnerLoop()(old(s)) == Model.SampleInnerLoop((a, v))(s)
-      {
-        Equivalence.SampleInnerLoopTailRecursiveEquivalence(s, (a, v));
-
-        a := BernoulliExpNegSample(Rationals.Int(1));
-        if a {
-          v := v + 1;
-        }
-      }
-
-      Equivalence.SampleInnerLoopLiftToEnsures(old(s), s, a, v);
-    }
-  }
-}
diff --git a/src/Distributions/DiscreteLaplace/Interface.dfy b/src/Distributions/DiscreteLaplace/Interface.dfy
deleted file mode 100644
index 60fdddfa..00000000
--- a/src/Distributions/DiscreteLaplace/Interface.dfy
+++ /dev/null
@@ -1,21 +0,0 @@
-/*******************************************************************************
- *  Copyright by the contributors to the Dafny Project
- *  SPDX-License-Identifier: MIT
- *******************************************************************************/
-
-module DiscreteLaplace.Interface {
-  import Rationals
-  import Coin
-  import Uniform
-  import BernoulliExpNeg
-
-  trait {:termination false} Trait extends Coin.Interface.Trait, Uniform.Interface.Trait, BernoulliExpNeg.Interface.Trait {
-
-    // Based on Algorithm 2 in https://arxiv.org/pdf/2004.00010.pdf; unverified
-    method DiscreteLaplaceSample(scale: Rationals.Rational) returns (z: int)
-      modifies `s
-      requires scale.numer >= 1
-      decreases *
-
-  }
-}
diff --git a/src/Distributions/DiscreteLaplace/Model.dfy b/src/Distributions/DiscreteLaplace/Model.dfy
deleted file mode 100644
index b69daf37..00000000
--- a/src/Distributions/DiscreteLaplace/Model.dfy
+++ /dev/null
@@ -1,96 +0,0 @@
-/*******************************************************************************
- *  Copyright by the contributors to the Dafny Project
- *  SPDX-License-Identifier: MIT
- *******************************************************************************/
-
-module DiscreteLaplace.Model {
-  import Monad
-  import Rand
-  import Rationals
-  import Uniform
-  import BernoulliExpNeg
-  import Coin
-  import Loops
-
-  /************
-   Definitions
-  ************/
-
-  ghost opaque function Sample(scale: Rationals.Rational): Monad.Hurd<int>
-    requires scale.numer >= 1
-  {
-    var f := (x: (bool, int)) => if x.0 then -x.1 else x.1;
-    Monad.Map(SampleLoop(scale), f)
-  }
-
-  ghost opaque function SampleLoop(scale: Rationals.Rational, x: (bool, int) := (true, 0)): Monad.Hurd<(bool, int)>
-    requires scale.numer >= 1
-  {
-    Loops.While(SampleLoopCondition, SampleLoopBody(scale))(x)
-  }
-
-  ghost opaque function SampleLoopBody(scale: Rationals.Rational): ((bool, int)) -> Monad.Hurd<(bool, int)>
-    requires scale.numer >= 1
-  {
-    (bY: (bool, int)) =>
-      Monad.Bind(
-        Uniform.Model.Sample(scale.numer),
-        (u: nat) =>
-          Monad.Bind(
-            BernoulliExpNeg.Model.Sample(Rationals.Rational(u, scale.numer)),
-            (d: bool) =>
-              if d then
-                Monad.Bind(
-                  SampleInnerLoopFull(),
-                  (v: int) =>
-                    Monad.Bind(
-                      Coin.Model.Sample,
-                      (b: bool) =>
-                        var x := u + scale.numer * v;
-                        var y := x / scale.denom;
-                        Monad.Return((b, y))
-                    )
-                )
-              else
-                Monad.Return(bY)
-          )
-      )
-  }
-
-  ghost opaque function SampleLoopCondition(x: (bool, int)): bool {
-    x.0 && (x.1 == 0)
-  }
-
-  ghost opaque function SampleInnerLoopFull(): Monad.Hurd<int> {
-    var f := (x: (bool, int)) => x.1;
-    Monad.Map(SampleInnerLoop(), f)
-  }
-
-  ghost opaque function SampleInnerLoop(x: (bool, int) := (true, 0)): Monad.Hurd<(bool, int)> {
-    Loops.While(SampleInnerLoopCondition, SampleInnerLoopBody)(x)
-  }
-
-  ghost opaque function SampleInnerLoopBody(x: (bool, int)): Monad.Hurd<(bool, int)> {
-    Monad.Bind(
-      BernoulliExpNeg.Model.Sample(Rationals.Int(1)),
-      (a: bool) => Monad.Return((a, if a then x.1 + 1 else x.1))
-    )
-  }
-
-  ghost opaque function SampleInnerLoopCondition(x: (bool, int)): bool {
-    x.0
-  }
-
-  /*******
-   Lemmas
-  *******/
-
-  // TODO: add later
-  lemma {:axiom} SampleInnerLoopTerminatesAlmostSurely()
-    ensures Loops.WhileTerminatesAlmostSurely(SampleInnerLoopCondition, SampleInnerLoopBody)
-
-  // TODO: add later
-  lemma {:axiom} SampleLoopTerminatesAlmostSurely(scale: Rationals.Rational)
-    requires scale.numer >= 1
-    ensures Loops.WhileTerminatesAlmostSurely(SampleLoopCondition, SampleLoopBody(scale))
-}
diff --git a/src/Distributions/Uniform/Equivalence.dfy b/src/Distributions/Uniform/Equivalence.dfy
deleted file mode 100644
index 7b7d36dc..00000000
--- a/src/Distributions/Uniform/Equivalence.dfy
+++ /dev/null
@@ -1,32 +0,0 @@
-/*******************************************************************************
- *  Copyright by the contributors to the Dafny Project
- *  SPDX-License-Identifier: MIT
- *******************************************************************************/
-
-module Uniform.Equivalence {
-  import Rand
-  import Model
-  import Monad
-  import Loops
-
-  lemma SampleUnroll(n: nat, s: Rand.Bitstream)
-    requires n > 0
-    ensures Model.Sample(n)(s) == Monad.Bind(Model.Proposal(n), (x: nat) => if Model.Accept(n)(x) then Monad.Return(x) else Model.Sample(n))(s)
-  {
-    Model.SampleTerminates(n);
-    reveal Model.Sample();
-    Loops.UntilUnroll(Model.Proposal(n), Model.Accept(n), s);
-  }
-
-  lemma SampleLifts(n: nat, u: nat, oldS: Rand.Bitstream, prevS: Rand.Bitstream, s: Rand.Bitstream)
-    requires n > 0
-    requires Model.Sample(n)(oldS) == Model.Sample(n)(prevS)
-    requires Monad.Result(u, s) == Model.Proposal(n)(prevS)
-    requires Model.Accept(n)(u)
-    ensures Model.Sample(n)(oldS) == Monad.Result(u, s)
-  {
-    reveal Model.Sample();
-    Model.SampleTerminates(n);
-    Loops.UntilUnroll(Model.Proposal(n), Model.Accept(n), prevS);
-  }
-}
diff --git a/src/Distributions/Uniform/Implementation.dfy b/src/Distributions/Uniform/Implementation.dfy
deleted file mode 100644
index a73ed78c..00000000
--- a/src/Distributions/Uniform/Implementation.dfy
+++ /dev/null
@@ -1,36 +0,0 @@
-/*******************************************************************************
- *  Copyright by the contributors to the Dafny Project
- *  SPDX-License-Identifier: MIT
- *******************************************************************************/
-
-module Uniform.Implementation {
-  import Monad
-  import UniformPowerOfTwo
-  import Model
-  import Interface
-  import Equivalence
-  import Loops
-
-  trait {:termination false} Trait extends Interface.Trait {
-    method UniformSample(n: nat) returns (u: nat)
-      modifies `s
-      decreases *
-      requires n > 0
-      ensures u < n
-      ensures Model.Sample(n)(old(s)) == Monad.Result(u, s)
-    {
-      ghost var prevS := s;
-      u := UniformPowerOfTwoSample(2 * n);
-      while u >= n
-        decreases *
-        invariant Model.Sample(n)(old(s)) == Model.Sample(n)(prevS)
-        invariant Monad.Result(u, s) == Model.Proposal(n)(prevS)
-      {
-        Equivalence.SampleUnroll(n, prevS);
-        prevS := s;
-        u := UniformPowerOfTwoSample(2 * n);
-      }
-      Equivalence.SampleLifts(n, u, old(s), prevS, s);
-    }
-  }
-}
diff --git a/src/Distributions/Uniform/Interface.dfy b/src/Distributions/Uniform/Interface.dfy
index 617175cd..af2306a0 100644
--- a/src/Distributions/Uniform/Interface.dfy
+++ b/src/Distributions/Uniform/Interface.dfy
@@ -8,10 +8,11 @@ module Uniform.Interface {
   import Coin
   import Model
   import UniformPowerOfTwo
+  import Pos
 
   trait {:termination false} Trait extends UniformPowerOfTwo.Interface.Trait {
 
-    method UniformSample(n: nat) returns (u: nat)
+    method UniformSample(n: Pos.pos) returns (u: nat)
       modifies `s
       decreases *
       requires n > 0
diff --git a/src/Math/Exponential.dfy b/src/Math/Exponential.dfy
deleted file mode 100644
index 248a9d4c..00000000
--- a/src/Math/Exponential.dfy
+++ /dev/null
@@ -1,106 +0,0 @@
-module Exponential {
-  import NatArith
-  import RealArith
-  import Series
-
-  ghost function Exp(x: real): real {
-    ExpSeriesConverges(x);
-    Series.Sum(ExpSeries(x))
-  }
-
-  // The sequence of summands of the exponential series
-  function ExpSeries(gamma: real): nat -> real {
-    (n: nat) => ExpTerm(gamma, n)
-  }
-
-  // ExpTerm(gamma, n) is gamma^n / n!, i.e. the n-th term in the power series of the exponential function.
-  // The start parameter is useful for inductive proofs.
-  opaque function ExpTerm(gamma: real, n: nat, start: nat := 1): real
-    requires 1 <= start
-  {
-    NatArith.FactoralPositive(n, start);
-    RealArith.Pow(gamma, n) / NatArith.Factorial(n, start) as real
-  }
-
-  // Decomposition of ExpTerm useful in inductive proofs.
-  lemma ExpTermStep(gamma: real, n: nat, start: nat := 1)
-    requires start >= 1
-    requires n >= 1
-    ensures ExpTerm(gamma, n, start) == ExpTerm(gamma, n - 1, start + 1) * (gamma / start as real)
-  {
-    NatArith.FactoralPositive(n, start);
-    var denom := NatArith.Factorial(n, start) as real;
-    var denom2 := (start * NatArith.Factorial(n - 1, start + 1)) as real;
-    var denom3 := (start as real) * (NatArith.Factorial(n - 1, start + 1) as real);
-    var numer := RealArith.Pow(gamma, n);
-    var numer2 := gamma * RealArith.Pow(gamma, n - 1);
-    calc {
-      ExpTerm(gamma, n, start);
-      { reveal ExpTerm(); }
-      numer / denom;
-      { reveal RealArith.Pow(); assert numer == numer2; }
-      numer2 / denom;
-      { reveal NatArith.Factorial(); assert denom == denom2; }
-      numer2 / denom2;
-      { assert denom2 == denom3; }
-      numer2 / denom3;
-      { RealArith.FractionOfMul(gamma, RealArith.Pow(gamma, n - 1), start as real, NatArith.Factorial(n - 1, start + 1) as real); }
-      (gamma / start as real) * (RealArith.Pow(gamma, n - 1) / NatArith.Factorial(n - 1, start + 1) as real);
-      { reveal ExpTerm(); }
-      (gamma / start as real) * ExpTerm(gamma, n - 1, start + 1);
-      ExpTerm(gamma, n - 1, start + 1) * (gamma / start as real);
-    }
-  }
-
-  lemma {:axiom} ExpSeriesConverges(x: real)
-    ensures Series.IsSummable(ExpSeries(x))
-
-  lemma {:axiom} FunctionalEquation(x: real, y: real)
-    ensures Exp(x + y) == Exp(x) * Exp(y)
-
-  lemma {:axiom} Increasing(x: real, y: real)
-    requires x < y
-    ensures Exp(x) < Exp(y)
-
-  lemma {:axiom} EvalOne()
-    ensures 2.718281828 <= Exp(1.0) <= 2.718281829
-
-  lemma Positive(x: real)
-    ensures Exp(x) > 0.0
-  {
-    assert Exp(x) >= 0.0 by {
-      var sqrt := Exp(x / 2.0);
-      calc {
-        Exp(x);
-        { FunctionalEquation(x / 2.0, x / 2.0); }
-        sqrt * sqrt;
-      >=
-        0.0;
-      }
-    }
-    if Exp(x) == 0.0 {
-      calc {
-        0.0;
-        Exp(x);
-      < { Increasing(x, x + 1.0); }
-        Exp(x + 1.0);
-        { FunctionalEquation(x, 1.0); }
-        Exp(x) * Exp(1.0);
-      ==
-        0.0;
-      }
-    }
-  }
-
-  lemma EvalZero()
-    ensures Exp(0.0) == 1.0
-  {
-    var one := Exp(0.0);
-    assert one > 0.0 by {
-      Positive(0.0);
-    }
-    assert one * one == one by {
-      FunctionalEquation(0.0, 0.0);
-    }
-  }
-}
diff --git a/src/ProbabilisticProgramming/Loops.dfy b/src/ProbabilisticProgramming/Loops.dfy
index 3ddb0bc3..95e6fc0f 100644
--- a/src/ProbabilisticProgramming/Loops.dfy
+++ b/src/ProbabilisticProgramming/Loops.dfy
@@ -6,7 +6,6 @@
 module Loops {
   import Helper
   import Measures
-  import Limits
   import Monad
   import Quantifier
   import Independence
@@ -364,40 +363,10 @@ module Loops {
     requires Quantifier.WithPosProb(WhileLoopGloballyExitsAfterOneIteration(condition, body))
     ensures WhileTerminatesAlmostSurely(condition, body)
 
-  lemma {:axiom} EnsureWhileTerminatesAlmostSurelyViaLimit<A>(condition: A -> bool, body: A -> Monad.Hurd<A>, init: A)
-    requires Limits.ConvergesTo(WhileCutDivergenceProbability(condition, body, init), 0.0)
-    ensures WhileTerminatesAlmostSurelyInit(condition, body, init)
-  /*
-    Proof strategy:
-
-    Prove that the event that WhileCut terminates grows with the fuel.
-    By monotonicity of probability, the probability is increasing with the fuel.
-    The event that while terminates is the union of WhileCut terminating over all possible values for fuel.
-    By standard measure theory results and the monotonicity of the events,
-    the probability that while terminates is the supremum of the probabilities that WhileCut terminates over all possible values for fuel.
-    Since the probability is increasing, the supremum is the same as the limit.
-  */
-
   ghost function WhileCutProbability<A(!new)>(condition: A -> bool, body: A -> Monad.Hurd<A>, init: A, resultSet: iset<A>): nat -> real {
     (fuel: nat) => Rand.prob(Monad.BitstreamsWithValueIn(WhileCut(condition, body, init, fuel), resultSet))
   }
 
-  lemma {:axiom} WhileProbabilityViaLimit<A(!new)>(condition: A -> bool, body: A -> Monad.Hurd<A>, init: A, resultSet: iset<A>, resultSetRestricted: iset<A>, limit: real)
-    // TODO: we should probably require measurability of condition and independence of body(a) for all a?
-    requires resultSetRestricted == iset a <- resultSet | !condition(a)
-    requires Limits.ConvergesTo(WhileCutProbability(condition, body, init, resultSetRestricted), limit)
-    ensures Rand.prob(Monad.BitstreamsWithValueIn(While(condition, body)(init), resultSet)) == limit
-  /*
-    Proof strategy (similar to EnsureWhileTerminates, can they be unified?):
-
-    Prove that the event that WhileCut yields a result in resultSet violating `condition` grows with the fuel.
-    By monotonicity of probability, the probability is increasing with the fuel.
-    The event that while yields a result in resultsSet (it always violates `condition`) is the union of the events of WhileCut yielding a value in `resultSetRestricted`, over all possible values for fuel.
-    By standard measure theory results and the monotonicity of the events,
-    the probability that while yields a result in resultSet is the supremum of the probabilities that WhileCut yields a result in resultSet violating `condition`, over all possible values for fuel.
-    Since the probability is increasing, the supremum is the same as the limit.
-  */
-
   // Theorem 45 (wrong!) / PROB_BERN_UNTIL (correct!)
   lemma {:axiom} UntilProbabilityFraction<A>(proposal: Monad.Hurd<A>, accept: A -> bool, d: A -> bool)
     requires Independence.IsIndep(proposal)
diff --git a/src/Util/Pos.dfy b/src/Util/Pos.dfy
new file mode 100644
index 00000000..15b8ee37
--- /dev/null
+++ b/src/Util/Pos.dfy
@@ -0,0 +1,6 @@
+
+module Pos {
+
+  type pos = x : nat | x > 0 witness 1
+
+}
\ No newline at end of file
diff --git a/src/interop/java/Full/Random.java b/src/interop/java/Full/Random.java
index 29862596..94b38f75 100644
--- a/src/interop/java/Full/Random.java
+++ b/src/interop/java/Full/Random.java
@@ -36,39 +36,67 @@ public java.math.BigInteger UniformIntervalSample(java.math.BigInteger a, java.m
   }
 
   public java.math.BigInteger UniformSample(java.math.BigInteger n) {
-    return Uniform.Implementation._Companion_Trait.UniformSample(this, n);
+    return DafnyVMCTrait._Companion_RandomTrait.UniformSample(this, n);
   }
 
-  public boolean CoinSample() {
-    return Coin.Implementation._Companion_Trait.CoinSample(this);
+  public boolean BernoulliSample(java.math.BigInteger num, java.math.BigInteger den) {
+    return DafnyVMCTrait._Companion_RandomTrait.BernoulliSample(this, num, den);
+  }  
+
+  public dafny.Tuple2<Boolean, java.math.BigInteger> BernoulliExpNegSampleUnitLoop(java.math.BigInteger num, java.math.BigInteger den, dafny.Tuple2<Boolean, java.math.BigInteger> state) {
+    return DafnyVMCTrait._Companion_RandomTrait.BernoulliExpNegSampleUnitLoop(this, num, den, state);
   }
 
-  public boolean BernoulliSample(Rationals.Rational p) {
-    return Bernoulli.Implementation._Companion_Trait.BernoulliSample(this, p);
-  }  
+  public java.math.BigInteger BernoulliExpNegSampleUnitAux(java.math.BigInteger num, java.math.BigInteger den) {
+    return DafnyVMCTrait._Companion_RandomTrait.BernoulliExpNegSampleUnitAux(this, num, den);
+  }
+
+
+  public boolean BernoulliExpNegSampleUnit(java.math.BigInteger num, java.math.BigInteger den) {
+    return DafnyVMCTrait._Companion_RandomTrait.BernoulliExpNegSampleUnit(this, num, den);
+  }
+
+
+  public boolean BernoulliExpNegSampleGenLoop(java.math.BigInteger iter) {
+    return DafnyVMCTrait._Companion_RandomTrait.BernoulliExpNegSampleGenLoop(this, iter);
+  }
+
+
+  public boolean BernoulliExpNegSample(java.math.BigInteger num, java.math.BigInteger den) {
+    return DafnyVMCTrait._Companion_RandomTrait.BernoulliExpNegSample(this, num, den);
+  }
+
+  public dafny.Tuple2<java.math.BigInteger, Boolean> DiscreteGaussianSampleLoop(java.math.BigInteger num, java.math.BigInteger den, java.math.BigInteger t) {
+    return DafnyVMCTrait._Companion_RandomTrait.DiscreteGaussianSampleLoop(this, num, den, t);
+  }
+
+
+  public java.math.BigInteger DiscreteGaussianSample(java.math.BigInteger num, java.math.BigInteger den) {
+    return DafnyVMCTrait._Companion_RandomTrait.DiscreteGaussianSample(this, num, den);
+  }
 
-  public boolean BernoulliExpNegSampleCaseLe1(Rationals.Rational gamma) {
-    return BernoulliExpNeg.Implementation._Companion_Trait.BernoulliExpNegSampleCaseLe1(this, gamma);
+  public dafny.Tuple2<java.math.BigInteger, Boolean> DiscreteLaplaceSampleLoopIn1Aux(java.math.BigInteger t) {
+    return DafnyVMCTrait._Companion_RandomTrait.DiscreteLaplaceSampleLoopIn1Aux(this ,t);
   }
 
-  public boolean BernoulliExpNegSample(Rationals.Rational gamma) {
-    return BernoulliExpNeg.Implementation._Companion_Trait.BernoulliExpNegSample(this, gamma);
+  public java.math.BigInteger DiscreteLaplaceSampleLoopIn1(java.math.BigInteger t) {
+    return DafnyVMCTrait._Companion_RandomTrait.DiscreteLaplaceSampleLoopIn1(this, t);
   }
 
-  public java.math.BigInteger DiscreteGaussianSample(Rationals.Rational sigma) {
-    return DiscreteGaussian.Implementation._Companion_Trait.DiscreteGaussianSample(this, sigma);
+  public dafny.Tuple2<Boolean, java.math.BigInteger> DiscreteLaplaceSampleLoopIn2Aux(java.math.BigInteger num, java.math.BigInteger den, dafny.Tuple2<Boolean, java.math.BigInteger> K) {
+    return DafnyVMCTrait._Companion_RandomTrait.DiscreteLaplaceSampleLoopIn2Aux(this, num, den, K);
   }
 
-  public dafny.Tuple2<Boolean, java.math.BigInteger> DiscreteLaplaceSampleLoop(Rationals.Rational scale) {
-    return DiscreteLaplace.Implementation._Companion_Trait.DiscreteLaplaceSampleLoop(this, scale);
+  public java.math.BigInteger DiscreteLaplaceSampleLoopIn2(java.math.BigInteger num, java.math.BigInteger den) {
+    return DafnyVMCTrait._Companion_RandomTrait.DiscreteLaplaceSampleLoopIn2(this, num, den);
   }
 
-  public java.math.BigInteger DisceteLaplaceSampleInnerLoop() {
-    return DiscreteLaplace.Implementation._Companion_Trait.DisceteLaplaceSampleInnerLoop(this);
+  public dafny.Tuple2<Boolean, java.math.BigInteger> DiscreteLaplaceSampleLoop(java.math.BigInteger num, java.math.BigInteger den) {
+    return DafnyVMCTrait._Companion_RandomTrait.DiscreteLaplaceSampleLoop(this, num, den);
   }
 
-  public java.math.BigInteger DiscreteLaplaceSample(Rationals.Rational scale) {
-    return DiscreteLaplace.Implementation._Companion_Trait.DiscreteLaplaceSample(this, scale);
+  public java.math.BigInteger DiscreteLaplaceSample(java.math.BigInteger num, java.math.BigInteger den) {
+    return DafnyVMCTrait._Companion_RandomTrait.DiscreteLaplaceSample(this, num, den);
   }
 
   public <__T> void Shuffle(dafny.TypeDescriptor<__T> _td___T, java.lang.Object a) {
diff --git a/tests/Tests.dfy b/tests/Tests.dfy
index 336e901d..6c131ef3 100644
--- a/tests/Tests.dfy
+++ b/tests/Tests.dfy
@@ -4,18 +4,10 @@
  *******************************************************************************/
 
 module Tests {
-  import Rationals
-  import Coin
-  import Uniform
-  import UniformPowerOfTwo
-  import Bernoulli
-  import BernoulliExpNeg
-  import DiscreteLaplace
-  import DiscreteGaussian
   import NatArith
   import RealArith
-  import FisherYates
   import Helper
+  import DafnyVMC
 
   method TestBernoulliIsWithin3SigmaOfTrueMean(
     n: nat,
@@ -51,21 +43,7 @@ module Tests {
     expect diff * diff <= threshold, "Empirical mean should be within 3 sigma of true mean. This individual test may fail with probability of about 6.3e-5.";
   }
 
-  method TestCoin(n: nat, r: Coin.Interface.Trait)
-    requires n > 0
-    modifies r
-  {
-    var t := 0;
-    for i := 0 to n {
-      var b := r.CoinSample();
-      if b {
-        t := t + 1;
-      }
-    }
-    TestBernoulliIsWithin3SigmaOfTrueMean(n, t as real, 0.5, "p(true)");
-  }
-
-  method TestUniformPowerOfTwo(n: nat, u: nat, r: UniformPowerOfTwo.Interface.Trait)
+  method TestUniformPowerOfTwo(n: nat, u: nat, r: DafnyVMC.Random)
     decreases *
     requires n > 0
     requires u > 0
@@ -87,7 +65,7 @@ module Tests {
     TestEmpiricalIsWithin3SigmaOfTrueMean(n, sum as real, (m - 1) as real / 2.0, (m * m - 1) as real / 12.0, "mean of UniformPowerOfTwo(" + Helper.NatToString(u) + ")");
   }
 
-  method TestUniformPowerOfTwoMean(n: nat, u: nat, r: UniformPowerOfTwo.Interface.Trait)
+  method TestUniformPowerOfTwoMean(n: nat, u: nat, r: DafnyVMC.Random)
     decreases *
     requires n > 0
     requires u > 0
@@ -104,7 +82,7 @@ module Tests {
     TestEmpiricalIsWithin3SigmaOfTrueMean(n, sum as real, (m - 1) as real / 2.0, (m * m - 1) as real / 12.0, "mean of UniformPowerOfTwo(" + Helper.NatToString(u) + ")");
   }
 
-  method TestUniform(n: nat, u: nat, r: Uniform.Interface.Trait)
+  method TestUniform(n: nat, u: nat, r: DafnyVMC.Random)
     decreases *
     requires n > 0
     requires u > 0
@@ -125,7 +103,7 @@ module Tests {
     TestEmpiricalIsWithin3SigmaOfTrueMean(n, sum as real, (u - 1) as real / 2.0, (u * u - 1) as real / 12.0, "mean of Uniform(" + Helper.NatToString(u) + ")");
   }
 
-  method TestUniformMean(n: nat, u: nat, r: Uniform.Interface.Trait)
+  method TestUniformMean(n: nat, u: nat, r: DafnyVMC.Random)
     decreases *
     requires n > 0
     requires u > 0
@@ -141,7 +119,7 @@ module Tests {
     TestEmpiricalIsWithin3SigmaOfTrueMean(n, sum as real, (u - 1) as real / 2.0, (u * u - 1) as real / 12.0, "mean of Uniform(" + Helper.NatToString(u) + ")");
   }
 
-  method TestUniformInterval(n: nat, r: Uniform.Interface.Trait)
+  method TestUniformInterval(n: nat, r: DafnyVMC.Random)
     decreases *
     requires n > 0
     modifies r
@@ -162,14 +140,14 @@ module Tests {
     TestBernoulliIsWithin3SigmaOfTrueMean(n, c as real, 1.0 / 3.0, "p(9)");
   }
 
-  method TestBernoulli(n: nat, r: Bernoulli.Interface.Trait)
+  method TestBernoulli(n: nat, r: DafnyVMC.Random)
     decreases *
     requires n > 0
     modifies r
   {
     var t := 0;
     for i := 0 to n {
-      var b := r.BernoulliSample(Rationals.Rational(1, 5));
+      var b := r.BernoulliSample(1, 5);
       if b {
         t := t + 1;
       }
@@ -177,13 +155,13 @@ module Tests {
     TestBernoulliIsWithin3SigmaOfTrueMean(n, t as real, 0.2, "p(true)");
   }
 
-  method TestBernoulli2(n: nat, r: Bernoulli.Interface.Trait)
+  method TestBernoulli2(n: nat, r: DafnyVMC.Random)
     decreases *
     modifies r
   {
     var t := 0;
     for i := 0 to n {
-      var b := r.BernoulliSample(Rationals.Rational(0, 5));
+      var b := r.BernoulliSample(0, 5);
       if b {
         t := t + 1;
       }
@@ -192,13 +170,13 @@ module Tests {
     expect t == 0;
   }
 
-  method TestBernoulli3(n: nat, r: Bernoulli.Interface.Trait)
+  method TestBernoulli3(n: nat, r: DafnyVMC.Random)
     decreases *
     modifies r
   {
     var t := 0;
     for i := 0 to n {
-      var b := r.BernoulliSample(Rationals.Rational(5, 5));
+      var b := r.BernoulliSample(5, 5);
       if b {
         t := t + 1;
       }
@@ -207,14 +185,29 @@ module Tests {
     expect t == n;
   }
 
-  method TestBernoulliExpNeg(n: nat, r: BernoulliExpNeg.Interface.Trait)
+  method TestBernoulliExpNegLe1(n: nat, r: DafnyVMC.Random)
+    decreases *
+    requires n > 0
+    modifies r
+  {
+    var t := 0;
+    for i := 0 to n {
+      var u := r.BernoulliExpNegSample(5108256, 10000000); // about -ln(0.6)
+      if u {
+        t := t + 1;
+      }
+    }
+    TestBernoulliIsWithin3SigmaOfTrueMean(n, t as real, 0.6, "p(true)");
+  }
+
+  method TestBernoulliExpNeg(n: nat, r: DafnyVMC.Random)
     decreases *
     requires n > 0
     modifies r
   {
     var t := 0;
     for i := 0 to n {
-      var u := r.BernoulliExpNegSample(Rationals.Rational(12381, 5377)); // about -ln(0.1)
+      var u := r.BernoulliExpNegSample(12381, 5377); // about -ln(0.1)
       if u {
         t := t + 1;
       }
@@ -222,7 +215,7 @@ module Tests {
     TestBernoulliIsWithin3SigmaOfTrueMean(n, t as real, 0.1, "p(true)");
   }
 
-  method TestDiscreteLaplace(n: nat, r: DiscreteLaplace.Interface.Trait)
+  method TestDiscreteLaplace(n: nat, r: DafnyVMC.Random)
     decreases *
     requires n > 0
     modifies r
@@ -232,7 +225,7 @@ module Tests {
     for i := 0 to n
       invariant -2 in counts && -1 in counts && 0 in counts && 1 in counts && 2 in counts
     {
-      var u := r.DiscreteLaplaceSample(Rationals.Rational(7, 5));
+      var u := r.DiscreteLaplaceSample(7, 5);
       sum := sum + u;
       if u !in counts {
         counts := counts[ u := 1 ];
@@ -253,7 +246,7 @@ module Tests {
     TestEmpiricalIsWithin3SigmaOfTrueMean(n, sum as real, 0.0, variance, "mean");
   }
 
-  method TestDiscreteGaussian(n: nat, r: DiscreteGaussian.Interface.Trait)
+  method TestDiscreteGaussian(n: nat, r: DafnyVMC.Random)
     decreases *
     requires n > 0
     modifies r
@@ -263,7 +256,7 @@ module Tests {
     for i := 0 to n
       invariant -2 in counts && -1 in counts && 0 in counts && 1 in counts && 2 in counts
     {
-      var u := r.DiscreteGaussianSample(Rationals.Rational(7, 5));
+      var u := r.DiscreteGaussianSample(7, 5);
       sum := sum + u;
       if u !in counts {
         counts := counts[ u := 1 ];
@@ -284,7 +277,7 @@ module Tests {
     TestEmpiricalIsWithin3SigmaOfTrueMean(n, sum as real, 0.0, varianceBound, "mean");
   }
 
-  method TestFisherYates<T(==, !new)>(n: nat, a: array<T>, r: FisherYates.Interface.Trait, printer: ((T, nat)) -> string) 
+  method TestFisherYates<T(==, !new)>(n: nat, a: array<T>, r: DafnyVMC.Random, printer: ((T, nat)) -> string) 
     decreases *
     modifies r
     modifies a
diff --git a/tests/TestsRandom.dfy b/tests/TestsRandom.dfy
index c1f45d32..eab4340b 100644
--- a/tests/TestsRandom.dfy
+++ b/tests/TestsRandom.dfy
@@ -9,11 +9,6 @@ module TestsRandom {
   import Helper
   import Tests
 
-  method {:test} TestCoin() {
-    var r := new DafnyVMC.Random();
-    Tests.TestCoin(1_000_000, r);
-  }
-
   method {:test} TestUniformPowerOfTwo_10()
     decreases *
   {
@@ -87,6 +82,13 @@ module TestsRandom {
     Tests.TestBernoulli3(1_000_000, r);
   }
 
+  method {:test} TestBernoulliExpNegLe1()
+    decreases *
+  {
+    var r := new DafnyVMC.Random();
+    Tests.TestBernoulliExpNegLe1(1_000_000, r);
+  }
+
   method {:test} TestBernoulliExpNeg()
     decreases *
   {