From 920dd7160aa4d49aa39ee517ae4f504b4366f4e9 Mon Sep 17 00:00:00 2001
From: sezna <alehansen@microsoft.com>
Date: Tue, 10 Sep 2024 14:24:00 -0700
Subject: [PATCH] migrate canon and undo unstable migration

---
 library/src/lib.rs                            |  20 +-
 library/std/qsharp.json                       |  12 +-
 library/std/src/Std/Arrays.qs                 |  44 --
 library/std/src/Std/Canon.qs                  | 588 +++++++++++++++
 library/std/src/Std/Convert.qs                |  13 -
 library/std/src/Std/Diagnostics.qs            |  11 -
 library/std/src/Std/Intrinsic.qs              |  29 -
 library/std/src/Std/Math.qs                   |  81 --
 library/std/src/Std/Measurement.qs            |   6 -
 library/std/src/Std/Random.qs                 |   2 -
 library/std/src/Std/ResourceEstimation.qs     |  12 -
 library/std/src/Std/Unstable/Arithmetic.qs    | 702 ------------------
 .../std/src/Std/Unstable/ArithmeticHelpers.qs | 586 ---------------
 .../std/src/Std/Unstable/StatePreparation.qs  | 390 ----------
 library/std/src/Std/Unstable/TableLookup.qs   | 281 -------
 library/std/src/canon.qs                      | 589 ---------------
 library/std/src/legacy_api.qs                 |   6 +-
 library/std/src/unstable_arithmetic.qs        | 701 +++++++++++++++++
 .../std/src/unstable_arithmetic_internal.qs   | 585 +++++++++++++++
 library/std/src/unstable_state_preparation.qs | 388 ++++++++++
 library/std/src/unstable_table_lookup.qs      | 278 +++++++
 21 files changed, 2557 insertions(+), 2767 deletions(-)
 create mode 100644 library/std/src/Std/Canon.qs
 delete mode 100644 library/std/src/Std/Unstable/Arithmetic.qs
 delete mode 100644 library/std/src/Std/Unstable/ArithmeticHelpers.qs
 delete mode 100644 library/std/src/Std/Unstable/StatePreparation.qs
 delete mode 100644 library/std/src/Std/Unstable/TableLookup.qs
 delete mode 100644 library/std/src/canon.qs
 create mode 100644 library/std/src/unstable_arithmetic.qs
 create mode 100644 library/std/src/unstable_arithmetic_internal.qs
 create mode 100644 library/std/src/unstable_state_preparation.qs
 create mode 100644 library/std/src/unstable_table_lookup.qs

diff --git a/library/src/lib.rs b/library/src/lib.rs
index f41b02b74e..c0f28dfeda 100644
--- a/library/src/lib.rs
+++ b/library/src/lib.rs
@@ -26,8 +26,8 @@ pub const STD_LIB: &[(&str, &str)] = &[
         include_str!("../std/src/Std/Arrays.qs"),
     ),
     (
-        "qsharp-library-source:canon.qs",
-        include_str!("../std/src/canon.qs"),
+        "qsharp-library-source:Std/Canon.qs",
+        include_str!("../std/src/Std/Canon.qs"),
     ),
     (
         "qsharp-library-source:Std/Convert.qs",
@@ -74,20 +74,20 @@ pub const STD_LIB: &[(&str, &str)] = &[
         include_str!("../std/src/Std/ResourceEstimation.qs"),
     ),
     (
-        "qsharp-library-source:Std/Unstable/Arithmetic.qs",
-        include_str!("../std/src/Std/Unstable/Arithmetic.qs"),
+        "qsharp-library-source:unstable_arithmetic_internal.qs",
+        include_str!("../std/src/unstable_arithmetic_internal.qs"),
     ),
     (
-        "qsharp-library-source:Std/Unstable/ArithmeticHelpers.qs",
-        include_str!("../std/src/Std/Unstable/ArithmeticHelpers.qs"),
+        "qsharp-library-source:unstable_arithmetic.qs",
+        include_str!("../std/src/unstable_arithmetic.qs"),
     ),
     (
-        "qsharp-library-source:Std/Unstable/StatePreparation.qs",
-        include_str!("../std/src/Std/Unstable/StatePreparation.qs"),
+        "qsharp-library-source:unstable_state_preparation.qs",
+        include_str!("../std/src/unstable_state_preparation.qs"),
     ),
     (
-        "qsharp-library-source:Std/Unstable/TableLookup.qs",
-        include_str!("../std/src/Std/Unstable/TableLookup.qs"),
+        "qsharp-library-source:unstable_table_lookup.qs",
+        include_str!("../std/src/unstable_table_lookup.qs"),
     ),
     (
         "qsharp-library-source:legacy_api.qs",
diff --git a/library/std/qsharp.json b/library/std/qsharp.json
index 4c27d0d140..e3a857383e 100644
--- a/library/std/qsharp.json
+++ b/library/std/qsharp.json
@@ -2,11 +2,15 @@
   "author": "Microsoft",
   "license": "MIT",
   "files": [
-    "src/canon.qs",
     "src/core.qs",
     "src/legacy_api.qs",
+    "src/unstable_arithmetic.qs",
+    "src/unstable_arithmetic_internal.qs",
+    "src/unstable_state_preparation.qs",
+    "src/unstable_table_lookup.qs",
     "src/QIR/Intrinsic.qs",
     "src/Std/Arrays.qs",
+    "src/Std/Canon.qs",
     "src/Std/Convert.qs",
     "src/Std/Diagnostics.qs",
     "src/Std/InternalHelpers.qs",
@@ -15,10 +19,6 @@
     "src/Std/Math.qs",
     "src/Std/Measurement.qs",
     "src/Std/Random.qs",
-    "src/Std/ResourceEstimation.qs",
-    "src/Std/Unstable/Arithmetic.qs",
-    "src/Std/Unstable/ArithmeticHelpers.qs",
-    "src/Std/Unstable/StatePreparation.qs",
-    "src/Std/Unstable/TableLookup.qs"
+    "src/Std/ResourceEstimation.qs"
   ]
 }
\ No newline at end of file
diff --git a/library/std/src/Std/Arrays.qs b/library/std/src/Std/Arrays.qs
index 3253f33a1d..f2a572b96b 100644
--- a/library/std/src/Std/Arrays.qs
+++ b/library/std/src/Std/Arrays.qs
@@ -37,7 +37,6 @@ function All<'T>(predicate : ('T -> Bool), array : 'T[]) : Bool {
     true
 }
 
-    /// # Summary
 /// Given an array and a predicate that is defined
 /// for the elements of the array, checks if at least one element of
 /// the array satisfies the predicate.
@@ -69,7 +68,6 @@ function Any<'T>(predicate : ('T -> Bool), array : 'T[]) : Bool {
     false
 }
 
-    /// # Summary
 /// Splits an array into multiple parts of equal length.
 ///
 /// # Input
@@ -97,7 +95,6 @@ function Chunks<'T>(chunkSize : Int, array : 'T[]) : 'T[][] {
     output
 }
 
-    /// # Summary
 /// Shift an array circularly left or right by a specific step size.
 ///
 /// # Type Parameters
@@ -145,7 +142,6 @@ function CircularlyShifted<'T>(stepCount : Int, array : 'T[]) : 'T[] {
     rightPart + leftPart
 }
 
-    /// # Summary
 /// Extracts a column from a matrix.
 ///
 /// # Description
@@ -182,7 +178,6 @@ function ColumnAt<'T>(column : Int, matrix : 'T[][]) : 'T[] {
     columnValues
 }
 
-    /// # Summary
 /// Given an array and a predicate that is defined
 /// for the elements of the array, returns the number of elements
 /// an array that consists of those elements that satisfy the predicate.
@@ -215,7 +210,6 @@ function Count<'T>(predicate : ('T -> Bool), array : 'T[]) : Int {
     count
 }
 
-    /// # Summary
 /// Returns an array of diagonal elements of a 2-dimensional array
 ///
 /// # Description
@@ -252,7 +246,6 @@ function Diagonal<'T>(matrix : 'T[][]) : 'T[] {
     diagonal
 }
 
-    /// # Summary
 /// Repeats an operation for a given number of samples, collecting its outputs
 /// in an array.
 ///
@@ -284,7 +277,6 @@ operation DrawMany<'TInput, 'TOutput>(op : ('TInput => 'TOutput), nSamples : Int
     outputs
 }
 
-    /// # Summary
 /// Given an array, returns a new array containing elements of the original
 /// array along with the indices of each element.
 ///
@@ -313,7 +305,6 @@ function Enumerated<'TElement>(array : 'TElement[]) : (Int, 'TElement)[] {
     MappedByIndex((index, element) -> (index, element), array)
 }
 
-    /// # Summary
 /// Returns an array containing the elements of another array,
 /// excluding elements at a given list of indices.
 ///
@@ -355,7 +346,6 @@ function Excluding<'T>(remove : Int[], array : 'T[]) : 'T[] {
     output
 }
 
-    /// # Summary
 /// Given an array and a predicate that is defined
 /// for the elements of the array, returns an array that consists of
 /// those elements that satisfy the predicate.
@@ -388,7 +378,6 @@ function Filtered<'T>(predicate : ('T -> Bool), array : 'T[]) : 'T[] {
     filtered
 }
 
-    /// # Summary
 /// Given an array and a function that maps an array element to some output
 /// array, returns the concatenated output arrays for each array element.
 ///
@@ -422,7 +411,6 @@ function FlatMapped<'TInput, 'TOutput>(mapper : ('TInput -> 'TOutput[]), array :
     output
 }
 
-    /// # Summary
 /// Given an array of arrays, returns the concatenation of all arrays.
 ///
 /// # Type Parameters
@@ -449,7 +437,6 @@ function Flattened<'T>(arrays : 'T[][]) : 'T[] {
     output
 }
 
-    /// # Summary
 /// Iterates a function `f` through an array `array`, returning
 /// `f(...f(f(initialState, array[0]), array[1]), ...)`.
 ///
@@ -484,7 +471,6 @@ function Fold<'State, 'T>(folder : (('State, 'T) -> 'State), state : 'State, arr
     current
 }
 
-    /// # Summary
 /// Given an array and an operation that is defined
 /// for the elements of the array, returns a new array that consists
 /// of the images of the original array under the operation.
@@ -514,7 +500,6 @@ operation ForEach<'T, 'U>(action : ('T => 'U), array : 'T[]) : 'U[] {
     output
 }
 
-    /// # Summary
 /// Returns the first element of the array.
 ///
 /// # Type Parameters
@@ -532,7 +517,6 @@ function Head<'A>(array : 'A[]) : 'A {
     array[0]
 }
 
-    /// # Summary
 /// Returns a tuple of first and all remaining elements of the array.
 ///
 /// # Type Parameters
@@ -549,7 +533,6 @@ function HeadAndRest<'A>(array : 'A[]) : ('A, 'A[]) {
     (Head(array), Rest(array))
 }
 
-    /// # Summary
 /// Returns the first index of the first element in an array that satisfies
 /// a given predicate. If no such element exists, returns -1.
 ///
@@ -578,7 +561,6 @@ function IndexOf<'T>(predicate : ('T -> Bool), array : 'T[]) : Int {
         -1
 }
 
-    /// # Summary
 /// Given an array, returns a range over the indices of that array, suitable
 /// for use in a for loop.
 ///
@@ -603,7 +585,6 @@ function IndexRange<'TElement>(array : 'TElement[]) : Range {
     0..Length(array) - 1
 }
 
-    /// # Summary
 /// Interleaves two arrays of (almost) same size.
 ///
 /// # Description
@@ -651,7 +632,6 @@ function Interleaved<'T>(first : 'T[], second : 'T[]) : 'T[] {
     interleaved
 }
 
-    /// # Summary
 /// Returns true if and only if an array is empty.
 ///
 /// # Input
@@ -664,7 +644,6 @@ function IsEmpty<'T>(array : 'T[]) : Bool {
     Length(array) == 0
 }
 
-    /// # Summary
 /// Returns whether a 2-dimensional array has a rectangular shape
 ///
 /// # Type Parameters
@@ -700,7 +679,6 @@ function IsRectangularArray<'T>(array : 'T[][]) : Bool {
     true
 }
 
-    /// # Summary
 /// Given an array, returns whether that array is sorted as defined by
 /// a given comparison function.
 ///
@@ -732,7 +710,6 @@ function IsSorted<'T>(comparison : (('T, 'T) -> Bool), array : 'T[]) : Bool {
     true
 }
 
-    /// # Summary
 /// Returns whether a 2-dimensional array has a square shape
 ///
 /// # Type Parameters
@@ -768,7 +745,6 @@ function IsSquareArray<'T>(array : 'T[][]) : Bool {
     true
 }
 
-    /// # Summary
 /// Given an array and a function that is defined
 /// for the elements of the array, returns a new array that consists
 /// of the images of the original array under the function.
@@ -798,7 +774,6 @@ function Mapped<'T, 'U>(mapper : ('T -> 'U), array : 'T[]) : 'U[] {
     mapped
 }
 
-    /// # Summary
 /// Given an array and a function that is defined
 /// for the indexed elements of the array, returns a new array that consists
 /// of the images of the original array under the function.
@@ -839,7 +814,6 @@ function MappedByIndex<'T, 'U>(mapper : ((Int, 'T) -> 'U), array : 'T[]) : 'U[]
     mapped
 }
 
-    /// # Summary
 /// Given a range and a function that takes an integer as input,
 /// returns a new array that consists
 /// of the images of the range values under the function.
@@ -874,7 +848,6 @@ function MappedOverRange<'T>(mapper : (Int -> 'T), range : Range) : 'T[] {
     output
 }
 
-    /// # Summary
 /// Creates an array that is equal to an input array except that the last array
 /// element is dropped.
 ///
@@ -892,7 +865,6 @@ function Most<'T>(array : 'T[]) : 'T[] {
     array[...Length(array) - 2]
 }
 
-    /// # Summary
 /// Returns a tuple of all but one and the last element of the array.
 ///
 /// # Type Parameters
@@ -909,7 +881,6 @@ function MostAndTail<'A>(array : 'A[]) : ('A[], 'A) {
     (Most(array), Tail(array))
 }
 
-    /// # Summary
 /// Returns an array padded at with specified values up to a
 /// specified length.
 ///
@@ -954,7 +925,6 @@ function Padded<'T>(paddedLength : Int, defaultElement : 'T, inputArray : 'T[])
     }
 }
 
-    /// # Summary
 /// Splits an array into multiple parts.
 ///
 /// # Input
@@ -994,7 +964,6 @@ function Partitioned<'T>(partitionSizes : Int[], array : 'T[]) : 'T[][] {
     output
 }
 
-    /// # Summary
 /// Creates an array that is equal to an input array except that the first array
 /// element is dropped.
 ///
@@ -1012,7 +981,6 @@ function Rest<'T>(array : 'T[]) : 'T[] {
     array[1...]
 }
 
-    /// # Summary
 /// Create an array that contains the same elements as an input array but in reversed
 /// order.
 ///
@@ -1030,7 +998,6 @@ function Reversed<'T>(array : 'T[]) : 'T[] {
     array[...-1...]
 }
 
-    /// # Summary
 /// Get an array of integers in a given interval.
 ///
 /// # Input
@@ -1061,7 +1028,6 @@ function SequenceI(from : Int, to : Int) : Int[] {
     array
 }
 
-    /// # Summary
 /// Get an array of integers in a given interval.
 ///
 /// # Input
@@ -1095,7 +1061,6 @@ function SequenceL(from : BigInt, to : BigInt) : BigInt[] {
     array
 }
 
-    /// # Summary
 /// Given an array, returns the elements of that array sorted by a given
 /// comparison function.
 ///
@@ -1145,7 +1110,6 @@ function Sorted<'T>(comparison : (('T, 'T) -> Bool), array : 'T[]) : 'T[] {
     )
 }
 
-    /// # Summary
 /// Given two sorted arrays, returns a single array containing the
 /// elements of both in sorted order. Used internally by `Sorted`.
 internal function SortedMerged<'T>(comparison : (('T, 'T) -> Bool), left : 'T[], right : 'T[]) : 'T[] {
@@ -1167,7 +1131,6 @@ internal function SortedMerged<'T>(comparison : (('T, 'T) -> Bool), left : 'T[],
     output + remainingLeft + remainingRight
 }
 
-    /// # Summary
 /// Takes an array and a list of locations and
 /// produces a new array formed from the elements of the original
 /// array that match the given locations.
@@ -1205,7 +1168,6 @@ function Subarray<'T>(locations : Int[], array : 'T[]) : 'T[] {
     subarray
 }
 
-    /// # Summary
 /// Applies a swap of two elements in an array.
 ///
 /// # Input
@@ -1232,7 +1194,6 @@ function Swapped<'T>(firstIndex : Int, secondIndex : Int, array : 'T[]) : 'T[] {
         w/ secondIndex <- array[firstIndex]
 }
 
-    /// # Summary
 /// Returns the transpose of a matrix represented as an array
 /// of arrays.
 ///
@@ -1276,7 +1237,6 @@ function Transposed<'T>(matrix : 'T[][]) : 'T[][] {
     transposed
 }
 
-    /// # Summary
 /// Returns the last element of the array.
 ///
 /// # Type Parameters
@@ -1295,7 +1255,6 @@ function Tail<'A>(array : 'A[]) : 'A {
     array[size - 1]
 }
 
-    /// # Summary
 /// Given an array of 2-tuples, returns a tuple of two arrays, each containing
 /// the elements of the tuples of the input array.
 ///
@@ -1332,7 +1291,6 @@ function Unzipped<'T, 'U>(array : ('T, 'U)[]) : ('T[], 'U[]) {
     return (first, second);
 }
 
-    /// # Summary
 /// Given a predicate and an array, returns the indices of that
 /// array where the predicate is true.
 ///
@@ -1358,7 +1316,6 @@ function Where<'T>(predicate : ('T -> Bool), array : 'T[]) : Int[] {
     indexes
 }
 
-    /// # Summary
 /// Returns all consecutive subarrays of length `size`.
 ///
 /// # Description
@@ -1400,7 +1357,6 @@ function Windows<'T>(size : Int, array : 'T[]) : 'T[][] {
     windows
 }
 
-    /// # Summary
 /// Given two arrays, returns a new array of pairs such that each pair
 /// contains an element from each original array.
 ///
diff --git a/library/std/src/Std/Canon.qs b/library/std/src/Std/Canon.qs
new file mode 100644
index 0000000000..be7def9b1b
--- /dev/null
+++ b/library/std/src/Std/Canon.qs
@@ -0,0 +1,588 @@
+// Copyright (c) Microsoft Corporation.
+// Licensed under the MIT License.
+
+
+import QIR.Intrinsic.*;
+open Microsoft.Quantum.Intrinsic;
+open Microsoft.Quantum.Diagnostics;
+open Microsoft.Quantum.Math;
+
+/// # Summary
+/// Applies an operation to each element in a register.
+///
+/// # Input
+/// ## singleElementOperation
+/// Operation to apply to each element.
+/// ## register
+/// Array of elements on which to apply the given operation.
+///
+/// # Type Parameters
+/// ## 'T
+/// The target on which the operation acts.
+///
+/// # Example
+/// Prepare a three-qubit |+⟩ state:
+/// ```qsharp
+/// use register = Qubit[3];
+/// ApplyToEach(H, register);
+/// ```
+operation ApplyToEach<'T>(singleElementOperation : ('T => Unit), register : 'T[]) : Unit {
+    for item in register {
+        singleElementOperation(item);
+    }
+}
+
+    /// # Summary
+/// Applies an operation to each element in a register.
+/// The modifier `A` indicates that the single-element operation is adjointable.
+///
+/// # Input
+/// ## singleElementOperation
+/// Operation to apply to each element.
+/// ## register
+/// Array of elements on which to apply the given operation.
+///
+/// # Type Parameters
+/// ## 'T
+/// The target on which the operation acts.
+///
+/// # Example
+/// Prepare a three-qubit |+⟩ state:
+/// ```qsharp
+/// use register = Qubit[3];
+/// ApplyToEach(H, register);
+/// ```
+///
+/// # See Also
+/// - Microsoft.Quantum.Canon.ApplyToEach
+operation ApplyToEachA<'T>(singleElementOperation : ('T => Unit is Adj), register : 'T[]) : Unit is Adj {
+    for item in register {
+        singleElementOperation(item);
+    }
+}
+
+    /// # Summary
+/// Applies an operation to each element in a register.
+/// The modifier `C` indicates that the single-element operation is controllable.
+///
+/// # Input
+/// ## singleElementOperation
+/// Operation to apply to each element.
+/// ## register
+/// Array of elements on which to apply the given operation.
+///
+/// # Type Parameters
+/// ## 'T
+/// The target on which the operation acts.
+///
+/// # Example
+/// Prepare a three-qubit |+⟩ state:
+/// ```qsharp
+/// use register = Qubit[3];
+/// ApplyToEach(H, register);
+/// ```
+///
+/// # See Also
+/// - Microsoft.Quantum.Canon.ApplyToEach
+operation ApplyToEachC<'T>(singleElementOperation : ('T => Unit is Ctl), register : 'T[]) : Unit is Ctl {
+    for item in register {
+        singleElementOperation(item);
+    }
+}
+
+    /// # Summary
+/// Applies an operation to each element in a register.
+/// The modifier `CA` indicates that the single-element operation is controllable and adjointable.
+///
+/// # Input
+/// ## singleElementOperation
+/// Operation to apply to each element.
+/// ## register
+/// Array of elements on which to apply the given operation.
+///
+/// # Type Parameters
+/// ## 'T
+/// The target on which the operation acts.
+///
+/// # Example
+/// Prepare a three-qubit |+⟩ state:
+/// ```qsharp
+/// use register = Qubit[3];
+/// ApplyToEach(H, register);
+/// ```
+///
+/// # See Also
+/// - Microsoft.Quantum.Canon.ApplyToEach
+operation ApplyToEachCA<'T>(singleElementOperation : ('T => Unit is Adj + Ctl), register : 'T[]) : Unit is Adj + Ctl {
+    for item in register {
+        singleElementOperation(item);
+    }
+}
+
+    /// # Summary
+/// Applies the controlled-X (CX) gate to a pair of qubits.
+///
+/// # Input
+/// ## control
+/// Control qubit for the CX gate.
+/// ## target
+/// Target qubit for the CX gate.
+///
+/// # Remarks
+/// This operation can be simulated by the unitary matrix
+/// $$
+/// \begin{align}
+///     \left(\begin{matrix}
+///         1 & 0 & 0 & 0 \\\\
+///         0 & 1 & 0 & 0 \\\\
+///         0 & 0 & 0 & 1 \\\\
+///         0 & 0 & 1 & 0
+///      \end{matrix}\right)
+/// \end{align},
+/// $$
+/// where rows and columns are organized as in the quantum concepts guide.
+///
+/// Equivalent to:
+/// ```qsharp
+/// Controlled X([control], target);
+/// ```
+/// and to:
+/// ```qsharp
+/// CNOT(control, target);
+/// ```
+operation CX(control : Qubit, target : Qubit) : Unit is Adj + Ctl {
+    body ... {
+        __quantum__qis__cx__body(control, target);
+    }
+    controlled (ctls, ...) {
+        Controlled X(ctls + [control], target);
+    }
+    adjoint self;
+}
+
+    /// # Summary
+/// Applies the controlled-Y (CY) gate to a pair of qubits.
+///
+/// # Input
+/// ## control
+/// Control qubit for the CY gate.
+/// ## target
+/// Target qubit for the CY gate.
+///
+/// # Remarks
+/// This operation can be simulated by the unitary matrix
+/// $$
+/// \begin{align}
+///     \left(\begin{matrix}
+///         1 & 0 & 0 & 0 \\\\
+///         0 & 1 & 0 & 0 \\\\
+///         0 & 0 & 0 & -i \\\\
+///         0 & 0 & i & 0
+///      \end{matrix}\right)
+/// \end{align},
+/// $$
+/// where rows and columns are organized as in the quantum concepts guide.
+///
+/// Equivalent to:
+/// ```qsharp
+/// Controlled Y([control], target);
+/// ```
+operation CY(control : Qubit, target : Qubit) : Unit is Adj + Ctl {
+    body ... {
+        __quantum__qis__cy__body(control, target);
+    }
+    controlled (ctls, ...) {
+        Controlled Y(ctls + [control], target);
+    }
+    adjoint self;
+}
+
+    /// # Summary
+/// Applies the controlled-Z (CZ) gate to a pair of qubits.
+///
+/// # Input
+/// ## control
+/// Control qubit for the CZ gate.
+/// ## target
+/// Target qubit for the CZ gate.
+///
+/// # Remarks
+/// This operation can be simulated by the unitary matrix
+/// $$
+/// \begin{align}
+///     \left(\begin{matrix}
+///         1 & 0 & 0 & 0 \\\\
+///         0 & 1 & 0 & 0 \\\\
+///         0 & 0 & 1 & 0 \\\\
+///         0 & 0 & 0 & -1
+///     \end{matrix}\right)
+/// \end{align},
+/// $$
+/// where rows and columns are organized as in the quantum concepts guide.
+///
+/// Equivalent to:
+/// ```qsharp
+/// Controlled Z([control], target);
+/// ```
+operation CZ(control : Qubit, target : Qubit) : Unit is Adj + Ctl {
+    body ... {
+        __quantum__qis__cz__body(control, target);
+    }
+    controlled (ctls, ...) {
+        Controlled Z(ctls + [control], target);
+    }
+    adjoint self;
+}
+
+    /// Given a pair, returns its first element.
+function Fst<'T, 'U>(pair : ('T, 'U)) : 'T {
+    let (fst, _) = pair;
+    return fst;
+}
+
+    /// Given a pair, returns its second element.
+function Snd<'T, 'U>(pair : ('T, 'U)) : 'U {
+    let (_, snd) = pair;
+    return snd;
+}
+
+    /// # Summary
+/// Computes the parity of a register of qubits in-place.
+///
+/// # Input
+/// ## qubits
+/// Array of qubits whose parity is to be computed and stored.
+///
+/// # Remarks
+/// This operation transforms the state of its input as
+/// $$
+/// \begin{align}
+///     \ket{q_0} \ket{q_1} \cdots \ket{q_{n - 1}} & \mapsto
+///     \ket{q_0} \ket{q_0 \oplus q_1} \ket{q_0 \oplus q_1 \oplus q_2} \cdots
+///         \ket{q_0 \oplus \cdots \oplus q_{n - 1}}.
+/// \end{align}
+/// $$
+operation ApplyCNOTChain(qubits : Qubit[]) : Unit is Adj + Ctl {
+    for i in 0..Length(qubits) - 2 {
+        CNOT(qubits[i], qubits[i + 1]);
+    }
+}
+
+    /// # Summary
+/// Given a single-qubit Pauli operator, applies the corresponding operation
+/// to a single qubit.
+///
+/// # Input
+/// ## pauli
+/// The Pauli operator to be applied.
+/// ## target
+/// The qubit to which `pauli` is to be applied as an operation.
+///
+/// # Example
+/// The following are equivalent:
+/// ```qsharp
+/// ApplyP(PauliX, q);
+/// ```
+/// and
+/// ```qsharp
+/// X(q);
+/// ```
+operation ApplyP(pauli : Pauli, target : Qubit) : Unit is Adj + Ctl {
+    if pauli == PauliX { X(target); } elif pauli == PauliY { Y(target); } elif pauli == PauliZ { Z(target); }
+}
+
+    /// # Summary
+/// Given a multi-qubit Pauli operator, applies the corresponding operation
+/// to a quantum register.
+///
+/// # Input
+/// ## pauli
+/// A multi-qubit Pauli operator represented as an array of single-qubit Pauli operators.
+/// ## target
+/// Register to apply the given Pauli operation on.
+///
+/// # Example
+/// The following are equivalent:
+/// ```qsharp
+/// ApplyPauli([PauliY, PauliZ, PauliX], target);
+/// ```
+/// and
+/// ```qsharp
+/// Y(target[0]);
+/// Z(target[1]);
+/// X(target[2]);
+/// ```
+operation ApplyPauli(pauli : Pauli[], target : Qubit[]) : Unit is Adj + Ctl {
+    Fact(Length(pauli) == Length(target), "`pauli` and `target` must be of the same length.");
+    for i in 0..Length(pauli) - 1 {
+        ApplyP(pauli[i], target[i]);
+    }
+}
+
+    /// # Summary
+/// Applies a Pauli operator on each qubit in an array if the corresponding
+/// bit of a Boolean array matches a given input.
+///
+/// # Input
+/// ## pauli
+/// Pauli operator to apply to `qubits[idx]` where `bitApply == bits[idx]`
+/// ## bitApply
+/// apply Pauli if bit is this value
+/// ## bits
+/// Boolean register specifying which corresponding qubit in `qubits` should be operated on
+/// ## qubits
+/// Quantum register on which to selectively apply the specified Pauli operator
+///
+/// # Remarks
+/// The Boolean array and the quantum register must be of equal length.
+///
+/// # Example
+/// The following applies an X operation on qubits 0 and 2, and a Z operation on qubits 1 and 3.
+/// ```qsharp
+/// use qubits = Qubit[4];
+/// let bits = [true, false, true, false];
+/// // Apply when index in `bits` is `true`.
+/// ApplyPauliFromBitString(PauliX, true, bits, qubits);
+/// // Apply when index in `bits` is `false`.
+/// ApplyPauliFromBitString(PauliZ, false, bits, qubits);
+/// ```
+operation ApplyPauliFromBitString(pauli : Pauli, bitApply : Bool, bits : Bool[], qubits : Qubit[]) : Unit is Adj + Ctl {
+    let nBits = Length(bits);
+    Fact(nBits == Length(qubits), "Number of bits must be equal to number of qubits.");
+    for i in 0..nBits - 1 {
+        if bits[i] == bitApply {
+            ApplyP(pauli, qubits[i]);
+        }
+    }
+}
+
+    /// # Summary
+/// Applies a Pauli operator on each qubit in an array if the corresponding
+/// bit of a Little-endian integer matches a given input.
+///
+/// # Input
+/// ## pauli
+/// Pauli operator to apply to `qubits[idx]` when bit of numberState
+/// in idx position is the same as bitApply.
+/// ## bitApply
+/// apply Pauli if bit is this value
+/// ## numberState
+/// Little-endian integer specifying which corresponding qubit in `qubits` should be operated on
+/// ## qubits
+/// Quantum register on which to selectively apply the specified Pauli operator
+///
+/// # Example
+/// The following applies an X operation on qubits 0 and 2, and a Z operation on qubits 1 and 3.
+/// ```qsharp
+/// use qubits = Qubit[4];
+/// let n = 5;
+/// // Apply when index in `bits` is `true`.
+/// ApplyPauliFromBitString(PauliX, true, n, qubits);
+/// // Apply when index in `bits` is `false`.
+/// ApplyPauliFromBitString(PauliZ, false, n, qubits);
+/// ```
+operation ApplyPauliFromInt(
+    pauli : Pauli,
+    bitApply : Bool,
+    numberState : Int,
+    qubits : Qubit[]
+) : Unit is Adj + Ctl {
+
+    let length = Length(qubits);
+    Fact(numberState >= 0, "number must be non-negative");
+    Fact(BitSizeI(numberState) <= length, "Bit size of numberState must not exceed qubits length");
+
+    for i in 0..length - 1 {
+        // If we assume loop unrolling, 2^i will be optimized to a constant.
+        if ((numberState &&& (1 <<< i)) != 0) == bitApply {
+            ApplyP(pauli, qubits[i]);
+        }
+    }
+}
+
+    /// # Summary
+/// Applies a unitary operation on the target if the control
+/// register state corresponds to a specified nonnegative integer.
+///
+/// # Input
+/// ## numberState
+/// A nonnegative integer on which the operation `oracle` should be
+/// controlled.
+/// ## oracle
+/// A unitary operation to be controlled.
+/// ## target
+/// A target on which to apply `oracle`.
+/// ## controlRegister
+/// A quantum register that controls application of `oracle`.
+///
+/// # Remarks
+/// The value of `numberState` is interpreted using a little-endian encoding.
+///
+/// `numberState` must be at most $2^\texttt{Length(controlRegister)} - 1$.
+/// For example, `numberState = 537` means that `oracle`
+/// is applied if and only if `controlRegister` is in the state $\ket{537}$.
+operation ApplyControlledOnInt<'T>(
+    numberState : Int,
+    oracle : ('T => Unit is Adj + Ctl),
+    controlRegister : Qubit[],
+    target : 'T
+) : Unit is Adj + Ctl {
+
+    within {
+        ApplyPauliFromInt(PauliX, false, numberState, controlRegister);
+    } apply {
+        Controlled oracle(controlRegister, target);
+    }
+}
+
+    /// # Summary
+/// Applies `oracle` on `target` when `controlRegister`
+/// is in the state specified by `bits`.
+///
+/// # Description
+/// Applies a unitary operation `oracle` on the `target`, controlled
+/// on a state specified by a given bit mask `bits`.
+/// The bit at `bits[i]` corresponds to qubit at `controlRegister[i]`.
+/// The pattern given by `bits` may be shorter than `controlRegister`,
+/// in which case additional control qubits are ignored (that is, neither
+/// controlled on |0⟩ nor |1⟩).
+/// If `bits` is longer than `controlRegister`, an error is raised.
+///
+/// # Input
+/// ## bits
+/// The bit string to control the given unitary operation on.
+/// ## oracle
+/// The unitary operation to be applied on the target.
+/// ## target
+/// The target to be passed to `oracle` as an input.
+/// ## controlRegister
+/// A quantum register that controls application of `oracle`.
+///
+/// # Example
+/// ```qsharp
+/// // When bits = [1,0,0] oracle is applied if and only if controlRegister
+/// // is in the state |100⟩.
+/// use t = Qubit();
+/// use c = Qubit[3];
+/// X(c[0]);
+/// ApplyControlledOnBitString([true, false, false], X, c, t);
+/// Message($"{M(t)}"); // Prints `One` since oracle `X` was applied.
+/// ```
+operation ApplyControlledOnBitString<'T>(
+    bits : Bool[],
+    oracle : ('T => Unit is Adj + Ctl),
+    controlRegister : Qubit[],
+    target : 'T
+) : Unit is Adj + Ctl {
+
+    // The control register must have enough bits to implement the requested control.
+    Fact(Length(bits) <= Length(controlRegister), "Control register shorter than control pattern.");
+
+    // Use a subregister of the controlled register when
+    // bits is shorter than controlRegister.
+    let controlSubregister = controlRegister[...Length(bits) - 1];
+    within {
+        ApplyPauliFromBitString(PauliX, false, bits, controlSubregister);
+    } apply {
+        Controlled oracle(controlSubregister, target);
+    }
+}
+
+    /// # Summary
+/// Applies the rotations of Quantum Fourier Transform (QFT) to a little-endian quantum register.
+///
+/// # Description
+/// Applies the rotations of QFT to a little-endian register `qs` of length n
+/// containing |x₁⟩⊗|x₂⟩⊗…⊗|xₙ⟩. The qs[0] initially contains the
+/// least significant bit xₙ. The state of qs[0] becomes
+/// (|0⟩+𝑒^(2π𝑖[0.xₙ])|1⟩)/sqrt(2) after the operation.
+///
+/// # Input
+/// ## qs
+/// Quantum register in a little-endian format to which the rotations are applied.
+///
+/// # Remarks
+/// Note that this operation applies only the rotations part of the QFT.
+/// To complete the transform, you need to reverse the order of qubits after this operation,
+/// for example, using the operation `SwapReverseRegister`.
+///
+/// # Reference
+///  - [Quantum Fourier transform](https://en.wikipedia.org/wiki/Quantum_Fourier_transform)
+operation ApplyQFT(qs : Qubit[]) : Unit is Adj + Ctl {
+    let length = Length(qs);
+    Fact(length >= 1, "ApplyQFT: Length(qs) must be at least 1.");
+    for i in length - 1..-1..0 {
+        H(qs[i]);
+        for j in 0..i - 1 {
+            Controlled R1Frac([qs[i]], (1, j + 1, qs[i - j - 1]));
+        }
+    }
+}
+
+    /// # Summary
+/// Uses SWAP gates to reverse the order of the qubits in a register.
+///
+/// # Input
+/// ## register
+/// The qubits order of which should be reversed using SWAP gates
+operation SwapReverseRegister(register : Qubit[]) : Unit is Adj + Ctl {
+    let length = Length(register);
+    for i in 0..length / 2 - 1 {
+        SWAP(register[i], register[(length - i) - 1]);
+    }
+}
+
+    /// # Summary
+/// Applies a bitwise-XOR operation between a classical integer and an
+/// integer represented by a register of qubits.
+///
+/// # Description
+/// Applies `X` operations to qubits in a little-endian register based on
+/// 1 bits in an integer.
+///
+/// Let us denote `value` by a and let y be an unsigned integer encoded in `target`,
+/// then `ApplyXorInPlace` performs an operation given by the following map:
+/// |y⟩ ↦ |y ⊕ a⟩, where ⊕ is the bitwise exclusive OR operator.
+operation ApplyXorInPlace(value : Int, target : Qubit[]) : Unit is Adj + Ctl {
+    body (...) {
+        Fact(value >= 0, "`value` must be non-negative.");
+        mutable runningValue = value;
+        for q in target {
+            if (runningValue &&& 1) != 0 {
+                X(q);
+            }
+            set runningValue >>>= 1;
+        }
+        Fact(runningValue == 0, "value is too large");
+    }
+    adjoint self;
+}
+
+    /// # Summary
+/// Applies a bitwise-XOR operation between a classical integer and an
+/// integer represented by a register of qubits.
+///
+/// # Description
+/// Applies `X` operations to qubits in a little-endian register based on
+/// 1 bits in an integer.
+///
+/// Let us denote `value` by a and let y be an unsigned integer encoded in `target`,
+/// then `ApplyXorInPlace` performs an operation given by the following map:
+/// |y⟩ ↦ |y ⊕ a⟩, where ⊕ is the bitwise exclusive OR operator.
+operation ApplyXorInPlaceL(value : BigInt, target : Qubit[]) : Unit is Adj + Ctl {
+    body (...) {
+        Fact(value >= 0L, "`value` must be non-negative.");
+        mutable runningValue = value;
+        for q in target {
+            if (runningValue &&& 1L) != 0L {
+                X(q);
+            }
+            set runningValue >>>= 1;
+        }
+        Fact(runningValue == 0L, "`value` is too large.");
+    }
+    adjoint self;
+}
+
+export ApplyToEach, ApplyToEachA, ApplyToEachC, ApplyToEachCA, CX, CY, CZ, Fst, Snd, ApplyCNOTChain, ApplyP, ApplyPauli, ApplyPauliFromBitString, ApplyPauliFromInt, ApplyControlledOnInt, ApplyControlledOnBitString, ApplyQFT, SwapReverseRegister, ApplyXorInPlace, ApplyXorInPlaceL;
+
diff --git a/library/std/src/Std/Convert.qs b/library/std/src/Std/Convert.qs
index 19a0ed6c7a..46ab3a3715 100644
--- a/library/std/src/Std/Convert.qs
+++ b/library/std/src/Std/Convert.qs
@@ -21,13 +21,11 @@ function IntAsDouble(number : Int) : Double {
     body intrinsic;
 }
 
-    /// # Summary
 /// Converts a given integer `number` to an equivalent big integer.
 function IntAsBigInt(number : Int) : BigInt {
     body intrinsic;
 }
 
-    /// # Summary
 /// Converts a `Result` type to a `Bool` type, where `One` is mapped to
 /// `true` and `Zero` is mapped to `false`.
 ///
@@ -41,7 +39,6 @@ function ResultAsBool(input : Result) : Bool {
     input == One
 }
 
-    /// # Summary
 /// Converts a `Bool` type to a `Result` type, where `true` is mapped to
 /// `One` and `false` is mapped to `Zero`.
 ///
@@ -55,7 +52,6 @@ function BoolAsResult(input : Bool) : Result {
     if input { One } else { Zero }
 }
 
-    /// # Summary
 /// Produces a non-negative integer from a string of bits in little-endian format.
 /// `bits[0]` represents the least significant bit.
 ///
@@ -76,7 +72,6 @@ function BoolArrayAsInt(bits : Bool[]) : Int {
     number
 }
 
-    /// # Summary
 /// Produces a binary representation of a non-negative integer, using the
 /// little-endian representation for the returned array.
 ///
@@ -106,7 +101,6 @@ function IntAsBoolArray(number : Int, bits : Int) : Bool[] {
     result
 }
 
-    /// # Summary
 /// Converts an array of Boolean values into a non-negative BigInt, interpreting the
 /// array as a binary representation in little-endian format.
 ///
@@ -132,7 +126,6 @@ function BoolArrayAsBigInt(boolArray : Bool[]) : BigInt {
     result
 }
 
-    /// # Summary
 /// Produces a binary representation of a non-negative BigInt, using the
 /// little-endian representation for the returned array.
 ///
@@ -162,7 +155,6 @@ function BigIntAsBoolArray(number : BigInt, bits : Int) : Bool[] {
     result
 }
 
-    /// # Summary
 /// Produces a non-negative integer from a string of Results in little-endian format.
 ///
 /// # Input
@@ -191,7 +183,6 @@ function ResultArrayAsInt(results : Result[]) : Int {
     number
 }
 
-    /// # Summary
 /// Converts a `Result[]` type to a `Bool[]` type, where `One`
 /// is mapped to `true` and `Zero` is mapped to `false`.
 ///
@@ -210,7 +201,6 @@ function ResultArrayAsBoolArray(input : Result[]) : Bool[] {
     output
 }
 
-    /// # Summary
 /// Converts a `Bool[]` type to a `Result[]` type, where `true`
 /// is mapped to `One` and `false` is mapped to `Zero`.
 ///
@@ -229,7 +219,6 @@ function BoolArrayAsResultArray(input : Bool[]) : Result[] {
     output
 }
 
-    /// # Summary
 /// Converts a complex number of type `Complex` to a complex
 /// number of type `ComplexPolar`.
 ///
@@ -243,7 +232,6 @@ function ComplexAsComplexPolar(input : Complex) : ComplexPolar {
     return ComplexPolar(AbsComplex(input), ArgComplex(input));
 }
 
-    /// # Summary
 /// Converts a complex number of type `ComplexPolar` to a complex
 /// number of type `Complex`.
 ///
@@ -260,7 +248,6 @@ function ComplexPolarAsComplex(input : ComplexPolar) : Complex {
     );
 }
 
-    /// # Summary
 /// Converts a given double-precision floating-point number to a string representation with desired precision, rounding if required.
 ///
 /// # Input
diff --git a/library/std/src/Std/Diagnostics.qs b/library/std/src/Std/Diagnostics.qs
index c250b29b48..127d05b4f2 100644
--- a/library/std/src/Std/Diagnostics.qs
+++ b/library/std/src/Std/Diagnostics.qs
@@ -29,7 +29,6 @@ function DumpMachine() : Unit {
     body intrinsic;
 }
 
-    /// # Summary
 /// Dumps the current target machine's status associated with the given qubits.
 ///
 /// # Input
@@ -65,7 +64,6 @@ function DumpRegister(register : Qubit[]) : Unit {
     body intrinsic;
 }
 
-    /// # Summary
 /// Checks whether a qubit is in the |0⟩ state, returning true if it is.
 ///
 /// # Description
@@ -87,7 +85,6 @@ operation CheckZero(qubit : Qubit) : Bool {
     body intrinsic;
 }
 
-    /// # Summary
 /// Checks whether all qubits in the provided array are in the |0⟩ state. Returns true if they are.
 ///
 /// # Description
@@ -115,7 +112,6 @@ operation CheckAllZero(qubits : Qubit[]) : Bool {
     return true;
 }
 
-    /// # Summary
 /// Checks whether a given condition is true, failing with a message if it is not.
 ///
 /// # Description
@@ -133,7 +129,6 @@ function Fact(actual : Bool, message : String) : Unit {
     }
 }
 
-    /// # Summary
 /// Given two operations, checks that they act identically for all input states.
 ///
 /// # Description
@@ -190,7 +185,6 @@ operation CheckOperationsAreEqual(
     areEqual
 }
 
-    /// # Summary
 /// Starts counting the number of times the given operation is called. Fails if the operation is already being counted.
 ///
 /// # Description
@@ -227,7 +221,6 @@ operation StartCountingOperation<'In, 'Out>(callable : 'In => 'Out) : Unit {
     body intrinsic;
 }
 
-    /// # Summary
 /// Stops counting the number of times the given operation is called and returns the count. Fails
 /// if the operation was not being counted.
 ///
@@ -245,7 +238,6 @@ operation StopCountingOperation<'In, 'Out>(callable : 'In => 'Out) : Int {
     body intrinsic;
 }
 
-    /// # Summary
 /// Starts counting the number of times the given function is called. Fails if the function is already being counted.
 ///
 /// # Description
@@ -273,7 +265,6 @@ operation StartCountingFunction<'In, 'Out>(callable : 'In -> 'Out) : Unit {
     body intrinsic;
 }
 
-    /// # Summary
 /// Stops counting the number of times the given function is called and returns the count. Fails
 /// if the function was not being counted.
 ///
@@ -291,7 +282,6 @@ operation StopCountingFunction<'In, 'Out>(callable : 'In -> 'Out) : Int {
     body intrinsic;
 }
 
-    /// # Summary
 /// Starts counting the number of qubits allocated. Fails if qubits are already being counted.
 ///
 /// # Description
@@ -314,7 +304,6 @@ operation StartCountingQubits() : Unit {
     body intrinsic;
 }
 
-    /// # Summary
 /// Stops counting the number of qubits allocated and returns the count. Fails if the qubits were not being counted.
 ///
 /// # Description
diff --git a/library/std/src/Std/Intrinsic.qs b/library/std/src/Std/Intrinsic.qs
index f7c8f01f7a..2e177663c9 100644
--- a/library/std/src/Std/Intrinsic.qs
+++ b/library/std/src/Std/Intrinsic.qs
@@ -36,7 +36,6 @@ operation AND(control1 : Qubit, control2 : Qubit, target : Qubit) : Unit is Adj
     }
 }
 
-    /// # Summary
 /// Applies the AND gate that is more efficient for use with decomposition of multi-controlled operations.
 /// Note that target qubit must be in |0⟩ state.
 ///
@@ -55,7 +54,6 @@ operation AND(control1 : Qubit, control2 : Qubit, target : Qubit) : Unit is Adj
     PhaseCCX(control1, control2, target);
 }
 
-    /// # Summary
 /// Applies the doubly controlled–NOT (CCNOT) gate to three qubits.
 ///
 /// # Input
@@ -81,7 +79,6 @@ operation CCNOT(control1 : Qubit, control2 : Qubit, target : Qubit) : Unit is Ad
     adjoint self;
 }
 
-    /// # Summary
 /// Applies the controlled-NOT (CNOT) gate to a pair of qubits.
 ///
 /// # Input
@@ -119,7 +116,6 @@ operation CNOT(control : Qubit, target : Qubit) : Unit is Adj + Ctl {
     adjoint self;
 }
 
-    /// # Summary
 /// Applies the exponential of a multi-qubit Pauli operator.
 ///
 /// # Input
@@ -183,7 +179,6 @@ operation Exp(paulis : Pauli[], theta : Double, qubits : Qubit[]) : Unit is Adj
     }
 }
 
-    /// # Summary
 /// Applies the Hadamard transformation to a single qubit.
 ///
 /// # Input
@@ -228,7 +223,6 @@ operation H(qubit : Qubit) : Unit is Adj + Ctl {
     adjoint self;
 }
 
-    /// # Summary
 /// Performs the identity operation (no-op) on a single qubit.
 ///
 /// # Remarks
@@ -239,7 +233,6 @@ operation I(target : Qubit) : Unit is Adj + Ctl {
     adjoint self;
 }
 
-    /// # Summary
 /// Performs a measurement of a single qubit in the
 /// Pauli _Z_ basis.
 ///
@@ -270,7 +263,6 @@ operation M(qubit : Qubit) : Result {
     __quantum__qis__m__body(qubit)
 }
 
-    /// # Summary
 /// Performs a measurement of a single qubit in the
 /// Pauli _Z_ basis.
 ///
@@ -301,7 +293,6 @@ operation M(qubit : Qubit) : Result {
     Measure([PauliZ], [qubit])
 }
 
-    /// # Summary
 /// Performs a joint measurement of one or more qubits in the
 /// specified Pauli bases.
 ///
@@ -350,7 +341,6 @@ operation Measure(bases : Pauli[], qubits : Qubit[]) : Result {
     }
 }
 
-    /// # Summary
 /// Performs a joint measurement of one or more qubits in the
 /// specified Pauli bases.
 ///
@@ -394,7 +384,6 @@ operation Measure(bases : Pauli[], qubits : Qubit[]) : Result {
     __quantum__qis__mresetz__body(aux)
 }
 
-    /// # Summary
 /// Applies a rotation about the given Pauli axis.
 ///
 /// # Input
@@ -430,7 +419,6 @@ operation R(pauli : Pauli, theta : Double, qubit : Qubit) : Unit is Adj + Ctl {
     }
 }
 
-    /// # Summary
 /// Applies a rotation about the |1⟩ state by a given angle.
 ///
 /// # Input
@@ -457,7 +445,6 @@ operation R1(theta : Double, qubit : Qubit) : Unit is Adj + Ctl {
     R(PauliI, -theta, qubit);
 }
 
-    /// # Summary
 /// Applies a rotation about the |1⟩ state by an angle specified
 /// as a dyadic fraction.
 ///
@@ -493,7 +480,6 @@ operation R1Frac(numerator : Int, power : Int, qubit : Qubit) : Unit is Adj + Ct
     RFrac(PauliI, numerator, power + 1, qubit);
 }
 
-    /// # Summary
 /// Given a single qubit, measures it and ensures it is in the |0⟩ state
 /// such that it can be safely released.
 ///
@@ -504,7 +490,6 @@ operation Reset(qubit : Qubit) : Unit {
     __quantum__qis__reset__body(qubit);
 }
 
-    /// # Summary
 /// Given an array of qubits, measure them and ensure they are in the |0⟩ state
 /// such that they can be safely released.
 ///
@@ -517,7 +502,6 @@ operation ResetAll(qubits : Qubit[]) : Unit {
     }
 }
 
-    /// # Summary
 /// Applies a rotation about the given Pauli axis by an angle specified
 /// as a dyadic fraction.
 ///
@@ -558,7 +542,6 @@ operation RFrac(pauli : Pauli, numerator : Int, power : Int, qubit : Qubit) : Un
     R(pauli, angle, qubit);
 }
 
-    /// # Summary
 /// Applies a rotation about the _x_-axis by a given angle.
 ///
 /// # Input
@@ -603,7 +586,6 @@ operation Rx(theta : Double, qubit : Qubit) : Unit is Adj + Ctl {
     }
 }
 
-    /// # Summary
 /// Applies the two qubit Ising _XX_ rotation gate.
 ///
 /// # Input
@@ -650,7 +632,6 @@ operation Rxx(theta : Double, qubit0 : Qubit, qubit1 : Qubit) : Unit is Adj + Ct
     }
 }
 
-    /// # Summary
 /// Applies a rotation about the _y_-axis by a given angle.
 ///
 /// # Input
@@ -695,7 +676,6 @@ operation Ry(theta : Double, qubit : Qubit) : Unit is Adj + Ctl {
     }
 }
 
-    /// # Summary
 /// Applies the two qubit Ising _YY_ rotation gate.
 ///
 /// # Input
@@ -742,7 +722,6 @@ operation Ryy(theta : Double, qubit0 : Qubit, qubit1 : Qubit) : Unit is Adj + Ct
     }
 }
 
-    /// # Summary
 /// Applies a rotation about the _z_-axis by a given angle.
 ///
 /// # Input
@@ -791,7 +770,6 @@ operation Rz(theta : Double, qubit : Qubit) : Unit is Adj + Ctl {
     }
 }
 
-    /// # Summary
 /// Applies the two qubit Ising _ZZ_ rotation gate.
 ///
 /// # Input
@@ -838,7 +816,6 @@ operation Rzz(theta : Double, qubit0 : Qubit, qubit1 : Qubit) : Unit is Adj + Ct
     }
 }
 
-    /// # Summary
 /// Applies the π/4 phase gate to a single qubit.
 ///
 /// # Input
@@ -904,7 +881,6 @@ operation S(qubit : Qubit) : Unit is Adj + Ctl {
     }
 }
 
-    /// # Summary
 /// Applies the SWAP gate to a pair of qubits.
 ///
 /// # Input
@@ -952,7 +928,6 @@ operation SWAP(qubit1 : Qubit, qubit2 : Qubit) : Unit is Adj + Ctl {
     }
 }
 
-    /// # Summary
 /// Applies the π/8 gate to a single qubit.
 ///
 /// # Input
@@ -1008,7 +983,6 @@ operation T(qubit : Qubit) : Unit is Adj + Ctl {
     }
 }
 
-    /// # Summary
 /// Applies the Pauli _X_ gate.
 ///
 /// # Input
@@ -1052,7 +1026,6 @@ operation X(qubit : Qubit) : Unit is Adj + Ctl {
     adjoint self;
 }
 
-    /// # Summary
 /// Applies the Pauli _Y_ gate.
 ///
 /// # Input
@@ -1096,7 +1069,6 @@ operation Y(qubit : Qubit) : Unit is Adj + Ctl {
     adjoint self;
 }
 
-    /// # Summary
 /// Applies the Pauli _Z_ gate.
 ///
 /// # Input
@@ -1140,7 +1112,6 @@ operation Z(qubit : Qubit) : Unit is Adj + Ctl {
     adjoint self;
 }
 
-    /// # Summary
 /// Logs a message.
 ///
 /// # Input
diff --git a/library/std/src/Std/Math.qs b/library/std/src/Std/Math.qs
index 52079e74fa..d046f91908 100644
--- a/library/std/src/Std/Math.qs
+++ b/library/std/src/Std/Math.qs
@@ -27,7 +27,6 @@ function PI() : Double {
     3.14159265358979323846
 }
 
-    /// # Summary
 /// Returns a double-precision approximation of the
 /// mathematical constant 𝒆 ≈ 2.7182818284590452354
 ///
@@ -44,7 +43,6 @@ function E() : Double {
     2.7182818284590452354
 }
 
-    /// # Summary
 /// Returns a double-precision approximation of the constant
 /// ㏑2 ≈ 0.6931471805599453
 ///
@@ -75,7 +73,6 @@ function IsNaN(d : Double) : Bool {
     return d != d;
 }
 
-    /// # Summary
 /// Returns whether a given floating-point value is either positive or
 /// negative infinity.
 ///
@@ -124,7 +121,6 @@ function SignI(a : Int) : Int {
     }
 }
 
-    /// # Summary
 /// Returns -1, 0 or +1 that indicates the sign of a number.
 function SignD(a : Double) : Int {
     if (a < 0.0) {
@@ -136,7 +132,6 @@ function SignD(a : Double) : Int {
     }
 }
 
-    /// # Summary
 /// Returns -1, 0 or +1 that indicates the sign of a number.
 function SignL(a : BigInt) : Int {
     if (a < 0L) {
@@ -148,60 +143,50 @@ function SignL(a : BigInt) : Int {
     }
 }
 
-    /// # Summary
 /// Returns the absolute value of an integer.
 function AbsI(a : Int) : Int {
     a < 0 ? -a | a
 }
 
-    /// # Summary
 /// Returns the absolute value of a double-precision floating-point number.
 function AbsD(a : Double) : Double {
     a < 0.0 ? -a | a
 }
 
-    /// # Summary
 function AbsL(a : BigInt) : BigInt {
     a < 0L ? -a | a
 }
 
-    /// # Summary
 /// Returns the larger of two specified numbers.
 function MaxI(a : Int, b : Int) : Int {
     a > b ? a | b
 }
 
-    /// # Summary
 /// Returns the larger of two specified numbers.
 function MaxD(a : Double, b : Double) : Double {
     a > b ? a | b
 }
 
-    /// # Summary
 /// Returns the larger of two specified numbers.
 function MaxL(a : BigInt, b : BigInt) : BigInt {
     a > b ? a | b
 }
 
-    /// # Summary
 /// Returns the smaller of two specified numbers.
 function MinI(a : Int, b : Int) : Int {
     a < b ? a | b
 }
 
-    /// # Summary
 /// Returns the smaller of two specified numbers.
 function MinD(a : Double, b : Double) : Double {
     a < b ? a | b
 }
 
-    /// # Summary
 /// Returns the smaller of two specified numbers.
 function MinL(a : BigInt, b : BigInt) : BigInt {
     a < b ? a | b
 }
 
-    /// # Summary
 /// Given an array of integers, returns the largest element.
 ///
 /// # Input
@@ -222,7 +207,6 @@ function Max(values : Int[]) : Int {
     max
 }
 
-    /// # Summary
 /// Given an array of integers, returns the smallest element.
 ///
 /// # Input
@@ -253,74 +237,62 @@ function ArcCos(x : Double) : Double {
     body intrinsic;
 }
 
-    /// # Summary
 /// Returns the angle whose sine is the specified number.
 function ArcSin(y : Double) : Double {
     body intrinsic;
 }
 
-    /// # Summary
 /// Returns the angle whose tangent is the specified number.
 function ArcTan(d : Double) : Double {
     body intrinsic;
 }
 
-    /// # Summary
 /// Returns the angle whose tangent is the quotient of two specified numbers.
 function ArcTan2(y : Double, x : Double) : Double {
     body intrinsic;
 }
 
-    /// # Summary
 /// Returns the cosine of the specified angle.
 function Cos(theta : Double) : Double {
     body intrinsic;
 }
 
-    /// # Summary
 /// Returns the hyperbolic cosine of the specified angle.
 function Cosh(d : Double) : Double {
     body intrinsic;
 }
 
-    /// # Summary
 /// Returns the sine of the specified angle.
 function Sin(theta : Double) : Double {
     body intrinsic;
 }
 
-    /// # Summary
 /// Returns the hyperbolic sine of the specified angle.
 function Sinh(d : Double) : Double {
     body intrinsic;
 }
 
-    /// # Summary
 /// Returns the tangent of the specified angle.
 function Tan(d : Double) : Double {
     body intrinsic;
 }
 
-    /// # Summary
 /// Returns the hyperbolic tangent of the specified angle.
 function Tanh(d : Double) : Double {
     body intrinsic;
 }
 
-    /// # Summary
 /// Computes the inverse hyperbolic cosine of a number.
 function ArcCosh(x : Double) : Double {
     Log(x + Sqrt(x * x - 1.0))
 }
 
-    /// # Summary
 /// Computes the inverse hyperbolic sine of a number.
 function ArcSinh(x : Double) : Double {
     Log(x + Sqrt(x * x + 1.0))
 }
 
 
-    /// # Summary
 /// Computes the inverse hyperbolic tangent of a number.
 function ArcTanh(x : Double) : Double {
     Log((1.0 + x) / (1.0 - x)) * 0.5
@@ -336,19 +308,16 @@ function Sqrt(d : Double) : Double {
     body intrinsic;
 }
 
-    /// # Summary
 /// Returns the natural (base _e_) logarithm of a specified number.
 function Log(input : Double) : Double {
     body intrinsic;
 }
 
-    /// # Summary
 /// Returns the base-10 logarithm of a specified number.
 function Log10(input : Double) : Double {
     Log(input) / Log(10.0)
 }
 
-    /// # Summary
 /// Computes the base-2 logarithm of a number.
 function Lg(input : Double) : Double {
     Log(input) / Log(2.0)
@@ -370,7 +339,6 @@ internal function ExtendedTruncation(value : Double) : (Int, Double, Bool) {
     (truncated, IntAsDouble(truncated) - value, value >= 0.0)
 }
 
-    /// # Summary
 /// Returns the smallest integer greater than or equal to the specified number.
 /// For example: Ceiling(3.1) = 4; Ceiling(-3.7) = -3
 function Ceiling(value : Double) : Int {
@@ -382,7 +350,6 @@ function Ceiling(value : Double) : Int {
     }
 }
 
-    /// # Summary
 /// Returns the largest integer less than or equal to the specified number.
 /// For example: Floor(3.7) = 3; Floor(-3.1) = -4
 function Floor(value : Double) : Int {
@@ -394,7 +361,6 @@ function Floor(value : Double) : Int {
     }
 }
 
-    /// # Summary
 /// Returns the nearest integer to the specified number.
 /// For example: Round(3.7) = 4; Round(-3.7) = -4
 function Round(value : Double) : Int {
@@ -417,13 +383,11 @@ function DivRemI(dividend : Int, divisor : Int) : (Int, Int) {
     (dividend / divisor, dividend % divisor)
 }
 
-    /// # Summary
 /// Divides one BigInteger value by another, returns the result and the remainder as a tuple.
 function DivRemL(dividend : BigInt, divisor : BigInt) : (BigInt, BigInt) {
     (dividend / divisor, dividend % divisor)
 }
 
-    /// # Summary
 /// Computes the canonical residue of `value` modulo `modulus`.
 /// The result is always in the range 0..modulus-1 even for negative numbers.
 function ModulusI(value : Int, modulus : Int) : Int {
@@ -432,7 +396,6 @@ function ModulusI(value : Int, modulus : Int) : Int {
     (r < 0) ? (r + modulus) | r
 }
 
-    /// # Summary
 /// Computes the canonical residue of `value` modulo `modulus`.
 /// The result is always in the range 0..modulus-1 even for negative numbers.
 function ModulusL(value : BigInt, modulus : BigInt) : BigInt {
@@ -441,7 +404,6 @@ function ModulusL(value : BigInt, modulus : BigInt) : BigInt {
     (r < 0L) ? (r + modulus) | r
 }
 
-    /// # Summary
 /// Returns an integer raised to a given power, with respect to a given
 /// modulus. I.e. (expBase^power) % modulus.
 function ExpModI(expBase : Int, power : Int, modulus : Int) : Int {
@@ -472,7 +434,6 @@ function ExpModI(expBase : Int, power : Int, modulus : Int) : Int {
     res
 }
 
-    /// # Summary
 /// Returns an integer raised to a given power, with respect to a given
 /// modulus. I.e. (expBase^power) % modulus.
 function ExpModL(expBase : BigInt, power : BigInt, modulus : BigInt) : BigInt {
@@ -503,7 +464,6 @@ function ExpModL(expBase : BigInt, power : BigInt, modulus : BigInt) : BigInt {
     res
 }
 
-    /// # Summary
 /// Returns the multiplicative inverse of a modular integer.
 ///
 /// # Description
@@ -516,7 +476,6 @@ function InverseModI(a : Int, modulus : Int) : Int {
     ModulusI(u, modulus)
 }
 
-    /// # Summary
 /// Returns the multiplicative inverse of a modular integer.
 ///
 /// # Description
@@ -547,7 +506,6 @@ function GreatestCommonDivisorI(a : Int, b : Int) : Int {
     aa
 }
 
-    /// # Summary
 /// Computes the greatest common divisor of two integers.
 /// Note: GCD is always positive except that GCD(0,0)=0.
 function GreatestCommonDivisorL(a : BigInt, b : BigInt) : BigInt {
@@ -561,7 +519,6 @@ function GreatestCommonDivisorL(a : BigInt, b : BigInt) : BigInt {
     aa
 }
 
-    /// # Summary
 /// Returns a tuple (u,v) such that u*a+v*b=GCD(a,b)
 /// Note: GCD is always positive except that GCD(0,0)=0.
 function ExtendedGreatestCommonDivisorI(a : Int, b : Int) : (Int, Int) {
@@ -581,7 +538,6 @@ function ExtendedGreatestCommonDivisorI(a : Int, b : Int) : (Int, Int) {
     (s1 * signA, t1 * signB)
 }
 
-    /// # Summary
 /// Returns a tuple (u,v) such that u*a+v*b=GCD(a,b)
 /// Note: GCD is always positive except that GCD(0,0)=0.
 function ExtendedGreatestCommonDivisorL(a : BigInt, b : BigInt) : (BigInt, BigInt) {
@@ -601,7 +557,6 @@ function ExtendedGreatestCommonDivisorL(a : BigInt, b : BigInt) : (BigInt, BigIn
     (s1 * signA, t1 * signB)
 }
 
-    /// # Summary
 /// Returns if two integers are co-prime.
 ///
 /// # Description
@@ -620,7 +575,6 @@ function IsCoprimeI(a : Int, b : Int) : Bool {
     GreatestCommonDivisorI(a, b) == 1
 }
 
-    /// # Summary
 /// Returns if two integers are co-prime.
 ///
 /// # Description
@@ -639,7 +593,6 @@ function IsCoprimeL(a : BigInt, b : BigInt) : Bool {
     GreatestCommonDivisorL(a, b) == 1L
 }
 
-    /// # Summary
 /// Finds the continued fraction convergent closest to `fraction`
 /// with the denominator less or equal to `denominatorBound`
 /// Using process similar to this: https://nrich.maths.org/1397
@@ -670,7 +623,6 @@ function ContinuedFractionConvergentI(
     }
 }
 
-    /// # Summary
 /// Finds the continued fraction convergent closest to `fraction`
 /// with the denominator less or equal to `denominatorBound`
 /// Using process similar to this: https://nrich.maths.org/1397
@@ -701,7 +653,6 @@ function ContinuedFractionConvergentL(
     }
 }
 
-    /// # Summary
 /// Computes the modulus between two real numbers.
 ///
 /// # Input
@@ -745,7 +696,6 @@ function BitSizeI(a : Int) : Int {
     size
 }
 
-    /// # Summary
 /// For a non-negative integer `a`, returns the number of bits required to represent `a`.
 /// NOTE: This function returns the smallest n such that a < 2^n.
 function BitSizeL(a : BigInt) : Int {
@@ -760,7 +710,6 @@ function BitSizeL(a : BigInt) : Int {
     size
 }
 
-    /// # Summary
 /// For a non-zero integer `a`, returns the number of trailing zero bits
 /// in the binary representation of `a`.
 function TrailingZeroCountI(a : Int) : Int {
@@ -776,7 +725,6 @@ function TrailingZeroCountI(a : Int) : Int {
     count
 }
 
-    /// # Summary
 /// For a non-zero integer `a`, returns the number of trailing zero bits
 /// in the binary representation of `a`.
 function TrailingZeroCountL(a : BigInt) : Int {
@@ -792,7 +740,6 @@ function TrailingZeroCountL(a : BigInt) : Int {
     count
 }
 
-    /// # Summary
 /// Returns the number of 1 bits in the binary representation of integer `n`.
 function HammingWeightI(n : Int) : Int {
     let i1 = n - ((n >>> 1) &&& 0x5555555555555555);
@@ -853,7 +800,6 @@ function FactorialI(n : Int) : Int {
     ][n]
 }
 
-    /// # Summary
 /// Returns the factorial of a given number.
 ///
 /// # Input
@@ -876,7 +822,6 @@ function FactorialL(n : Int) : BigInt {
     result
 }
 
-    /// # Summary
 /// Returns an approximate factorial of a given number.
 ///
 /// # Description
@@ -914,7 +859,6 @@ function ApproximateFactorial(n : Int) : Double {
     a * b * c
 }
 
-    /// # Summary
 /// Returns the natural logarithm of the gamma function (aka the log-gamma
 /// function).
 ///
@@ -966,7 +910,6 @@ function LogGammaD(x : Double) : Double {
     Log(2.506628274631000 * acc / x) + ((x + 0.5) * Log(tmp) - tmp)
 }
 
-    /// # Summary
 /// Returns the approximate natural logarithm of the factorial of a given
 /// integer.
 ///
@@ -985,7 +928,6 @@ function LogFactorialD(n : Int) : Double {
     LogGammaD(IntAsDouble(n) + 1.0)
 }
 
-    /// # Summary
 /// Returns the approximate binomial coefficient of two integers.
 ///
 /// # Description
@@ -1035,7 +977,6 @@ function SquaredNorm(array : Double[]) : Double {
     sum
 }
 
-    /// # Summary
 /// Returns the `L(p)` norm of a vector of `Double`s.
 ///
 /// That is, given an array x of type `Double[]`, this returns the p-norm
@@ -1060,7 +1001,6 @@ function PNorm(p : Double, array : Double[]) : Double {
     sum^(1.0 / p)
 }
 
-    /// # Summary
 /// Normalizes a vector of `Double`s in the `L(p)` norm.
 ///
 /// That is, given an array x of type `Double[]`, this returns an array where
@@ -1106,7 +1046,6 @@ function PNormalized(p : Double, array : Double[]) : Double[] {
 /// ```
 struct Complex { Real : Double, Imag : Double }
 
-    /// # Summary
 /// Represents a complex number in polar form.
 /// The polar representation of a complex number is c = r⋅𝑒^(t𝑖).
 ///
@@ -1117,7 +1056,6 @@ struct Complex { Real : Double, Imag : Double }
 /// The phase t ∈ ℝ of c.
 struct ComplexPolar { Magnitude : Double, Argument : Double }
 
-    /// # Summary
 /// Returns the squared absolute value of a complex number of type
 /// `Complex`.
 ///
@@ -1131,7 +1069,6 @@ function AbsSquaredComplex(input : Complex) : Double {
     input.Real * input.Real + input.Imag * input.Imag
 }
 
-    /// # Summary
 /// Returns the absolute value of a complex number of type
 /// `Complex`.
 ///
@@ -1145,7 +1082,6 @@ function AbsComplex(input : Complex) : Double {
     Sqrt(AbsSquaredComplex(input))
 }
 
-    /// # Summary
 /// Returns the phase of a complex number of type
 /// `Complex`.
 ///
@@ -1159,7 +1095,6 @@ function ArgComplex(input : Complex) : Double {
     ArcTan2(input.Imag, input.Real)
 }
 
-    /// # Summary
 /// Returns the squared absolute value of a complex number of type
 /// `ComplexPolar`.
 ///
@@ -1173,7 +1108,6 @@ function AbsSquaredComplexPolar(input : ComplexPolar) : Double {
     input.Magnitude * input.Magnitude
 }
 
-    /// # Summary
 /// Returns the absolute value of a complex number of type
 /// `ComplexPolar`.
 ///
@@ -1185,7 +1119,6 @@ function AbsSquaredComplexPolar(input : ComplexPolar) : Double {
 /// Absolute value |c| = r.
 function AbsComplexPolar(input : ComplexPolar) : Double { input.Magnitude }
 
-    /// # Summary
 /// Returns the phase of a complex number of type `ComplexPolar`.
 ///
 /// # Input
@@ -1196,7 +1129,6 @@ function AbsComplexPolar(input : ComplexPolar) : Double { input.Magnitude }
 /// Phase Arg(c) = t.
 function ArgComplexPolar(input : ComplexPolar) : Double { input.Argument }
 
-    /// # Summary
 /// Returns the unary negation of an input of type `Complex`.
 ///
 /// # Input
@@ -1209,7 +1141,6 @@ function NegationC(input : Complex) : Complex {
     Complex(-input.Real, -input.Imag)
 }
 
-    /// # Summary
 /// Returns the unary negation of an input of type `ComplexPolar`
 ///
 /// # Input
@@ -1222,7 +1153,6 @@ function NegationCP(input : ComplexPolar) : ComplexPolar {
     ComplexPolar(input.Magnitude, input.Argument + PI())
 }
 
-    /// # Summary
 /// Returns the sum of two inputs of type `Complex`.
 ///
 /// # Input
@@ -1237,7 +1167,6 @@ function PlusC(a : Complex, b : Complex) : Complex {
     Complex(a.Real + b.Real, a.Imag + b.Imag)
 }
 
-    /// # Summary
 /// Returns the sum of two inputs of type `ComplexPolar`.
 ///
 /// # Input
@@ -1257,7 +1186,6 @@ function PlusCP(a : ComplexPolar, b : ComplexPolar) : ComplexPolar {
     )
 }
 
-    /// # Summary
 /// Returns the difference between two inputs of type `Complex`.
 ///
 /// # Input
@@ -1272,7 +1200,6 @@ function MinusC(a : Complex, b : Complex) : Complex {
     Complex(a.Real - b.Real, a.Imag - b.Imag)
 }
 
-    /// # Summary
 /// Returns the difference between two inputs of type `ComplexPolar`.
 ///
 /// # Input
@@ -1287,7 +1214,6 @@ function MinusCP(a : ComplexPolar, b : ComplexPolar) : ComplexPolar {
     PlusCP(a, NegationCP(b))
 }
 
-    /// # Summary
 /// Returns the product of two inputs of type `Complex`.
 ///
 /// # Input
@@ -1305,7 +1231,6 @@ function TimesC(a : Complex, b : Complex) : Complex {
     )
 }
 
-    /// # Summary
 /// Returns the product of two inputs of type `ComplexPolar`.
 ///
 /// # Input
@@ -1323,7 +1248,6 @@ function TimesCP(a : ComplexPolar, b : ComplexPolar) : ComplexPolar {
     )
 }
 
-    /// # Summary
 /// Internal. Since it is easiest to define the power of two complex numbers
 /// in Cartesian form as returning in polar form, we define that here, then
 /// convert as needed.
@@ -1350,7 +1274,6 @@ internal function PowCAsCP(base : Complex, power : Complex) : ComplexPolar {
     ComplexPolar(magnitude, angle)
 }
 
-    /// # Summary
 /// Returns a number raised to a given power of type `Complex`.
 /// Note that this is a multi-valued function, but only one value is returned.
 ///
@@ -1366,7 +1289,6 @@ function PowC(a : Complex, power : Complex) : Complex {
     ComplexPolarAsComplex(PowCAsCP(a, power))
 }
 
-    /// # Summary
 /// Returns a number raised to a given power of type `ComplexPolar`.
 /// Note that this is a multi-valued function, but only one value is returned.
 ///
@@ -1382,7 +1304,6 @@ function PowCP(a : ComplexPolar, power : ComplexPolar) : ComplexPolar {
     PowCAsCP(ComplexPolarAsComplex(a), ComplexPolarAsComplex(power))
 }
 
-    /// # Summary
 /// Returns the quotient of two inputs of type `Complex`.
 ///
 /// # Input
@@ -1401,7 +1322,6 @@ function DividedByC(a : Complex, b : Complex) : Complex {
     )
 }
 
-    /// # Summary
 /// Returns the quotient of two inputs of type `ComplexPolar`.
 ///
 /// # Input
@@ -1435,7 +1355,6 @@ function SmallestFixedPoint(integerBits : Int, fractionalBits : Int) : Double {
         -(2.0^IntAsDouble(integerBits - 1))
 }
 
-    /// # Summary
 /// Returns the largest representable number for specific fixed point dimensions.
 ///
 /// # Input
diff --git a/library/std/src/Std/Measurement.qs b/library/std/src/Std/Measurement.qs
index 81009c68a6..a37580b552 100644
--- a/library/std/src/Std/Measurement.qs
+++ b/library/std/src/Std/Measurement.qs
@@ -29,7 +29,6 @@ operation MeasureAllZ(register : Qubit[]) : Result {
     Measure(Repeated(PauliZ, Length(register)), register)
 }
 
-    /// # Summary
 /// Measures each qubit in a given array in the standard basis.
 ///
 /// # Description
@@ -64,7 +63,6 @@ operation MeasureEachZ(register : Qubit[]) : Result[] {
     results
 }
 
-    /// # Summary
 /// Measures each qubit in a given array in the Z basis
 /// and resets them to a fixed initial state.
 ///
@@ -85,7 +83,6 @@ operation MResetEachZ(register : Qubit[]) : Result[] {
     results
 }
 
-    /// # Summary
 /// Measures a single qubit in the X basis,
 /// and resets it to a fixed initial state
 /// following the measurement.
@@ -108,7 +105,6 @@ operation MResetX(target : Qubit) : Result {
     MResetZ(target)
 }
 
-    /// # Summary
 /// Measures a single qubit in the Y basis,
 /// and resets it to a fixed initial state
 /// following the measurement.
@@ -133,7 +129,6 @@ operation MResetY(target : Qubit) : Result {
     MResetZ(target)
 }
 
-    /// # Summary
 /// Measures a single qubit in the Z basis,
 /// and resets it to a fixed initial state
 /// following the measurement.
@@ -153,7 +148,6 @@ operation MResetZ(target : Qubit) : Result {
     __quantum__qis__mresetz__body(target)
 }
 
-    /// # Summary
 /// Measures the content of a quantum register and converts
 /// it to an integer. The measurement is performed with respect
 /// to the standard computational basis, i.e., the eigenbasis of `PauliZ`.
diff --git a/library/std/src/Std/Random.qs b/library/std/src/Std/Random.qs
index 7709d6d2e1..e1007dbe5d 100644
--- a/library/std/src/Std/Random.qs
+++ b/library/std/src/Std/Random.qs
@@ -27,7 +27,6 @@ operation DrawRandomInt(min : Int, max : Int) : Int {
     body intrinsic;
 }
 
-    /// # Summary
 /// Draws a random real number from a uniform distribution
 /// in a given inclusive interval. Fails if `max < min`.
 ///
@@ -51,7 +50,6 @@ operation DrawRandomDouble(min : Double, max : Double) : Double {
     body intrinsic;
 }
 
-    /// # Summary
 /// Given a success probability, returns a single Bernoulli trial
 /// that is true with the given probability.
 ///
diff --git a/library/std/src/Std/ResourceEstimation.qs b/library/std/src/Std/ResourceEstimation.qs
index a834565d66..6dfd5540fd 100644
--- a/library/std/src/Std/ResourceEstimation.qs
+++ b/library/std/src/Std/ResourceEstimation.qs
@@ -15,7 +15,6 @@ function SingleVariant() : Int {
     return 0;
 }
 
-    /// # Summary
 /// Informs the resource estimator of the start of the code fragment
 /// for which estimates caching can be done. This function
 /// is only available when using resource estimator execution target.
@@ -37,7 +36,6 @@ function BeginEstimateCaching(name : String, variant : Int) : Bool {
     body intrinsic;
 }
 
-    /// # Summary
 /// Instructs the resource estimator to stop estimates caching
 /// because the code fragment in consideration is over. This function
 /// is only available when using resource estimator execution target.
@@ -58,42 +56,36 @@ function AuxQubitCount(amount : Int) : (Int, Int) {
     return (0, amount);
 }
 
-    /// # Summary
 /// Returns a tuple that can be passed to the `AccountForEstimates` operation
 /// to specify that the number of the T gates is equal to the `amount`.
 function TCount(amount : Int) : (Int, Int) {
     return (1, amount);
 }
 
-    /// # Summary
 /// Returns a tuple that can be passed to the `AccountForEstimates` operation
 /// to specify that the number of rotations is equal to the `amount`.
 function RotationCount(amount : Int) : (Int, Int) {
     return (2, amount);
 }
 
-    /// # Summary
 /// Returns a tuple that can be passed to the `AccountForEstimates` operation
 /// to specify that the rotation depth is equal to the `amount`.
 function RotationDepth(amount : Int) : (Int, Int) {
     return (3, amount);
 }
 
-    /// # Summary
 /// Returns a tuple that can be passed to the `AccountForEstimates` operation
 /// to specify that the number of the CCZ gates is equal to the `amount`.
 function CczCount(amount : Int) : (Int, Int) {
     return (4, amount);
 }
 
-    /// # Summary
 /// Returns a tuple that can be passed to the `AccountForEstimates` operation
 /// to specify that the number Measurements is equal to the `amount`.
 function MeasurementCount(amount : Int) : (Int, Int) {
     return (5, amount);
 }
 
-    /// # Summary
 /// Pass the value returned by the function to the `AccountForEstimates` operation
 /// to indicate Parallel Synthesis Sequential Pauli Computation (PSSPC) layout.
 /// See https://arxiv.org/pdf/2211.07629.pdf for details.
@@ -101,7 +93,6 @@ function PSSPCLayout() : Int {
     return 1;
 }
 
-    /// # Summary
 /// Account for the resource estimates of an unimplemented operation,
 /// which were obtained separately. This operation is only available
 /// when using resource estimator execution target.
@@ -126,7 +117,6 @@ internal operation AccountForEstimatesInternal(estimates : (Int, Int)[], layout
     body intrinsic;
 }
 
-    /// # Summary
 /// Instructs the resource estimator to assume that the resources from the
 /// call of this operation until a call to `EndRepeatEstimates` are
 /// accounted for `count` times, without the need to execute the code that many
@@ -148,7 +138,6 @@ internal operation BeginRepeatEstimatesInternal(count : Int) : Unit {
     body intrinsic;
 }
 
-    /// # Summary
 /// Companion operation to `BeginRepeatEstimates`.
 operation EndRepeatEstimates() : Unit {
     body ... {
@@ -161,7 +150,6 @@ internal operation EndRepeatEstimatesInternal() : Unit {
     body intrinsic;
 }
 
-    /// # Summary
 /// Instructs the resource estimator to assume that the resources from the
 /// call of this operation until a call to `Adjoint RepeatEstimates` are
 /// accounted for `count` times, without the need to execute the code that many
diff --git a/library/std/src/Std/Unstable/Arithmetic.qs b/library/std/src/Std/Unstable/Arithmetic.qs
deleted file mode 100644
index 8cd6f370b3..0000000000
--- a/library/std/src/Std/Unstable/Arithmetic.qs
+++ /dev/null
@@ -1,702 +0,0 @@
-// Copyright (c) Microsoft Corporation.
-// Licensed under the MIT License.
-
-
-import Std.Arrays.*;
-import QIR.Intrinsic.*;
-open Microsoft.Quantum.Diagnostics;
-open Microsoft.Quantum.Math;
-open Microsoft.Quantum.Convert;
-import Std.Unstable.ArithmeticHelpers.*;
-
-/// # Summary
-/// This applies the in-place majority operation to 3 qubits.
-///
-/// # Description
-/// Assuming the state of the input qubits are |x⟩, |y⟩ and |z⟩, then
-/// this operation performs the following transformation:
-/// |x⟩|y⟩|z⟩ ↦ |x ⊕ z⟩|y ⊕ z⟩MAJ(x, y, z).
-///
-/// # Input
-/// ## x
-/// The first input qubit.
-/// ## y
-/// The second input qubit.
-/// ## z
-/// A qubit onto which the majority function will be applied.
-operation MAJ(x : Qubit, y : Qubit, z : Qubit) : Unit is Adj + Ctl {
-    CNOT(z, y);
-    CNOT(z, x);
-    CCNOT(y, x, z);
-}
-
-    /// # Summary
-/// Reflects a quantum register about a given classical integer.
-///
-/// # Description
-/// Given a quantum register initially in the state ∑ᵢ(αᵢ|i⟩),
-/// where each |i⟩ is a basis state representing an integer i,
-/// reflects the state of the register about the basis state |j⟩
-/// for a given integer j: ∑ᵢ(-1)^(δᵢⱼ)(αᵢ|i⟩)
-/// This operation is implemented in-place, without explicit allocation of
-/// additional auxiliary qubits.
-///
-/// # Input
-/// ## index
-/// The classical integer j indexing the basis state about which to reflect.
-/// ## reg
-/// Little-endian quantum register to reflect.
-operation ReflectAboutInteger(index : Int, reg : Qubit[]) : Unit is Adj + Ctl {
-    within {
-        // Evaluation optimization for case index == 0
-        if index == 0 {
-            ApplyToEachA(X, reg);
-        } else {
-            // We want to reduce to the problem of reflecting about the all-ones
-            // state. To do that, we apply our reflection within an application
-            // of X instructions that flip all the zeros in our index.
-            ApplyPauliFromInt(PauliX, false, index, reg);
-        }
-    } apply {
-        Controlled ApplyAsSinglyControlled(Most(reg), (Z, Tail(reg)));
-    }
-}
-
-//
-// Add, Increment      |   Operation    | Description
-// ____________________|________________|_______________________________________________________________
-// y += 5              | IncByI, IncByL | Increment LE register in-place by integer
-// y += x              | IncByLE        | Increment LE register in-place by LE register
-// z = x + 5 (z was 0) |                | Add integer to LE register creating result out-of-place
-// z = x + y (z was 0) | AddLE          | Add two LE register creating result out-of-place
-// z += x + 5          |                | Increment LE register by the sum of integer and LE register
-// z += x + y          |                | Increment LE register by the sum of two LE registers
-//
-// IncByLE implementations:
-//     RippleCarryTTKIncByLE (default)
-//     RippleCarryCGIncByLE
-//     FourierTDIncByLE
-//     via IncByLEUsingAddLE and any out-of-place addition
-// IncByI implementations:
-//     via IncByIUsingIncByLE and any in-place LE adder
-// IncByL implementations:
-//     via IncByLUsingIncByLE and any in-place LE adder
-// AddLE implementations:
-//     RippleCarryCGAddLE (default)
-//     LookAheadDKRSAddLE
-//
-
-/// # Summary
-/// Increments a little-endian register ys by an integer number c
-///
-/// # Description
-/// Computes ys += c modulo 2ⁿ, where ys is a little-endian register,
-/// Length(ys) = n > 0, c is a Int number, 0 ≤ c < 2ⁿ.
-/// NOTE: Use IncByIUsingIncByLE directly if the choice of implementation
-/// is important.
-operation IncByI(c : Int, ys : Qubit[]) : Unit is Adj + Ctl {
-    IncByIUsingIncByLE(RippleCarryTTKIncByLE, c, ys);
-}
-
-    /// # Summary
-/// Increments a little-endian register ys by a BigInt number c
-///
-/// # Description
-/// Computes ys += c modulo 2ⁿ, where ys is a little-endian register,
-/// Length(ys) = n > 0, c is a BigInt number, 0 ≤ c < 2ⁿ.
-/// NOTE: Use IncByLUsingIncByLE directly if the choice of implementation
-/// is important.
-operation IncByL(c : BigInt, ys : Qubit[]) : Unit is Adj + Ctl {
-    IncByLUsingIncByLE(RippleCarryTTKIncByLE, c, ys);
-}
-
-    /// # Summary
-/// Increments a little-endian register ys by a little-endian register xs
-///
-/// # Description
-/// Computes ys += xs modulo 2ⁿ, where xs and ys are little-endian registers,
-/// and Length(xs) ≤ Length(ys) = n.
-/// NOTE: Use operations like RippleCarryCGIncByLE directly if
-/// the choice of implementation is important.
-operation IncByLE(xs : Qubit[], ys : Qubit[]) : Unit is Adj + Ctl {
-    RippleCarryTTKIncByLE(xs, ys);
-}
-
-    /// # Summary
-/// Sets a zero-initialized little-endian register zs to the sum of
-/// little-endian registers xs and ys
-///
-/// # Description
-/// Computes zs := xs + ys modulo 2ⁿ, where xs, ys, and zs are little-endian registers,
-/// Length(xs) = Length(ys) ≤ Length(zs) = n, assuming zs is 0-initialized.
-/// NOTE: Use operations like RippleCarryCGAddLE directly if
-/// the choice of implementation is important.
-operation AddLE(xs : Qubit[], ys : Qubit[], zs : Qubit[]) : Unit is Adj {
-    RippleCarryCGAddLE(xs, ys, zs);
-}
-
-    /// # Summary
-/// Reversible, in-place ripple-carry addition of two integers.
-///
-/// # Description
-/// Computes ys += xs modulo 2ⁿ, where xs and ys are little-endian registers,
-/// and Length(xs) ≤ Length(ys) = n.
-/// This operation uses the ripple-carry algorithm.
-/// Note that if Length(ys) >= Length(xs)+2, xs is padded with 0-initialized
-/// qubits to match ys's length. The operation doesn't use any auxiliary
-/// qubits otherwise.
-///
-/// # References
-/// - [arXiv:0910.2530](https://arxiv.org/abs/0910.2530)
-///   "Quantum Addition Circuits and Unbounded Fan-Out",
-///   Yasuhiro Takahashi, Seiichiro Tani, Noboru Kunihiro
-operation RippleCarryTTKIncByLE(xs : Qubit[], ys : Qubit[]) : Unit is Adj + Ctl {
-    let xsLen = Length(xs);
-    let ysLen = Length(ys);
-
-    Fact(ysLen >= xsLen, "Register `ys` must be longer than register `xs`.");
-    Fact(xsLen >= 1, "Registers `xs` and `ys` must contain at least one qubit.");
-
-    if xsLen == ysLen {
-        if xsLen > 1 {
-            within {
-                ApplyOuterTTKAdder(xs, ys);
-            } apply {
-                ApplyInnerTTKAdderNoCarry(xs, ys);
-            }
-        }
-        CNOT(xs[0], ys[0]);
-    } elif xsLen + 1 == ysLen {
-        if xsLen > 1 {
-            CNOT(xs[xsLen - 1], ys[ysLen - 1]);
-            within {
-                ApplyOuterTTKAdder(xs, ys);
-            } apply {
-                ApplyInnerTTKAdderWithCarry(xs, ys);
-            }
-        } else {
-            CCNOT(xs[0], ys[0], ys[1]);
-        }
-        CNOT(xs[0], ys[0]);
-    } elif xsLen + 2 <= ysLen {
-        // Pad xs so that its length is one qubit shorter than ys.
-        use padding = Qubit[ysLen - xsLen - 1];
-        RippleCarryTTKIncByLE(xs + padding, ys);
-    }
-}
-
-    /// # Summary
-/// Increments a little-endian register ys by a little-endian register xs
-/// using the ripple-carry algorithm.
-///
-/// # Description
-/// Computes ys += xs modulo 2ⁿ, where xs and ys are little-endian registers,
-/// and Length(xs) ≤ Length(ys) = n.
-/// Note that if Length(xs) != Length(ys), xs is padded with 0-initialized
-/// qubits to match ys's length.
-/// This operation uses the ripple-carry algorithm.
-///
-/// # Reference
-/// - [arXiv:1709.06648](https://arxiv.org/pdf/1709.06648.pdf)
-///   "Halving the cost of quantum addition", Craig Gidney.
-operation RippleCarryCGIncByLE(xs : Qubit[], ys : Qubit[]) : Unit is Adj + Ctl {
-    let xsLen = Length(xs);
-    let ysLen = Length(ys);
-
-    Fact(ysLen >= xsLen, "Register `ys` must be longer than register `xs`.");
-    Fact(xsLen >= 1, "Registers `xs` and `ys` must contain at least one qubit.");
-
-    if ysLen - xsLen >= 2 {
-        // Pad xs so that its length is one qubit shorter than ys.
-        use padding = Qubit[ysLen - xsLen - 1];
-        RippleCarryCGIncByLE(xs + padding, ys);
-    } elif xsLen == 1 {
-        if ysLen == 1 {
-            CNOT(xs[0], ys[0]);
-        } elif ysLen == 2 {
-            HalfAdderForInc(xs[0], ys[0], ys[1]);
-        }
-    } else {
-        use carries = Qubit[xsLen];
-        within {
-            ApplyAndAssuming0Target(xs[0], ys[0], carries[0]);
-        } apply {
-            for i in 1..xsLen - 2 {
-                CarryForInc(carries[i - 1], xs[i], ys[i], carries[i]);
-            }
-            if xsLen == ysLen {
-                within {
-                    CNOT(carries[xsLen - 2], xs[xsLen - 1]);
-                } apply {
-                    CNOT(xs[xsLen - 1], ys[xsLen - 1]);
-                }
-            } else {
-                FullAdderForInc(carries[xsLen - 2], xs[xsLen - 1], ys[xsLen - 1], ys[xsLen]);
-            }
-            for i in xsLen - 2..-1..1 {
-                UncarryForInc(carries[i - 1], xs[i], ys[i], carries[i]);
-            }
-        }
-        CNOT(xs[0], ys[0]);
-    }
-}
-
-    /// # Summary
-/// Sets a zero-initialized little-endian register zs to the sum of
-/// little-endian registers xs and ys using the ripple-carry algorithm.
-///
-/// # Description
-/// Computes zs := xs + ys + zs[0] modulo 2ⁿ, where xs, ys, and zs are
-/// little-endian registers, Length(xs) = Length(ys) ≤ Length(zs) = n,
-/// assuming zs is 0-initialized, except for maybe zs[0], which can be
-// in |0> or |1> state and can be used as carry-in.
-/// This operation uses the ripple-carry algorithm.
-/// NOTE: `zs[Length(xs)]` can be used as carry-out, if `zs` is longer than `xs`.
-///
-/// # Reference
-/// - [arXiv:1709.06648](https://arxiv.org/pdf/1709.06648.pdf)
-///   "Halving the cost of quantum addition", Craig Gidney.
-operation RippleCarryCGAddLE(xs : Qubit[], ys : Qubit[], zs : Qubit[]) : Unit is Adj {
-    let xsLen = Length(xs);
-    let zsLen = Length(zs);
-    Fact(Length(ys) == xsLen, "Registers `xs` and `ys` must be of same length.");
-    Fact(zsLen >= xsLen, "Register `zs` must be no shorter than register `xs`.");
-
-    // Since zs is zero-initialized, its bits at indexes higher than
-    // xsLen remain unused as there will be no carry into them.
-    let top = MinI(zsLen - 2, xsLen - 1);
-    for k in 0..top {
-        FullAdder(zs[k], xs[k], ys[k], zs[k + 1]);
-    }
-
-    if xsLen > 0 and xsLen == zsLen {
-        CNOT(Tail(xs), Tail(zs));
-        CNOT(Tail(ys), Tail(zs));
-    }
-}
-
-    /// # Summary
-/// Sets a zero-initialized little-endian register zs to the sum of
-/// little-endian registers xs and ys using the carry-lookahead algorithm.
-///
-/// # Description
-/// Computes zs := xs + ys + zs[0] modulo 2ⁿ, where xs, ys, and zs are
-/// little-endian registers, Length(xs) = Length(ys) ≤ Length(zs) = n,
-/// assuming zs is 0-initialized, except for maybe zs[0], which can be
-/// in |0> or |1> state and can be used as carry-in.
-/// NOTE: `zs[Length(xs)]` can be used as carry-out, if `zs` is longer than `xs`.
-/// This operation uses the carry-lookahead algorithm.
-///
-/// # Reference
-/// - [arXiv:quant-ph/0406142](https://arxiv.org/abs/quant-ph/0406142)
-///   "A logarithmic-depth quantum carry-lookahead adder",
-///   Thomas G. Draper, Samuel A. Kutin, Eric M. Rains, Krysta M. Svore
-operation LookAheadDKRSAddLE(xs : Qubit[], ys : Qubit[], zs : Qubit[]) : Unit is Adj {
-    let xsLen = Length(xs);
-    let zsLen = Length(zs);
-    Fact(Length(ys) == xsLen, "Registers `xs` and `ys` must be of same length.");
-    Fact(zsLen >= xsLen, "Register `zs` must be no shorter than register `xs`.");
-
-    if zsLen > xsLen {
-        // with carry-out
-        // compute initial generate values
-        for k in 0..xsLen - 1 {
-            ApplyAndAssuming0Target(xs[k], ys[k], zs[k + 1]);
-        }
-
-        within {
-            // compute initial propagate values
-            for i in IndexRange(xs) {
-                CNOT(xs[i], ys[i]);
-            }
-        } apply {
-            if xsLen > 1 {
-                ComputeCarries(Rest(ys), zs[1..xsLen]);
-            }
-
-            // compute sum into carries
-            for k in 0..xsLen - 1 {
-                CNOT(ys[k], zs[k]);
-            }
-        }
-    } else {
-        // xsLen == zsLen, so without carry-out
-        LookAheadDKRSAddLE(Most(xs), Most(ys), zs);
-        CNOT(Tail(xs), Tail(zs));
-        CNOT(Tail(ys), Tail(zs));
-    }
-}
-
-    /// # Summary
-/// Increments a little-endian register ys by a little-endian register xs
-/// using Quantum Fourier Transform.
-///
-/// # Description
-/// Computes ys += xs modulo 2ⁿ, where xs and ys are little-endian registers,
-/// and Length(xs) = Length(ys) = n.
-/// This operation uses Quantum Fourier Transform.
-///
-/// # Reference
-/// - [arXiv:quant-ph/0008033](https://arxiv.org/abs/quant-ph/0008033)
-///   "Addition on a Quantum Computer", Thomas G. Draper
-operation FourierTDIncByLE(xs : Qubit[], ys : Qubit[]) : Unit is Adj + Ctl {
-    within {
-        ApplyQFT(ys);
-    } apply {
-        for i in IndexRange(xs) {
-            Controlled PhaseGradient([xs[i]], ys[i...]);
-        }
-    }
-}
-
-    /// # Summary
-/// Increments a little-endian register ys by a BigInt number c
-/// using provided adder.
-///
-/// # Description
-/// Computes ys += c modulo 2ⁿ, where ys is a little-endian register
-/// Length(ys) = n > 0, c is a BigInt number, 0 ≤ c < 2ⁿ.
-operation IncByLUsingIncByLE(
-    adder : (Qubit[], Qubit[]) => Unit is Adj + Ctl,
-    c : BigInt,
-    ys : Qubit[]
-) : Unit is Adj + Ctl {
-
-    let ysLen = Length(ys);
-    Fact(ysLen > 0, "Length of `ys` must be at least 1.");
-    Fact(c >= 0L, "Constant `c` must be non-negative.");
-    Fact(c < 2L^ysLen, "Constant `c` must be smaller than 2^Length(ys).");
-
-    if c != 0L {
-        // If c has j trailing zeros, then the j least significant
-        // bits of y won't be affected by the addition and can
-        // therefore be ignored by applying the addition only to
-        // the other qubits and shifting c accordingly.
-        let j = TrailingZeroCountL(c);
-        use x = Qubit[ysLen - j];
-        within {
-            ApplyXorInPlaceL(c >>> j, x);
-        } apply {
-            adder(x, ys[j...]);
-        }
-    }
-}
-
-    /// # Summary
-/// Increments a little-endian register ys by an Int number c
-/// using provided adder.
-///
-/// # Description
-/// Computes ys += c modulo 2ⁿ, where ys is a little-endian register
-/// Length(ys) = n > 0, c is an Int number, 0 ≤ c < 2ⁿ.
-operation IncByIUsingIncByLE(
-    adder : (Qubit[], Qubit[]) => Unit is Adj + Ctl,
-    c : Int,
-    ys : Qubit[]
-) : Unit is Adj + Ctl {
-
-    let ysLen = Length(ys);
-    Fact(ysLen > 0, "Length of `ys` must be at least 1.");
-    Fact(c >= 0, "Constant `c` must be non-negative.");
-    Fact(c < 2^ysLen, "Constant `c` must be smaller than 2^Length(ys).");
-
-    if c != 0 {
-        // If c has j trailing zeros than the j least significant
-        // bits of y won't be affected by the addition and can
-        // therefore be ignored by applying the addition only to
-        // the other qubits and shifting c accordingly.
-        let j = TrailingZeroCountI(c);
-        use x = Qubit[ysLen - j];
-        within {
-            ApplyXorInPlace(c >>> j, x);
-        } apply {
-            adder(x, ys[j...]);
-        }
-    }
-}
-
-    /// # Summary
-/// Generic operation to turn two out-place adders into one in-place adder
-///
-/// # Description
-/// This implementation allows to specify two distinct adders for forward
-/// and backward direction.  The forward adder is always applied in its
-/// body variant, whereas the backward adder is always applied in its adjoint
-/// variant.  Therefore, it's possible to, for example, use the ripple-carry
-/// out-of-place adder in backwards direction to require no T gates.
-///
-/// The controlled variant is also optimized in a way that everything but
-/// the adders is controlled,
-///
-/// # Reference
-/// - [arXiv:2012.01624](https://arxiv.org/abs/2012.01624)
-///   "Quantum block lookahead adders and the wait for magic states",
-///   Craig Gidney.
-operation IncByLEUsingAddLE(
-    forwardAdder : (Qubit[], Qubit[], Qubit[]) => Unit is Adj,
-    backwardAdder : (Qubit[], Qubit[], Qubit[]) => Unit is Adj,
-    xs : Qubit[],
-    ys : Qubit[]
-) : Unit is Adj + Ctl {
-
-    body (...) {
-        let n = Length(xs);
-
-        Fact(Length(ys) == n, "Registers xs and ys must be of same length");
-
-        use qs = Qubit[n];
-
-        forwardAdder(xs, ys, qs);
-        for i in IndexRange(ys) {
-            SWAP(ys[i], qs[i]);
-        }
-        ApplyToEachA(X, qs);
-        within {
-            ApplyToEachA(X, ys);
-        } apply {
-            Adjoint backwardAdder(xs, ys, qs);
-        }
-    }
-    adjoint (...) {
-        let n = Length(xs);
-
-        Fact(Length(ys) == n, "Registers xs and ys must be of same length");
-
-        use qs = Qubit[n];
-
-        within {
-            ApplyToEachA(X, ys);
-        } apply {
-            forwardAdder(xs, ys, qs);
-        }
-        ApplyToEachA(X, qs);
-        for i in IndexRange(ys) {
-            SWAP(ys[i], qs[i]);
-        }
-        Adjoint backwardAdder(xs, ys, qs);
-    }
-    controlled (ctls, ...) {
-        // When we control everything except the adders, the adders will
-        // cancel themselves.
-        let n = Length(xs);
-
-        Fact(Length(ys) == n, "Registers xs and ys must be of same length");
-
-        use qs = Qubit[n];
-
-        forwardAdder(xs, ys, qs);
-        for i in IndexRange(ys) {
-            Controlled SWAP(ctls, (ys[i], qs[i]))
-        }
-        ApplyToEachA(tgt => Controlled X(ctls, tgt), qs);
-        within {
-            ApplyToEachA(tgt => Controlled X(ctls, tgt), ys);
-        } apply {
-            Adjoint backwardAdder(xs, ys, qs);
-        }
-    }
-    controlled adjoint (ctls, ...) {
-        // When we control everything except the adders, the adders will
-        // cancel themselves.
-        let n = Length(xs);
-
-        Fact(Length(ys) == n, "Registers xs and ys must be of same length");
-
-        use qs = Qubit[n];
-
-        within {
-            ApplyToEachA(tgt => Controlled X(ctls, tgt), ys);
-        } apply {
-            forwardAdder(xs, ys, qs);
-        }
-        ApplyToEachA(tgt => Controlled X(ctls, tgt), qs);
-        for i in IndexRange(ys) {
-            Controlled SWAP(ctls, (ys[i], qs[i]))
-        }
-        Adjoint backwardAdder(xs, ys, qs);
-    }
-}
-
-//
-// Comparisons
-//
-// Compare BigInt and qubit register in a little-endian format and apply action
-//   if c < x { action(target) }  | ApplyIfLessL
-//   if c <= x { action(target) } | ApplyIfLessOrEqualL
-//   if c == x { action(target) } | ApplyIfEqualL
-//   if c >= x { action(target) } | ApplyIfGreaterOrEqualL
-//   if c > x { action(target) }  | ApplyIfGreaterL
-//
-// Compare two qubit registers in a little-endian format and apply action
-//   if x < y { action(target) }  | ApplyIfLessLE
-//   if x <= y { action(target) } | ApplyIfLessOrEqualLE
-//   if x == y { action(target) } | ApplyIfEqualLE
-//   if x >= y { action(target) } | ApplyIfGreaterOrEqualLE
-//   if x > y { action(target) }  | ApplyIfGreaterLE
-//
-
-/// # Summary
-/// Computes `if (c < x) { action(target) }`, that is, applies `action` to `target`
-/// if a BigInt value `c` is less than the little-endian qubit register `x`
-operation ApplyIfLessL<'T>(
-    action : 'T => Unit is Adj + Ctl,
-    c : BigInt,
-    x : Qubit[],
-    target : 'T
-) : Unit is Adj + Ctl {
-
-    ApplyActionIfGreaterThanOrEqualConstant(false, action, c + 1L, x, target);
-}
-
-    /// # Summary
-/// Computes `if (c <= x) { action(target) }`, that is, applies `action` to `target`
-/// if a BigInt value `c` is less or equal to the little-endian qubit register `x`
-operation ApplyIfLessOrEqualL<'T>(
-    action : 'T => Unit is Adj + Ctl,
-    c : BigInt,
-    x : Qubit[],
-    target : 'T
-) : Unit is Adj + Ctl {
-
-    ApplyActionIfGreaterThanOrEqualConstant(false, action, c, x, target);
-}
-
-    /// # Summary
-/// Computes `if (c == x) { action(target) }`, that is, applies `action` to `target`
-/// if a BigInt value `c` is equal to the little-endian qubit register `x`
-operation ApplyIfEqualL<'T>(
-    action : 'T => Unit is Adj + Ctl,
-    c : BigInt,
-    xs : Qubit[],
-    target : 'T
-) : Unit is Adj + Ctl {
-
-    let cBitSize = BitSizeL(c);
-    let xLen = Length(xs);
-    if (cBitSize <= xLen) {
-        let bits = BigIntAsBoolArray(c, Length(xs));
-        within {
-            ApplyPauliFromBitString(PauliX, false, bits, xs);
-        } apply {
-            Controlled ApplyAsSinglyControlled(xs, (a => action(a), target));
-        }
-    }
-}
-
-    /// # Summary
-/// Computes `if (c >= x) { action(target) }`, that is, applies `action` to `target`
-/// if a BigInt value `c` is greater or equal to the little-endian qubit register `x`
-operation ApplyIfGreaterOrEqualL<'T>(
-    action : 'T => Unit is Adj + Ctl,
-    c : BigInt,
-    x : Qubit[],
-    target : 'T
-) : Unit is Adj + Ctl {
-
-    ApplyActionIfGreaterThanOrEqualConstant(true, action, c + 1L, x, target);
-}
-
-    /// # Summary
-/// Computes `if (c > x) { action(target) }`, that is, applies `action` to `target`
-/// if a BigInt value `c` is greater than the little-endian qubit register `x`
-operation ApplyIfGreaterL<'T>(
-    action : 'T => Unit is Adj + Ctl,
-    c : BigInt,
-    x : Qubit[],
-    target : 'T
-) : Unit is Adj + Ctl {
-
-    ApplyActionIfGreaterThanOrEqualConstant(true, action, c, x, target);
-}
-
-    /// # Summary
-/// Computes `if x < y { action(target) }`, that is, applies `action` to `target`
-/// if register `x` is less than the register `y`.
-/// Both qubit registers should be in a little-endian format.
-operation ApplyIfLessLE<'T>(
-    action : 'T => Unit is Adj + Ctl,
-    x : Qubit[],
-    y : Qubit[],
-    target : 'T
-) : Unit is Adj + Ctl {
-
-    ApplyIfGreaterLE(action, y, x, target);
-}
-
-    /// # Summary
-/// Computes `if x <= y { action(target) }`, that is, applies `action` to `target`
-/// if register `x` is less or equal to the register `y`.
-/// Both qubit registers should be in a little-endian format.
-operation ApplyIfLessOrEqualLE<'T>(
-    action : 'T => Unit is Adj + Ctl,
-    x : Qubit[],
-    y : Qubit[],
-    target : 'T
-) : Unit is Adj + Ctl {
-
-    Fact(Length(x) > 0, "Bitwidth must be at least 1");
-    within {
-        ApplyToEachA(X, x);
-    } apply {
-        // control is not inverted
-        ApplyActionIfSumOverflows(action, x, y, false, target);
-    }
-}
-
-    /// # Summary
-/// Computes `if x == y { action(target) }`, that is, applies `action` to `target`
-/// if register `x` is equal to the register `y`.
-/// Both qubit registers should be in a little-endian format.
-operation ApplyIfEqualLE<'T>(
-    action : 'T => Unit is Adj + Ctl,
-    x : Qubit[],
-    y : Qubit[],
-    target : 'T
-) : Unit is Adj + Ctl {
-
-    Fact(Length(x) == Length(y), "x and y must be of same length");
-    within {
-        for i in IndexRange(x) {
-            CNOT(x[i], y[i]);
-            X(y[i]);
-        }
-    } apply {
-        Controlled ApplyAsSinglyControlled(y, (a => action(a), target))
-    }
-}
-
-    /// # Summary
-/// Computes `if x >= y { action(target) }`, that is, applies `action` to `target`
-/// if register `x` is greater or equal to the register `y`.
-/// Both qubit registers should be in a little-endian format.
-operation ApplyIfGreaterOrEqualLE<'T>(
-    action : 'T => Unit is Adj + Ctl,
-    x : Qubit[],
-    y : Qubit[],
-    target : 'T
-) : Unit is Adj + Ctl {
-
-    ApplyIfLessOrEqualLE(action, y, x, target);
-}
-
-    /// # Summary
-/// Computes `if x > y { action(target) }`, that is, applies `action` to `target`
-/// if register `x` is greater than the register `y`.
-/// Both qubit registers should be in a little-endian format.
-operation ApplyIfGreaterLE<'T>(
-    action : 'T => Unit is Adj + Ctl,
-    x : Qubit[],
-    y : Qubit[],
-    target : 'T
-) : Unit is Adj + Ctl {
-
-    Fact(Length(x) > 0, "Bitwidth must be at least 1");
-    within {
-        ApplyToEachA(X, x);
-    } apply {
-        // control is inverted
-        ApplyActionIfSumOverflows(action, x, y, true, target);
-    }
-}
-
-export AddLE, ApplyIfEqualLE, ApplyIfEqualL, ApplyIfGreaterLE, ApplyIfGreaterL, ApplyIfGreaterOrEqualLE, ApplyIfGreaterOrEqualL, ApplyIfLessLE, ApplyIfLessL, ApplyIfLessOrEqualLE, ApplyIfLessOrEqualL, IncByI, IncByIUsingIncByLE, IncByL, IncByLUsingIncByLE, IncByLE, IncByLEUsingAddLE, LookAheadDKRSAddLE, MAJ, ReflectAboutInteger, RippleCarryCGAddLE, RippleCarryCGIncByLE, RippleCarryTTKIncByLE, FourierTDIncByLE;
diff --git a/library/std/src/Std/Unstable/ArithmeticHelpers.qs b/library/std/src/Std/Unstable/ArithmeticHelpers.qs
deleted file mode 100644
index bb28c85c4f..0000000000
--- a/library/std/src/Std/Unstable/ArithmeticHelpers.qs
+++ /dev/null
@@ -1,586 +0,0 @@
-// Copyright (c) Microsoft Corporation.
-// Licensed under the MIT License.
-
-
-open Microsoft.Quantum.Diagnostics;
-import Std.Arrays.*;
-open Microsoft.Quantum.Math;
-open Microsoft.Quantum.Convert;
-import QIR.Intrinsic.*;
-
-
-/// # Summary
-/// Implements the outer operation for RippleCarryTTKIncByLE to conjugate
-/// the inner operation to construct the full adder. Only Length(xs)
-/// qubits are processed.
-///
-/// # Input
-/// ## xs
-/// Qubit register in a little-endian format containing the first summand
-/// input to RippleCarryTTKIncByLE.
-/// ## ys
-/// Qubit register in a little-endian format containing the second summand
-/// input to RippleCarryTTKIncByLE.
-///
-/// # References
-/// - Yasuhiro Takahashi, Seiichiro Tani, Noboru Kunihiro: "Quantum
-///   Addition Circuits and Unbounded Fan-Out", Quantum Information and
-///   Computation, Vol. 10, 2010.
-///   https://arxiv.org/abs/0910.2530
-internal operation ApplyOuterTTKAdder(xs : Qubit[], ys : Qubit[]) : Unit is Adj + Ctl {
-    Fact(Length(xs) <= Length(ys), "Input register ys must be at lease as long as xs.");
-    for i in 1..Length(xs) - 1 {
-        CNOT(xs[i], ys[i]);
-    }
-    for i in Length(xs) - 2..-1..1 {
-        CNOT(xs[i], xs[i + 1]);
-    }
-}
-
-    /// # Summary
-/// Implements the inner addition function for the operation
-/// RippleCarryTTKIncByLE. This is the inner operation that is conjugated
-/// with the outer operation to construct the full adder.
-///
-/// # Input
-/// ## xs
-/// Qubit register in a little-endian format containing the first summand
-/// input to RippleCarryTTKIncByLE.
-/// ## ys
-/// Qubit register in a little-endian format containing the second summand
-/// input to RippleCarryTTKIncByLE.
-///
-/// # References
-/// - Yasuhiro Takahashi, Seiichiro Tani, Noboru Kunihiro: "Quantum
-///   Addition Circuits and Unbounded Fan-Out", Quantum Information and
-///   Computation, Vol. 10, 2010.
-///   https://arxiv.org/abs/0910.2530
-///
-/// # Remarks
-/// The specified controlled operation makes use of symmetry and mutual
-/// cancellation of operations to improve on the default implementation
-/// that adds a control to every operation.
-internal operation ApplyInnerTTKAdderNoCarry(xs : Qubit[], ys : Qubit[]) : Unit is Adj + Ctl {
-    body (...) {
-        (Controlled ApplyInnerTTKAdderNoCarry)([], (xs, ys));
-    }
-    controlled (controls, ...) {
-        Fact(Length(xs) == Length(ys), "Input registers must have the same number of qubits.");
-
-        for idx in 0..Length(xs) - 2 {
-            CCNOT(xs[idx], ys[idx], xs[idx + 1]);
-        }
-        for idx in Length(xs) - 1..-1..1 {
-            Controlled CNOT(controls, (xs[idx], ys[idx]));
-            CCNOT(xs[idx - 1], ys[idx - 1], xs[idx]);
-        }
-    }
-}
-
-    /// # Summary
-/// Implements the inner addition function for the operation
-/// RippleCarryTTKIncByLE. This is the inner operation that is conjugated
-/// with the outer operation to construct the full adder.
-///
-/// # Input
-/// ## xs
-/// Qubit register in a little-endian format containing the first summand
-/// input to RippleCarryTTKIncByLE.
-/// ## ys
-/// Qubit register in a little-endian format containing the second summand
-/// input to RippleCarryTTKIncByLE.
-///
-/// # References
-/// - Yasuhiro Takahashi, Seiichiro Tani, Noboru Kunihiro: "Quantum
-///   Addition Circuits and Unbounded Fan-Out", Quantum Information and
-///   Computation, Vol. 10, 2010.
-///   https://arxiv.org/abs/0910.2530
-///
-/// # Remarks
-/// The specified controlled operation makes use of symmetry and mutual
-/// cancellation of operations to improve on the default implementation
-/// that adds a control to every operation.
-internal operation ApplyInnerTTKAdderWithCarry(xs : Qubit[], ys : Qubit[]) : Unit is Adj + Ctl {
-    body (...) {
-        (Controlled ApplyInnerTTKAdderWithCarry)([], (xs, ys));
-    }
-    controlled (controls, ...) {
-        Fact(Length(xs) + 1 == Length(ys), "ys must be one qubit longer then xs.");
-        Fact(Length(xs) > 0, "Array should not be empty.");
-
-
-        let nQubits = Length(xs);
-        for idx in 0..nQubits - 2 {
-            CCNOT(xs[idx], ys[idx], xs[idx + 1]);
-        }
-        (Controlled CCNOT)(controls, (xs[nQubits - 1], ys[nQubits - 1], ys[nQubits]));
-        for idx in nQubits - 1..-1..1 {
-            Controlled CNOT(controls, (xs[idx], ys[idx]));
-            CCNOT(xs[idx - 1], ys[idx - 1], xs[idx]);
-        }
-    }
-}
-
-    /// # Summary
-/// Implements Half-adder. Adds qubit x to qubit y and sets carryOut appropriately
-internal operation HalfAdderForInc(x : Qubit, y : Qubit, carryOut : Qubit) : Unit is Adj + Ctl {
-    body (...) {
-        CCNOT(x, y, carryOut);
-        CNOT(x, y);
-    }
-    adjoint auto;
-
-    controlled (ctls, ...) {
-        Fact(Length(ctls) == 1, "HalfAdderForInc should be controlled by exactly one control qubit.");
-
-        let ctl = ctls[0];
-        use helper = Qubit();
-
-        within {
-            ApplyAndAssuming0Target(x, y, helper);
-        } apply {
-            ApplyAndAssuming0Target(ctl, helper, carryOut);
-        }
-        CCNOT(ctl, x, y);
-    }
-    controlled adjoint auto;
-}
-
-    /// # Summary
-/// Implements Full-adder. Adds qubit carryIn and x to qubit y and sets carryOut appropriately.
-internal operation FullAdderForInc(carryIn : Qubit, x : Qubit, y : Qubit, carryOut : Qubit) : Unit is Adj + Ctl {
-    body (...) {
-        // TODO: cannot use `Carry` operation here
-        CNOT(carryIn, x);
-        CNOT(carryIn, y);
-        CCNOT(x, y, carryOut);
-        CNOT(carryIn, carryOut);
-        CNOT(carryIn, x);
-        CNOT(x, y);
-    }
-    adjoint auto;
-
-    controlled (ctls, ...) {
-        Fact(Length(ctls) == 1, "FullAdderForInc should be controlled by exactly one control qubit.");
-
-        let ctl = ctls[0];
-        use helper = Qubit();
-
-        CarryForInc(carryIn, x, y, helper);
-        CCNOT(ctl, helper, carryOut);
-        Controlled UncarryForInc(ctls, (carryIn, x, y, helper));
-    }
-    controlled adjoint auto;
-}
-
-// Computes carryOut := carryIn + x + y
-internal operation FullAdder(carryIn : Qubit, x : Qubit, y : Qubit, carryOut : Qubit) : Unit is Adj {
-    CNOT(x, y);
-    CNOT(x, carryIn);
-    ApplyAndAssuming0Target(y, carryIn, carryOut);
-    CNOT(x, y);
-    CNOT(x, carryOut);
-    CNOT(y, carryIn);
-}
-
-    /// # Summary
-/// Computes carry bit for a full adder.
-internal operation CarryForInc(carryIn : Qubit, x : Qubit, y : Qubit, carryOut : Qubit) : Unit is Adj + Ctl {
-    body (...) {
-        CNOT(carryIn, x);
-        CNOT(carryIn, y);
-        ApplyAndAssuming0Target(x, y, carryOut);
-        CNOT(carryIn, carryOut);
-    }
-    adjoint auto;
-    controlled (ctls, ...) {
-        // This CarryForInc is intended to be used only in an in-place
-        // ripple-carry implementation. Only such particular use case allows
-        // for this simple implementation where controlled version
-        // is the same as uncontrolled body.
-        CarryForInc(carryIn, x, y, carryOut);
-    }
-    controlled adjoint auto;
-}
-
-    /// # Summary
-/// Uncomputes carry bit for a full adder.
-internal operation UncarryForInc(carryIn : Qubit, x : Qubit, y : Qubit, carryOut : Qubit) : Unit is Adj + Ctl {
-    body (...) {
-        CNOT(carryIn, carryOut);
-        Adjoint ApplyAndAssuming0Target(x, y, carryOut);
-        CNOT(carryIn, x);
-        CNOT(x, y);
-    }
-    adjoint auto;
-    controlled (ctls, ...) {
-        Fact(Length(ctls) == 1, "UncarryForInc should be controlled by exactly one control qubit.");
-
-        let ctl = ctls[0];
-
-        CNOT(carryIn, carryOut);
-        Adjoint ApplyAndAssuming0Target(x, y, carryOut);
-        CCNOT(ctl, x, y); // Controlled X(ctls + [x], y);
-        CNOT(carryIn, x);
-        CNOT(carryIn, y);
-    }
-    controlled adjoint auto;
-}
-
-    /// # Summary
-/// Applies AND gate between `control1` and `control2` and stores the result
-/// in `target` assuming `target` is in |0> state.
-///
-/// # Description
-/// Inverts `target` if and only if both controls are 1, but assumes that
-/// `target` is in state 0. The operation has T-count 4, T-depth 2 and
-/// requires no helper qubit, and may therefore be preferable to a CCNOT
-/// operation, if `target` is known to be 0.
-/// The adjoint of this operation is measurement based and requires no T
-/// gates (but requires target to support branching on measurements).
-/// Although the Toffoli gate (CCNOT) will perform faster in simulations,
-/// this version has lower T gate requirements.
-/// # References
-/// - Cody Jones: "Novel constructions for the fault-tolerant Toffoli gate",
-///   Phys. Rev. A 87, 022328, 2013
-///   [arXiv:1212.5069](https://arxiv.org/abs/1212.5069)
-///   doi:10.1103/PhysRevA.87.022328
-@Config(Adaptive)
-internal operation ApplyAndAssuming0Target(control1 : Qubit, control2 : Qubit, target : Qubit) : Unit is Adj {
-    // NOTE: Eventually this operation will be public and intrinsic.
-    body (...) {
-        CCNOT(control1, control2, target);
-    }
-    adjoint (...) {
-        H(target);
-        if M(target) == One {
-            Reset(target);
-            CZ(control1, control2);
-        }
-    }
-}
-
-internal operation ApplyOrAssuming0Target(control1 : Qubit, control2 : Qubit, target : Qubit) : Unit is Adj {
-    within {
-        X(control1);
-        X(control2);
-    } apply {
-        ApplyAndAssuming0Target(control1, control2, target);
-        X(target);
-    }
-}
-
-    /// # Summary
-/// Applies AND gate between `control1` and `control2` and stores the result
-/// in `target` assuming `target` is in |0> state.
-///
-/// # Description
-/// Inverts `target` if and only if both controls are 1, but assumes that
-/// `target` is in state 0. The operation has T-count 4, T-depth 2 and
-/// requires no helper qubit, and may therefore be preferable to a CCNOT
-/// operation, if `target` is known to be 0.
-/// This version is suitable for Base profile.
-/// Although the Toffoli gate (CCNOT) will perform faster in simulations,
-/// this version has lower T gate requirements.
-/// # References
-/// - Cody Jones: "Novel constructions for the fault-tolerant Toffoli gate",
-///   Phys. Rev. A 87, 022328, 2013
-///   [arXiv:1212.5069](https://arxiv.org/abs/1212.5069)
-///   doi:10.1103/PhysRevA.87.022328
-@Config(not Adaptive)
-internal operation ApplyAndAssuming0Target(control1 : Qubit, control2 : Qubit, target : Qubit) : Unit is Adj {
-    H(target);
-    T(target);
-    CNOT(control1, target);
-    CNOT(control2, target);
-    within {
-        CNOT(target, control1);
-        CNOT(target, control2);
-    } apply {
-        Adjoint T(control1);
-        Adjoint T(control2);
-        T(target);
-    }
-    H(target);
-    S(target);
-}
-
-    /// # Summary
-/// Computes carries for the look-ahead adder
-internal operation ComputeCarries(ps : Qubit[], gs : Qubit[]) : Unit is Adj {
-    let n = Length(gs);
-    Fact(Length(ps) + 1 == n, "Register gs must be one qubit longer than register gs.");
-
-    let T = Floor(Lg(IntAsDouble(n)));
-    use qs = Qubit[n - HammingWeightI(n) - T];
-
-    let registerPartition = MappedOverRange(t -> Floor(IntAsDouble(n) / IntAsDouble(2^t)) - 1, 1..T - 1);
-    let pWorkspace = [ps] + Partitioned(registerPartition, qs);
-
-    within {
-        PRounds(pWorkspace);
-    } apply {
-        // U_G
-        GRounds(pWorkspace, gs);
-
-        // U_C
-        CRounds(pWorkspace, gs);
-    }
-}
-
-    /// # Summary
-/// Computes all p[i, j] values in workspace for the look-ahead adder.
-///
-/// The register array `pWorkspace` has T entries, where T = ⌊log₂ n⌋.
-///
-/// The first entry `pWorkspace[0]` is initialized with `P_0` which is
-/// computed before `ComputeCarries` is called.  The other registers are
-/// 0-initialized and will be computed in successive rounds t = 1, ..., T - 1.
-///
-/// In each round t we compute
-///
-/// p[i, j] = p[2ᵗ × m, 2ᵗ × (m + 1)] = p[i, k] ∧ p[k, j]
-///
-/// in `pWorkspace[t][m - 1]` and use that for k = 2ᵗ × m + 2ᵗ⁻¹, p[i, k] and p[k, j]
-/// have already been computed in round t - 1 in `pWorkspace[t - 1][2 * m - 1]` and
-/// `pWorkspace[t - 1][2 * m]`, respectively.
-internal operation PRounds(pWorkspace : Qubit[][]) : Unit is Adj {
-    for ws in Windows(2, pWorkspace) {
-        // note that we are using Rest, since pWorkspace[t - 1][0] is never
-        // accessed in round t.
-        let (current, next) = (Rest(ws[0]), ws[1]);
-
-        for m in IndexRange(next) {
-            ApplyAndAssuming0Target(current[2 * m], current[2 * m + 1], next[m]);
-        }
-    }
-}
-
-    /// # Summary
-/// Computes g[i ∧ (i + 1), i + 1] into gs[i] for the look-ahead adder.
-///
-/// The register gs has n entries initialized to gs[i] = g[i, i + 1].
-///
-/// After successive rounds t = 1, ..., T, the register is updated to
-/// gs[i] = g[i ∧ (i + 1), i + 1], from which we can compute the carries
-/// in the C-rounds.
-internal operation GRounds(pWorkspace : Qubit[][], gs : Qubit[]) : Unit is Adj {
-    let T = Length(pWorkspace);
-    let n = Length(gs);
-
-    for t in 1..T {
-        let length = Floor(IntAsDouble(n) / IntAsDouble(2^t)) - 1;
-        let ps = pWorkspace[t - 1][0..2...];
-
-        for m in 0..length {
-            CCNOT(gs[2^t * m + 2^(t - 1) - 1], ps[m], gs[2^t * m + 2^t - 1]);
-        }
-    }
-}
-
-    /// # Summary
-/// Computes carries into gs for the look-ahead adder.
-internal operation CRounds(pWorkspace : Qubit[][], gs : Qubit[]) : Unit is Adj {
-    let n = Length(gs);
-
-    let start = Floor(Lg(IntAsDouble(2 * n) / 3.0));
-    for t in start..-1..1 {
-        let length = Floor(IntAsDouble(n - 2^(t - 1)) / IntAsDouble(2^t));
-        let ps = pWorkspace[t - 1][1..2...];
-
-        for m in 1..length {
-            CCNOT(gs[2^t * m - 1], ps[m - 1], gs[2^t * m + 2^(t - 1) - 1]);
-        }
-    }
-}
-
-internal operation PhaseGradient(qs : Qubit[]) : Unit is Adj + Ctl {
-    for i in IndexRange(qs) {
-        R1Frac(1, i, qs[i]);
-    }
-}
-
-//
-// Internal operations for comparisons
-//
-
-/// # Summary
-/// Applies `action` to `target` if register `x` is greater or equal to BigInt `c`
-/// (if `invertControl` is false). If `invertControl` is true, the `action`
-/// is applied in the opposite situation.
-internal operation ApplyActionIfGreaterThanOrEqualConstant<'T>(
-    invertControl : Bool,
-    action : 'T => Unit is Adj + Ctl,
-    c : BigInt,
-    x : Qubit[],
-    target : 'T
-) : Unit is Adj + Ctl {
-
-    let bitWidth = Length(x);
-    if c == 0L {
-        if not invertControl {
-            action(target);
-        }
-    } elif c >= (2L^bitWidth) {
-        if invertControl {
-            action(target);
-        }
-    } else {
-        // normalize constant
-        let l = TrailingZeroCountL(c);
-
-        let cNormalized = c >>> l;
-        let xNormalized = x[l...];
-        let bitWidthNormalized = Length(xNormalized);
-
-        // If c == 2L^(bitwidth - 1), then bitWidthNormalized will be 1,
-        // and qs will be empty.  In that case, we do not need to compute
-        // any temporary values, and some optimizations are apply, which
-        // are considered in the remainder.
-        use qs = Qubit[bitWidthNormalized - 1];
-        let cs1 = IsEmpty(qs) ? [] | [Head(xNormalized)] + Most(qs);
-
-        Fact(Length(cs1) == Length(qs), "Arrays should be of the same length.");
-
-        within {
-            for i in 0..Length(cs1) - 1 {
-                let op = cNormalized &&& (1L <<< (i + 1)) != 0L ? ApplyAndAssuming0Target | ApplyOrAssuming0Target;
-                op(cs1[i], xNormalized[i + 1], qs[i]);
-            }
-        } apply {
-            let control = IsEmpty(qs) ? Tail(x) | Tail(qs);
-            within {
-                if invertControl {
-                    X(control);
-                }
-            } apply {
-                Controlled action([control], target);
-            }
-        }
-    }
-}
-
-    /// # Summary
-/// Applies `action` to `target` if the sum of `x` and `y` registers
-/// overflows, i.e. there's a carry out (if `invertControl` is false).
-/// If `invertControl` is true, the `action` is applied when there's no carry out.
-internal operation ApplyActionIfSumOverflows<'T>(
-    action : 'T => Unit is Adj + Ctl,
-    x : Qubit[],
-    y : Qubit[],
-    invertControl : Bool,
-    target : 'T
-) : Unit is Adj + Ctl {
-
-    let n = Length(x);
-    Fact(n >= 1, "Registers must contain at least one qubit.");
-    Fact(Length(y) == n, "Registers must be of the same length.");
-
-    use carries = Qubit[n];
-
-    within {
-        CarryWith1CarryIn(x[0], y[0], carries[0]);
-        for i in 1..n - 1 {
-            CarryForInc(carries[i - 1], x[i], y[i], carries[i]);
-        }
-    } apply {
-        within {
-            if invertControl {
-                X(carries[n - 1]);
-            }
-        } apply {
-            Controlled action([carries[n - 1]], target);
-        }
-    }
-}
-
-    /// # Summary
-/// Computes carry out assuming carry in is 1.
-/// Simplified version that is only applicable for scenarios
-/// where controlled version is the same as non-controlled.
-internal operation CarryWith1CarryIn(
-    x : Qubit,
-    y : Qubit,
-    carryOut : Qubit
-) : Unit is Adj + Ctl {
-
-    body (...) {
-        X(x);
-        X(y);
-        ApplyAndAssuming0Target(x, y, carryOut);
-        X(carryOut);
-    }
-
-    adjoint auto;
-
-    controlled (ctls, ...) {
-        Fact(Length(ctls) <= 1, "Number of control lines must be at most 1");
-        CarryWith1CarryIn(x, y, carryOut);
-    }
-
-    controlled adjoint auto;
-}
-
-    /// # Summary
-/// This wrapper allows operations that support only one control
-/// qubit to be used in a multi-controlled scenarios. It provides
-/// controlled version that collects controls into one qubit
-/// by applying AND chain using auxiliary qubit array.
-internal operation ApplyAsSinglyControlled<'TIn>(
-    op : ('TIn => Unit is Adj + Ctl),
-    input : 'TIn
-) : Unit is Adj + Ctl {
-
-    body (...) {
-        op(input);
-    }
-
-    controlled (ctls, ...) {
-        let n = Length(ctls);
-        if n == 0 {
-            op(input);
-        } elif n == 1 {
-            Controlled op(ctls, input);
-        } else {
-            use aux = Qubit[n - 1];
-            within {
-                LogDepthAndChain(ctls, aux);
-            } apply {
-                Controlled op([Tail(aux)], input);
-            }
-        }
-    }
-}
-
-    /// # Summary
-/// This helper function computes the AND of all control bits in `ctls` into
-/// the last qubit of `tgts`, using the other qubits in `tgts` as helper
-/// qubits for the AND of subsets of control bits.  The operation has a
-/// logarithmic depth of AND gates by aligning them using a balanced binary
-/// tree.
-internal operation LogDepthAndChain(ctls : Qubit[], tgts : Qubit[]) : Unit is Adj {
-    let lc = Length(ctls);
-    let lt = Length(tgts);
-
-    Fact(lc == lt + 1, $"There must be exactly one more control qubit than target qubits (got {lc}, {lt})");
-
-    if lt == 1 {
-        ApplyAndAssuming0Target(ctls[0], ctls[1], tgts[0]);
-    } elif lt == 2 {
-        ApplyAndAssuming0Target(ctls[0], ctls[1], tgts[0]);
-        ApplyAndAssuming0Target(ctls[2], tgts[0], tgts[1]);
-    } else {
-        let left = lc / 2;
-        let right = lc - left;
-
-        let ctlsLeft = ctls[...left - 1];
-        let tgtsLeft = tgts[...left - 2];
-
-        let ctlsRight = ctls[left..left + right - 1];
-        let tgtsRight = tgts[left - 1..left + right - 3];
-
-        LogDepthAndChain(ctlsLeft, tgtsLeft);
-        LogDepthAndChain(ctlsRight, tgtsRight);
-        ApplyAndAssuming0Target(Tail(tgtsLeft), Tail(tgtsRight), Tail(tgts));
-    }
-}
diff --git a/library/std/src/Std/Unstable/StatePreparation.qs b/library/std/src/Std/Unstable/StatePreparation.qs
deleted file mode 100644
index 1d3cb742ef..0000000000
--- a/library/std/src/Std/Unstable/StatePreparation.qs
+++ /dev/null
@@ -1,390 +0,0 @@
-// Copyright (c) Microsoft Corporation.
-// Licensed under the MIT License.
-
-
-open Microsoft.Quantum.Convert;
-open Microsoft.Quantum.Diagnostics;
-import Std.Arrays.*;
-open Microsoft.Quantum.Math;
-import QIR.Intrinsic.*;
-import Std.Unstable.ArithmeticHelpers.*;
-
-
-/// # Summary
-/// Given a set of coefficients and a big-endian quantum register,
-/// prepares a state on that register described by the given coefficients.
-///
-/// # Description
-/// This operation prepares an arbitrary quantum
-/// state |𝜓⟩ with coefficients 𝑎ⱼ from
-/// the n-qubit computational basis state |0...0⟩.
-///
-/// The action of U on the all-zeros state is given by
-/// $$
-/// \begin{align}
-///     U \ket{0\cdots 0} = \ket{\psi} = \frac{\sum_{j=0}^{2^n-1}\alpha_j \ket{j}}{\sqrt{\sum_{j=0}^{2^n-1}|\alpha_j|^2}}.
-/// \end{align}
-/// $$
-///
-/// # Input
-/// ## coefficients
-/// Array of up to 2ⁿ real coefficients. The j-th coefficient
-/// indexes the number state |j⟩ encoded in big-endian format.
-///
-/// ## qubits
-/// Qubit register encoding number states in a big-endian format. This is
-/// expected to be initialized in the computational basis state |0...0⟩.
-///
-/// # Remarks
-/// `coefficients` will be normalized and padded with
-/// elements 𝑎ⱼ = 0.0 if fewer than 2ⁿ are specified.
-///
-/// # Example
-/// The following snippet prepares the quantum state |𝜓⟩=√(1/8)|0⟩+√(7/8)|2⟩=√(1/8)|00⟩+√(7/8)|10⟩
-/// in the qubit register `qubits`.
-/// ```qsharp
-/// let amplitudes = [Sqrt(0.125), 0.0, Sqrt(0.875), 0.0];
-/// use qubits = Qubit[2];
-/// PreparePureStateD(amplitudes, qubits);
-/// ```
-///
-/// # References
-/// - [arXiv:quant-ph/0406176](https://arxiv.org/abs/quant-ph/0406176)
-///   "Synthesis of Quantum Logic Circuits",
-///   Vivek V. Shende, Stephen S. Bullock, Igor L. Markov
-///
-/// # See Also
-/// - Microsoft.Quantum.Unstable.StatePreparation.ApproximatelyPreparePureStateCP
-operation PreparePureStateD(coefficients : Double[], qubits : Qubit[]) : Unit is Adj + Ctl {
-    let coefficientsAsComplexPolar = Mapped(a -> ComplexAsComplexPolar(Complex(a, 0.0)), coefficients);
-    ApproximatelyPreparePureStateCP(0.0, coefficientsAsComplexPolar, qubits);
-}
-
-    /// # Summary
-/// Given a set of coefficients and a big-endian quantum register,
-/// prepares a state on that register described by the given coefficients,
-/// up to a given approximation tolerance.
-///
-/// # Description
-/// This operation prepares an arbitrary quantum
-/// state |𝜓⟩ with complex coefficients rⱼ·𝒆^(𝒊·tⱼ) from
-/// the n-qubit computational basis state |0...0⟩.
-/// In particular, the action of this operation can be simulated by the
-/// a unitary transformation U which acts on the all-zeros state as
-///
-/// $$
-/// \begin{align}
-///     U\ket{0...0}
-///         & = \ket{\psi} \\\\
-///         & = \frac{
-///                 \sum_{j=0}^{2^n-1} r_j e^{i t_j} \ket{j}
-///             }{
-///                 \sqrt{\sum_{j=0}^{2^n-1} |r_j|^2}
-///             }.
-/// \end{align}
-/// $$
-///
-/// # Input
-/// ## tolerance
-/// The approximation tolerance to be used when preparing the given state.
-///
-/// ## coefficients
-/// Array of up to 2ⁿ complex coefficients represented by their
-/// absolute value and phase (rⱼ, tⱼ). The j-th coefficient
-/// indexes the number state |j⟩ encoded in a big-endian format.
-///
-/// ## qubits
-/// Qubit register encoding number states in a big-endian format. This is
-/// expected to be initialized in the computational basis state
-/// |0...0⟩.
-///
-/// # Remarks
-/// `coefficients` will be padded with
-/// elements (rⱼ, tⱼ) = (0.0, 0.0) if fewer than 2ⁿ are
-/// specified.
-///
-/// # References
-/// - [arXiv:quant-ph/0406176](https://arxiv.org/abs/quant-ph/0406176)
-///   "Synthesis of Quantum Logic Circuits",
-///   Vivek V. Shende, Stephen S. Bullock, Igor L. Markov
-operation ApproximatelyPreparePureStateCP(
-    tolerance : Double,
-    coefficients : ComplexPolar[],
-    qubits : Qubit[]
-) : Unit is Adj + Ctl {
-
-    let nQubits = Length(qubits);
-    // pad coefficients at tail length to a power of 2.
-    let coefficientsPadded = Padded(-2^nQubits, ComplexPolar(0.0, 0.0), coefficients);
-    let idxTarget = 0;
-    // Determine what controls to apply
-    let rngControl = nQubits > 1 ? (1..(nQubits - 1)) | (1..0);
-    // Note we use the reversed qubits array to get the endianness ordering that we expect
-    // when corresponding qubit state to state vector index.
-    Adjoint ApproximatelyUnprepareArbitraryState(
-        tolerance,
-        coefficientsPadded,
-        rngControl,
-        idxTarget,
-        Reversed(qubits)
-    );
-}
-
-    /// # Summary
-/// Implementation step of arbitrary state preparation procedure.
-internal operation ApproximatelyUnprepareArbitraryState(
-    tolerance : Double,
-    coefficients : ComplexPolar[],
-    rngControl : Range,
-    idxTarget : Int,
-    register : Qubit[]
-) : Unit is Adj + Ctl {
-
-    // For each 2D block, compute disentangling single-qubit rotation parameters
-    let (disentanglingY, disentanglingZ, newCoefficients) = StatePreparationSBMComputeCoefficients(coefficients);
-    if (AnyOutsideToleranceD(tolerance, disentanglingZ)) {
-        ApproximatelyMultiplexPauli(tolerance, disentanglingZ, PauliZ, register[rngControl], register[idxTarget]);
-
-    }
-    if (AnyOutsideToleranceD(tolerance, disentanglingY)) {
-        ApproximatelyMultiplexPauli(tolerance, disentanglingY, PauliY, register[rngControl], register[idxTarget]);
-    }
-    // target is now in |0> state up to the phase given by arg of newCoefficients.
-
-    // Continue recursion while there are control qubits.
-    if (IsRangeEmpty(rngControl)) {
-        let (abs, arg) = (newCoefficients[0].Magnitude, newCoefficients[0].Argument);
-        if (AbsD(arg) > tolerance) {
-            Exp([PauliI], -1.0 * arg, [register[idxTarget]]);
-        }
-    } elif (Any(c -> AbsComplexPolar(c) > tolerance, newCoefficients)) {
-        // Some coefficients are outside tolerance
-        let newControl = (RangeStart(rngControl) + 1)..RangeStep(rngControl)..RangeEnd(rngControl);
-        let newTarget = RangeStart(rngControl);
-        ApproximatelyUnprepareArbitraryState(tolerance, newCoefficients, newControl, newTarget, register);
-    }
-}
-
-    /// # Summary
-/// Applies a Pauli rotation conditioned on an array of qubits, truncating
-/// small rotation angles according to a given tolerance.
-///
-/// # Description
-/// This applies a multiply controlled unitary operation that performs
-/// rotations by angle $\theta_j$ about single-qubit Pauli operator $P$
-/// when controlled by the $n$-qubit number state $\ket{j}$.
-/// In particular, the action of this operation is represented by the
-/// unitary
-///
-/// $$
-/// \begin{align}
-///     U = \sum^{2^n - 1}_{j=0} \ket{j}\bra{j} \otimes e^{i P \theta_j}.
-/// \end{align}
-/// ##
-///
-/// # Input
-/// ## tolerance
-/// A tolerance below which small coefficients are truncated.
-///
-/// ## coefficients
-/// Array of up to $2^n$ coefficients $\theta_j$. The $j$th coefficient
-/// indexes the number state $\ket{j}$ encoded in little-endian format.
-///
-/// ## pauli
-/// Pauli operator $P$ that determines axis of rotation.
-///
-/// ## control
-/// $n$-qubit control register that encodes number states $\ket{j}$ in
-/// little-endian format.
-///
-/// ## target
-/// Single qubit register that is rotated by $e^{i P \theta_j}$.
-///
-/// # Remarks
-/// `coefficients` will be padded with elements $\theta_j = 0.0$ if
-/// fewer than $2^n$ are specified.
-internal operation ApproximatelyMultiplexPauli(
-    tolerance : Double,
-    coefficients : Double[],
-    pauli : Pauli,
-    control : Qubit[],
-    target : Qubit
-) : Unit is Adj + Ctl {
-
-    if pauli == PauliZ {
-        ApproximatelyMultiplexZ(tolerance, coefficients, control, target);
-    } elif pauli == PauliX {
-        within {
-            H(target);
-        } apply {
-            ApproximatelyMultiplexPauli(tolerance, coefficients, PauliZ, control, target);
-        }
-    } elif pauli == PauliY {
-        within {
-            Adjoint S(target);
-        } apply {
-            ApproximatelyMultiplexPauli(tolerance, coefficients, PauliX, control, target);
-        }
-    } else {
-        fail $"MultiplexPauli failed. Invalid pauli {pauli}.";
-    }
-}
-
-    /// # Summary
-/// Implementation step of arbitrary state preparation procedure.
-internal function StatePreparationSBMComputeCoefficients(
-    coefficients : ComplexPolar[]
-) : (Double[], Double[], ComplexPolar[]) {
-
-    mutable disentanglingZ = [];
-    mutable disentanglingY = [];
-    mutable newCoefficients = [];
-
-    for idxCoeff in 0..2..Length(coefficients) - 1 {
-        let (rt, phi, theta) = BlochSphereCoordinates(coefficients[idxCoeff], coefficients[idxCoeff + 1]);
-        set disentanglingZ += [0.5 * phi];
-        set disentanglingY += [0.5 * theta];
-        set newCoefficients += [rt];
-    }
-
-    return (disentanglingY, disentanglingZ, newCoefficients);
-}
-
-    /// # Summary
-/// Computes the Bloch sphere coordinates for a single-qubit state.
-///
-/// Given two complex numbers $a0, a1$ that represent the qubit state, computes coordinates
-/// on the Bloch sphere such that
-/// $a0 \ket{0} + a1 \ket{1} = r e^{it}(e^{-i \phi /2}\cos{(\theta/2)}\ket{0}+e^{i \phi /2}\sin{(\theta/2)}\ket{1})$.
-///
-/// # Input
-/// ## a0
-/// Complex coefficient of state $\ket{0}$.
-/// ## a1
-/// Complex coefficient of state $\ket{1}$.
-///
-/// # Output
-/// A tuple containing `(ComplexPolar(r, t), phi, theta)`.
-internal function BlochSphereCoordinates(
-    a0 : ComplexPolar,
-    a1 : ComplexPolar
-) : (ComplexPolar, Double, Double) {
-
-    let abs0 = AbsComplexPolar(a0);
-    let abs1 = AbsComplexPolar(a1);
-    let arg0 = ArgComplexPolar(a0);
-    let arg1 = ArgComplexPolar(a1);
-    let r = Sqrt(abs0 * abs0 + abs1 * abs1);
-    let t = 0.5 * (arg0 + arg1);
-    let phi = arg1 - arg0;
-    let theta = 2.0 * ArcTan2(abs1, abs0);
-    return (ComplexPolar(r, t), phi, theta);
-}
-
-    /// # Summary
-/// Applies a Pauli Z rotation conditioned on an array of qubits, truncating
-/// small rotation angles according to a given tolerance.
-///
-/// # Description
-/// This applies the multiply controlled unitary operation that performs
-/// rotations by angle $\theta_j$ about single-qubit Pauli operator $Z$
-/// when controlled by the $n$-qubit number state $\ket{j}$.
-/// In particular, this operation can be represented by the unitary
-///
-/// $$
-/// \begin{align}
-///     U = \sum^{2^n-1}_{j=0} \ket{j}\bra{j} \otimes e^{i Z \theta_j}.
-/// \end{align}
-/// $$
-///
-/// # Input
-/// ## tolerance
-/// A tolerance below which small coefficients are truncated.
-///
-/// ## coefficients
-/// Array of up to $2^n$ coefficients $\theta_j$. The $j$th coefficient
-/// indexes the number state $\ket{j}$ encoded in little-endian format.
-///
-/// ## control
-/// $n$-qubit control register that encodes number states $\ket{j}$ in
-/// little-endian format.
-///
-/// ## target
-/// Single qubit register that is rotated by $e^{i P \theta_j}$.
-///
-/// # Remarks
-/// `coefficients` will be padded with elements $\theta_j = 0.0$ if
-/// fewer than $2^n$ are specified.
-///
-/// # References
-/// - [arXiv:quant-ph/0406176](https://arxiv.org/abs/quant-ph/0406176)
-///   "Synthesis of Quantum Logic Circuits",
-///   Vivek V. Shende, Stephen S. Bullock, Igor L. Markov
-internal operation ApproximatelyMultiplexZ(
-    tolerance : Double,
-    coefficients : Double[],
-    control : Qubit[],
-    target : Qubit
-) : Unit is Adj + Ctl {
-
-    body (...) {
-        // pad coefficients length at tail to a power of 2.
-        let coefficientsPadded = Padded(-2^Length(control), 0.0, coefficients);
-
-        if Length(coefficientsPadded) == 1 {
-            // Termination case
-            if AbsD(coefficientsPadded[0]) > tolerance {
-                Exp([PauliZ], coefficientsPadded[0], [target]);
-            }
-        } else {
-            // Compute new coefficients.
-            let (coefficients0, coefficients1) = MultiplexZCoefficients(coefficientsPadded);
-            ApproximatelyMultiplexZ(tolerance, coefficients0, Most(control), target);
-            if AnyOutsideToleranceD(tolerance, coefficients1) {
-                within {
-                    CNOT(Tail(control), target);
-                } apply {
-                    ApproximatelyMultiplexZ(tolerance, coefficients1, Most(control), target);
-                }
-            }
-        }
-    }
-
-    controlled (controlRegister, ...) {
-        // pad coefficients length to a power of 2.
-        let coefficientsPadded = Padded(2^(Length(control) + 1), 0.0, Padded(-2^Length(control), 0.0, coefficients));
-        let (coefficients0, coefficients1) = MultiplexZCoefficients(coefficientsPadded);
-        ApproximatelyMultiplexZ(tolerance, coefficients0, control, target);
-        if AnyOutsideToleranceD(tolerance, coefficients1) {
-            within {
-                Controlled X(controlRegister, target);
-            } apply {
-                ApproximatelyMultiplexZ(tolerance, coefficients1, control, target);
-            }
-        }
-    }
-}
-
-    /// # Summary
-/// Implementation step of multiply-controlled Z rotations.
-internal function MultiplexZCoefficients(coefficients : Double[]) : (Double[], Double[]) {
-    let newCoefficientsLength = Length(coefficients) / 2;
-    mutable coefficients0 = [];
-    mutable coefficients1 = [];
-
-    for idxCoeff in 0..newCoefficientsLength - 1 {
-        set coefficients0 += [0.5 * (coefficients[idxCoeff] + coefficients[idxCoeff + newCoefficientsLength])];
-        set coefficients1 += [0.5 * (coefficients[idxCoeff] - coefficients[idxCoeff + newCoefficientsLength])];
-    }
-
-    return (coefficients0, coefficients1);
-}
-
-internal function AnyOutsideToleranceD(tolerance : Double, coefficients : Double[]) : Bool {
-    // NOTE: This function is not used as the only recursion termination condition
-    // only to determine if the multiplex step needs to be applied.
-    // For tolerance 0.0 it is always applied due to >= comparison.
-    Any(coefficient -> AbsD(coefficient) >= tolerance, coefficients)
-}
-
-export PreparePureStateD, ApproximatelyPreparePureStateCP;
diff --git a/library/std/src/Std/Unstable/TableLookup.qs b/library/std/src/Std/Unstable/TableLookup.qs
deleted file mode 100644
index 89a2d45da0..0000000000
--- a/library/std/src/Std/Unstable/TableLookup.qs
+++ /dev/null
@@ -1,281 +0,0 @@
-// Copyright (c) Microsoft Corporation.
-// Licensed under the MIT License.
-
-
-import Std.Arrays.*;
-open Microsoft.Quantum.Convert;
-open Microsoft.Quantum.Diagnostics;
-open Microsoft.Quantum.Math;
-open Microsoft.Quantum.ResourceEstimation;
-open Microsoft.Quantum.Unstable.Arithmetic;
-import Std.Unstable.ArithmeticHelpers.*;
-
-import QIR.Intrinsic.*;
-
-
-/// # Summary
-/// Performs table lookup using a SELECT network
-///
-/// # Description
-/// Assuming a zero-initialized `target` register, this operation will
-/// initialize it with the bitstrings in `data` at indices according to the
-/// computational values of the `address` register.
-///
-/// # Input
-/// ## data
-/// The classical table lookup data which is prepared in `target` with
-/// respect to the state in `address`. The length of data must be less than
-/// 2ⁿ, where 𝑛 is the length of `address`. Each entry in data must have
-/// the same length that must be equal to the length of `target`.
-/// ## address
-/// Address register
-/// ## target
-/// Zero-initialized target register
-///
-/// # Remarks
-/// The implementation of the SELECT network is based on unary encoding as
-/// presented in [1].  The recursive implementation of that algorithm is
-/// presented in [3].  The adjoint variant is optimized using a
-/// measurement-based unlookup operation [3]. The controlled adjoint variant
-/// is not optimized using this technique.
-///
-/// # References
-/// 1. [arXiv:1805.03662](https://arxiv.org/abs/1805.03662)
-///    "Encoding Electronic Spectra in Quantum Circuits with Linear T
-///    Complexity"
-/// 2. [arXiv:1905.07682](https://arxiv.org/abs/1905.07682)
-///    "Windowed arithmetic"
-/// 3. [arXiv:2211.01133](https://arxiv.org/abs/2211.01133)
-///    "Space-time optimized table lookup"
-operation Select(
-    data : Bool[][],
-    address : Qubit[],
-    target : Qubit[]
-) : Unit is Adj + Ctl {
-    body (...) {
-        let (N, n) = DimensionsForSelect(data, address);
-
-        if N == 1 {
-            // base case
-            WriteMemoryContents(Head(data), target);
-        } else {
-            let (most, tail) = MostAndTail(address[...n - 1]);
-            let parts = Partitioned([2^(n - 1)], data);
-
-            within {
-                X(tail);
-            } apply {
-                SinglyControlledSelect(tail, parts[0], most, target);
-            }
-
-            SinglyControlledSelect(tail, parts[1], most, target);
-        }
-    }
-    adjoint (...) {
-        Unlookup(Select, data, address, target);
-    }
-
-    controlled (ctls, ...) {
-        let numCtls = Length(ctls);
-
-        if numCtls == 0 {
-            Select(data, address, target);
-        } elif numCtls == 1 {
-            SinglyControlledSelect(ctls[0], data, address, target);
-        } else {
-            use andChainTarget = Qubit();
-            let andChain = MakeAndChain(ctls, andChainTarget);
-            use helper = Qubit[andChain.NGarbageQubits];
-
-            within {
-                andChain.Apply(helper);
-            } apply {
-                SinglyControlledSelect(andChainTarget, data, address, target);
-            }
-        }
-    }
-
-    controlled adjoint (ctls, ...) {
-        Controlled Select(ctls, (data, address, target));
-    }
-}
-
-internal operation SinglyControlledSelect(
-    ctl : Qubit,
-    data : Bool[][],
-    address : Qubit[],
-    target : Qubit[]
-) : Unit {
-    let (N, n) = DimensionsForSelect(data, address);
-
-    if BeginEstimateCaching("Microsoft.Quantum.Unstable.TableLookup.SinglyControlledSelect", N) {
-        if N == 1 {
-            // base case
-            Controlled WriteMemoryContents([ctl], (Head(data), target));
-        } else {
-            use helper = Qubit();
-
-            let (most, tail) = MostAndTail(address[...n - 1]);
-            let parts = Partitioned([2^(n - 1)], data);
-
-            within {
-                X(tail);
-            } apply {
-                ApplyAndAssuming0Target(ctl, tail, helper);
-            }
-
-            SinglyControlledSelect(helper, parts[0], most, target);
-
-            CNOT(ctl, helper);
-
-            SinglyControlledSelect(helper, parts[1], most, target);
-
-            Adjoint ApplyAndAssuming0Target(ctl, tail, helper);
-        }
-
-        EndEstimateCaching();
-    }
-}
-
-internal function DimensionsForSelect(
-    data : Bool[][],
-    address : Qubit[]
-) : (Int, Int) {
-    let N = Length(data);
-    Fact(N > 0, "data cannot be empty");
-
-    let n = Ceiling(Lg(IntAsDouble(N)));
-    Fact(
-        Length(address) >= n,
-        $"address register is too small, requires at least {n} qubits"
-    );
-
-    return (N, n);
-}
-
-internal operation WriteMemoryContents(
-    value : Bool[],
-    target : Qubit[]
-) : Unit is Adj + Ctl {
-    Fact(
-        Length(value) == Length(target),
-        "number of data bits must equal number of target qubits"
-    );
-
-    ApplyPauliFromBitString(PauliX, true, value, target);
-}
-
-    /// # References
-/// - [arXiv:1905.07682](https://arxiv.org/abs/1905.07682)
-///   "Windowed arithmetic"
-internal operation Unlookup(
-    lookup : (Bool[][], Qubit[], Qubit[]) => Unit,
-    data : Bool[][],
-    select : Qubit[],
-    target : Qubit[]
-) : Unit {
-    let numBits = Length(target);
-    let numAddressBits = Length(select);
-
-    let l = MinI(Floor(Lg(IntAsDouble(numBits))), numAddressBits - 1);
-    Fact(
-        l < numAddressBits,
-        $"l = {l} must be smaller than {numAddressBits}"
-    );
-
-    let res = Mapped(r -> r == One, ForEach(MResetX, target));
-
-    let dataFixup = Chunks(2^l, Padded(-2^numAddressBits, false, Mapped(MustBeFixed(res, _), data)));
-
-    let numAddressBitsFixup = numAddressBits - l;
-
-    let selectParts = Partitioned([l], select);
-    let targetFixup = target[...2^l - 1];
-
-    within {
-        EncodeUnary(selectParts[0], targetFixup);
-        ApplyToEachA(H, targetFixup);
-    } apply {
-        lookup(dataFixup, selectParts[1], targetFixup);
-    }
-}
-
-// Checks whether specific bit string `data` must be fixed for a given
-// measurement result `result`.
-//
-// Returns true if the number of indices for which both result and data are
-// `true` is odd.
-internal function MustBeFixed(result : Bool[], data : Bool[]) : Bool {
-    mutable state = false;
-    for i in IndexRange(result) {
-        set state = state != (result[i] and data[i]);
-    }
-    state
-}
-
-// Computes unary encoding of value in `input` into `target`
-//
-// Assumptions:
-//    - `target` is zero-initialized
-//    - length of `input` is n
-//    - length of `target` is 2^n
-internal operation EncodeUnary(
-    input : Qubit[],
-    target : Qubit[]
-) : Unit is Adj {
-    Fact(
-        Length(target) == 2^Length(input),
-        $"target register should be of length {2^Length(input)}, but is {Length(target)}"
-    );
-
-    X(Head(target));
-
-    for i in IndexRange(input) {
-        if i == 0 {
-            CNOT(input[i], target[1]);
-            CNOT(target[1], target[0]);
-        } else {
-            // targets are the first and second 2^i qubits of the target register
-            let split = Partitioned([2^i, 2^i], target);
-            for j in IndexRange(split[0]) {
-                ApplyAndAssuming0Target(input[i], split[0][j], split[1][j]);
-                CNOT(split[1][j], split[0][j]);
-            }
-        }
-    }
-
-}
-
-internal newtype AndChain = (
-    NGarbageQubits : Int,
-    Apply : Qubit[] => Unit is Adj
-);
-
-internal function MakeAndChain(ctls : Qubit[], target : Qubit) : AndChain {
-    AndChain(
-        MaxI(Length(ctls) - 2, 0),
-        helper => AndChainOperation(ctls, helper, target)
-    )
-}
-
-internal operation AndChainOperation(ctls : Qubit[], helper : Qubit[], target : Qubit) : Unit is Adj {
-    let n = Length(ctls);
-
-    Fact(Length(helper) == MaxI(n - 2, 0), "Invalid number of helper qubits");
-
-    if n == 0 {
-        X(target);
-    } elif n == 1 {
-        CNOT(ctls[0], target);
-    } else {
-        let ctls1 = ctls[0..0] + helper;
-        let ctls2 = ctls[1...];
-        let tgts = helper + [target];
-
-        for idx in IndexRange(tgts) {
-            ApplyAndAssuming0Target(ctls1[idx], ctls2[idx], tgts[idx]);
-        }
-    }
-}
-
-export Select;
diff --git a/library/std/src/canon.qs b/library/std/src/canon.qs
deleted file mode 100644
index a8dd54f66b..0000000000
--- a/library/std/src/canon.qs
+++ /dev/null
@@ -1,589 +0,0 @@
-// Copyright (c) Microsoft Corporation.
-// Licensed under the MIT License.
-
-namespace Microsoft.Quantum.Canon {
-    import QIR.Intrinsic.*;
-    open Microsoft.Quantum.Intrinsic;
-    open Microsoft.Quantum.Diagnostics;
-    open Microsoft.Quantum.Math;
-
-    /// # Summary
-    /// Applies an operation to each element in a register.
-    ///
-    /// # Input
-    /// ## singleElementOperation
-    /// Operation to apply to each element.
-    /// ## register
-    /// Array of elements on which to apply the given operation.
-    ///
-    /// # Type Parameters
-    /// ## 'T
-    /// The target on which the operation acts.
-    ///
-    /// # Example
-    /// Prepare a three-qubit |+⟩ state:
-    /// ```qsharp
-    /// use register = Qubit[3];
-    /// ApplyToEach(H, register);
-    /// ```
-    operation ApplyToEach<'T>(singleElementOperation : ('T => Unit), register : 'T[]) : Unit {
-        for item in register {
-            singleElementOperation(item);
-        }
-    }
-
-    /// # Summary
-    /// Applies an operation to each element in a register.
-    /// The modifier `A` indicates that the single-element operation is adjointable.
-    ///
-    /// # Input
-    /// ## singleElementOperation
-    /// Operation to apply to each element.
-    /// ## register
-    /// Array of elements on which to apply the given operation.
-    ///
-    /// # Type Parameters
-    /// ## 'T
-    /// The target on which the operation acts.
-    ///
-    /// # Example
-    /// Prepare a three-qubit |+⟩ state:
-    /// ```qsharp
-    /// use register = Qubit[3];
-    /// ApplyToEach(H, register);
-    /// ```
-    ///
-    /// # See Also
-    /// - Microsoft.Quantum.Canon.ApplyToEach
-    operation ApplyToEachA<'T>(singleElementOperation : ('T => Unit is Adj), register : 'T[]) : Unit is Adj {
-        for item in register {
-            singleElementOperation(item);
-        }
-    }
-
-    /// # Summary
-    /// Applies an operation to each element in a register.
-    /// The modifier `C` indicates that the single-element operation is controllable.
-    ///
-    /// # Input
-    /// ## singleElementOperation
-    /// Operation to apply to each element.
-    /// ## register
-    /// Array of elements on which to apply the given operation.
-    ///
-    /// # Type Parameters
-    /// ## 'T
-    /// The target on which the operation acts.
-    ///
-    /// # Example
-    /// Prepare a three-qubit |+⟩ state:
-    /// ```qsharp
-    /// use register = Qubit[3];
-    /// ApplyToEach(H, register);
-    /// ```
-    ///
-    /// # See Also
-    /// - Microsoft.Quantum.Canon.ApplyToEach
-    operation ApplyToEachC<'T>(singleElementOperation : ('T => Unit is Ctl), register : 'T[]) : Unit is Ctl {
-        for item in register {
-            singleElementOperation(item);
-        }
-    }
-
-    /// # Summary
-    /// Applies an operation to each element in a register.
-    /// The modifier `CA` indicates that the single-element operation is controllable and adjointable.
-    ///
-    /// # Input
-    /// ## singleElementOperation
-    /// Operation to apply to each element.
-    /// ## register
-    /// Array of elements on which to apply the given operation.
-    ///
-    /// # Type Parameters
-    /// ## 'T
-    /// The target on which the operation acts.
-    ///
-    /// # Example
-    /// Prepare a three-qubit |+⟩ state:
-    /// ```qsharp
-    /// use register = Qubit[3];
-    /// ApplyToEach(H, register);
-    /// ```
-    ///
-    /// # See Also
-    /// - Microsoft.Quantum.Canon.ApplyToEach
-    operation ApplyToEachCA<'T>(singleElementOperation : ('T => Unit is Adj + Ctl), register : 'T[]) : Unit is Adj + Ctl {
-        for item in register {
-            singleElementOperation(item);
-        }
-    }
-
-    /// # Summary
-    /// Applies the controlled-X (CX) gate to a pair of qubits.
-    ///
-    /// # Input
-    /// ## control
-    /// Control qubit for the CX gate.
-    /// ## target
-    /// Target qubit for the CX gate.
-    ///
-    /// # Remarks
-    /// This operation can be simulated by the unitary matrix
-    /// $$
-    /// \begin{align}
-    ///     \left(\begin{matrix}
-    ///         1 & 0 & 0 & 0 \\\\
-    ///         0 & 1 & 0 & 0 \\\\
-    ///         0 & 0 & 0 & 1 \\\\
-    ///         0 & 0 & 1 & 0
-    ///      \end{matrix}\right)
-    /// \end{align},
-    /// $$
-    /// where rows and columns are organized as in the quantum concepts guide.
-    ///
-    /// Equivalent to:
-    /// ```qsharp
-    /// Controlled X([control], target);
-    /// ```
-    /// and to:
-    /// ```qsharp
-    /// CNOT(control, target);
-    /// ```
-    operation CX(control : Qubit, target : Qubit) : Unit is Adj + Ctl {
-        body ... {
-            __quantum__qis__cx__body(control, target);
-        }
-        controlled (ctls, ...) {
-            Controlled X(ctls + [control], target);
-        }
-        adjoint self;
-    }
-
-    /// # Summary
-    /// Applies the controlled-Y (CY) gate to a pair of qubits.
-    ///
-    /// # Input
-    /// ## control
-    /// Control qubit for the CY gate.
-    /// ## target
-    /// Target qubit for the CY gate.
-    ///
-    /// # Remarks
-    /// This operation can be simulated by the unitary matrix
-    /// $$
-    /// \begin{align}
-    ///     \left(\begin{matrix}
-    ///         1 & 0 & 0 & 0 \\\\
-    ///         0 & 1 & 0 & 0 \\\\
-    ///         0 & 0 & 0 & -i \\\\
-    ///         0 & 0 & i & 0
-    ///      \end{matrix}\right)
-    /// \end{align},
-    /// $$
-    /// where rows and columns are organized as in the quantum concepts guide.
-    ///
-    /// Equivalent to:
-    /// ```qsharp
-    /// Controlled Y([control], target);
-    /// ```
-    operation CY(control : Qubit, target : Qubit) : Unit is Adj + Ctl {
-        body ... {
-            __quantum__qis__cy__body(control, target);
-        }
-        controlled (ctls, ...) {
-            Controlled Y(ctls + [control], target);
-        }
-        adjoint self;
-    }
-
-    /// # Summary
-    /// Applies the controlled-Z (CZ) gate to a pair of qubits.
-    ///
-    /// # Input
-    /// ## control
-    /// Control qubit for the CZ gate.
-    /// ## target
-    /// Target qubit for the CZ gate.
-    ///
-    /// # Remarks
-    /// This operation can be simulated by the unitary matrix
-    /// $$
-    /// \begin{align}
-    ///     \left(\begin{matrix}
-    ///         1 & 0 & 0 & 0 \\\\
-    ///         0 & 1 & 0 & 0 \\\\
-    ///         0 & 0 & 1 & 0 \\\\
-    ///         0 & 0 & 0 & -1
-    ///     \end{matrix}\right)
-    /// \end{align},
-    /// $$
-    /// where rows and columns are organized as in the quantum concepts guide.
-    ///
-    /// Equivalent to:
-    /// ```qsharp
-    /// Controlled Z([control], target);
-    /// ```
-    operation CZ(control : Qubit, target : Qubit) : Unit is Adj + Ctl {
-        body ... {
-            __quantum__qis__cz__body(control, target);
-        }
-        controlled (ctls, ...) {
-            Controlled Z(ctls + [control], target);
-        }
-        adjoint self;
-    }
-
-    /// Given a pair, returns its first element.
-    function Fst<'T, 'U>(pair : ('T, 'U)) : 'T {
-        let (fst, _) = pair;
-        return fst;
-    }
-
-    /// Given a pair, returns its second element.
-    function Snd<'T, 'U>(pair : ('T, 'U)) : 'U {
-        let (_, snd) = pair;
-        return snd;
-    }
-
-    /// # Summary
-    /// Computes the parity of a register of qubits in-place.
-    ///
-    /// # Input
-    /// ## qubits
-    /// Array of qubits whose parity is to be computed and stored.
-    ///
-    /// # Remarks
-    /// This operation transforms the state of its input as
-    /// $$
-    /// \begin{align}
-    ///     \ket{q_0} \ket{q_1} \cdots \ket{q_{n - 1}} & \mapsto
-    ///     \ket{q_0} \ket{q_0 \oplus q_1} \ket{q_0 \oplus q_1 \oplus q_2} \cdots
-    ///         \ket{q_0 \oplus \cdots \oplus q_{n - 1}}.
-    /// \end{align}
-    /// $$
-    operation ApplyCNOTChain(qubits : Qubit[]) : Unit is Adj + Ctl {
-        for i in 0..Length(qubits) - 2 {
-            CNOT(qubits[i], qubits[i + 1]);
-        }
-    }
-
-    /// # Summary
-    /// Given a single-qubit Pauli operator, applies the corresponding operation
-    /// to a single qubit.
-    ///
-    /// # Input
-    /// ## pauli
-    /// The Pauli operator to be applied.
-    /// ## target
-    /// The qubit to which `pauli` is to be applied as an operation.
-    ///
-    /// # Example
-    /// The following are equivalent:
-    /// ```qsharp
-    /// ApplyP(PauliX, q);
-    /// ```
-    /// and
-    /// ```qsharp
-    /// X(q);
-    /// ```
-    operation ApplyP(pauli : Pauli, target : Qubit) : Unit is Adj + Ctl {
-        if pauli == PauliX { X(target); } elif pauli == PauliY { Y(target); } elif pauli == PauliZ { Z(target); }
-    }
-
-    /// # Summary
-    /// Given a multi-qubit Pauli operator, applies the corresponding operation
-    /// to a quantum register.
-    ///
-    /// # Input
-    /// ## pauli
-    /// A multi-qubit Pauli operator represented as an array of single-qubit Pauli operators.
-    /// ## target
-    /// Register to apply the given Pauli operation on.
-    ///
-    /// # Example
-    /// The following are equivalent:
-    /// ```qsharp
-    /// ApplyPauli([PauliY, PauliZ, PauliX], target);
-    /// ```
-    /// and
-    /// ```qsharp
-    /// Y(target[0]);
-    /// Z(target[1]);
-    /// X(target[2]);
-    /// ```
-    operation ApplyPauli(pauli : Pauli[], target : Qubit[]) : Unit is Adj + Ctl {
-        Fact(Length(pauli) == Length(target), "`pauli` and `target` must be of the same length.");
-        for i in 0..Length(pauli) - 1 {
-            ApplyP(pauli[i], target[i]);
-        }
-    }
-
-    /// # Summary
-    /// Applies a Pauli operator on each qubit in an array if the corresponding
-    /// bit of a Boolean array matches a given input.
-    ///
-    /// # Input
-    /// ## pauli
-    /// Pauli operator to apply to `qubits[idx]` where `bitApply == bits[idx]`
-    /// ## bitApply
-    /// apply Pauli if bit is this value
-    /// ## bits
-    /// Boolean register specifying which corresponding qubit in `qubits` should be operated on
-    /// ## qubits
-    /// Quantum register on which to selectively apply the specified Pauli operator
-    ///
-    /// # Remarks
-    /// The Boolean array and the quantum register must be of equal length.
-    ///
-    /// # Example
-    /// The following applies an X operation on qubits 0 and 2, and a Z operation on qubits 1 and 3.
-    /// ```qsharp
-    /// use qubits = Qubit[4];
-    /// let bits = [true, false, true, false];
-    /// // Apply when index in `bits` is `true`.
-    /// ApplyPauliFromBitString(PauliX, true, bits, qubits);
-    /// // Apply when index in `bits` is `false`.
-    /// ApplyPauliFromBitString(PauliZ, false, bits, qubits);
-    /// ```
-    operation ApplyPauliFromBitString(pauli : Pauli, bitApply : Bool, bits : Bool[], qubits : Qubit[]) : Unit is Adj + Ctl {
-        let nBits = Length(bits);
-        Fact(nBits == Length(qubits), "Number of bits must be equal to number of qubits.");
-        for i in 0..nBits - 1 {
-            if bits[i] == bitApply {
-                ApplyP(pauli, qubits[i]);
-            }
-        }
-    }
-
-    /// # Summary
-    /// Applies a Pauli operator on each qubit in an array if the corresponding
-    /// bit of a Little-endian integer matches a given input.
-    ///
-    /// # Input
-    /// ## pauli
-    /// Pauli operator to apply to `qubits[idx]` when bit of numberState
-    /// in idx position is the same as bitApply.
-    /// ## bitApply
-    /// apply Pauli if bit is this value
-    /// ## numberState
-    /// Little-endian integer specifying which corresponding qubit in `qubits` should be operated on
-    /// ## qubits
-    /// Quantum register on which to selectively apply the specified Pauli operator
-    ///
-    /// # Example
-    /// The following applies an X operation on qubits 0 and 2, and a Z operation on qubits 1 and 3.
-    /// ```qsharp
-    /// use qubits = Qubit[4];
-    /// let n = 5;
-    /// // Apply when index in `bits` is `true`.
-    /// ApplyPauliFromBitString(PauliX, true, n, qubits);
-    /// // Apply when index in `bits` is `false`.
-    /// ApplyPauliFromBitString(PauliZ, false, n, qubits);
-    /// ```
-    operation ApplyPauliFromInt(
-        pauli : Pauli,
-        bitApply : Bool,
-        numberState : Int,
-        qubits : Qubit[]
-    ) : Unit is Adj + Ctl {
-
-        let length = Length(qubits);
-        Fact(numberState >= 0, "number must be non-negative");
-        Fact(BitSizeI(numberState) <= length, "Bit size of numberState must not exceed qubits length");
-
-        for i in 0..length - 1 {
-            // If we assume loop unrolling, 2^i will be optimized to a constant.
-            if ((numberState &&& (1 <<< i)) != 0) == bitApply {
-                ApplyP(pauli, qubits[i]);
-            }
-        }
-    }
-
-    /// # Summary
-    /// Applies a unitary operation on the target if the control
-    /// register state corresponds to a specified nonnegative integer.
-    ///
-    /// # Input
-    /// ## numberState
-    /// A nonnegative integer on which the operation `oracle` should be
-    /// controlled.
-    /// ## oracle
-    /// A unitary operation to be controlled.
-    /// ## target
-    /// A target on which to apply `oracle`.
-    /// ## controlRegister
-    /// A quantum register that controls application of `oracle`.
-    ///
-    /// # Remarks
-    /// The value of `numberState` is interpreted using a little-endian encoding.
-    ///
-    /// `numberState` must be at most $2^\texttt{Length(controlRegister)} - 1$.
-    /// For example, `numberState = 537` means that `oracle`
-    /// is applied if and only if `controlRegister` is in the state $\ket{537}$.
-    operation ApplyControlledOnInt<'T>(
-        numberState : Int,
-        oracle : ('T => Unit is Adj + Ctl),
-        controlRegister : Qubit[],
-        target : 'T
-    ) : Unit is Adj + Ctl {
-
-        within {
-            ApplyPauliFromInt(PauliX, false, numberState, controlRegister);
-        } apply {
-            Controlled oracle(controlRegister, target);
-        }
-    }
-
-    /// # Summary
-    /// Applies `oracle` on `target` when `controlRegister`
-    /// is in the state specified by `bits`.
-    ///
-    /// # Description
-    /// Applies a unitary operation `oracle` on the `target`, controlled
-    /// on a state specified by a given bit mask `bits`.
-    /// The bit at `bits[i]` corresponds to qubit at `controlRegister[i]`.
-    /// The pattern given by `bits` may be shorter than `controlRegister`,
-    /// in which case additional control qubits are ignored (that is, neither
-    /// controlled on |0⟩ nor |1⟩).
-    /// If `bits` is longer than `controlRegister`, an error is raised.
-    ///
-    /// # Input
-    /// ## bits
-    /// The bit string to control the given unitary operation on.
-    /// ## oracle
-    /// The unitary operation to be applied on the target.
-    /// ## target
-    /// The target to be passed to `oracle` as an input.
-    /// ## controlRegister
-    /// A quantum register that controls application of `oracle`.
-    ///
-    /// # Example
-    /// ```qsharp
-    /// // When bits = [1,0,0] oracle is applied if and only if controlRegister
-    /// // is in the state |100⟩.
-    /// use t = Qubit();
-    /// use c = Qubit[3];
-    /// X(c[0]);
-    /// ApplyControlledOnBitString([true, false, false], X, c, t);
-    /// Message($"{M(t)}"); // Prints `One` since oracle `X` was applied.
-    /// ```
-    operation ApplyControlledOnBitString<'T>(
-        bits : Bool[],
-        oracle : ('T => Unit is Adj + Ctl),
-        controlRegister : Qubit[],
-        target : 'T
-    ) : Unit is Adj + Ctl {
-
-        // The control register must have enough bits to implement the requested control.
-        Fact(Length(bits) <= Length(controlRegister), "Control register shorter than control pattern.");
-
-        // Use a subregister of the controlled register when
-        // bits is shorter than controlRegister.
-        let controlSubregister = controlRegister[...Length(bits) - 1];
-        within {
-            ApplyPauliFromBitString(PauliX, false, bits, controlSubregister);
-        } apply {
-            Controlled oracle(controlSubregister, target);
-        }
-    }
-
-    /// # Summary
-    /// Applies the rotations of Quantum Fourier Transform (QFT) to a little-endian quantum register.
-    ///
-    /// # Description
-    /// Applies the rotations of QFT to a little-endian register `qs` of length n
-    /// containing |x₁⟩⊗|x₂⟩⊗…⊗|xₙ⟩. The qs[0] initially contains the
-    /// least significant bit xₙ. The state of qs[0] becomes
-    /// (|0⟩+𝑒^(2π𝑖[0.xₙ])|1⟩)/sqrt(2) after the operation.
-    ///
-    /// # Input
-    /// ## qs
-    /// Quantum register in a little-endian format to which the rotations are applied.
-    ///
-    /// # Remarks
-    /// Note that this operation applies only the rotations part of the QFT.
-    /// To complete the transform, you need to reverse the order of qubits after this operation,
-    /// for example, using the operation `SwapReverseRegister`.
-    ///
-    /// # Reference
-    ///  - [Quantum Fourier transform](https://en.wikipedia.org/wiki/Quantum_Fourier_transform)
-    operation ApplyQFT(qs : Qubit[]) : Unit is Adj + Ctl {
-        let length = Length(qs);
-        Fact(length >= 1, "ApplyQFT: Length(qs) must be at least 1.");
-        for i in length - 1..-1..0 {
-            H(qs[i]);
-            for j in 0..i - 1 {
-                Controlled R1Frac([qs[i]], (1, j + 1, qs[i - j - 1]));
-            }
-        }
-    }
-
-    /// # Summary
-    /// Uses SWAP gates to reverse the order of the qubits in a register.
-    ///
-    /// # Input
-    /// ## register
-    /// The qubits order of which should be reversed using SWAP gates
-    operation SwapReverseRegister(register : Qubit[]) : Unit is Adj + Ctl {
-        let length = Length(register);
-        for i in 0..length / 2 - 1 {
-            SWAP(register[i], register[(length - i) - 1]);
-        }
-    }
-
-    /// # Summary
-    /// Applies a bitwise-XOR operation between a classical integer and an
-    /// integer represented by a register of qubits.
-    ///
-    /// # Description
-    /// Applies `X` operations to qubits in a little-endian register based on
-    /// 1 bits in an integer.
-    ///
-    /// Let us denote `value` by a and let y be an unsigned integer encoded in `target`,
-    /// then `ApplyXorInPlace` performs an operation given by the following map:
-    /// |y⟩ ↦ |y ⊕ a⟩, where ⊕ is the bitwise exclusive OR operator.
-    operation ApplyXorInPlace(value : Int, target : Qubit[]) : Unit is Adj + Ctl {
-        body (...) {
-            Fact(value >= 0, "`value` must be non-negative.");
-            mutable runningValue = value;
-            for q in target {
-                if (runningValue &&& 1) != 0 {
-                    X(q);
-                }
-                set runningValue >>>= 1;
-            }
-            Fact(runningValue == 0, "value is too large");
-        }
-        adjoint self;
-    }
-
-    /// # Summary
-    /// Applies a bitwise-XOR operation between a classical integer and an
-    /// integer represented by a register of qubits.
-    ///
-    /// # Description
-    /// Applies `X` operations to qubits in a little-endian register based on
-    /// 1 bits in an integer.
-    ///
-    /// Let us denote `value` by a and let y be an unsigned integer encoded in `target`,
-    /// then `ApplyXorInPlace` performs an operation given by the following map:
-    /// |y⟩ ↦ |y ⊕ a⟩, where ⊕ is the bitwise exclusive OR operator.
-    operation ApplyXorInPlaceL(value : BigInt, target : Qubit[]) : Unit is Adj + Ctl {
-        body (...) {
-            Fact(value >= 0L, "`value` must be non-negative.");
-            mutable runningValue = value;
-            for q in target {
-                if (runningValue &&& 1L) != 0L {
-                    X(q);
-                }
-                set runningValue >>>= 1;
-            }
-            Fact(runningValue == 0L, "`value` is too large.");
-        }
-        adjoint self;
-    }
-
-    export ApplyToEach, ApplyToEachA, ApplyToEachC, ApplyToEachCA, CX, CY, CZ, Fst, Snd, ApplyCNOTChain, ApplyP, ApplyPauli, ApplyPauliFromBitString, ApplyPauliFromInt, ApplyControlledOnInt, ApplyControlledOnBitString, ApplyQFT, SwapReverseRegister, ApplyXorInPlace, ApplyXorInPlaceL;
-
-}
diff --git a/library/std/src/legacy_api.qs b/library/std/src/legacy_api.qs
index b4de6e7ee1..d351c868de 100644
--- a/library/std/src/legacy_api.qs
+++ b/library/std/src/legacy_api.qs
@@ -5,11 +5,7 @@
 // This file re-exports the standard library under the name `Std`, which will be the preferred standard library API going forward.
 
 namespace Microsoft.Quantum {
-    export Std.Arrays, Std.Convert, Std.Diagnostics, Std.Logical, Std.Math, Std.Measurement, Std.Intrinsic, Std.Random, Std.ResourceEstimation;
-}
-
-namespace Microsoft.Quantum.Unstable {
-    export Std.Unstable.Arithmetic, Std.Unstable.TableLookup, Std.Unstable.StatePreparation;
+    export Std.Arrays, Std.Convert, Std.Diagnostics, Std.Logical, Std.Math, Std.Measurement, Std.Intrinsic, Std.Random, Std.ResourceEstimation, Std.Canon;
 }
 
 namespace Std {
diff --git a/library/std/src/unstable_arithmetic.qs b/library/std/src/unstable_arithmetic.qs
new file mode 100644
index 0000000000..22ae84ae13
--- /dev/null
+++ b/library/std/src/unstable_arithmetic.qs
@@ -0,0 +1,701 @@
+// Copyright (c) Microsoft Corporation.
+// Licensed under the MIT License.
+
+namespace Microsoft.Quantum.Unstable.Arithmetic {
+    open Microsoft.Quantum.Arrays;
+    open Microsoft.Quantum.Diagnostics;
+    open Microsoft.Quantum.Math;
+    open Microsoft.Quantum.Convert;
+
+    /// # Summary
+    /// This applies the in-place majority operation to 3 qubits.
+    ///
+    /// # Description
+    /// Assuming the state of the input qubits are |x⟩, |y⟩ and |z⟩, then
+    /// this operation performs the following transformation:
+    /// |x⟩|y⟩|z⟩ ↦ |x ⊕ z⟩|y ⊕ z⟩MAJ(x, y, z).
+    ///
+    /// # Input
+    /// ## x
+    /// The first input qubit.
+    /// ## y
+    /// The second input qubit.
+    /// ## z
+    /// A qubit onto which the majority function will be applied.
+    operation MAJ(x : Qubit, y : Qubit, z : Qubit) : Unit is Adj + Ctl {
+        CNOT(z, y);
+        CNOT(z, x);
+        CCNOT(y, x, z);
+    }
+
+    /// # Summary
+    /// Reflects a quantum register about a given classical integer.
+    ///
+    /// # Description
+    /// Given a quantum register initially in the state ∑ᵢ(αᵢ|i⟩),
+    /// where each |i⟩ is a basis state representing an integer i,
+    /// reflects the state of the register about the basis state |j⟩
+    /// for a given integer j: ∑ᵢ(-1)^(δᵢⱼ)(αᵢ|i⟩)
+    /// This operation is implemented in-place, without explicit allocation of
+    /// additional auxiliary qubits.
+    ///
+    /// # Input
+    /// ## index
+    /// The classical integer j indexing the basis state about which to reflect.
+    /// ## reg
+    /// Little-endian quantum register to reflect.
+    operation ReflectAboutInteger(index : Int, reg : Qubit[]) : Unit is Adj + Ctl {
+        within {
+            // Evaluation optimization for case index == 0
+            if index == 0 {
+                ApplyToEachA(X, reg);
+            } else {
+                // We want to reduce to the problem of reflecting about the all-ones
+                // state. To do that, we apply our reflection within an application
+                // of X instructions that flip all the zeros in our index.
+                ApplyPauliFromInt(PauliX, false, index, reg);
+            }
+        } apply {
+            Controlled ApplyAsSinglyControlled(Most(reg), (Z, Tail(reg)));
+        }
+    }
+
+    //
+    // Add, Increment      |   Operation    | Description
+    // ____________________|________________|_______________________________________________________________
+    // y += 5              | IncByI, IncByL | Increment LE register in-place by integer
+    // y += x              | IncByLE        | Increment LE register in-place by LE register
+    // z = x + 5 (z was 0) |                | Add integer to LE register creating result out-of-place
+    // z = x + y (z was 0) | AddLE          | Add two LE register creating result out-of-place
+    // z += x + 5          |                | Increment LE register by the sum of integer and LE register
+    // z += x + y          |                | Increment LE register by the sum of two LE registers
+    //
+    // IncByLE implementations:
+    //     RippleCarryTTKIncByLE (default)
+    //     RippleCarryCGIncByLE
+    //     FourierTDIncByLE
+    //     via IncByLEUsingAddLE and any out-of-place addition
+    // IncByI implementations:
+    //     via IncByIUsingIncByLE and any in-place LE adder
+    // IncByL implementations:
+    //     via IncByLUsingIncByLE and any in-place LE adder
+    // AddLE implementations:
+    //     RippleCarryCGAddLE (default)
+    //     LookAheadDKRSAddLE
+    //
+
+    /// # Summary
+    /// Increments a little-endian register ys by an integer number c
+    ///
+    /// # Description
+    /// Computes ys += c modulo 2ⁿ, where ys is a little-endian register,
+    /// Length(ys) = n > 0, c is a Int number, 0 ≤ c < 2ⁿ.
+    /// NOTE: Use IncByIUsingIncByLE directly if the choice of implementation
+    /// is important.
+    operation IncByI(c : Int, ys : Qubit[]) : Unit is Adj + Ctl {
+        IncByIUsingIncByLE(RippleCarryTTKIncByLE, c, ys);
+    }
+
+    /// # Summary
+    /// Increments a little-endian register ys by a BigInt number c
+    ///
+    /// # Description
+    /// Computes ys += c modulo 2ⁿ, where ys is a little-endian register,
+    /// Length(ys) = n > 0, c is a BigInt number, 0 ≤ c < 2ⁿ.
+    /// NOTE: Use IncByLUsingIncByLE directly if the choice of implementation
+    /// is important.
+    operation IncByL(c : BigInt, ys : Qubit[]) : Unit is Adj + Ctl {
+        IncByLUsingIncByLE(RippleCarryTTKIncByLE, c, ys);
+    }
+
+    /// # Summary
+    /// Increments a little-endian register ys by a little-endian register xs
+    ///
+    /// # Description
+    /// Computes ys += xs modulo 2ⁿ, where xs and ys are little-endian registers,
+    /// and Length(xs) ≤ Length(ys) = n.
+    /// NOTE: Use operations like RippleCarryCGIncByLE directly if
+    /// the choice of implementation is important.
+    operation IncByLE(xs : Qubit[], ys : Qubit[]) : Unit is Adj + Ctl {
+        RippleCarryTTKIncByLE(xs, ys);
+    }
+
+    /// # Summary
+    /// Sets a zero-initialized little-endian register zs to the sum of
+    /// little-endian registers xs and ys
+    ///
+    /// # Description
+    /// Computes zs := xs + ys modulo 2ⁿ, where xs, ys, and zs are little-endian registers,
+    /// Length(xs) = Length(ys) ≤ Length(zs) = n, assuming zs is 0-initialized.
+    /// NOTE: Use operations like RippleCarryCGAddLE directly if
+    /// the choice of implementation is important.
+    operation AddLE(xs : Qubit[], ys : Qubit[], zs : Qubit[]) : Unit is Adj {
+        RippleCarryCGAddLE(xs, ys, zs);
+    }
+
+    /// # Summary
+    /// Reversible, in-place ripple-carry addition of two integers.
+    ///
+    /// # Description
+    /// Computes ys += xs modulo 2ⁿ, where xs and ys are little-endian registers,
+    /// and Length(xs) ≤ Length(ys) = n.
+    /// This operation uses the ripple-carry algorithm.
+    /// Note that if Length(ys) >= Length(xs)+2, xs is padded with 0-initialized
+    /// qubits to match ys's length. The operation doesn't use any auxiliary
+    /// qubits otherwise.
+    ///
+    /// # References
+    /// - [arXiv:0910.2530](https://arxiv.org/abs/0910.2530)
+    ///   "Quantum Addition Circuits and Unbounded Fan-Out",
+    ///   Yasuhiro Takahashi, Seiichiro Tani, Noboru Kunihiro
+    operation RippleCarryTTKIncByLE(xs : Qubit[], ys : Qubit[]) : Unit is Adj + Ctl {
+        let xsLen = Length(xs);
+        let ysLen = Length(ys);
+
+        Fact(ysLen >= xsLen, "Register `ys` must be longer than register `xs`.");
+        Fact(xsLen >= 1, "Registers `xs` and `ys` must contain at least one qubit.");
+
+        if xsLen == ysLen {
+            if xsLen > 1 {
+                within {
+                    ApplyOuterTTKAdder(xs, ys);
+                } apply {
+                    ApplyInnerTTKAdderNoCarry(xs, ys);
+                }
+            }
+            CNOT(xs[0], ys[0]);
+        } elif xsLen + 1 == ysLen {
+            if xsLen > 1 {
+                CNOT(xs[xsLen - 1], ys[ysLen - 1]);
+                within {
+                    ApplyOuterTTKAdder(xs, ys);
+                } apply {
+                    ApplyInnerTTKAdderWithCarry(xs, ys);
+                }
+            } else {
+                CCNOT(xs[0], ys[0], ys[1]);
+            }
+            CNOT(xs[0], ys[0]);
+        } elif xsLen + 2 <= ysLen {
+            // Pad xs so that its length is one qubit shorter than ys.
+            use padding = Qubit[ysLen - xsLen - 1];
+            RippleCarryTTKIncByLE(xs + padding, ys);
+        }
+    }
+
+    /// # Summary
+    /// Increments a little-endian register ys by a little-endian register xs
+    /// using the ripple-carry algorithm.
+    ///
+    /// # Description
+    /// Computes ys += xs modulo 2ⁿ, where xs and ys are little-endian registers,
+    /// and Length(xs) ≤ Length(ys) = n.
+    /// Note that if Length(xs) != Length(ys), xs is padded with 0-initialized
+    /// qubits to match ys's length.
+    /// This operation uses the ripple-carry algorithm.
+    ///
+    /// # Reference
+    /// - [arXiv:1709.06648](https://arxiv.org/pdf/1709.06648.pdf)
+    ///   "Halving the cost of quantum addition", Craig Gidney.
+    operation RippleCarryCGIncByLE(xs : Qubit[], ys : Qubit[]) : Unit is Adj + Ctl {
+        let xsLen = Length(xs);
+        let ysLen = Length(ys);
+
+        Fact(ysLen >= xsLen, "Register `ys` must be longer than register `xs`.");
+        Fact(xsLen >= 1, "Registers `xs` and `ys` must contain at least one qubit.");
+
+        if ysLen - xsLen >= 2 {
+            // Pad xs so that its length is one qubit shorter than ys.
+            use padding = Qubit[ysLen - xsLen - 1];
+            RippleCarryCGIncByLE(xs + padding, ys);
+        } elif xsLen == 1 {
+            if ysLen == 1 {
+                CNOT(xs[0], ys[0]);
+            } elif ysLen == 2 {
+                HalfAdderForInc(xs[0], ys[0], ys[1]);
+            }
+        } else {
+            use carries = Qubit[xsLen];
+            within {
+                ApplyAndAssuming0Target(xs[0], ys[0], carries[0]);
+            } apply {
+                for i in 1..xsLen - 2 {
+                    CarryForInc(carries[i - 1], xs[i], ys[i], carries[i]);
+                }
+                if xsLen == ysLen {
+                    within {
+                        CNOT(carries[xsLen - 2], xs[xsLen - 1]);
+                    } apply {
+                        CNOT(xs[xsLen - 1], ys[xsLen - 1]);
+                    }
+                } else {
+                    FullAdderForInc(carries[xsLen - 2], xs[xsLen - 1], ys[xsLen - 1], ys[xsLen]);
+                }
+                for i in xsLen - 2..-1..1 {
+                    UncarryForInc(carries[i - 1], xs[i], ys[i], carries[i]);
+                }
+            }
+            CNOT(xs[0], ys[0]);
+        }
+    }
+
+    /// # Summary
+    /// Sets a zero-initialized little-endian register zs to the sum of
+    /// little-endian registers xs and ys using the ripple-carry algorithm.
+    ///
+    /// # Description
+    /// Computes zs := xs + ys + zs[0] modulo 2ⁿ, where xs, ys, and zs are
+    /// little-endian registers, Length(xs) = Length(ys) ≤ Length(zs) = n,
+    /// assuming zs is 0-initialized, except for maybe zs[0], which can be
+    // in |0> or |1> state and can be used as carry-in.
+    /// This operation uses the ripple-carry algorithm.
+    /// NOTE: `zs[Length(xs)]` can be used as carry-out, if `zs` is longer than `xs`.
+    ///
+    /// # Reference
+    /// - [arXiv:1709.06648](https://arxiv.org/pdf/1709.06648.pdf)
+    ///   "Halving the cost of quantum addition", Craig Gidney.
+    operation RippleCarryCGAddLE(xs : Qubit[], ys : Qubit[], zs : Qubit[]) : Unit is Adj {
+        let xsLen = Length(xs);
+        let zsLen = Length(zs);
+        Fact(Length(ys) == xsLen, "Registers `xs` and `ys` must be of same length.");
+        Fact(zsLen >= xsLen, "Register `zs` must be no shorter than register `xs`.");
+
+        // Since zs is zero-initialized, its bits at indexes higher than
+        // xsLen remain unused as there will be no carry into them.
+        let top = MinI(zsLen - 2, xsLen - 1);
+        for k in 0..top {
+            FullAdder(zs[k], xs[k], ys[k], zs[k + 1]);
+        }
+
+        if xsLen > 0 and xsLen == zsLen {
+            CNOT(Tail(xs), Tail(zs));
+            CNOT(Tail(ys), Tail(zs));
+        }
+    }
+
+    /// # Summary
+    /// Sets a zero-initialized little-endian register zs to the sum of
+    /// little-endian registers xs and ys using the carry-lookahead algorithm.
+    ///
+    /// # Description
+    /// Computes zs := xs + ys + zs[0] modulo 2ⁿ, where xs, ys, and zs are
+    /// little-endian registers, Length(xs) = Length(ys) ≤ Length(zs) = n,
+    /// assuming zs is 0-initialized, except for maybe zs[0], which can be
+    /// in |0> or |1> state and can be used as carry-in.
+    /// NOTE: `zs[Length(xs)]` can be used as carry-out, if `zs` is longer than `xs`.
+    /// This operation uses the carry-lookahead algorithm.
+    ///
+    /// # Reference
+    /// - [arXiv:quant-ph/0406142](https://arxiv.org/abs/quant-ph/0406142)
+    ///   "A logarithmic-depth quantum carry-lookahead adder",
+    ///   Thomas G. Draper, Samuel A. Kutin, Eric M. Rains, Krysta M. Svore
+    operation LookAheadDKRSAddLE(xs : Qubit[], ys : Qubit[], zs : Qubit[]) : Unit is Adj {
+        let xsLen = Length(xs);
+        let zsLen = Length(zs);
+        Fact(Length(ys) == xsLen, "Registers `xs` and `ys` must be of same length.");
+        Fact(zsLen >= xsLen, "Register `zs` must be no shorter than register `xs`.");
+
+        if zsLen > xsLen {
+            // with carry-out
+            // compute initial generate values
+            for k in 0..xsLen - 1 {
+                ApplyAndAssuming0Target(xs[k], ys[k], zs[k + 1]);
+            }
+
+            within {
+                // compute initial propagate values
+                for i in IndexRange(xs) {
+                    CNOT(xs[i], ys[i]);
+                }
+            } apply {
+                if xsLen > 1 {
+                    ComputeCarries(Rest(ys), zs[1..xsLen]);
+                }
+
+                // compute sum into carries
+                for k in 0..xsLen - 1 {
+                    CNOT(ys[k], zs[k]);
+                }
+            }
+        } else {
+            // xsLen == zsLen, so without carry-out
+            LookAheadDKRSAddLE(Most(xs), Most(ys), zs);
+            CNOT(Tail(xs), Tail(zs));
+            CNOT(Tail(ys), Tail(zs));
+        }
+    }
+
+    /// # Summary
+    /// Increments a little-endian register ys by a little-endian register xs
+    /// using Quantum Fourier Transform.
+    ///
+    /// # Description
+    /// Computes ys += xs modulo 2ⁿ, where xs and ys are little-endian registers,
+    /// and Length(xs) = Length(ys) = n.
+    /// This operation uses Quantum Fourier Transform.
+    ///
+    /// # Reference
+    /// - [arXiv:quant-ph/0008033](https://arxiv.org/abs/quant-ph/0008033)
+    ///   "Addition on a Quantum Computer", Thomas G. Draper
+    operation FourierTDIncByLE(xs : Qubit[], ys : Qubit[]) : Unit is Adj + Ctl {
+        within {
+            ApplyQFT(ys);
+        } apply {
+            for i in IndexRange(xs) {
+                Controlled PhaseGradient([xs[i]], ys[i...]);
+            }
+        }
+    }
+
+    /// # Summary
+    /// Increments a little-endian register ys by a BigInt number c
+    /// using provided adder.
+    ///
+    /// # Description
+    /// Computes ys += c modulo 2ⁿ, where ys is a little-endian register
+    /// Length(ys) = n > 0, c is a BigInt number, 0 ≤ c < 2ⁿ.
+    operation IncByLUsingIncByLE(
+        adder : (Qubit[], Qubit[]) => Unit is Adj + Ctl,
+        c : BigInt,
+        ys : Qubit[]
+    ) : Unit is Adj + Ctl {
+
+        let ysLen = Length(ys);
+        Fact(ysLen > 0, "Length of `ys` must be at least 1.");
+        Fact(c >= 0L, "Constant `c` must be non-negative.");
+        Fact(c < 2L^ysLen, "Constant `c` must be smaller than 2^Length(ys).");
+
+        if c != 0L {
+            // If c has j trailing zeros, then the j least significant
+            // bits of y won't be affected by the addition and can
+            // therefore be ignored by applying the addition only to
+            // the other qubits and shifting c accordingly.
+            let j = TrailingZeroCountL(c);
+            use x = Qubit[ysLen - j];
+            within {
+                ApplyXorInPlaceL(c >>> j, x);
+            } apply {
+                adder(x, ys[j...]);
+            }
+        }
+    }
+
+    /// # Summary
+    /// Increments a little-endian register ys by an Int number c
+    /// using provided adder.
+    ///
+    /// # Description
+    /// Computes ys += c modulo 2ⁿ, where ys is a little-endian register
+    /// Length(ys) = n > 0, c is an Int number, 0 ≤ c < 2ⁿ.
+    operation IncByIUsingIncByLE(
+        adder : (Qubit[], Qubit[]) => Unit is Adj + Ctl,
+        c : Int,
+        ys : Qubit[]
+    ) : Unit is Adj + Ctl {
+
+        let ysLen = Length(ys);
+        Fact(ysLen > 0, "Length of `ys` must be at least 1.");
+        Fact(c >= 0, "Constant `c` must be non-negative.");
+        Fact(c < 2^ysLen, "Constant `c` must be smaller than 2^Length(ys).");
+
+        if c != 0 {
+            // If c has j trailing zeros than the j least significant
+            // bits of y won't be affected by the addition and can
+            // therefore be ignored by applying the addition only to
+            // the other qubits and shifting c accordingly.
+            let j = TrailingZeroCountI(c);
+            use x = Qubit[ysLen - j];
+            within {
+                ApplyXorInPlace(c >>> j, x);
+            } apply {
+                adder(x, ys[j...]);
+            }
+        }
+    }
+
+    /// # Summary
+    /// Generic operation to turn two out-place adders into one in-place adder
+    ///
+    /// # Description
+    /// This implementation allows to specify two distinct adders for forward
+    /// and backward direction.  The forward adder is always applied in its
+    /// body variant, whereas the backward adder is always applied in its adjoint
+    /// variant.  Therefore, it's possible to, for example, use the ripple-carry
+    /// out-of-place adder in backwards direction to require no T gates.
+    ///
+    /// The controlled variant is also optimized in a way that everything but
+    /// the adders is controlled,
+    ///
+    /// # Reference
+    /// - [arXiv:2012.01624](https://arxiv.org/abs/2012.01624)
+    ///   "Quantum block lookahead adders and the wait for magic states",
+    ///   Craig Gidney.
+    operation IncByLEUsingAddLE(
+        forwardAdder : (Qubit[], Qubit[], Qubit[]) => Unit is Adj,
+        backwardAdder : (Qubit[], Qubit[], Qubit[]) => Unit is Adj,
+        xs : Qubit[],
+        ys : Qubit[]
+    ) : Unit is Adj + Ctl {
+
+        body (...) {
+            let n = Length(xs);
+
+            Fact(Length(ys) == n, "Registers xs and ys must be of same length");
+
+            use qs = Qubit[n];
+
+            forwardAdder(xs, ys, qs);
+            for i in IndexRange(ys) {
+                SWAP(ys[i], qs[i]);
+            }
+            ApplyToEachA(X, qs);
+            within {
+                ApplyToEachA(X, ys);
+            } apply {
+                Adjoint backwardAdder(xs, ys, qs);
+            }
+        }
+        adjoint (...) {
+            let n = Length(xs);
+
+            Fact(Length(ys) == n, "Registers xs and ys must be of same length");
+
+            use qs = Qubit[n];
+
+            within {
+                ApplyToEachA(X, ys);
+            } apply {
+                forwardAdder(xs, ys, qs);
+            }
+            ApplyToEachA(X, qs);
+            for i in IndexRange(ys) {
+                SWAP(ys[i], qs[i]);
+            }
+            Adjoint backwardAdder(xs, ys, qs);
+        }
+        controlled (ctls, ...) {
+            // When we control everything except the adders, the adders will
+            // cancel themselves.
+            let n = Length(xs);
+
+            Fact(Length(ys) == n, "Registers xs and ys must be of same length");
+
+            use qs = Qubit[n];
+
+            forwardAdder(xs, ys, qs);
+            for i in IndexRange(ys) {
+                Controlled SWAP(ctls, (ys[i], qs[i]))
+            }
+            ApplyToEachA(tgt => Controlled X(ctls, tgt), qs);
+            within {
+                ApplyToEachA(tgt => Controlled X(ctls, tgt), ys);
+            } apply {
+                Adjoint backwardAdder(xs, ys, qs);
+            }
+        }
+        controlled adjoint (ctls, ...) {
+            // When we control everything except the adders, the adders will
+            // cancel themselves.
+            let n = Length(xs);
+
+            Fact(Length(ys) == n, "Registers xs and ys must be of same length");
+
+            use qs = Qubit[n];
+
+            within {
+                ApplyToEachA(tgt => Controlled X(ctls, tgt), ys);
+            } apply {
+                forwardAdder(xs, ys, qs);
+            }
+            ApplyToEachA(tgt => Controlled X(ctls, tgt), qs);
+            for i in IndexRange(ys) {
+                Controlled SWAP(ctls, (ys[i], qs[i]))
+            }
+            Adjoint backwardAdder(xs, ys, qs);
+        }
+    }
+
+    //
+    // Comparisons
+    //
+    // Compare BigInt and qubit register in a little-endian format and apply action
+    //   if c < x { action(target) }  | ApplyIfLessL
+    //   if c <= x { action(target) } | ApplyIfLessOrEqualL
+    //   if c == x { action(target) } | ApplyIfEqualL
+    //   if c >= x { action(target) } | ApplyIfGreaterOrEqualL
+    //   if c > x { action(target) }  | ApplyIfGreaterL
+    //
+    // Compare two qubit registers in a little-endian format and apply action
+    //   if x < y { action(target) }  | ApplyIfLessLE
+    //   if x <= y { action(target) } | ApplyIfLessOrEqualLE
+    //   if x == y { action(target) } | ApplyIfEqualLE
+    //   if x >= y { action(target) } | ApplyIfGreaterOrEqualLE
+    //   if x > y { action(target) }  | ApplyIfGreaterLE
+    //
+
+    /// # Summary
+    /// Computes `if (c < x) { action(target) }`, that is, applies `action` to `target`
+    /// if a BigInt value `c` is less than the little-endian qubit register `x`
+    operation ApplyIfLessL<'T>(
+        action : 'T => Unit is Adj + Ctl,
+        c : BigInt,
+        x : Qubit[],
+        target : 'T
+    ) : Unit is Adj + Ctl {
+
+        ApplyActionIfGreaterThanOrEqualConstant(false, action, c + 1L, x, target);
+    }
+
+    /// # Summary
+    /// Computes `if (c <= x) { action(target) }`, that is, applies `action` to `target`
+    /// if a BigInt value `c` is less or equal to the little-endian qubit register `x`
+    operation ApplyIfLessOrEqualL<'T>(
+        action : 'T => Unit is Adj + Ctl,
+        c : BigInt,
+        x : Qubit[],
+        target : 'T
+    ) : Unit is Adj + Ctl {
+
+        ApplyActionIfGreaterThanOrEqualConstant(false, action, c, x, target);
+    }
+
+    /// # Summary
+    /// Computes `if (c == x) { action(target) }`, that is, applies `action` to `target`
+    /// if a BigInt value `c` is equal to the little-endian qubit register `x`
+    operation ApplyIfEqualL<'T>(
+        action : 'T => Unit is Adj + Ctl,
+        c : BigInt,
+        xs : Qubit[],
+        target : 'T
+    ) : Unit is Adj + Ctl {
+
+        let cBitSize = BitSizeL(c);
+        let xLen = Length(xs);
+        if (cBitSize <= xLen) {
+            let bits = BigIntAsBoolArray(c, Length(xs));
+            within {
+                ApplyPauliFromBitString(PauliX, false, bits, xs);
+            } apply {
+                Controlled ApplyAsSinglyControlled(xs, (a => action(a), target));
+            }
+        }
+    }
+
+    /// # Summary
+    /// Computes `if (c >= x) { action(target) }`, that is, applies `action` to `target`
+    /// if a BigInt value `c` is greater or equal to the little-endian qubit register `x`
+    operation ApplyIfGreaterOrEqualL<'T>(
+        action : 'T => Unit is Adj + Ctl,
+        c : BigInt,
+        x : Qubit[],
+        target : 'T
+    ) : Unit is Adj + Ctl {
+
+        ApplyActionIfGreaterThanOrEqualConstant(true, action, c + 1L, x, target);
+    }
+
+    /// # Summary
+    /// Computes `if (c > x) { action(target) }`, that is, applies `action` to `target`
+    /// if a BigInt value `c` is greater than the little-endian qubit register `x`
+    operation ApplyIfGreaterL<'T>(
+        action : 'T => Unit is Adj + Ctl,
+        c : BigInt,
+        x : Qubit[],
+        target : 'T
+    ) : Unit is Adj + Ctl {
+
+        ApplyActionIfGreaterThanOrEqualConstant(true, action, c, x, target);
+    }
+
+    /// # Summary
+    /// Computes `if x < y { action(target) }`, that is, applies `action` to `target`
+    /// if register `x` is less than the register `y`.
+    /// Both qubit registers should be in a little-endian format.
+    operation ApplyIfLessLE<'T>(
+        action : 'T => Unit is Adj + Ctl,
+        x : Qubit[],
+        y : Qubit[],
+        target : 'T
+    ) : Unit is Adj + Ctl {
+
+        ApplyIfGreaterLE(action, y, x, target);
+    }
+
+    /// # Summary
+    /// Computes `if x <= y { action(target) }`, that is, applies `action` to `target`
+    /// if register `x` is less or equal to the register `y`.
+    /// Both qubit registers should be in a little-endian format.
+    operation ApplyIfLessOrEqualLE<'T>(
+        action : 'T => Unit is Adj + Ctl,
+        x : Qubit[],
+        y : Qubit[],
+        target : 'T
+    ) : Unit is Adj + Ctl {
+
+        Fact(Length(x) > 0, "Bitwidth must be at least 1");
+        within {
+            ApplyToEachA(X, x);
+        } apply {
+            // control is not inverted
+            ApplyActionIfSumOverflows(action, x, y, false, target);
+        }
+    }
+
+    /// # Summary
+    /// Computes `if x == y { action(target) }`, that is, applies `action` to `target`
+    /// if register `x` is equal to the register `y`.
+    /// Both qubit registers should be in a little-endian format.
+    operation ApplyIfEqualLE<'T>(
+        action : 'T => Unit is Adj + Ctl,
+        x : Qubit[],
+        y : Qubit[],
+        target : 'T
+    ) : Unit is Adj + Ctl {
+
+        Fact(Length(x) == Length(y), "x and y must be of same length");
+        within {
+            for i in IndexRange(x) {
+                CNOT(x[i], y[i]);
+                X(y[i]);
+            }
+        } apply {
+            Controlled ApplyAsSinglyControlled(y, (a => action(a), target))
+        }
+    }
+
+    /// # Summary
+    /// Computes `if x >= y { action(target) }`, that is, applies `action` to `target`
+    /// if register `x` is greater or equal to the register `y`.
+    /// Both qubit registers should be in a little-endian format.
+    operation ApplyIfGreaterOrEqualLE<'T>(
+        action : 'T => Unit is Adj + Ctl,
+        x : Qubit[],
+        y : Qubit[],
+        target : 'T
+    ) : Unit is Adj + Ctl {
+
+        ApplyIfLessOrEqualLE(action, y, x, target);
+    }
+
+    /// # Summary
+    /// Computes `if x > y { action(target) }`, that is, applies `action` to `target`
+    /// if register `x` is greater than the register `y`.
+    /// Both qubit registers should be in a little-endian format.
+    operation ApplyIfGreaterLE<'T>(
+        action : 'T => Unit is Adj + Ctl,
+        x : Qubit[],
+        y : Qubit[],
+        target : 'T
+    ) : Unit is Adj + Ctl {
+
+        Fact(Length(x) > 0, "Bitwidth must be at least 1");
+        within {
+            ApplyToEachA(X, x);
+        } apply {
+            // control is inverted
+            ApplyActionIfSumOverflows(action, x, y, true, target);
+        }
+    }
+
+    export AddLE, ApplyIfEqualLE, ApplyIfEqualL, ApplyIfGreaterLE, ApplyIfGreaterL, ApplyIfGreaterOrEqualLE, ApplyIfGreaterOrEqualL, ApplyIfLessLE, ApplyIfLessL, ApplyIfLessOrEqualLE, ApplyIfLessOrEqualL, IncByI, IncByIUsingIncByLE, IncByL, IncByLUsingIncByLE, IncByLE, IncByLEUsingAddLE, LookAheadDKRSAddLE, MAJ, ReflectAboutInteger, RippleCarryCGAddLE, RippleCarryCGIncByLE, RippleCarryTTKIncByLE, FourierTDIncByLE;
+}
diff --git a/library/std/src/unstable_arithmetic_internal.qs b/library/std/src/unstable_arithmetic_internal.qs
new file mode 100644
index 0000000000..5125859fa8
--- /dev/null
+++ b/library/std/src/unstable_arithmetic_internal.qs
@@ -0,0 +1,585 @@
+// Copyright (c) Microsoft Corporation.
+// Licensed under the MIT License.
+
+namespace Microsoft.Quantum.Unstable.Arithmetic {
+    open Microsoft.Quantum.Diagnostics;
+    open Microsoft.Quantum.Arrays;
+    open Microsoft.Quantum.Math;
+    open Microsoft.Quantum.Convert;
+
+    /// # Summary
+    /// Implements the outer operation for RippleCarryTTKIncByLE to conjugate
+    /// the inner operation to construct the full adder. Only Length(xs)
+    /// qubits are processed.
+    ///
+    /// # Input
+    /// ## xs
+    /// Qubit register in a little-endian format containing the first summand
+    /// input to RippleCarryTTKIncByLE.
+    /// ## ys
+    /// Qubit register in a little-endian format containing the second summand
+    /// input to RippleCarryTTKIncByLE.
+    ///
+    /// # References
+    /// - Yasuhiro Takahashi, Seiichiro Tani, Noboru Kunihiro: "Quantum
+    ///   Addition Circuits and Unbounded Fan-Out", Quantum Information and
+    ///   Computation, Vol. 10, 2010.
+    ///   https://arxiv.org/abs/0910.2530
+    internal operation ApplyOuterTTKAdder(xs : Qubit[], ys : Qubit[]) : Unit is Adj + Ctl {
+        Fact(Length(xs) <= Length(ys), "Input register ys must be at lease as long as xs.");
+        for i in 1..Length(xs) - 1 {
+            CNOT(xs[i], ys[i]);
+        }
+        for i in Length(xs) - 2..-1..1 {
+            CNOT(xs[i], xs[i + 1]);
+        }
+    }
+
+    /// # Summary
+    /// Implements the inner addition function for the operation
+    /// RippleCarryTTKIncByLE. This is the inner operation that is conjugated
+    /// with the outer operation to construct the full adder.
+    ///
+    /// # Input
+    /// ## xs
+    /// Qubit register in a little-endian format containing the first summand
+    /// input to RippleCarryTTKIncByLE.
+    /// ## ys
+    /// Qubit register in a little-endian format containing the second summand
+    /// input to RippleCarryTTKIncByLE.
+    ///
+    /// # References
+    /// - Yasuhiro Takahashi, Seiichiro Tani, Noboru Kunihiro: "Quantum
+    ///   Addition Circuits and Unbounded Fan-Out", Quantum Information and
+    ///   Computation, Vol. 10, 2010.
+    ///   https://arxiv.org/abs/0910.2530
+    ///
+    /// # Remarks
+    /// The specified controlled operation makes use of symmetry and mutual
+    /// cancellation of operations to improve on the default implementation
+    /// that adds a control to every operation.
+    internal operation ApplyInnerTTKAdderNoCarry(xs : Qubit[], ys : Qubit[]) : Unit is Adj + Ctl {
+        body (...) {
+            (Controlled ApplyInnerTTKAdderNoCarry)([], (xs, ys));
+        }
+        controlled (controls, ...) {
+            Fact(Length(xs) == Length(ys), "Input registers must have the same number of qubits.");
+
+            for idx in 0..Length(xs) - 2 {
+                CCNOT(xs[idx], ys[idx], xs[idx + 1]);
+            }
+            for idx in Length(xs) - 1..-1..1 {
+                Controlled CNOT(controls, (xs[idx], ys[idx]));
+                CCNOT(xs[idx - 1], ys[idx - 1], xs[idx]);
+            }
+        }
+    }
+
+    /// # Summary
+    /// Implements the inner addition function for the operation
+    /// RippleCarryTTKIncByLE. This is the inner operation that is conjugated
+    /// with the outer operation to construct the full adder.
+    ///
+    /// # Input
+    /// ## xs
+    /// Qubit register in a little-endian format containing the first summand
+    /// input to RippleCarryTTKIncByLE.
+    /// ## ys
+    /// Qubit register in a little-endian format containing the second summand
+    /// input to RippleCarryTTKIncByLE.
+    ///
+    /// # References
+    /// - Yasuhiro Takahashi, Seiichiro Tani, Noboru Kunihiro: "Quantum
+    ///   Addition Circuits and Unbounded Fan-Out", Quantum Information and
+    ///   Computation, Vol. 10, 2010.
+    ///   https://arxiv.org/abs/0910.2530
+    ///
+    /// # Remarks
+    /// The specified controlled operation makes use of symmetry and mutual
+    /// cancellation of operations to improve on the default implementation
+    /// that adds a control to every operation.
+    internal operation ApplyInnerTTKAdderWithCarry(xs : Qubit[], ys : Qubit[]) : Unit is Adj + Ctl {
+        body (...) {
+            (Controlled ApplyInnerTTKAdderWithCarry)([], (xs, ys));
+        }
+        controlled (controls, ...) {
+            Fact(Length(xs) + 1 == Length(ys), "ys must be one qubit longer then xs.");
+            Fact(Length(xs) > 0, "Array should not be empty.");
+
+
+            let nQubits = Length(xs);
+            for idx in 0..nQubits - 2 {
+                CCNOT(xs[idx], ys[idx], xs[idx + 1]);
+            }
+            (Controlled CCNOT)(controls, (xs[nQubits - 1], ys[nQubits - 1], ys[nQubits]));
+            for idx in nQubits - 1..-1..1 {
+                Controlled CNOT(controls, (xs[idx], ys[idx]));
+                CCNOT(xs[idx - 1], ys[idx - 1], xs[idx]);
+            }
+        }
+    }
+
+    /// # Summary
+    /// Implements Half-adder. Adds qubit x to qubit y and sets carryOut appropriately
+    internal operation HalfAdderForInc(x : Qubit, y : Qubit, carryOut : Qubit) : Unit is Adj + Ctl {
+        body (...) {
+            CCNOT(x, y, carryOut);
+            CNOT(x, y);
+        }
+        adjoint auto;
+
+        controlled (ctls, ...) {
+            Fact(Length(ctls) == 1, "HalfAdderForInc should be controlled by exactly one control qubit.");
+
+            let ctl = ctls[0];
+            use helper = Qubit();
+
+            within {
+                ApplyAndAssuming0Target(x, y, helper);
+            } apply {
+                ApplyAndAssuming0Target(ctl, helper, carryOut);
+            }
+            CCNOT(ctl, x, y);
+        }
+        controlled adjoint auto;
+    }
+
+    /// # Summary
+    /// Implements Full-adder. Adds qubit carryIn and x to qubit y and sets carryOut appropriately.
+    internal operation FullAdderForInc(carryIn : Qubit, x : Qubit, y : Qubit, carryOut : Qubit) : Unit is Adj + Ctl {
+        body (...) {
+            // TODO: cannot use `Carry` operation here
+            CNOT(carryIn, x);
+            CNOT(carryIn, y);
+            CCNOT(x, y, carryOut);
+            CNOT(carryIn, carryOut);
+            CNOT(carryIn, x);
+            CNOT(x, y);
+        }
+        adjoint auto;
+
+        controlled (ctls, ...) {
+            Fact(Length(ctls) == 1, "FullAdderForInc should be controlled by exactly one control qubit.");
+
+            let ctl = ctls[0];
+            use helper = Qubit();
+
+            CarryForInc(carryIn, x, y, helper);
+            CCNOT(ctl, helper, carryOut);
+            Controlled UncarryForInc(ctls, (carryIn, x, y, helper));
+        }
+        controlled adjoint auto;
+    }
+
+    // Computes carryOut := carryIn + x + y
+    internal operation FullAdder(carryIn : Qubit, x : Qubit, y : Qubit, carryOut : Qubit) : Unit is Adj {
+        CNOT(x, y);
+        CNOT(x, carryIn);
+        ApplyAndAssuming0Target(y, carryIn, carryOut);
+        CNOT(x, y);
+        CNOT(x, carryOut);
+        CNOT(y, carryIn);
+    }
+
+    /// # Summary
+    /// Computes carry bit for a full adder.
+    internal operation CarryForInc(carryIn : Qubit, x : Qubit, y : Qubit, carryOut : Qubit) : Unit is Adj + Ctl {
+        body (...) {
+            CNOT(carryIn, x);
+            CNOT(carryIn, y);
+            ApplyAndAssuming0Target(x, y, carryOut);
+            CNOT(carryIn, carryOut);
+        }
+        adjoint auto;
+        controlled (ctls, ...) {
+            // This CarryForInc is intended to be used only in an in-place
+            // ripple-carry implementation. Only such particular use case allows
+            // for this simple implementation where controlled version
+            // is the same as uncontrolled body.
+            CarryForInc(carryIn, x, y, carryOut);
+        }
+        controlled adjoint auto;
+    }
+
+    /// # Summary
+    /// Uncomputes carry bit for a full adder.
+    internal operation UncarryForInc(carryIn : Qubit, x : Qubit, y : Qubit, carryOut : Qubit) : Unit is Adj + Ctl {
+        body (...) {
+            CNOT(carryIn, carryOut);
+            Adjoint ApplyAndAssuming0Target(x, y, carryOut);
+            CNOT(carryIn, x);
+            CNOT(x, y);
+        }
+        adjoint auto;
+        controlled (ctls, ...) {
+            Fact(Length(ctls) == 1, "UncarryForInc should be controlled by exactly one control qubit.");
+
+            let ctl = ctls[0];
+
+            CNOT(carryIn, carryOut);
+            Adjoint ApplyAndAssuming0Target(x, y, carryOut);
+            CCNOT(ctl, x, y); // Controlled X(ctls + [x], y);
+            CNOT(carryIn, x);
+            CNOT(carryIn, y);
+        }
+        controlled adjoint auto;
+    }
+
+    /// # Summary
+    /// Applies AND gate between `control1` and `control2` and stores the result
+    /// in `target` assuming `target` is in |0> state.
+    ///
+    /// # Description
+    /// Inverts `target` if and only if both controls are 1, but assumes that
+    /// `target` is in state 0. The operation has T-count 4, T-depth 2 and
+    /// requires no helper qubit, and may therefore be preferable to a CCNOT
+    /// operation, if `target` is known to be 0.
+    /// The adjoint of this operation is measurement based and requires no T
+    /// gates (but requires target to support branching on measurements).
+    /// Although the Toffoli gate (CCNOT) will perform faster in simulations,
+    /// this version has lower T gate requirements.
+    /// # References
+    /// - Cody Jones: "Novel constructions for the fault-tolerant Toffoli gate",
+    ///   Phys. Rev. A 87, 022328, 2013
+    ///   [arXiv:1212.5069](https://arxiv.org/abs/1212.5069)
+    ///   doi:10.1103/PhysRevA.87.022328
+    @Config(Adaptive)
+    internal operation ApplyAndAssuming0Target(control1 : Qubit, control2 : Qubit, target : Qubit) : Unit is Adj {
+        // NOTE: Eventually this operation will be public and intrinsic.
+        body (...) {
+            CCNOT(control1, control2, target);
+        }
+        adjoint (...) {
+            H(target);
+            if M(target) == One {
+                Reset(target);
+                CZ(control1, control2);
+            }
+        }
+    }
+
+    internal operation ApplyOrAssuming0Target(control1 : Qubit, control2 : Qubit, target : Qubit) : Unit is Adj {
+        within {
+            X(control1);
+            X(control2);
+        } apply {
+            ApplyAndAssuming0Target(control1, control2, target);
+            X(target);
+        }
+    }
+
+    /// # Summary
+    /// Applies AND gate between `control1` and `control2` and stores the result
+    /// in `target` assuming `target` is in |0> state.
+    ///
+    /// # Description
+    /// Inverts `target` if and only if both controls are 1, but assumes that
+    /// `target` is in state 0. The operation has T-count 4, T-depth 2 and
+    /// requires no helper qubit, and may therefore be preferable to a CCNOT
+    /// operation, if `target` is known to be 0.
+    /// This version is suitable for Base profile.
+    /// Although the Toffoli gate (CCNOT) will perform faster in simulations,
+    /// this version has lower T gate requirements.
+    /// # References
+    /// - Cody Jones: "Novel constructions for the fault-tolerant Toffoli gate",
+    ///   Phys. Rev. A 87, 022328, 2013
+    ///   [arXiv:1212.5069](https://arxiv.org/abs/1212.5069)
+    ///   doi:10.1103/PhysRevA.87.022328
+    @Config(not Adaptive)
+    internal operation ApplyAndAssuming0Target(control1 : Qubit, control2 : Qubit, target : Qubit) : Unit is Adj {
+        H(target);
+        T(target);
+        CNOT(control1, target);
+        CNOT(control2, target);
+        within {
+            CNOT(target, control1);
+            CNOT(target, control2);
+        } apply {
+            Adjoint T(control1);
+            Adjoint T(control2);
+            T(target);
+        }
+        H(target);
+        S(target);
+    }
+
+    /// # Summary
+    /// Computes carries for the look-ahead adder
+    internal operation ComputeCarries(ps : Qubit[], gs : Qubit[]) : Unit is Adj {
+        let n = Length(gs);
+        Fact(Length(ps) + 1 == n, "Register gs must be one qubit longer than register gs.");
+
+        let T = Floor(Lg(IntAsDouble(n)));
+        use qs = Qubit[n - HammingWeightI(n) - T];
+
+        let registerPartition = MappedOverRange(t -> Floor(IntAsDouble(n) / IntAsDouble(2^t)) - 1, 1..T - 1);
+        let pWorkspace = [ps] + Partitioned(registerPartition, qs);
+
+        within {
+            PRounds(pWorkspace);
+        } apply {
+            // U_G
+            GRounds(pWorkspace, gs);
+
+            // U_C
+            CRounds(pWorkspace, gs);
+        }
+    }
+
+    /// # Summary
+    /// Computes all p[i, j] values in workspace for the look-ahead adder.
+    ///
+    /// The register array `pWorkspace` has T entries, where T = ⌊log₂ n⌋.
+    ///
+    /// The first entry `pWorkspace[0]` is initialized with `P_0` which is
+    /// computed before `ComputeCarries` is called.  The other registers are
+    /// 0-initialized and will be computed in successive rounds t = 1, ..., T - 1.
+    ///
+    /// In each round t we compute
+    ///
+    /// p[i, j] = p[2ᵗ × m, 2ᵗ × (m + 1)] = p[i, k] ∧ p[k, j]
+    ///
+    /// in `pWorkspace[t][m - 1]` and use that for k = 2ᵗ × m + 2ᵗ⁻¹, p[i, k] and p[k, j]
+    /// have already been computed in round t - 1 in `pWorkspace[t - 1][2 * m - 1]` and
+    /// `pWorkspace[t - 1][2 * m]`, respectively.
+    internal operation PRounds(pWorkspace : Qubit[][]) : Unit is Adj {
+        for ws in Windows(2, pWorkspace) {
+            // note that we are using Rest, since pWorkspace[t - 1][0] is never
+            // accessed in round t.
+            let (current, next) = (Rest(ws[0]), ws[1]);
+
+            for m in IndexRange(next) {
+                ApplyAndAssuming0Target(current[2 * m], current[2 * m + 1], next[m]);
+            }
+        }
+    }
+
+    /// # Summary
+    /// Computes g[i ∧ (i + 1), i + 1] into gs[i] for the look-ahead adder.
+    ///
+    /// The register gs has n entries initialized to gs[i] = g[i, i + 1].
+    ///
+    /// After successive rounds t = 1, ..., T, the register is updated to
+    /// gs[i] = g[i ∧ (i + 1), i + 1], from which we can compute the carries
+    /// in the C-rounds.
+    internal operation GRounds(pWorkspace : Qubit[][], gs : Qubit[]) : Unit is Adj {
+        let T = Length(pWorkspace);
+        let n = Length(gs);
+
+        for t in 1..T {
+            let length = Floor(IntAsDouble(n) / IntAsDouble(2^t)) - 1;
+            let ps = pWorkspace[t - 1][0..2...];
+
+            for m in 0..length {
+                CCNOT(gs[2^t * m + 2^(t - 1) - 1], ps[m], gs[2^t * m + 2^t - 1]);
+            }
+        }
+    }
+
+    /// # Summary
+    /// Computes carries into gs for the look-ahead adder.
+    internal operation CRounds(pWorkspace : Qubit[][], gs : Qubit[]) : Unit is Adj {
+        let n = Length(gs);
+
+        let start = Floor(Lg(IntAsDouble(2 * n) / 3.0));
+        for t in start..-1..1 {
+            let length = Floor(IntAsDouble(n - 2^(t - 1)) / IntAsDouble(2^t));
+            let ps = pWorkspace[t - 1][1..2...];
+
+            for m in 1..length {
+                CCNOT(gs[2^t * m - 1], ps[m - 1], gs[2^t * m + 2^(t - 1) - 1]);
+            }
+        }
+    }
+
+    internal operation PhaseGradient(qs : Qubit[]) : Unit is Adj + Ctl {
+        for i in IndexRange(qs) {
+            R1Frac(1, i, qs[i]);
+        }
+    }
+
+    //
+    // Internal operations for comparisons
+    //
+
+    /// # Summary
+    /// Applies `action` to `target` if register `x` is greater or equal to BigInt `c`
+    /// (if `invertControl` is false). If `invertControl` is true, the `action`
+    /// is applied in the opposite situation.
+    internal operation ApplyActionIfGreaterThanOrEqualConstant<'T>(
+        invertControl : Bool,
+        action : 'T => Unit is Adj + Ctl,
+        c : BigInt,
+        x : Qubit[],
+        target : 'T
+    ) : Unit is Adj + Ctl {
+
+        let bitWidth = Length(x);
+        if c == 0L {
+            if not invertControl {
+                action(target);
+            }
+        } elif c >= (2L^bitWidth) {
+            if invertControl {
+                action(target);
+            }
+        } else {
+            // normalize constant
+            let l = TrailingZeroCountL(c);
+
+            let cNormalized = c >>> l;
+            let xNormalized = x[l...];
+            let bitWidthNormalized = Length(xNormalized);
+
+            // If c == 2L^(bitwidth - 1), then bitWidthNormalized will be 1,
+            // and qs will be empty.  In that case, we do not need to compute
+            // any temporary values, and some optimizations are apply, which
+            // are considered in the remainder.
+            use qs = Qubit[bitWidthNormalized - 1];
+            let cs1 = IsEmpty(qs) ? [] | [Head(xNormalized)] + Most(qs);
+
+            Fact(Length(cs1) == Length(qs), "Arrays should be of the same length.");
+
+            within {
+                for i in 0..Length(cs1) - 1 {
+                    let op = cNormalized &&& (1L <<< (i + 1)) != 0L ? ApplyAndAssuming0Target | ApplyOrAssuming0Target;
+                    op(cs1[i], xNormalized[i + 1], qs[i]);
+                }
+            } apply {
+                let control = IsEmpty(qs) ? Tail(x) | Tail(qs);
+                within {
+                    if invertControl {
+                        X(control);
+                    }
+                } apply {
+                    Controlled action([control], target);
+                }
+            }
+        }
+    }
+
+    /// # Summary
+    /// Applies `action` to `target` if the sum of `x` and `y` registers
+    /// overflows, i.e. there's a carry out (if `invertControl` is false).
+    /// If `invertControl` is true, the `action` is applied when there's no carry out.
+    internal operation ApplyActionIfSumOverflows<'T>(
+        action : 'T => Unit is Adj + Ctl,
+        x : Qubit[],
+        y : Qubit[],
+        invertControl : Bool,
+        target : 'T
+    ) : Unit is Adj + Ctl {
+
+        let n = Length(x);
+        Fact(n >= 1, "Registers must contain at least one qubit.");
+        Fact(Length(y) == n, "Registers must be of the same length.");
+
+        use carries = Qubit[n];
+
+        within {
+            CarryWith1CarryIn(x[0], y[0], carries[0]);
+            for i in 1..n - 1 {
+                CarryForInc(carries[i - 1], x[i], y[i], carries[i]);
+            }
+        } apply {
+            within {
+                if invertControl {
+                    X(carries[n - 1]);
+                }
+            } apply {
+                Controlled action([carries[n - 1]], target);
+            }
+        }
+    }
+
+    /// # Summary
+    /// Computes carry out assuming carry in is 1.
+    /// Simplified version that is only applicable for scenarios
+    /// where controlled version is the same as non-controlled.
+    internal operation CarryWith1CarryIn(
+        x : Qubit,
+        y : Qubit,
+        carryOut : Qubit
+    ) : Unit is Adj + Ctl {
+
+        body (...) {
+            X(x);
+            X(y);
+            ApplyAndAssuming0Target(x, y, carryOut);
+            X(carryOut);
+        }
+
+        adjoint auto;
+
+        controlled (ctls, ...) {
+            Fact(Length(ctls) <= 1, "Number of control lines must be at most 1");
+            CarryWith1CarryIn(x, y, carryOut);
+        }
+
+        controlled adjoint auto;
+    }
+
+    /// # Summary
+    /// This wrapper allows operations that support only one control
+    /// qubit to be used in a multi-controlled scenarios. It provides
+    /// controlled version that collects controls into one qubit
+    /// by applying AND chain using auxiliary qubit array.
+    internal operation ApplyAsSinglyControlled<'TIn>(
+        op : ('TIn => Unit is Adj + Ctl),
+        input : 'TIn
+    ) : Unit is Adj + Ctl {
+
+        body (...) {
+            op(input);
+        }
+
+        controlled (ctls, ...) {
+            let n = Length(ctls);
+            if n == 0 {
+                op(input);
+            } elif n == 1 {
+                Controlled op(ctls, input);
+            } else {
+                use aux = Qubit[n - 1];
+                within {
+                    LogDepthAndChain(ctls, aux);
+                } apply {
+                    Controlled op([Tail(aux)], input);
+                }
+            }
+        }
+    }
+
+    /// # Summary
+    /// This helper function computes the AND of all control bits in `ctls` into
+    /// the last qubit of `tgts`, using the other qubits in `tgts` as helper
+    /// qubits for the AND of subsets of control bits.  The operation has a
+    /// logarithmic depth of AND gates by aligning them using a balanced binary
+    /// tree.
+    internal operation LogDepthAndChain(ctls : Qubit[], tgts : Qubit[]) : Unit is Adj {
+        let lc = Length(ctls);
+        let lt = Length(tgts);
+
+        Fact(lc == lt + 1, $"There must be exactly one more control qubit than target qubits (got {lc}, {lt})");
+
+        if lt == 1 {
+            ApplyAndAssuming0Target(ctls[0], ctls[1], tgts[0]);
+        } elif lt == 2 {
+            ApplyAndAssuming0Target(ctls[0], ctls[1], tgts[0]);
+            ApplyAndAssuming0Target(ctls[2], tgts[0], tgts[1]);
+        } else {
+            let left = lc / 2;
+            let right = lc - left;
+
+            let ctlsLeft = ctls[...left - 1];
+            let tgtsLeft = tgts[...left - 2];
+
+            let ctlsRight = ctls[left..left + right - 1];
+            let tgtsRight = tgts[left - 1..left + right - 3];
+
+            LogDepthAndChain(ctlsLeft, tgtsLeft);
+            LogDepthAndChain(ctlsRight, tgtsRight);
+            ApplyAndAssuming0Target(Tail(tgtsLeft), Tail(tgtsRight), Tail(tgts));
+        }
+    }
+}
diff --git a/library/std/src/unstable_state_preparation.qs b/library/std/src/unstable_state_preparation.qs
new file mode 100644
index 0000000000..dfc3d5ebb7
--- /dev/null
+++ b/library/std/src/unstable_state_preparation.qs
@@ -0,0 +1,388 @@
+// Copyright (c) Microsoft Corporation.
+// Licensed under the MIT License.
+
+namespace Microsoft.Quantum.Unstable.StatePreparation {
+    open Microsoft.Quantum.Convert;
+    open Microsoft.Quantum.Diagnostics;
+    open Microsoft.Quantum.Arrays;
+    open Microsoft.Quantum.Math;
+
+    /// # Summary
+    /// Given a set of coefficients and a big-endian quantum register,
+    /// prepares a state on that register described by the given coefficients.
+    ///
+    /// # Description
+    /// This operation prepares an arbitrary quantum
+    /// state |𝜓⟩ with coefficients 𝑎ⱼ from
+    /// the n-qubit computational basis state |0...0⟩.
+    ///
+    /// The action of U on the all-zeros state is given by
+    /// $$
+    /// \begin{align}
+    ///     U \ket{0\cdots 0} = \ket{\psi} = \frac{\sum_{j=0}^{2^n-1}\alpha_j \ket{j}}{\sqrt{\sum_{j=0}^{2^n-1}|\alpha_j|^2}}.
+    /// \end{align}
+    /// $$
+    ///
+    /// # Input
+    /// ## coefficients
+    /// Array of up to 2ⁿ real coefficients. The j-th coefficient
+    /// indexes the number state |j⟩ encoded in big-endian format.
+    ///
+    /// ## qubits
+    /// Qubit register encoding number states in a big-endian format. This is
+    /// expected to be initialized in the computational basis state |0...0⟩.
+    ///
+    /// # Remarks
+    /// `coefficients` will be normalized and padded with
+    /// elements 𝑎ⱼ = 0.0 if fewer than 2ⁿ are specified.
+    ///
+    /// # Example
+    /// The following snippet prepares the quantum state |𝜓⟩=√(1/8)|0⟩+√(7/8)|2⟩=√(1/8)|00⟩+√(7/8)|10⟩
+    /// in the qubit register `qubits`.
+    /// ```qsharp
+    /// let amplitudes = [Sqrt(0.125), 0.0, Sqrt(0.875), 0.0];
+    /// use qubits = Qubit[2];
+    /// PreparePureStateD(amplitudes, qubits);
+    /// ```
+    ///
+    /// # References
+    /// - [arXiv:quant-ph/0406176](https://arxiv.org/abs/quant-ph/0406176)
+    ///   "Synthesis of Quantum Logic Circuits",
+    ///   Vivek V. Shende, Stephen S. Bullock, Igor L. Markov
+    ///
+    /// # See Also
+    /// - Microsoft.Quantum.Unstable.StatePreparation.ApproximatelyPreparePureStateCP
+    operation PreparePureStateD(coefficients : Double[], qubits : Qubit[]) : Unit is Adj + Ctl {
+        let coefficientsAsComplexPolar = Mapped(a -> ComplexAsComplexPolar(Complex(a, 0.0)), coefficients);
+        ApproximatelyPreparePureStateCP(0.0, coefficientsAsComplexPolar, qubits);
+    }
+
+    /// # Summary
+    /// Given a set of coefficients and a big-endian quantum register,
+    /// prepares a state on that register described by the given coefficients,
+    /// up to a given approximation tolerance.
+    ///
+    /// # Description
+    /// This operation prepares an arbitrary quantum
+    /// state |𝜓⟩ with complex coefficients rⱼ·𝒆^(𝒊·tⱼ) from
+    /// the n-qubit computational basis state |0...0⟩.
+    /// In particular, the action of this operation can be simulated by the
+    /// a unitary transformation U which acts on the all-zeros state as
+    ///
+    /// $$
+    /// \begin{align}
+    ///     U\ket{0...0}
+    ///         & = \ket{\psi} \\\\
+    ///         & = \frac{
+    ///                 \sum_{j=0}^{2^n-1} r_j e^{i t_j} \ket{j}
+    ///             }{
+    ///                 \sqrt{\sum_{j=0}^{2^n-1} |r_j|^2}
+    ///             }.
+    /// \end{align}
+    /// $$
+    ///
+    /// # Input
+    /// ## tolerance
+    /// The approximation tolerance to be used when preparing the given state.
+    ///
+    /// ## coefficients
+    /// Array of up to 2ⁿ complex coefficients represented by their
+    /// absolute value and phase (rⱼ, tⱼ). The j-th coefficient
+    /// indexes the number state |j⟩ encoded in a big-endian format.
+    ///
+    /// ## qubits
+    /// Qubit register encoding number states in a big-endian format. This is
+    /// expected to be initialized in the computational basis state
+    /// |0...0⟩.
+    ///
+    /// # Remarks
+    /// `coefficients` will be padded with
+    /// elements (rⱼ, tⱼ) = (0.0, 0.0) if fewer than 2ⁿ are
+    /// specified.
+    ///
+    /// # References
+    /// - [arXiv:quant-ph/0406176](https://arxiv.org/abs/quant-ph/0406176)
+    ///   "Synthesis of Quantum Logic Circuits",
+    ///   Vivek V. Shende, Stephen S. Bullock, Igor L. Markov
+    operation ApproximatelyPreparePureStateCP(
+        tolerance : Double,
+        coefficients : ComplexPolar[],
+        qubits : Qubit[]
+    ) : Unit is Adj + Ctl {
+
+        let nQubits = Length(qubits);
+        // pad coefficients at tail length to a power of 2.
+        let coefficientsPadded = Padded(-2^nQubits, ComplexPolar(0.0, 0.0), coefficients);
+        let idxTarget = 0;
+        // Determine what controls to apply
+        let rngControl = nQubits > 1 ? (1..(nQubits - 1)) | (1..0);
+        // Note we use the reversed qubits array to get the endianness ordering that we expect
+        // when corresponding qubit state to state vector index.
+        Adjoint ApproximatelyUnprepareArbitraryState(
+            tolerance,
+            coefficientsPadded,
+            rngControl,
+            idxTarget,
+            Reversed(qubits)
+        );
+    }
+
+    /// # Summary
+    /// Implementation step of arbitrary state preparation procedure.
+    internal operation ApproximatelyUnprepareArbitraryState(
+        tolerance : Double,
+        coefficients : ComplexPolar[],
+        rngControl : Range,
+        idxTarget : Int,
+        register : Qubit[]
+    ) : Unit is Adj + Ctl {
+
+        // For each 2D block, compute disentangling single-qubit rotation parameters
+        let (disentanglingY, disentanglingZ, newCoefficients) = StatePreparationSBMComputeCoefficients(coefficients);
+        if (AnyOutsideToleranceD(tolerance, disentanglingZ)) {
+            ApproximatelyMultiplexPauli(tolerance, disentanglingZ, PauliZ, register[rngControl], register[idxTarget]);
+
+        }
+        if (AnyOutsideToleranceD(tolerance, disentanglingY)) {
+            ApproximatelyMultiplexPauli(tolerance, disentanglingY, PauliY, register[rngControl], register[idxTarget]);
+        }
+        // target is now in |0> state up to the phase given by arg of newCoefficients.
+
+        // Continue recursion while there are control qubits.
+        if (IsRangeEmpty(rngControl)) {
+            let (abs, arg) = (newCoefficients[0].Magnitude, newCoefficients[0].Argument);
+            if (AbsD(arg) > tolerance) {
+                Exp([PauliI], -1.0 * arg, [register[idxTarget]]);
+            }
+        } elif (Any(c -> AbsComplexPolar(c) > tolerance, newCoefficients)) {
+            // Some coefficients are outside tolerance
+            let newControl = (RangeStart(rngControl) + 1)..RangeStep(rngControl)..RangeEnd(rngControl);
+            let newTarget = RangeStart(rngControl);
+            ApproximatelyUnprepareArbitraryState(tolerance, newCoefficients, newControl, newTarget, register);
+        }
+    }
+
+    /// # Summary
+    /// Applies a Pauli rotation conditioned on an array of qubits, truncating
+    /// small rotation angles according to a given tolerance.
+    ///
+    /// # Description
+    /// This applies a multiply controlled unitary operation that performs
+    /// rotations by angle $\theta_j$ about single-qubit Pauli operator $P$
+    /// when controlled by the $n$-qubit number state $\ket{j}$.
+    /// In particular, the action of this operation is represented by the
+    /// unitary
+    ///
+    /// $$
+    /// \begin{align}
+    ///     U = \sum^{2^n - 1}_{j=0} \ket{j}\bra{j} \otimes e^{i P \theta_j}.
+    /// \end{align}
+    /// ##
+    ///
+    /// # Input
+    /// ## tolerance
+    /// A tolerance below which small coefficients are truncated.
+    ///
+    /// ## coefficients
+    /// Array of up to $2^n$ coefficients $\theta_j$. The $j$th coefficient
+    /// indexes the number state $\ket{j}$ encoded in little-endian format.
+    ///
+    /// ## pauli
+    /// Pauli operator $P$ that determines axis of rotation.
+    ///
+    /// ## control
+    /// $n$-qubit control register that encodes number states $\ket{j}$ in
+    /// little-endian format.
+    ///
+    /// ## target
+    /// Single qubit register that is rotated by $e^{i P \theta_j}$.
+    ///
+    /// # Remarks
+    /// `coefficients` will be padded with elements $\theta_j = 0.0$ if
+    /// fewer than $2^n$ are specified.
+    internal operation ApproximatelyMultiplexPauli(
+        tolerance : Double,
+        coefficients : Double[],
+        pauli : Pauli,
+        control : Qubit[],
+        target : Qubit
+    ) : Unit is Adj + Ctl {
+
+        if pauli == PauliZ {
+            ApproximatelyMultiplexZ(tolerance, coefficients, control, target);
+        } elif pauli == PauliX {
+            within {
+                H(target);
+            } apply {
+                ApproximatelyMultiplexPauli(tolerance, coefficients, PauliZ, control, target);
+            }
+        } elif pauli == PauliY {
+            within {
+                Adjoint S(target);
+            } apply {
+                ApproximatelyMultiplexPauli(tolerance, coefficients, PauliX, control, target);
+            }
+        } else {
+            fail $"MultiplexPauli failed. Invalid pauli {pauli}.";
+        }
+    }
+
+    /// # Summary
+    /// Implementation step of arbitrary state preparation procedure.
+    internal function StatePreparationSBMComputeCoefficients(
+        coefficients : ComplexPolar[]
+    ) : (Double[], Double[], ComplexPolar[]) {
+
+        mutable disentanglingZ = [];
+        mutable disentanglingY = [];
+        mutable newCoefficients = [];
+
+        for idxCoeff in 0..2..Length(coefficients) - 1 {
+            let (rt, phi, theta) = BlochSphereCoordinates(coefficients[idxCoeff], coefficients[idxCoeff + 1]);
+            set disentanglingZ += [0.5 * phi];
+            set disentanglingY += [0.5 * theta];
+            set newCoefficients += [rt];
+        }
+
+        return (disentanglingY, disentanglingZ, newCoefficients);
+    }
+
+    /// # Summary
+    /// Computes the Bloch sphere coordinates for a single-qubit state.
+    ///
+    /// Given two complex numbers $a0, a1$ that represent the qubit state, computes coordinates
+    /// on the Bloch sphere such that
+    /// $a0 \ket{0} + a1 \ket{1} = r e^{it}(e^{-i \phi /2}\cos{(\theta/2)}\ket{0}+e^{i \phi /2}\sin{(\theta/2)}\ket{1})$.
+    ///
+    /// # Input
+    /// ## a0
+    /// Complex coefficient of state $\ket{0}$.
+    /// ## a1
+    /// Complex coefficient of state $\ket{1}$.
+    ///
+    /// # Output
+    /// A tuple containing `(ComplexPolar(r, t), phi, theta)`.
+    internal function BlochSphereCoordinates(
+        a0 : ComplexPolar,
+        a1 : ComplexPolar
+    ) : (ComplexPolar, Double, Double) {
+
+        let abs0 = AbsComplexPolar(a0);
+        let abs1 = AbsComplexPolar(a1);
+        let arg0 = ArgComplexPolar(a0);
+        let arg1 = ArgComplexPolar(a1);
+        let r = Sqrt(abs0 * abs0 + abs1 * abs1);
+        let t = 0.5 * (arg0 + arg1);
+        let phi = arg1 - arg0;
+        let theta = 2.0 * ArcTan2(abs1, abs0);
+        return (ComplexPolar(r, t), phi, theta);
+    }
+
+    /// # Summary
+    /// Applies a Pauli Z rotation conditioned on an array of qubits, truncating
+    /// small rotation angles according to a given tolerance.
+    ///
+    /// # Description
+    /// This applies the multiply controlled unitary operation that performs
+    /// rotations by angle $\theta_j$ about single-qubit Pauli operator $Z$
+    /// when controlled by the $n$-qubit number state $\ket{j}$.
+    /// In particular, this operation can be represented by the unitary
+    ///
+    /// $$
+    /// \begin{align}
+    ///     U = \sum^{2^n-1}_{j=0} \ket{j}\bra{j} \otimes e^{i Z \theta_j}.
+    /// \end{align}
+    /// $$
+    ///
+    /// # Input
+    /// ## tolerance
+    /// A tolerance below which small coefficients are truncated.
+    ///
+    /// ## coefficients
+    /// Array of up to $2^n$ coefficients $\theta_j$. The $j$th coefficient
+    /// indexes the number state $\ket{j}$ encoded in little-endian format.
+    ///
+    /// ## control
+    /// $n$-qubit control register that encodes number states $\ket{j}$ in
+    /// little-endian format.
+    ///
+    /// ## target
+    /// Single qubit register that is rotated by $e^{i P \theta_j}$.
+    ///
+    /// # Remarks
+    /// `coefficients` will be padded with elements $\theta_j = 0.0$ if
+    /// fewer than $2^n$ are specified.
+    ///
+    /// # References
+    /// - [arXiv:quant-ph/0406176](https://arxiv.org/abs/quant-ph/0406176)
+    ///   "Synthesis of Quantum Logic Circuits",
+    ///   Vivek V. Shende, Stephen S. Bullock, Igor L. Markov
+    internal operation ApproximatelyMultiplexZ(
+        tolerance : Double,
+        coefficients : Double[],
+        control : Qubit[],
+        target : Qubit
+    ) : Unit is Adj + Ctl {
+
+        body (...) {
+            // pad coefficients length at tail to a power of 2.
+            let coefficientsPadded = Padded(-2^Length(control), 0.0, coefficients);
+
+            if Length(coefficientsPadded) == 1 {
+                // Termination case
+                if AbsD(coefficientsPadded[0]) > tolerance {
+                    Exp([PauliZ], coefficientsPadded[0], [target]);
+                }
+            } else {
+                // Compute new coefficients.
+                let (coefficients0, coefficients1) = MultiplexZCoefficients(coefficientsPadded);
+                ApproximatelyMultiplexZ(tolerance, coefficients0, Most(control), target);
+                if AnyOutsideToleranceD(tolerance, coefficients1) {
+                    within {
+                        CNOT(Tail(control), target);
+                    } apply {
+                        ApproximatelyMultiplexZ(tolerance, coefficients1, Most(control), target);
+                    }
+                }
+            }
+        }
+
+        controlled (controlRegister, ...) {
+            // pad coefficients length to a power of 2.
+            let coefficientsPadded = Padded(2^(Length(control) + 1), 0.0, Padded(-2^Length(control), 0.0, coefficients));
+            let (coefficients0, coefficients1) = MultiplexZCoefficients(coefficientsPadded);
+            ApproximatelyMultiplexZ(tolerance, coefficients0, control, target);
+            if AnyOutsideToleranceD(tolerance, coefficients1) {
+                within {
+                    Controlled X(controlRegister, target);
+                } apply {
+                    ApproximatelyMultiplexZ(tolerance, coefficients1, control, target);
+                }
+            }
+        }
+    }
+
+    /// # Summary
+    /// Implementation step of multiply-controlled Z rotations.
+    internal function MultiplexZCoefficients(coefficients : Double[]) : (Double[], Double[]) {
+        let newCoefficientsLength = Length(coefficients) / 2;
+        mutable coefficients0 = [];
+        mutable coefficients1 = [];
+
+        for idxCoeff in 0..newCoefficientsLength - 1 {
+            set coefficients0 += [0.5 * (coefficients[idxCoeff] + coefficients[idxCoeff + newCoefficientsLength])];
+            set coefficients1 += [0.5 * (coefficients[idxCoeff] - coefficients[idxCoeff + newCoefficientsLength])];
+        }
+
+        return (coefficients0, coefficients1);
+    }
+
+    internal function AnyOutsideToleranceD(tolerance : Double, coefficients : Double[]) : Bool {
+        // NOTE: This function is not used as the only recursion termination condition
+        // only to determine if the multiplex step needs to be applied.
+        // For tolerance 0.0 it is always applied due to >= comparison.
+        Any(coefficient -> AbsD(coefficient) >= tolerance, coefficients)
+    }
+
+    export PreparePureStateD, ApproximatelyPreparePureStateCP;
+}
diff --git a/library/std/src/unstable_table_lookup.qs b/library/std/src/unstable_table_lookup.qs
new file mode 100644
index 0000000000..5ff35395b6
--- /dev/null
+++ b/library/std/src/unstable_table_lookup.qs
@@ -0,0 +1,278 @@
+// Copyright (c) Microsoft Corporation.
+// Licensed under the MIT License.
+
+namespace Microsoft.Quantum.Unstable.TableLookup {
+    open Microsoft.Quantum.Arrays;
+    open Microsoft.Quantum.Convert;
+    open Microsoft.Quantum.Diagnostics;
+    open Microsoft.Quantum.Math;
+    open Microsoft.Quantum.ResourceEstimation;
+    open Microsoft.Quantum.Unstable.Arithmetic;
+
+    /// # Summary
+    /// Performs table lookup using a SELECT network
+    ///
+    /// # Description
+    /// Assuming a zero-initialized `target` register, this operation will
+    /// initialize it with the bitstrings in `data` at indices according to the
+    /// computational values of the `address` register.
+    ///
+    /// # Input
+    /// ## data
+    /// The classical table lookup data which is prepared in `target` with
+    /// respect to the state in `address`. The length of data must be less than
+    /// 2ⁿ, where 𝑛 is the length of `address`. Each entry in data must have
+    /// the same length that must be equal to the length of `target`.
+    /// ## address
+    /// Address register
+    /// ## target
+    /// Zero-initialized target register
+    ///
+    /// # Remarks
+    /// The implementation of the SELECT network is based on unary encoding as
+    /// presented in [1].  The recursive implementation of that algorithm is
+    /// presented in [3].  The adjoint variant is optimized using a
+    /// measurement-based unlookup operation [3]. The controlled adjoint variant
+    /// is not optimized using this technique.
+    ///
+    /// # References
+    /// 1. [arXiv:1805.03662](https://arxiv.org/abs/1805.03662)
+    ///    "Encoding Electronic Spectra in Quantum Circuits with Linear T
+    ///    Complexity"
+    /// 2. [arXiv:1905.07682](https://arxiv.org/abs/1905.07682)
+    ///    "Windowed arithmetic"
+    /// 3. [arXiv:2211.01133](https://arxiv.org/abs/2211.01133)
+    ///    "Space-time optimized table lookup"
+    operation Select(
+        data : Bool[][],
+        address : Qubit[],
+        target : Qubit[]
+    ) : Unit is Adj + Ctl {
+        body (...) {
+            let (N, n) = DimensionsForSelect(data, address);
+
+            if N == 1 {
+                // base case
+                WriteMemoryContents(Head(data), target);
+            } else {
+                let (most, tail) = MostAndTail(address[...n - 1]);
+                let parts = Partitioned([2^(n - 1)], data);
+
+                within {
+                    X(tail);
+                } apply {
+                    SinglyControlledSelect(tail, parts[0], most, target);
+                }
+
+                SinglyControlledSelect(tail, parts[1], most, target);
+            }
+        }
+        adjoint (...) {
+            Unlookup(Select, data, address, target);
+        }
+
+        controlled (ctls, ...) {
+            let numCtls = Length(ctls);
+
+            if numCtls == 0 {
+                Select(data, address, target);
+            } elif numCtls == 1 {
+                SinglyControlledSelect(ctls[0], data, address, target);
+            } else {
+                use andChainTarget = Qubit();
+                let andChain = MakeAndChain(ctls, andChainTarget);
+                use helper = Qubit[andChain.NGarbageQubits];
+
+                within {
+                    andChain.Apply(helper);
+                } apply {
+                    SinglyControlledSelect(andChainTarget, data, address, target);
+                }
+            }
+        }
+
+        controlled adjoint (ctls, ...) {
+            Controlled Select(ctls, (data, address, target));
+        }
+    }
+
+    internal operation SinglyControlledSelect(
+        ctl : Qubit,
+        data : Bool[][],
+        address : Qubit[],
+        target : Qubit[]
+    ) : Unit {
+        let (N, n) = DimensionsForSelect(data, address);
+
+        if BeginEstimateCaching("Microsoft.Quantum.Unstable.TableLookup.SinglyControlledSelect", N) {
+            if N == 1 {
+                // base case
+                Controlled WriteMemoryContents([ctl], (Head(data), target));
+            } else {
+                use helper = Qubit();
+
+                let (most, tail) = MostAndTail(address[...n - 1]);
+                let parts = Partitioned([2^(n - 1)], data);
+
+                within {
+                    X(tail);
+                } apply {
+                    ApplyAndAssuming0Target(ctl, tail, helper);
+                }
+
+                SinglyControlledSelect(helper, parts[0], most, target);
+
+                CNOT(ctl, helper);
+
+                SinglyControlledSelect(helper, parts[1], most, target);
+
+                Adjoint ApplyAndAssuming0Target(ctl, tail, helper);
+            }
+
+            EndEstimateCaching();
+        }
+    }
+
+    internal function DimensionsForSelect(
+        data : Bool[][],
+        address : Qubit[]
+    ) : (Int, Int) {
+        let N = Length(data);
+        Fact(N > 0, "data cannot be empty");
+
+        let n = Ceiling(Lg(IntAsDouble(N)));
+        Fact(
+            Length(address) >= n,
+            $"address register is too small, requires at least {n} qubits"
+        );
+
+        return (N, n);
+    }
+
+    internal operation WriteMemoryContents(
+        value : Bool[],
+        target : Qubit[]
+    ) : Unit is Adj + Ctl {
+        Fact(
+            Length(value) == Length(target),
+            "number of data bits must equal number of target qubits"
+        );
+
+        ApplyPauliFromBitString(PauliX, true, value, target);
+    }
+
+    /// # References
+    /// - [arXiv:1905.07682](https://arxiv.org/abs/1905.07682)
+    ///   "Windowed arithmetic"
+    internal operation Unlookup(
+        lookup : (Bool[][], Qubit[], Qubit[]) => Unit,
+        data : Bool[][],
+        select : Qubit[],
+        target : Qubit[]
+    ) : Unit {
+        let numBits = Length(target);
+        let numAddressBits = Length(select);
+
+        let l = MinI(Floor(Lg(IntAsDouble(numBits))), numAddressBits - 1);
+        Fact(
+            l < numAddressBits,
+            $"l = {l} must be smaller than {numAddressBits}"
+        );
+
+        let res = Mapped(r -> r == One, ForEach(MResetX, target));
+
+        let dataFixup = Chunks(2^l, Padded(-2^numAddressBits, false, Mapped(MustBeFixed(res, _), data)));
+
+        let numAddressBitsFixup = numAddressBits - l;
+
+        let selectParts = Partitioned([l], select);
+        let targetFixup = target[...2^l - 1];
+
+        within {
+            EncodeUnary(selectParts[0], targetFixup);
+            ApplyToEachA(H, targetFixup);
+        } apply {
+            lookup(dataFixup, selectParts[1], targetFixup);
+        }
+    }
+
+    // Checks whether specific bit string `data` must be fixed for a given
+    // measurement result `result`.
+    //
+    // Returns true if the number of indices for which both result and data are
+    // `true` is odd.
+    internal function MustBeFixed(result : Bool[], data : Bool[]) : Bool {
+        mutable state = false;
+        for i in IndexRange(result) {
+            set state = state != (result[i] and data[i]);
+        }
+        state
+    }
+
+    // Computes unary encoding of value in `input` into `target`
+    //
+    // Assumptions:
+    //    - `target` is zero-initialized
+    //    - length of `input` is n
+    //    - length of `target` is 2^n
+    internal operation EncodeUnary(
+        input : Qubit[],
+        target : Qubit[]
+    ) : Unit is Adj {
+        Fact(
+            Length(target) == 2^Length(input),
+            $"target register should be of length {2^Length(input)}, but is {Length(target)}"
+        );
+
+        X(Head(target));
+
+        for i in IndexRange(input) {
+            if i == 0 {
+                CNOT(input[i], target[1]);
+                CNOT(target[1], target[0]);
+            } else {
+                // targets are the first and second 2^i qubits of the target register
+                let split = Partitioned([2^i, 2^i], target);
+                for j in IndexRange(split[0]) {
+                    ApplyAndAssuming0Target(input[i], split[0][j], split[1][j]);
+                    CNOT(split[1][j], split[0][j]);
+                }
+            }
+        }
+
+    }
+
+    internal newtype AndChain = (
+        NGarbageQubits : Int,
+        Apply : Qubit[] => Unit is Adj
+    );
+
+    internal function MakeAndChain(ctls : Qubit[], target : Qubit) : AndChain {
+        AndChain(
+            MaxI(Length(ctls) - 2, 0),
+            helper => AndChainOperation(ctls, helper, target)
+        )
+    }
+
+    internal operation AndChainOperation(ctls : Qubit[], helper : Qubit[], target : Qubit) : Unit is Adj {
+        let n = Length(ctls);
+
+        Fact(Length(helper) == MaxI(n - 2, 0), "Invalid number of helper qubits");
+
+        if n == 0 {
+            X(target);
+        } elif n == 1 {
+            CNOT(ctls[0], target);
+        } else {
+            let ctls1 = ctls[0..0] + helper;
+            let ctls2 = ctls[1...];
+            let tgts = helper + [target];
+
+            for idx in IndexRange(tgts) {
+                ApplyAndAssuming0Target(ctls1[idx], ctls2[idx], tgts[idx]);
+            }
+        }
+    }
+
+    export Select;
+}