diff --git a/library/src/tests/intrinsic.rs b/library/src/tests/intrinsic.rs
index 01ad9cdb48..944676e72c 100644
--- a/library/src/tests/intrinsic.rs
+++ b/library/src/tests/intrinsic.rs
@@ -1,6 +1,8 @@
 // Copyright (c) Microsoft Corporation.
 // Licensed under the MIT License.
 
+#![allow(clippy::needless_raw_string_hashes, clippy::too_many_lines)]
+
 use expect_test::expect;
 use indoc::indoc;
 use qsc::{interpret::Value, target::Profile, SparseSim};
@@ -2825,3 +2827,184 @@ fn global_phase_correct_for_triply_controlled_adjoint_r1() {
     "#]]
     .assert_eq(&dump);
 }
+
+#[test]
+fn test_exp() {
+    let dump = test_expression(
+        indoc! {r#"
+        {
+            open Microsoft.Quantum.Math;
+            open Microsoft.Quantum.Diagnostics;
+            for p in [PauliX, PauliY, PauliZ, PauliI] {
+                for i in 1 .. 4 {
+                    Message($"Exp with {p} on {i} qubits:");
+                    use qs = Qubit[i];
+                    for q in qs {
+                        H(q);
+                    }
+                    Exp(Repeated(p, i), PI() / 7.0, qs);
+                    DumpMachine();
+                    ResetAll(qs);
+                }
+            }
+        }
+        "#},
+        &Value::unit(),
+    );
+
+    expect![[r#"
+        Exp with PauliX on 1 qubits:
+        STATE:
+        |0⟩: 0.6371+0.3068𝑖
+        |1⟩: 0.6371+0.3068𝑖
+        Exp with PauliX on 2 qubits:
+        STATE:
+        |00⟩: 0.4505+0.2169𝑖
+        |01⟩: 0.4505+0.2169𝑖
+        |10⟩: 0.4505+0.2169𝑖
+        |11⟩: 0.4505+0.2169𝑖
+        Exp with PauliX on 3 qubits:
+        STATE:
+        |000⟩: 0.3185+0.1534𝑖
+        |001⟩: 0.3185+0.1534𝑖
+        |010⟩: 0.3185+0.1534𝑖
+        |011⟩: 0.3185+0.1534𝑖
+        |100⟩: 0.3185+0.1534𝑖
+        |101⟩: 0.3185+0.1534𝑖
+        |110⟩: 0.3185+0.1534𝑖
+        |111⟩: 0.3185+0.1534𝑖
+        Exp with PauliX on 4 qubits:
+        STATE:
+        |0000⟩: 0.2252+0.1085𝑖
+        |0001⟩: 0.2252+0.1085𝑖
+        |0010⟩: 0.2252+0.1085𝑖
+        |0011⟩: 0.2252+0.1085𝑖
+        |0100⟩: 0.2252+0.1085𝑖
+        |0101⟩: 0.2252+0.1085𝑖
+        |0110⟩: 0.2252+0.1085𝑖
+        |0111⟩: 0.2252+0.1085𝑖
+        |1000⟩: 0.2252+0.1085𝑖
+        |1001⟩: 0.2252+0.1085𝑖
+        |1010⟩: 0.2252+0.1085𝑖
+        |1011⟩: 0.2252+0.1085𝑖
+        |1100⟩: 0.2252+0.1085𝑖
+        |1101⟩: 0.2252+0.1085𝑖
+        |1110⟩: 0.2252+0.1085𝑖
+        |1111⟩: 0.2252+0.1085𝑖
+        Exp with PauliY on 1 qubits:
+        STATE:
+        |0⟩: 0.9439+0.0000𝑖
+        |1⟩: 0.3303+0.0000𝑖
+        Exp with PauliY on 2 qubits:
+        STATE:
+        |00⟩: 0.4505−0.2169𝑖
+        |01⟩: 0.4505+0.2169𝑖
+        |10⟩: 0.4505+0.2169𝑖
+        |11⟩: 0.4505−0.2169𝑖
+        Exp with PauliY on 3 qubits:
+        STATE:
+        |000⟩: 0.1651+0.0000𝑖
+        |001⟩: 0.4719+0.0000𝑖
+        |010⟩: 0.4719+0.0000𝑖
+        |011⟩: 0.1651+0.0000𝑖
+        |100⟩: 0.4719+0.0000𝑖
+        |101⟩: 0.1651+0.0000𝑖
+        |110⟩: 0.1651+0.0000𝑖
+        |111⟩: 0.4719+0.0000𝑖
+        Exp with PauliY on 4 qubits:
+        STATE:
+        |0000⟩: 0.2252+0.1085𝑖
+        |0001⟩: 0.2252−0.1085𝑖
+        |0010⟩: 0.2252−0.1085𝑖
+        |0011⟩: 0.2252+0.1085𝑖
+        |0100⟩: 0.2252−0.1085𝑖
+        |0101⟩: 0.2252+0.1085𝑖
+        |0110⟩: 0.2252+0.1085𝑖
+        |0111⟩: 0.2252−0.1085𝑖
+        |1000⟩: 0.2252−0.1085𝑖
+        |1001⟩: 0.2252+0.1085𝑖
+        |1010⟩: 0.2252+0.1085𝑖
+        |1011⟩: 0.2252−0.1085𝑖
+        |1100⟩: 0.2252+0.1085𝑖
+        |1101⟩: 0.2252−0.1085𝑖
+        |1110⟩: 0.2252−0.1085𝑖
+        |1111⟩: 0.2252+0.1085𝑖
+        Exp with PauliZ on 1 qubits:
+        STATE:
+        |0⟩: 0.6371+0.3068𝑖
+        |1⟩: 0.6371−0.3068𝑖
+        Exp with PauliZ on 2 qubits:
+        STATE:
+        |00⟩: 0.4505+0.2169𝑖
+        |01⟩: 0.4505−0.2169𝑖
+        |10⟩: 0.4505−0.2169𝑖
+        |11⟩: 0.4505+0.2169𝑖
+        Exp with PauliZ on 3 qubits:
+        STATE:
+        |000⟩: 0.3185+0.1534𝑖
+        |001⟩: 0.3185−0.1534𝑖
+        |010⟩: 0.3185−0.1534𝑖
+        |011⟩: 0.3185+0.1534𝑖
+        |100⟩: 0.3185−0.1534𝑖
+        |101⟩: 0.3185+0.1534𝑖
+        |110⟩: 0.3185+0.1534𝑖
+        |111⟩: 0.3185−0.1534𝑖
+        Exp with PauliZ on 4 qubits:
+        STATE:
+        |0000⟩: 0.2252+0.1085𝑖
+        |0001⟩: 0.2252−0.1085𝑖
+        |0010⟩: 0.2252−0.1085𝑖
+        |0011⟩: 0.2252+0.1085𝑖
+        |0100⟩: 0.2252−0.1085𝑖
+        |0101⟩: 0.2252+0.1085𝑖
+        |0110⟩: 0.2252+0.1085𝑖
+        |0111⟩: 0.2252−0.1085𝑖
+        |1000⟩: 0.2252−0.1085𝑖
+        |1001⟩: 0.2252+0.1085𝑖
+        |1010⟩: 0.2252+0.1085𝑖
+        |1011⟩: 0.2252−0.1085𝑖
+        |1100⟩: 0.2252+0.1085𝑖
+        |1101⟩: 0.2252−0.1085𝑖
+        |1110⟩: 0.2252−0.1085𝑖
+        |1111⟩: 0.2252+0.1085𝑖
+        Exp with PauliI on 1 qubits:
+        STATE:
+        |0⟩: 0.6371+0.3068𝑖
+        |1⟩: 0.6371+0.3068𝑖
+        Exp with PauliI on 2 qubits:
+        STATE:
+        |00⟩: 0.4505+0.2169𝑖
+        |01⟩: 0.4505+0.2169𝑖
+        |10⟩: 0.4505+0.2169𝑖
+        |11⟩: 0.4505+0.2169𝑖
+        Exp with PauliI on 3 qubits:
+        STATE:
+        |000⟩: 0.3185+0.1534𝑖
+        |001⟩: 0.3185+0.1534𝑖
+        |010⟩: 0.3185+0.1534𝑖
+        |011⟩: 0.3185+0.1534𝑖
+        |100⟩: 0.3185+0.1534𝑖
+        |101⟩: 0.3185+0.1534𝑖
+        |110⟩: 0.3185+0.1534𝑖
+        |111⟩: 0.3185+0.1534𝑖
+        Exp with PauliI on 4 qubits:
+        STATE:
+        |0000⟩: 0.2252+0.1085𝑖
+        |0001⟩: 0.2252+0.1085𝑖
+        |0010⟩: 0.2252+0.1085𝑖
+        |0011⟩: 0.2252+0.1085𝑖
+        |0100⟩: 0.2252+0.1085𝑖
+        |0101⟩: 0.2252+0.1085𝑖
+        |0110⟩: 0.2252+0.1085𝑖
+        |0111⟩: 0.2252+0.1085𝑖
+        |1000⟩: 0.2252+0.1085𝑖
+        |1001⟩: 0.2252+0.1085𝑖
+        |1010⟩: 0.2252+0.1085𝑖
+        |1011⟩: 0.2252+0.1085𝑖
+        |1100⟩: 0.2252+0.1085𝑖
+        |1101⟩: 0.2252+0.1085𝑖
+        |1110⟩: 0.2252+0.1085𝑖
+        |1111⟩: 0.2252+0.1085𝑖
+    "#]]
+    .assert_eq(&dump);
+}
diff --git a/library/std/internal.qs b/library/std/internal.qs
index 2b80f00ef3..2c4f9430fd 100644
--- a/library/std/internal.qs
+++ b/library/std/internal.qs
@@ -241,13 +241,28 @@ namespace Microsoft.Quantum.Intrinsic {
     }
 
     internal operation SpreadZ(from : Qubit, to : Qubit[]) : Unit is Adj {
-        if (Length(to) > 0) {
-            if (Length(to) > 1) {
-                let half = Length(to) / 2;
-                SpreadZ(to[0], to[half + 1..Length(to) - 1]);
-                SpreadZ(from, to[1..half]);
+        let targets = GetSpread(from, to);
+        for (ctl, tgt) in targets {
+            CNOT(ctl, tgt);
+        }
+    }
+
+    internal function GetSpread(from : Qubit, to : Qubit[]) : (Qubit, Qubit)[] {
+        mutable queue = [(from, to)];
+        mutable targets = [];
+        while Length(queue) > 0 {
+            mutable (next, rest) = (queue[0], queue[1...]);
+            set queue = rest;
+            let (next_from, next_to) = next;
+            if Length(next_to) > 0 {
+                set targets = [(next_to[0], next_from)] + targets;
+                if Length(next_to) > 1 {
+                    let half = Length(next_to) / 2;
+                    set queue = [(next_from, next_to[1..half]), (next_to[0], next_to[(half + 1)...])] + rest;
+                }
             }
-            CNOT(to[0], from);
         }
+
+        targets
     }
 }