diff --git a/src/Distributions/Bernoulli/Implementation.dfy b/src/Distributions/Bernoulli/Implementation.dfy
index 79817f65..4d213504 100644
--- a/src/Distributions/Bernoulli/Implementation.dfy
+++ b/src/Distributions/Bernoulli/Implementation.dfy
@@ -12,7 +12,7 @@ module Bernoulli.Implementation {
   trait {:termination false} Trait extends Interface.Trait {
 
     method BernoulliSample(p: Rationals.Rational) returns (c: bool)
-      modifies this
+      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/Interface.dfy b/src/Distributions/Bernoulli/Interface.dfy
index fefeefaa..c640b6e5 100644
--- a/src/Distributions/Bernoulli/Interface.dfy
+++ b/src/Distributions/Bernoulli/Interface.dfy
@@ -12,7 +12,7 @@ module Bernoulli.Interface {
   trait {:termination false} Trait extends Uniform.Interface.Trait {
 
     method BernoulliSample(p: Rationals.Rational) returns (c: bool)
-      modifies this
+      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/BernoulliExpNeg/Implementation.dfy b/src/Distributions/BernoulliExpNeg/Implementation.dfy
index e8233dc4..ea90cf75 100644
--- a/src/Distributions/BernoulliExpNeg/Implementation.dfy
+++ b/src/Distributions/BernoulliExpNeg/Implementation.dfy
@@ -17,7 +17,7 @@ module BernoulliExpNeg.Implementation {
 
     // Based on Algorithm 1 in https://arxiv.org/pdf/2004.00010.pdf; unverified
     method BernoulliExpNegSample(gamma: Rationals.Rational) returns (c: bool)
-      modifies this
+      modifies `s
       requires gamma.numer >= 0
       decreases *
       ensures Monad.Result(c, s) == Model.Sample(gamma)(old(s))
@@ -42,7 +42,7 @@ module BernoulliExpNeg.Implementation {
     }
 
     method BernoulliExpNegSampleCaseLe1(gamma: Rationals.Rational) returns (c: bool)
-      modifies this
+      modifies `s
       requires 0 <= gamma.numer <= gamma.denom
       decreases *
       ensures Monad.Result(c, s) == Model.SampleLe1(gamma)(old(s))
diff --git a/src/Distributions/BernoulliExpNeg/Interface.dfy b/src/Distributions/BernoulliExpNeg/Interface.dfy
index 79e175bd..9d5ab259 100644
--- a/src/Distributions/BernoulliExpNeg/Interface.dfy
+++ b/src/Distributions/BernoulliExpNeg/Interface.dfy
@@ -12,7 +12,7 @@ module BernoulliExpNeg.Interface {
   trait {:termination false} Trait extends Bernoulli.Interface.Trait {
 
     method BernoulliExpNegSample(gamma: Rationals.Rational) returns (c: bool)
-      modifies this
+      modifies `s
       decreases *
       requires gamma.numer >= 0
       ensures Monad.Result(c, s) == Model.Sample(gamma)(old(s))
diff --git a/src/Distributions/Coin/Implementation.dfy b/src/Distributions/Coin/Implementation.dfy
index d0ec31e4..3bec2089 100644
--- a/src/Distributions/Coin/Implementation.dfy
+++ b/src/Distributions/Coin/Implementation.dfy
@@ -12,7 +12,7 @@ module Coin.Implementation {
   trait {:termination false} Trait extends Interface.Trait {
 
     method CoinSample() returns (b: bool)
-      modifies this
+      modifies `s
       ensures Model.Sample(old(s)) == Monad.Result(b, s)
     {
       var x := UniformPowerOfTwoSample(2);
diff --git a/src/Distributions/Coin/Interface.dfy b/src/Distributions/Coin/Interface.dfy
index eca0c61b..c506c519 100644
--- a/src/Distributions/Coin/Interface.dfy
+++ b/src/Distributions/Coin/Interface.dfy
@@ -12,7 +12,7 @@ module Coin.Interface {
   trait {:termination false} Trait extends UniformPowerOfTwo.Interface.Trait {
 
     method CoinSample() returns (b: bool)
-      modifies this
+      modifies `s
       ensures Model.Sample(old(s)) == Monad.Result(b, s)
 
   }
diff --git a/src/Distributions/DiscreteGaussian/Implementation.dfy b/src/Distributions/DiscreteGaussian/Implementation.dfy
index fa9cd248..503df92b 100644
--- a/src/Distributions/DiscreteGaussian/Implementation.dfy
+++ b/src/Distributions/DiscreteGaussian/Implementation.dfy
@@ -12,7 +12,7 @@ module DiscreteGaussian.Implementation {
     // 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 this
+      modifies `s
       requires sigma.numer >= 1
       decreases *
     {
diff --git a/src/Distributions/DiscreteGaussian/Interface.dfy b/src/Distributions/DiscreteGaussian/Interface.dfy
index f6938556..1f937adf 100644
--- a/src/Distributions/DiscreteGaussian/Interface.dfy
+++ b/src/Distributions/DiscreteGaussian/Interface.dfy
@@ -11,7 +11,7 @@ module DiscreteGaussian.Interface {
   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 this
+      modifies `s
       requires sigma.numer >= 1
       decreases *
   }
diff --git a/src/Distributions/DiscreteLaplace/Implementation.dfy b/src/Distributions/DiscreteLaplace/Implementation.dfy
index 92b1d30b..a44cad09 100644
--- a/src/Distributions/DiscreteLaplace/Implementation.dfy
+++ b/src/Distributions/DiscreteLaplace/Implementation.dfy
@@ -19,7 +19,7 @@ module DiscreteLaplace.Implementation {
 
     // Based on Algorithm 2 in https://arxiv.org/pdf/2004.00010.pdf; unverified
     method DiscreteLaplaceSample(scale: Rationals.Rational) returns (z: int)
-      modifies this
+      modifies `s
       requires scale.numer >= 1
       decreases *
       ensures Model.Sample(scale)(old(s)) == Monad.Result(z, s)
@@ -30,7 +30,7 @@ module DiscreteLaplace.Implementation {
     }
 
     method {:rlimit 100000} DiscreteLaplaceSampleLoop(scale: Rationals.Rational) returns (bY: (bool, int))
-      modifies this
+      modifies `s
       requires scale.numer >= 1
       decreases *
       ensures Model.SampleLoop(scale)(old(s)) == Monad.Result(bY, s)
@@ -61,7 +61,7 @@ module DiscreteLaplace.Implementation {
     }
 
     method DisceteLaplaceSampleInnerLoop() returns (v: int)
-      modifies this
+      modifies `s
       decreases *
       ensures Model.SampleInnerLoopFull()(old(s)) == Monad.Result(v, s)
     {
diff --git a/src/Distributions/DiscreteLaplace/Interface.dfy b/src/Distributions/DiscreteLaplace/Interface.dfy
index 7b2be999..60fdddfa 100644
--- a/src/Distributions/DiscreteLaplace/Interface.dfy
+++ b/src/Distributions/DiscreteLaplace/Interface.dfy
@@ -13,7 +13,7 @@ module DiscreteLaplace.Interface {
 
     // Based on Algorithm 2 in https://arxiv.org/pdf/2004.00010.pdf; unverified
     method DiscreteLaplaceSample(scale: Rationals.Rational) returns (z: int)
-      modifies this
+      modifies `s
       requires scale.numer >= 1
       decreases *
 
diff --git a/src/Distributions/Uniform/Implementation.dfy b/src/Distributions/Uniform/Implementation.dfy
index 2538c9e7..7d3bf8b4 100644
--- a/src/Distributions/Uniform/Implementation.dfy
+++ b/src/Distributions/Uniform/Implementation.dfy
@@ -12,7 +12,7 @@ module Uniform.Implementation {
 
   trait {:termination false} Trait extends Interface.Trait {
     method UniformSample(n: nat) returns (u: nat)
-      modifies this
+      modifies `s
       decreases *
       requires n > 0
       ensures u < n
diff --git a/src/Distributions/Uniform/Interface.dfy b/src/Distributions/Uniform/Interface.dfy
index 1095581b..617175cd 100644
--- a/src/Distributions/Uniform/Interface.dfy
+++ b/src/Distributions/Uniform/Interface.dfy
@@ -12,14 +12,14 @@ module Uniform.Interface {
   trait {:termination false} Trait extends UniformPowerOfTwo.Interface.Trait {
 
     method UniformSample(n: nat) returns (u: nat)
-      modifies this
+      modifies `s
       decreases *
       requires n > 0
       ensures u < n
       ensures Model.Sample(n)(old(s)) == Monad.Result(u, s)
 
     method UniformIntervalSample(a: int, b: int) returns (u: int)
-      modifies this
+      modifies `s
       decreases *
       requires a < b
       ensures a <= u < b
diff --git a/src/Distributions/UniformPowerOfTwo/Implementation.dfy b/src/Distributions/UniformPowerOfTwo/Implementation.dfy
index ce81bfaf..063cdcfd 100644
--- a/src/Distributions/UniformPowerOfTwo/Implementation.dfy
+++ b/src/Distributions/UniformPowerOfTwo/Implementation.dfy
@@ -11,7 +11,7 @@ module UniformPowerOfTwo.Implementation {
   trait {:termination false} Trait extends Interface.Trait {
     method UniformPowerOfTwoSample(n: nat) returns (u: nat)
       requires n >= 1
-      modifies this
+      modifies `s
       ensures Model.Sample(n)(old(s)) == Monad.Result(u, s)
 
   }
diff --git a/src/Distributions/UniformPowerOfTwo/Interface.dfy b/src/Distributions/UniformPowerOfTwo/Interface.dfy
index 17f78c62..77c6f356 100644
--- a/src/Distributions/UniformPowerOfTwo/Interface.dfy
+++ b/src/Distributions/UniformPowerOfTwo/Interface.dfy
@@ -14,7 +14,7 @@ module UniformPowerOfTwo.Interface {
     // The return value u is uniformly distributed between 0 <= u < 2^k where 2^k <= n < 2^(k + 1).
     method UniformPowerOfTwoSample(n: nat) returns (u: nat)
       requires n >= 1
-      modifies this
+      modifies `s
       ensures Model.Sample(n)(old(s)) == Monad.Result(u, s)
 
   }
diff --git a/src/Util/FisherYates/Equivalence.dfy b/src/Util/FisherYates/Equivalence.dfy
new file mode 100644
index 00000000..cf20a1cd
--- /dev/null
+++ b/src/Util/FisherYates/Equivalence.dfy
@@ -0,0 +1,19 @@
+/*******************************************************************************
+ *  Copyright by the contributors to the Dafny Project
+ *  SPDX-License-Identifier: MIT
+ *******************************************************************************/
+
+module FisherYates.Equivalence {
+  import Model
+  import Rand
+
+  ghost predicate LoopInvariant<T>(prevI: nat, i: nat, a: array<T>, prevASeq: seq<T>, oldASeq: seq<T>, oldS: Rand.Bitstream, prevS: Rand.Bitstream, s: Rand.Bitstream)
+    reads a
+  {
+    && prevI <= |prevASeq|
+    && i <= a.Length - 1
+    && Model.Shuffle(oldASeq)(oldS) == Model.Shuffle(prevASeq, prevI)(prevS)
+    && Model.Shuffle(prevASeq, prevI)(prevS) == Model.Shuffle(a[..], i)(s)
+  }
+
+}
\ No newline at end of file
diff --git a/src/Util/FisherYates/Implementation.dfy b/src/Util/FisherYates/Implementation.dfy
index e58bdcb6..20a1b5b1 100644
--- a/src/Util/FisherYates/Implementation.dfy
+++ b/src/Util/FisherYates/Implementation.dfy
@@ -5,20 +5,41 @@
 
 module FisherYates.Implementation {
   import Interface
+  import Monad
+  import Model
+  import Uniform
+  import Equivalence
 
   trait {:termination false} Trait extends Interface.Trait {
 
     method Shuffle<T>(a: array<T>)
       decreases *
-      modifies this, a
+      modifies `s, a
+      ensures Model.Shuffle(old(a[..]))(old(s)) == Monad.Result(a[..], s)
     {
+      ghost var prevI, prevASeq, prevS := 0, a[..], s; // ghost
       if a.Length > 1 {
-        for i := 0 to a.Length - 1 {
+        for i := 0 to a.Length - 1
+          invariant Equivalence.LoopInvariant(prevI, i, a, prevASeq, old(a[..]), old(s), prevS, s)
+        {
+          prevI, prevASeq, prevS := i, a[..], s; // ghost
           var j := UniformIntervalSample(i, a.Length);
-          a[i], a[j] := a[j], a[i];
+          assert prevASeq == a[..]; // ghost
+          Swap(a, i, j);
         }
       }
     }
 
+    method Swap<T>(a: array<T>, i: nat, j: nat)
+      modifies a
+      requires i <= j
+      requires 0 <= i < a.Length
+      requires 0 <= j < a.Length
+      ensures Model.Swap(old(a[..]), i, j) == a[..]
+      ensures old(s) == s
+    {
+      a[i], a[j] := a[j], a[i];
+    }
+
   }
 }
\ No newline at end of file
diff --git a/src/Util/FisherYates/Interface.dfy b/src/Util/FisherYates/Interface.dfy
index 7375474d..1571870d 100644
--- a/src/Util/FisherYates/Interface.dfy
+++ b/src/Util/FisherYates/Interface.dfy
@@ -5,12 +5,15 @@
 
 module FisherYates.Interface {
   import Uniform
+  import Model
+  import Monad
 
   trait {:termination false} Trait extends Uniform.Interface.Trait {
 
     method Shuffle<T>(a: array<T>)
       decreases *
-      modifies this, a
+      modifies `s, a
+      ensures Model.Shuffle(old(a[..]))(old(s)) == Monad.Result(a[..], s)
 
   }
 }
\ No newline at end of file
diff --git a/src/Util/FisherYates/Model.dfy b/src/Util/FisherYates/Model.dfy
index a43682f7..af72424d 100644
--- a/src/Util/FisherYates/Model.dfy
+++ b/src/Util/FisherYates/Model.dfy
@@ -21,8 +21,8 @@ module FisherYates.Model {
     (s: Rand.Bitstream) =>
       if |xs[i..]| > 1 then
         var (j, s) :- Uniform.Model.IntervalSample(i, |xs|)(s);
-        var xs := Swap(xs, i, j);
-        Shuffle(xs, i + 1)(s)
+        var ys := Swap(xs, i, j);
+        Shuffle(ys, i + 1)(s)
       else
         Monad.Return(xs)(s)
   }
diff --git a/src/interop/java/Full/Random.java b/src/interop/java/Full/Random.java
index 18eec74e..9652b1e3 100644
--- a/src/interop/java/Full/Random.java
+++ b/src/interop/java/Full/Random.java
@@ -101,4 +101,9 @@ public void Shuffle(long[] arr) {
   public void Shuffle(short[] arr) {
     FisherYates.Implementation._Companion_Trait.Shuffle(TypeDescriptor.SHORT, this, arr);
   }
+
+  public <__T> void Swap(dafny.TypeDescriptor<__T> _td___T, java.lang.Object a, java.math.BigInteger i, java.math.BigInteger j) {
+    FisherYates.Implementation._Companion_Trait.Swap(_td___T, this, a, i, j);
+  }
+    
 }
\ No newline at end of file