diff --git a/library/table_lookup/README.md b/library/table_lookup/README.md
new file mode 100644
index 0000000000..db269f5f32
--- /dev/null
+++ b/library/table_lookup/README.md
@@ -0,0 +1,31 @@
+# table_lookup library
+
+The `table_lookup` library defines various primitives useful to perform computation and uncomputation of table lookup. It also defines wrapper function
+which uses one of the approaches depending on options.
+
+## Lookup
+
+`Lookup` is the main operation implementing various table lookup algorithms and options. Note, that most unlookup algorithms are measurement-based and return target register to zero state.
+
+### Options for lookup
+
+* `DoStdLookup` - Use lookup algorithm defined in the Q# standard library.
+* `DoMCXLookup` - Use naive lookup algorithm via multicontrolled X gates. See [arXiv:1805.03662](https://arxiv.org/abs/1805.03662), Section A.
+* `DoRecursiveSelectLookup` - Use select network algorithm via recursion. See [arXiv:2211.01133](https://arxiv.org/abs/2211.01133), Section 2.
+* `DoPPLookup` - Use lookup algorithm via power products without address split. See [arXiv:2505.15917](https://arxiv.org/abs/2505.15917), Section A.4.
+* `DoSplitPPLookup` - Use lookup algorithm via power products with address split. See [arXiv:2505.15917](https://arxiv.org/abs/2505.15917), Section A.4.
+
+### Options for unlookup
+
+* `DoStdUnlookup` - Use unlookup algorithm defined in the Q# standard library.
+* `DoUnlookupViaLookup` - Perform unlookup via the same algorithm as lookup as it is self-adjoint.
+* `DoMCXUnlookup` - Perform measurement-based unlookup with corrections via multicontrolled X gates. See [arXiv:2211.01133](https://arxiv.org/abs/2211.01133), Section 2.
+* `DoPPUnlookup` - Perform measurement-based unlookup with corrections via power products without address split (Phase lookup). See [arXiv:2505.15917](https://arxiv.org/abs/2505.15917), Section A.3.
+* `DoSplitPPUnlookup` - Perform measurement-based unlookup with corrections via power products with address split (Phase lookup). See [arXiv:2505.15917](https://arxiv.org/abs/2505.15917), Section A.3.
+
+## Potential future work
+
+* Add more control how uncomputation of AND gate is performed.
+* Add resource estimation hints.
+* If gate set includes multi-target gates, code can be optimized to use those.
+* Implement delayed combined corrections.
diff --git a/library/table_lookup/qsharp.json b/library/table_lookup/qsharp.json
new file mode 100644
index 0000000000..4f03cf1bdd
--- /dev/null
+++ b/library/table_lookup/qsharp.json
@@ -0,0 +1,13 @@
+{
+  "files": [
+    "src/Lookup.qs",
+    "src/LookupViaPP.qs",
+    "src/Main.qs",
+    "src/Multicontrolled.qs",
+    "src/PhaseLookup.qs",
+    "src/PowerProducts.qs",
+    "src/RecursiveSelect.qs",
+    "src/Tests.qs",
+    "src/Utils.qs"
+  ]
+}
diff --git a/library/table_lookup/src/Lookup.qs b/library/table_lookup/src/Lookup.qs
new file mode 100644
index 0000000000..2b8cfad06c
--- /dev/null
+++ b/library/table_lookup/src/Lookup.qs
@@ -0,0 +1,292 @@
+// Copyright (c) Microsoft Corporation.
+// Licensed under the MIT License.
+
+import Std.Arrays.*;
+
+import Utils.*;
+import Multicontrolled.*;
+import RecursiveSelect.*;
+import LookupViaPP.*;
+import PhaseLookup.*;
+
+// ----------------------------------------------
+// Lookup algorithm options.
+
+/// Use lookup algorithm defined in the standard library.
+function DoStdLookup() : Int {
+    0
+}
+
+/// Use basic lookup algorithm with multicontrolled X gates.
+function DoMCXLookup() : Int {
+    1
+}
+
+/// Use recursive SELECT network as lookup algorithm.
+function DoRecursiveSelectLookup() : Int {
+    2
+}
+
+/// Use lookup algorithm via power products without address split.
+function DoPPLookup() : Int {
+    3
+}
+
+/// Use lookup algorithm via power products with address split.
+function DoSplitPPLookup() : Int {
+    4
+}
+
+// ----------------------------------------------
+// Unlookup algorithm options.
+
+/// Use unlookup algorithm defined in the standard library.
+function DoStdUnlookup() : Int {
+    0
+}
+
+/// Use the same unlookup algorithm as lookup.
+/// This is always possible as lookup is self-adjoint.
+function DoUnlookupViaLookup() : Int {
+    1
+}
+
+/// Use unlookup algorithm with multicontrolled X gates.
+/// This options is measurement based and returns target to zero state.
+function DoMCXUnlookup() : Int {
+    2
+}
+
+/// Use unlookup algorithm via power products without address split (Phase lookup).
+/// This options is measurement based and returns target to zero state.
+function DoPPUnlookup() : Int {
+    3
+}
+
+/// Use unlookup algorithm via power products with address split (Phase lookup).
+/// This options is measurement based and returns target to zero state.
+function DoSplitPPUnlookup() : Int {
+    4
+}
+
+/// # Summary
+/// Options for table lookup and unlookup operations.
+struct LookupOptions {
+    // Specifies lookup algorithm. Options:
+    // `DoStdLookup`, `DoMCXLookup`, `DoRecursiveSelectLookup`, `DoPPLookup`, `DoSplitPPLookup`.
+    lookupAlgorithm : Int,
+
+    // Specifies unlookup algorithm. Options:
+    // `DoStdUnlookup`, `DoUnlookupViaLookup`, `DoMCXUnlookup`, `DoPPUnlookup`, `DoSplitPPUnlookup`.
+    unlookupAlgorithm : Int,
+    // Suggests using measurement-based uncomputation where applicable.
+    // Note that some algorithms are measurement-based only and some cannot use measurements.
+    // If `true`, use measurement-based uncomputations. Example: prefer adjoint AND.
+    // If `false`, avoid measurement-based uncomputations. Example: prefer adjoint CCNOT.
+    preferMeasurementBasedUncomputation : Bool,
+
+    // If `true`, an error is raised if data is longer than addressable space.
+    // If `false`, longer data beyond addressable space is ignored.
+    failOnLongData : Bool,
+
+    // If `true`, an error is raised if data is shorter than addressable space.
+    // If `false`, shorter data is tolerated according to respectExcessiveAddress.
+    failOnShortData : Bool,
+
+    // If `true`, all address qubits are respected and used.
+    // Addressing beyond data length yields the same result as if the data was padded with `false` values.
+    // If `false`, addressing beyond data length yields undefined results.
+    // As one consequence, when data is shorter than addressable space, higher address qubits are ignored.
+    respectExcessiveAddress : Bool,
+}
+
+/// # Summary
+/// Default lookup options. Use power products with register split for lookup and unlookup.
+function DefaultLookupOptions() : LookupOptions {
+    new LookupOptions {
+        lookupAlgorithm = DoSplitPPLookup(),
+        unlookupAlgorithm = DoSplitPPUnlookup(),
+        failOnLongData = false,
+        failOnShortData = false,
+        respectExcessiveAddress = false,
+        preferMeasurementBasedUncomputation = true,
+    }
+}
+
+/// # Summary
+/// Performs table lookup/unlookup operation using specified algorithm and other options.
+///
+/// # Input
+/// ## options
+/// LookupOptions defining lookup and unlookup algorithms and other parameters.
+/// ## data
+/// The data table to be looked up. Each entry is a Bool array the size of the target register.
+/// ## address
+/// Qubit register representing the address in little-endian format.
+/// If the state of this register is one of the basis states |i⟩, and the target register is in |0⟩,
+/// the target register will be set to the value data[i] during lookup. Address can also be in superposition.
+/// ## target
+/// Qubit register to accept the target data. Must be in the |0⟩ state for some algorithms.
+/// For specifics see the corresponding algorithm implementation.
+operation Lookup(
+    options : LookupOptions,
+    data : Bool[][],
+    address : Qubit[],
+    target : Qubit[]
+) : Unit is Adj + Ctl {
+    body (...) {
+        if (options.lookupAlgorithm == DoStdLookup()) {
+            // Don't do anthing beyond standard library select.
+            Std.TableLookup.Select(data, address, target);
+            return ();
+        }
+
+        let input = PrepareAddressAndData(options, address, data);
+
+        if options.lookupAlgorithm == DoMCXLookup() {
+            // Basic lookup via multicontrolled X gates.
+            LookupViaMCX(input.fitData, input.fitAddress, target);
+            return ();
+        }
+
+        if options.lookupAlgorithm == DoRecursiveSelectLookup() {
+            // Recursive select implementation.
+            if (options.respectExcessiveAddress) {
+                RecursiveLookup(options.preferMeasurementBasedUncomputation, input.fitData, input.fitAddress, target);
+            } else {
+                RecursiveLookupOpt(options.preferMeasurementBasedUncomputation, input.fitData, input.fitAddress, target);
+            }
+            return ();
+        }
+
+        if options.lookupAlgorithm == DoPPLookup() {
+            // Lookup via power products without address split.
+            LookupViaPP(input.fitData, input.fitAddress, target);
+            return ();
+        }
+
+        if options.lookupAlgorithm == DoSplitPPLookup() {
+            LookupViaSplitPP(input.fitData, input.fitAddress, target);
+            return ();
+        }
+
+        fail $"Unknown lookup algorithm specified ({options.lookupAlgorithm}).";
+    }
+
+    controlled (controls, ...) {
+        let control_size = Length(controls);
+        if control_size == 0 {
+            Lookup(options, data, address, target);
+            return ();
+        }
+
+        if options.lookupAlgorithm == DoStdLookup() {
+            // Don't do anthing beyond standard library select.
+            Controlled Std.TableLookup.Select(controls, (data, address, target));
+            return ();
+        }
+
+        let input = PrepareAddressAndData(options, address, data);
+
+        if options.lookupAlgorithm == DoMCXLookup() {
+            // This is already a multicontrolled approach. Just add more controls.
+            Controlled LookupViaMCX(controls, (data, address, target));
+            return ();
+        }
+
+        // Combine multiple controls into one.
+        use aux = Qubit[control_size - 1];
+        within {
+            CombineControls(controls, aux);
+        } apply {
+            let single_control = GetCombinedControl(controls, aux);
+
+            if options.lookupAlgorithm == DoRecursiveSelectLookup() {
+                // Recursive select implementation.
+                if (options.respectExcessiveAddress) {
+                    ControlledRecursiveSelect(
+                        options.preferMeasurementBasedUncomputation,
+                        single_control,
+                        input.fitData,
+                        input.fitAddress,
+                        target
+                    );
+                } else {
+                    ControlledRecursiveSelectOpt(
+                        options.preferMeasurementBasedUncomputation,
+                        single_control,
+                        input.fitData,
+                        input.fitAddress,
+                        target
+                    );
+                }
+            } else {
+                // To use control qubit as an extra address qubit we need to respect entire address.
+                // Power products implementation already does that.
+                within {
+                    // Invert control so that data is selected when control is |1⟩.
+                    X(single_control);
+                } apply {
+                    // Add control as the most significant address qubit.
+                    if options.lookupAlgorithm == DoPPLookup() {
+                        LookupViaPP(input.fitData, input.fitAddress + [single_control], target);
+                    } elif options.lookupAlgorithm == DoSplitPPLookup() {
+                        LookupViaSplitPP(input.fitData, input.fitAddress + [single_control], target);
+                    } else {
+                        fail $"Unknown lookup algorithm specified ({options.lookupAlgorithm}).";
+                    }
+                }
+            }
+        }
+    }
+
+    adjoint (...) {
+        if (options.unlookupAlgorithm == DoStdUnlookup()) {
+            // Don't do anthing beyond standard library select.
+            Adjoint Std.TableLookup.Select(data, address, target);
+            return ();
+        }
+        if (options.unlookupAlgorithm == DoUnlookupViaLookup()) {
+            // Perform same lookup operation (as it is self-adjoint).
+            Lookup(options, data, address, target);
+            return ();
+        }
+
+        // Perform measurement-based uncomputation.
+        let input = PrepareAddressAndData(options, address, data);
+        let phaseData = MeasureAndComputePhaseData(target, input.fitData, Length(input.fitAddress));
+        // Now apply phase corrections after measurement-based uncomputation.
+
+        if options.unlookupAlgorithm == DoMCXUnlookup() {
+            // Phase lookup via multicontrolled X gates.
+            PhaseLookupViaMCX(phaseData, input.fitAddress);
+            return ();
+        }
+
+        if options.unlookupAlgorithm == DoPPUnlookup() {
+            // Phase lookup via power products without address split.
+            PhaseLookupViaPP(input.fitAddress, phaseData);
+            return ();
+        }
+
+        if options.unlookupAlgorithm == DoSplitPPUnlookup() {
+            // Phase lookup via power products with address split.
+            PhaseLookupViaSplitPP(input.fitAddress, phaseData);
+            return ();
+        }
+
+        fail $"Unknown unlookup algorithm specified ({options.unlookupAlgorithm}).";
+    }
+
+    controlled adjoint (controls, ...) {
+        if options.unlookupAlgorithm == DoStdUnlookup() {
+            // Don't do anthing beyond standard library select.
+            Controlled Adjoint Std.TableLookup.Select(controls, (data, address, target));
+            return ();
+        }
+
+        // In all other cases we perform controlled lookup as
+        // we cannot do controlled measurement-based uncomputation.
+        Controlled Lookup(controls, (options, data, address, target));
+    }
+}
diff --git a/library/table_lookup/src/LookupViaPP.qs b/library/table_lookup/src/LookupViaPP.qs
new file mode 100644
index 0000000000..f1321ef8a5
--- /dev/null
+++ b/library/table_lookup/src/LookupViaPP.qs
@@ -0,0 +1,320 @@
+// Copyright (c) Microsoft Corporation.
+// Licensed under the MIT License.
+
+import Std.Arrays.*;
+import Std.Diagnostics.*;
+
+import Utils.*;
+import PowerProducts.*;
+
+/// # Summary
+/// Performs table lookup using power products without register split.
+///
+/// # Description
+/// Table lookup is preformed using power products constructed from the address qubits.
+/// Data is processed using Fast Mobius Transform to fit power products structure.
+/// Longer data is ignored, shorter data is padded with false values.
+/// Little-endian format is used throughout.
+/// This version uses O(2^n) auxiliary qubits for n address qubits.
+///
+/// # Reference
+/// 1. [arXiv:2505.15917](https://arxiv.org/abs/2505.15917)
+///    "How to factor 2048 bit RSA integers with less than a million noisy qubits"
+///    by Craig Gidney, May 2025. Section A.4.
+operation LookupViaPP(
+    data : Bool[][],
+    address : Qubit[],
+    target : Qubit[]
+) : Unit {
+    let data_length = Length(data);
+    let address_size = Length(address);
+    let addressable_space = 1 <<< address_size;
+    let data = if (addressable_space < data_length) {
+        data[...addressable_space-1]
+    } elif (addressable_space > data_length) {
+        Padded(-addressable_space, Repeated(false, Length(target)), data)
+    } else {
+        data
+    };
+
+    // Allocate auxilliary qubits.
+    use aux_qubits = Qubit[GetAuxCountForPP(address_size)];
+
+    // Construct power products.
+    let products = ConstructPowerProducts(address, aux_qubits);
+
+    ApplyFlips(data, products, [], target);
+
+    // Undo power products.
+    DestructPowerProducts(products);
+
+}
+
+/// # Summary
+/// Performs table lookup using power products and register split.
+///
+/// # Description
+/// Table lookup is preformed using power products constructed from the address qubits.
+/// Data is processed using Fast Mobius Transform to fit power products structure.
+/// Longer data is ignored, shorter data is padded with false values.
+/// Little-endian format is used throughout. Address register is split into two halves.
+/// This version uses O(2^(n/2)) auxiliary qubits for n address qubits.
+///
+/// # Reference
+/// 1. [arXiv:2505.15917](https://arxiv.org/abs/2505.15917)
+///    "How to factor 2048 bit RSA integers with less than a million noisy qubits"
+///    by Craig Gidney, May 2025. Section A.4.
+operation LookupViaSplitPP(
+    data : Bool[][],
+    address : Qubit[],
+    target : Qubit[]
+) : Unit {
+    let data_length = Length(data);
+    let address_size = Length(address);
+    let addressable_space = 1 <<< address_size;
+    let data = if (addressable_space < data_length) {
+        data[...addressable_space-1]
+    } elif (addressable_space > data_length) {
+        Padded(-addressable_space, Repeated(false, Length(target)), data)
+    } else {
+        data
+    };
+
+    let m = 2^address_size;
+    Fact(address_size >= 1, "Qubit register must be at least 1.");
+    Fact(Length(data) == m, "Data length must match 2^Length(qs).");
+    let n1 = address_size >>> 1; // Number of qubits in the first half
+    let n2 = address_size - n1; // Number of qubits in the second half
+    let h1 = address[...n1-1]; // Note that h1 will be empty if address_size == 1
+    let h2 = address[n1...];
+    let m1 = 1 <<< n1;
+    let m2 = 1 <<< n2;
+    Fact(m1 * m2 == m, "Length of halves must match total length.");
+
+    // Allocate auxilliary qubits.
+    use aux_qubits1 = Qubit[2^n1 - n1 - 1];
+    use aux_qubits2 = Qubit[2^n2 - n2 - 1];
+
+    // Construct power products for both halves.
+    let products1 = ConstructPowerProducts(h1, aux_qubits1);
+    let products2 = ConstructPowerProducts(h2, aux_qubits2);
+
+    ApplyFlips(data, products1, products2, target);
+
+    // Undo power products of both halves.
+    DestructPowerProducts(products1);
+    DestructPowerProducts(products2);
+}
+
+/// # Summary
+/// Applies flips to the target register based on the data and power products.
+operation ApplyFlips(
+    data : Bool[][],
+    products1 : Qubit[],
+    products2 : Qubit[],
+    target : Qubit[]
+) : Unit {
+    let m1 = Length(products1) + 1;
+    let m2 = Length(products2) + 1;
+
+    for bit_index in IndexRange(target) {
+        let sourceData = Mapped(a -> a[bit_index], data);
+        let flipData = FastMobiusTransform(sourceData);
+        let mask_as_matrix = Chunks(m1, flipData);
+
+        // Apply X to target[bit_index] if the empty product (index 0) is set.
+        if mask_as_matrix[0][0] {
+            X(target[bit_index]);
+        }
+
+        for row in 0..m2-2 {
+            if (mask_as_matrix[row + 1][0]) {
+                CX(products2[row], target[bit_index]);
+            }
+        }
+
+        for col in 0..m1-2 {
+            if (mask_as_matrix[0][col + 1]) {
+                CX(products1[col], target[bit_index]);
+            }
+        }
+
+        for row in 0..m2-2 {
+            for col in 0..m1-2 {
+                if mask_as_matrix[row + 1][col + 1] {
+                    CCNOT(products2[row], products1[col], target[bit_index]);
+                }
+            }
+        }
+
+    }
+}
+
+// =============================
+// Tests
+
+@Test()
+operation CheckLookupViaPP() : Unit {
+    let n = 3;
+    let data = [[true, false, false], [false, true, false], [false, false, true], [false, false, false], [true, true, false], [false, true, true], [true, false, true], [true, true, true]];
+
+    use addr = Qubit[n];
+    use target = Qubit[3];
+
+    // Check that data at all indices is looked up correctly.
+    for i in IndexRange(data) {
+        ApplyXorInPlace(i, addr);
+        LookupViaPP(data, addr, target);
+
+        ApplyPauliFromBitString(PauliX, true, data[i], target);
+        let zero = CheckAllZero(target);
+        Fact(zero, $"Target should match {data[i]} at index {i}.");
+        ResetAll(addr);
+    }
+}
+
+@Test()
+operation CheckLookupViaPPShorterData() : Unit {
+    let n = 3;
+    let width = 3;
+    let data = [[true, false, false], [false, true, false], [false, false, true]];
+
+    use addr = Qubit[n];
+    use target = Qubit[width];
+
+    // Check that shorter data at all indices is looked up correctly.
+    for i in 0..2^n-1 {
+        ApplyXorInPlace(i, addr);
+        LookupViaPP(data, addr, target);
+
+        mutable expected_data = [false, false, false];
+        if i < Length(data) {
+            ApplyPauliFromBitString(PauliX, true, data[i], target);
+            set expected_data = data[i];
+        } else {
+            // For out-of-bounds indices, target should remain |0...0⟩.
+        }
+        let zero = CheckAllZero(target);
+        Fact(zero, $"Target should match {expected_data} at index {i}.");
+        ResetAll(addr);
+    }
+}
+
+@Test()
+operation CheckLookupViaPPLongerData() : Unit {
+    let n = 2;
+    let width = 3;
+    let data = [[true, false, false], [false, true, false], [false, false, true], [false, false, false], [true, true, false], [false, true, true], [true, true, true]];
+
+    use addr = Qubit[n];
+    use target = Qubit[width];
+
+    // Check that longer data at all available indices is looked up correctly.
+    for i in 0..2^n-1 {
+        ApplyXorInPlace(i, addr);
+        LookupViaPP(data, addr, target);
+
+        ApplyPauliFromBitString(PauliX, true, data[i], target);
+        let zero = CheckAllZero(target);
+        Fact(zero, $"Target should match {data[i]} at index {i}.");
+        ResetAll(addr);
+    }
+}
+
+@Test()
+operation TestLookupViaPPMatchesStd() : Unit {
+    let n = 3;
+    let width = 4;
+    let data = [[true, false, false, false], [false, true, false, false], [false, false, true, false], [false, false, false, false], [true, true, false, false], [false, true, true, false], [true, false, true, true], [true, true, true, true]];
+
+    // Use adjoint Std.TableLookup.Select because this check takes adjoint of that.
+    let equal = CheckOperationsAreEqual(
+        n + width,
+        qs => LookupViaPP(data, qs[0..n-1], qs[n...]),
+        qs => Adjoint Std.TableLookup.Select(data, qs[0..n-1], qs[n...])
+    );
+    Fact(equal, "LookupViaPP should match Std.TableLookup.Select.");
+}
+
+@Test()
+operation CheckLookupViaSplitPP() : Unit {
+    let n = 3;
+    let data = [[true, false, false], [false, true, false], [false, false, true], [false, false, false], [true, true, false], [false, true, true], [true, false, true], [true, true, true]];
+
+    use addr = Qubit[n];
+    use target = Qubit[3];
+
+    // Check that data at all indices is looked up correctly.
+    for i in IndexRange(data) {
+        ApplyXorInPlace(i, addr);
+        LookupViaSplitPP(data, addr, target);
+
+        ApplyPauliFromBitString(PauliX, true, data[i], target);
+        let zero = CheckAllZero(target);
+        Fact(zero, $"Target should match {data[i]} at index {i}.");
+        ResetAll(addr);
+    }
+}
+
+@Test()
+operation CheckLookupViaSplitPPShorterData() : Unit {
+    let n = 3;
+    let width = 3;
+    let data = [[true, false, false], [false, true, false], [false, false, true]];
+
+    use addr = Qubit[n];
+    use target = Qubit[width];
+
+    // Check that shorter data at all indices is looked up correctly.
+    for i in 0..2^n-1 {
+        ApplyXorInPlace(i, addr);
+        LookupViaSplitPP(data, addr, target);
+
+        mutable expected_data = [false, false, false];
+        if i < Length(data) {
+            ApplyPauliFromBitString(PauliX, true, data[i], target);
+            set expected_data = data[i];
+        } else {
+            // For out-of-bounds indices, target should remain |0...0⟩.
+        }
+        let zero = CheckAllZero(target);
+        Fact(zero, $"Target should match {expected_data} at index {i}.");
+        ResetAll(addr);
+    }
+}
+
+@Test()
+operation CheckLookupViaSplitPPLongerData() : Unit {
+    let n = 2;
+    let width = 3;
+    let data = [[true, false, false], [false, true, false], [false, false, true], [false, false, false], [true, true, false], [false, true, true], [true, true, true]];
+
+    use addr = Qubit[n];
+    use target = Qubit[width];
+
+    // Check that longer data at all available indices is looked up correctly.
+    for i in 0..2^n-1 {
+        ApplyXorInPlace(i, addr);
+        LookupViaSplitPP(data, addr, target);
+
+        ApplyPauliFromBitString(PauliX, true, data[i], target);
+        let zero = CheckAllZero(target);
+        Fact(zero, $"Target should match {data[i]} at index {i}.");
+        ResetAll(addr);
+    }
+}
+
+@Test()
+operation TestLookupViaSplitPPMatchesStd() : Unit {
+    let n = 3;
+    let width = 4;
+    let data = [[true, false, false, false], [false, true, false, false], [false, false, true, false], [false, false, false, false], [true, true, false, false], [false, true, true, false], [true, false, true, true], [true, true, true, true]];
+
+    // Use adjoint Std.TableLookup.Select because this check takes adjoint of that.
+    let equal = CheckOperationsAreEqual(
+        n + width,
+        qs => LookupViaSplitPP(data, qs[0..n-1], qs[n...]),
+        qs => Adjoint Std.TableLookup.Select(data, qs[0..n-1], qs[n...])
+    );
+    Fact(equal, "LookupViaSplitPP should match Std.TableLookup.Select.");
+}
diff --git a/library/table_lookup/src/Main.qs b/library/table_lookup/src/Main.qs
new file mode 100644
index 0000000000..575f5057b7
--- /dev/null
+++ b/library/table_lookup/src/Main.qs
@@ -0,0 +1,55 @@
+// Copyright (c) Microsoft Corporation.
+// Licensed under the MIT License.
+
+// Main wrapper operation for lookup.
+
+export Lookup.Lookup;
+
+// Options and available algorithms.
+
+export Lookup.LookupOptions;
+export Lookup.DefaultLookupOptions;
+
+export Lookup.DoStdLookup;
+export Lookup.DoMCXLookup;
+export Lookup.DoRecursiveSelectLookup;
+export Lookup.DoPPLookup;
+export Lookup.DoSplitPPLookup;
+
+export Lookup.DoStdUnlookup;
+export Lookup.DoUnlookupViaLookup;
+export Lookup.DoMCXUnlookup;
+export Lookup.DoPPUnlookup;
+export Lookup.DoSplitPPUnlookup;
+
+// Lookup implementations via multicontrolled X gates.
+
+export Multicontrolled.LookupViaMCX;
+export Multicontrolled.BitLookupViaMCX;
+export Multicontrolled.PhaseLookupViaMCX;
+
+// Lookup implementations via recursive SELECT network.
+
+export RecursiveSelect.RecursiveLookup;
+export RecursiveSelect.RecursiveLookupOpt;
+export RecursiveSelect.ControlledRecursiveSelect;
+export RecursiveSelect.ControlledRecursiveSelectOpt;
+
+// Lookup implementations via power products.
+
+export PowerProducts.GetAuxCountForPP;
+export PowerProducts.ConstructPowerProducts;
+export PowerProducts.DestructPowerProducts;
+
+export LookupViaPP.LookupViaPP;
+export LookupViaPP.LookupViaSplitPP;
+export PhaseLookup.PhaseLookupViaPP;
+export PhaseLookup.PhaseLookupViaSplitPP;
+
+// Utility functions.
+
+export Utils.FastMobiusTransform;
+export Utils.MeasureAndComputePhaseData;
+export Utils.BinaryInnerProduct;
+export Utils.CombineControls;
+export Utils.GetCombinedControl;
diff --git a/library/table_lookup/src/Multicontrolled.qs b/library/table_lookup/src/Multicontrolled.qs
new file mode 100644
index 0000000000..e01fea1aa5
--- /dev/null
+++ b/library/table_lookup/src/Multicontrolled.qs
@@ -0,0 +1,198 @@
+// Copyright (c) Microsoft Corporation.
+// Licensed under the MIT License.
+
+// Simple implementations of lookup operations using multicontrolled X gates.
+// Data shorter or longer than addressable space is allowed:
+// * Longer data is ignored.
+// * Shorter data is treated as if it is padded with false values.
+// Little-endian format is used throughout.
+
+import Std.Arrays.IndexRange;
+import Std.Diagnostics.*;
+import Std.Math.*;
+
+// Lookup of a single bit using multicontrolled X gates.
+operation BitLookupViaMCX(data : Bool[], address : Qubit[], target : Qubit) : Unit is Adj + Ctl {
+    let address_size = Length(address);
+    let addressable_space = 2^address_size;
+    let data_length = Length(data);
+    for basis_vector in 0..MinI(data_length, addressable_space)-1 {
+        if data[basis_vector] {
+            within {
+                // Invert address qubits for 0-es in basis_vector.
+                ApplyXorInPlace(addressable_space-1-basis_vector, address)
+            } apply {
+                Controlled X(address, target);
+            }
+        }
+    }
+}
+
+// Phase lookup of a single bit using multicontrolled X gates.
+operation PhaseLookupViaMCX(data : Bool[], address : Qubit[]) : Unit is Adj + Ctl {
+    use aux = Qubit();
+    within {
+        X(aux);
+        H(aux);
+    } apply {
+        BitLookupViaMCX(data, address, aux);
+    }
+}
+
+// Lookup of mult-bit register using multicontrolled X gates.
+operation LookupViaMCX(data : Bool[][], address : Qubit[], target : Qubit[]) : Unit is Adj + Ctl {
+    let address_size = Length(address);
+    let addressable_space = 2^address_size;
+    let data_length = Length(data);
+    let target_size = Length(target);
+    for basis_vector in 0..MinI(data_length, addressable_space)-1 {
+        let data_vector = data[basis_vector];
+        Fact(Length(data_vector) == target_size, $"Data vector length {Length(data_vector)} must match target size {target_size}.");
+        within {
+            // Invert address qubits for 0-es in basis_vector.
+            ApplyXorInPlace(addressable_space-1-basis_vector, address)
+        } apply {
+            Controlled ApplyPauliFromBitString(address, (PauliX, true, data_vector, target));
+        }
+    }
+}
+
+// =============================
+// Tests
+
+@Test()
+operation CheckLookupViaMCX() : Unit {
+    let n = 3;
+    let data = [[true, false, false], [false, true, false], [false, false, true], [false, false, false], [true, true, false], [false, true, true], [true, false, true], [true, true, true]];
+
+    use addr = Qubit[n];
+    use target = Qubit[3];
+
+    // Check that data at all indices is looked up correctly.
+    for i in IndexRange(data) {
+        ApplyXorInPlace(i, addr);
+        LookupViaMCX(data, addr, target);
+
+        ApplyPauliFromBitString(PauliX, true, data[i], target);
+        let zero = CheckAllZero(target);
+        Fact(zero, $"Target should match {data[i]} at index {i}.");
+        ResetAll(addr);
+    }
+}
+
+@Test()
+operation CheckLookupViaMCXShorterData() : Unit {
+    let n = 3;
+    let width = 3;
+    let data = [[true, false, false], [false, true, false], [false, false, true]];
+
+    use addr = Qubit[n];
+    use target = Qubit[width];
+
+    // Check that shorter data at all indices is looked up correctly.
+    for i in 0..2^n-1 {
+        ApplyXorInPlace(i, addr);
+        LookupViaMCX(data, addr, target);
+
+        mutable expected_data = [false, false, false];
+        if i < Length(data) {
+            ApplyPauliFromBitString(PauliX, true, data[i], target);
+            set expected_data = data[i];
+        } else {
+            // For out-of-bounds indices, target should remain |0...0⟩.
+        }
+        let zero = CheckAllZero(target);
+        Fact(zero, $"Target should match {expected_data} at index {i}.");
+        ResetAll(addr);
+    }
+}
+
+@Test()
+operation CheckLookupViaMCXLongerData() : Unit {
+    let n = 2;
+    let width = 3;
+    let data = [[true, false, false], [false, true, false], [false, false, true], [false, false, false], [true, true, false], [false, true, true], [true, true, true]];
+
+    use addr = Qubit[n];
+    use target = Qubit[width];
+
+    // Check that longer data at all available indices is looked up correctly.
+    for i in 0..2^n-1 {
+        ApplyXorInPlace(i, addr);
+        LookupViaMCX(data, addr, target);
+
+        ApplyPauliFromBitString(PauliX, true, data[i], target);
+        let zero = CheckAllZero(target);
+        Fact(zero, $"Target should match {data[i]} at index {i}.");
+        ResetAll(addr);
+    }
+}
+
+@Test()
+operation CheckBitLookupViaMCX() : Unit {
+    let n = 4;
+    let data = [true, false, true, false, false, false, false, false, false, false, false, false, false, false, true, true];
+
+    use addr = Qubit[n];
+    use target = Qubit();
+
+    // Check that data at all indices is looked up correctly.
+    for i in IndexRange(data) {
+        ApplyXorInPlace(i, addr);
+        BitLookupViaMCX(data, addr, target);
+
+        let value = Std.Convert.ResultAsBool(MResetZ(target));
+        ResetAll(addr);
+
+        Fact(value == data[i], $"Target qubit measurement mismatch at index {i}.");
+    }
+}
+
+@Test()
+operation TestPhaseLookupViaMCX() : Unit {
+    let n = 4;
+    let data = [true, false, true, false, false, false, false, false, false, false, false, false, false, false, true, true];
+    let coeffs = [-0.25, 0.25, -0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, -0.25, -0.25];
+
+    use qs = Qubit[n];
+    ApplyToEach(H, qs);
+
+    // `Reversed` to match big-endian state preparation coefficients order.
+    PhaseLookupViaMCX(data, Std.Arrays.Reversed(qs));
+    Adjoint Std.StatePreparation.PreparePureStateD(coeffs, qs);
+
+    Fact(CheckAllZero(qs), "All qubits should be back to |0⟩ state.");
+}
+
+@Test()
+operation TestBitLookupViaMCXMatchesStd() : Unit {
+    let n = 4;
+    let data = [true, false, true, false, false, false, false, false, true, false, false, true, false, false, true, true];
+    let select_data = Std.Arrays.Mapped(x -> [x], data);
+
+    // Use adjoint Std.TableLookup.Select because this check takes adjoint of that.
+    let equal = CheckOperationsAreEqual(
+        n + 1,
+        qs => BitLookupViaMCX(data, qs[0..n-1], qs[n]),
+        qs => Adjoint Std.TableLookup.Select(select_data, qs[0..n-1], qs[n..n])
+    );
+    Fact(equal, "BitLookupViaMCX should match Std.TableLookup.Select.");
+}
+
+@Test()
+operation TestLookupViaMCXMatchesStd() : Unit {
+    let n = 3;
+    let width = 4;
+    let data = [[true, false, false, false], [false, true, false, false], [false, false, true, false], [false, false, false, false], [true, true, false, false], [false, true, true, false], [true, false, true, true], [true, true, true, true]];
+
+    use addr = Qubit[n];
+    use target = Qubit[width];
+
+    // Use adjoint Std.TableLookup.Select because this check takes adjoint of that.
+    let equal = CheckOperationsAreEqual(
+        n + width,
+        qs => LookupViaMCX(data, qs[0..n-1], qs[n...]),
+        qs => Adjoint Std.TableLookup.Select(data, qs[0..n-1], qs[n...])
+    );
+    Fact(equal, "LookupViaMCX should match Std.TableLookup.Select.");
+}
diff --git a/library/table_lookup/src/PhaseLookup.qs b/library/table_lookup/src/PhaseLookup.qs
new file mode 100644
index 0000000000..723fdca058
--- /dev/null
+++ b/library/table_lookup/src/PhaseLookup.qs
@@ -0,0 +1,193 @@
+// Copyright (c) Microsoft Corporation.
+// Licensed under the MIT License.
+
+import Std.Diagnostics.*;
+import Std.Arrays.*;
+import Std.Math.*;
+import Std.Convert.*;
+
+import PowerProducts.*;
+import Utils.*;
+
+/// # Summary
+/// Implements phaseup operation using power products without address split.
+operation PhaseLookupViaPP(address : Qubit[], data : Bool[]) : Unit {
+    let data_length = Length(data);
+    let address_size = Length(address);
+    let addressable_space = 1 <<< address_size;
+    let data = if (addressable_space < data_length) {
+        data[...addressable_space-1]
+    } elif (addressable_space > data_length) {
+        Padded(-addressable_space, false, data)
+    } else {
+        data
+    };
+    use aux_qubits = Qubit[GetAuxCountForPP(address_size)];
+    // Transform data from minterm coefficients to polynomial coefficients.
+    let corrections = FastMobiusTransform(data);
+    let products = ConstructPowerProducts(address, aux_qubits);
+    ApplyPhasingViaZ(products, corrections);
+    DestructPowerProducts(products);
+}
+
+operation ApplyPhasingViaZ(qs : Qubit[], mask : Bool[]) : Unit {
+    Fact(Length(mask) > 0, "Mask must be a non-empty array.");
+    Fact(Length(mask) == Length(qs) + 1, "Mask row count must match qs length.");
+
+    // Ignore the first element of mask, it affects the global phase.
+    ApplyPauliFromBitString(PauliZ, true, Std.Arrays.Rest(mask), qs);
+}
+
+/// # Summary
+/// Invert phases of `qs` basis states according to the provided boolean array.
+/// If `data[i]` is `true`, the phase of |i⟩ gets is inverted (multiplied by -1).
+/// Qubit register `qs` is expected to be in little-endian order.
+///
+/// # Description
+/// This operation implements phase lookup using power products and address split.
+/// It is a Q# implementation of the "phaseup" operation from the referenced paper.
+/// This operation assumes that `Length(data)` matches `2^Length(qs)`.
+///
+/// # Input
+/// ## qs
+/// Qubit register whose basis states will have their phases inverted.
+///
+/// ## data
+/// Boolean array indicating which basis states to invert. If `data[i]` is `true`,
+/// the phase of |i⟩ gets inverted (multiplied by -1).
+///
+/// # Reference
+/// 1. [arXiv:2505.15917](https://arxiv.org/abs/2505.15917)
+///    "How to factor 2048 bit RSA integers with less than a million noisy qubits"
+///    by Craig Gidney, May 2025. Section A.3.
+operation PhaseLookupViaSplitPP(address : Qubit[], data : Bool[]) : Unit {
+    let data_length = Length(data);
+    let address_size = Length(address);
+    let addressable_space = 1 <<< address_size;
+    let data = if (addressable_space < data_length) {
+        data[...addressable_space-1]
+    } elif (addressable_space > data_length) {
+        Padded(-addressable_space, false, data)
+    } else {
+        data
+    };
+
+    Fact(address_size >= 1, "Qubit register must be at least 1.");
+    Fact(Length(data) == addressable_space, "Data length must match 2^Length(qs).");
+    let n1 = address_size >>> 1; // Number of qubits in the first half
+    let n2 = address_size - n1; // Number of qubits in the second half
+    let address_low = address[...n1-1]; // Note that address_low will be empty if n == 1
+    let address_high = address[n1...];
+    let m1 = 1 <<< n1;
+    let m2 = 1 <<< n2;
+    Fact(m1 * m2 == addressable_space, "Length of halves must match total length.");
+
+    // Allocate auxilliary qubits.
+    use aux_qubits1 = Qubit[GetAuxCountForPP(n1)];
+    use aux_qubits2 = Qubit[GetAuxCountForPP(n2)];
+
+    // Construct power products for both halves.
+    let products1 = ConstructPowerProducts(address_low, aux_qubits1);
+    let products2 = ConstructPowerProducts(address_high, aux_qubits2);
+
+    // Convert data from minterm to monomial basis using Fast Möbius Transform.
+    // and chunk it into a matrix
+    let mask_as_matrix = Chunks(m1, FastMobiusTransform(data));
+
+    // Apply phasing within each half and between halves.
+    ApplyPhasingViaZandCZ(products1, products2, mask_as_matrix);
+
+    // Undo power products of both halves.
+    DestructPowerProducts(products1);
+    DestructPowerProducts(products2);
+}
+
+/// # Summary
+/// Applies phase corrections using Z and CZ gates based on power product coefficients.
+/// This is the core quantum operation in the address-split phase lookup algorithm.
+///
+/// # Description
+/// This operation applies conditional phase flips based on a 2D mask that represents
+/// power product coefficients after Fast Möbius Transform. The algorithm treats the
+/// input qubits as split into two halves, with separate power products for each half.
+///
+/// The phase correction is applied as follows:
+/// 1. Apply Z gates to products2 based on products1[0] (for products from first half only).
+/// 2. Apply Z gates to products1 based on products2[0] (for products from second half only).
+/// 3. Apply CZ gates between corresponding products from both halves.
+///
+/// # Input
+/// ## products1
+/// Power product qubits from the first half of the address register.
+///
+/// ## products2
+/// Power product qubits from the second half of the address register.
+///
+/// ## mask
+/// 2D boolean array containing power product coefficients.
+/// - `mask[i][j]` indicates whether to apply phase correction for the product
+///   of subset i from second half and subset j from first half.
+///
+/// # Remarks
+/// The mask is obtained by applying Fast Möbius Transform to phase data
+/// and reshaping into a 2D matrix. This allows efficient quantum evaluation of
+/// the phase function using O(2^(n/2)) quantum resources instead of O(2^n).
+operation ApplyPhasingViaZandCZ(products1 : Qubit[], products2 : Qubit[], mask : Bool[][]) : Unit {
+    Fact(Length(mask) > 0, "Mask must be a non-empty array.");
+    Fact(Length(mask) == Length(products2) + 1, "Mask row count must match products2 length.");
+    Fact(Length(mask[0]) == Length(products1) + 1, "Mask column count must match products1 length.");
+
+    // ColumnAt(0, mask) doesn't correspond to any qubits from the first half,
+    // so we can apply Z (rather than CZ) based on mask values.
+    ApplyPauliFromBitString(PauliZ, true, Rest(ColumnAt(0, mask)), products2);
+
+    // mask[0] row doesn't correspond to any qubits from the second half,
+    // so we can apply Z (rather than CZ) based on mask values.
+    ApplyPauliFromBitString(PauliZ, true, Rest(mask[0]), products1);
+
+    // From the second row on, take control from the first half and apply
+    // masked multi-target CZ gates via Controlled ApplyPauliFromBitString.
+    for row in IndexRange(products1) {
+        Controlled ApplyPauliFromBitString(
+            [products1[row]],
+            (PauliZ, true, Rest(ColumnAt(row + 1, mask)), products2)
+        );
+    }
+}
+
+// =============================
+// Tests
+
+@Test()
+operation TestPhaseLookupViaPPandZ() : Unit {
+    let address_size = 3;
+    let data_length = 2^address_size;
+    let data_value_count = 2^data_length;
+
+    for i in 0..data_value_count-1 {
+        let data = Std.Convert.IntAsBoolArray(i, data_length);
+        let same = CheckOperationsAreEqual(
+            address_size,
+            PhaseLookupViaPP(_, data),
+            Multicontrolled.PhaseLookupViaMCX(data, _)
+        );
+        Fact(same, $"PhaseLookupViaPPandZ must be the same as PhaseLookupViaMCX for {data}.");
+    }
+}
+
+@Test()
+operation TestPhaseLookupViaPPandCZ() : Unit {
+    let address_size = 3;
+    let data_length = 2^address_size;
+    let data_value_count = 2^data_length;
+
+    for i in 0..data_value_count-1 {
+        let data = Std.Convert.IntAsBoolArray(i, data_length);
+        let same = CheckOperationsAreEqual(
+            address_size,
+            PhaseLookupViaSplitPP(_, data),
+            Multicontrolled.PhaseLookupViaMCX(data, _)
+        );
+        Fact(same, $"PhaseLookupViaPPandCZ must be the same as PhaseLookupViaMCX for {data}.");
+    }
+}
diff --git a/library/table_lookup/src/PowerProducts.qs b/library/table_lookup/src/PowerProducts.qs
new file mode 100644
index 0000000000..fcb632faf9
--- /dev/null
+++ b/library/table_lookup/src/PowerProducts.qs
@@ -0,0 +1,181 @@
+// Copyright (c) Microsoft Corporation.
+// Licensed under the MIT License.
+
+import Std.Arrays.IndexRange;
+import Std.Diagnostics.*;
+
+/// # Summary
+/// Constructs power products - AND-ed subsets of qubits from the input register `qs`.
+/// `2^Length(qs) - 1` qubits corresponding to non-empty subsets of `qs` are placed into the result array.
+///
+/// # Description
+/// Resulting subsets correspond to an integer index that runs from `1` to `(2^Length(qs))-1`.
+/// (Since the empty set (index 0) is not included in the result, actual array indexes should be shifted.)
+/// Indexes are treated as bitmasks indicating if a particular qubit is included.
+/// Bitmasks `2^i` includes only qubit `qs[i]`, which is placed into the resulting array at index 2^i - 1.
+/// Bitmasks with more than one bit set correspond to subsets with multiple qubits from `qs`.
+/// Qubits for these masks are taken from aux_qubits register and their value is set using AND gates.
+/// Note:
+///     1. Empty set is not included in the result.
+///     2. For sets that only contain one qubit, the input qubits are reused.
+///
+/// # Alt summary
+/// Takes a register of qubits and returns "power products" - qubits corresponding to all non-empty subsets
+/// of the qubits from the input register: each power product qubit state is a result of AND operation
+/// for the qubits in corresponding subset.
+operation ConstructPowerProducts(qubits : Qubit[], aux_qubits : Qubit[]) : Qubit[] {
+    // Start with empty array - no dummy qubit for empty set.
+    mutable power_products = [];
+    // Index to take next free qubit from aux_qubits array.
+    mutable next_available = 0;
+    // Consider every index in the input qubit register.
+    for qubit_index in IndexRange(qubits) {
+        // First, add the set that consists of only one qubit at index qubit_index.
+        power_products += qubits[qubit_index..qubit_index];
+        // Then, construct and add sets that include this new qubit as the last one.
+        for existing_set_index in 0..Length(power_products)-2 {
+            // Take the next qubit for the new set.
+            let next_power_product = aux_qubits[next_available];
+            next_available += 1;
+            // Create appropriate set and add it to the result.
+            AND(power_products[existing_set_index], qubits[qubit_index], next_power_product);
+            power_products += [next_power_product];
+        }
+    }
+    Fact(next_available == Length(aux_qubits), "ConstructPowerProducts: All auxilliary qubits should be used.");
+    return power_products;
+}
+
+/// # Summary
+/// Uncomputes construction of power products done by `ConstructPowerProducts`.
+/// Pass array returned by `ConstructPowerProducts` to this function
+/// to reset auxiliary qubits used to hold power products back to |0⟩ state.
+///
+/// # Description
+/// `products` array has no qubit that corresponds to an empty product (≡1).
+/// All entries at indexes `2^i - 1` contain original qubits.
+/// Qubits from `2^i - 1` to `2^(i+1) - 2` represent power products that
+/// end in original qubit at `2^i - 1`.
+/// To undo power products this function goes over original qubits backwards.
+/// Then measures out qubits from `2^i - 1` to `2^(i+1) - 2` in X basis,
+/// targeting corresponding qubits from 0 to `2^i - 2` in CZ gates if necessary.
+operation DestructPowerProducts(products : Qubit[]) : Unit {
+    let length = Length(products);
+    if length <= 1 {
+        // Nothing to undo - this was one of the source qubits.
+        return ();
+    }
+    // Adjust for empty set.
+    let extended_len = length + 1;
+    Fact((extended_len &&& length) == 0, "DestructPowerProducts: Length + 1 of a qubit register should be a power of 2");
+
+    // At index h-1 a source qubit is located (shifted by 1 to account for the lack of empty set).
+    // To the right are all power products ending in it.
+    // We are going backwards over all original qubits.
+    mutable h = extended_len / 2;
+    // If h is 1 we have nothing else to undo.
+    while h > 1 {
+        // Go over all sets that end in original qubit currently at index h-1.
+        // NOTE: The order of targets here doesn't matter.
+        for k in 0..h-2 {
+            // Measure and reset the qubit that represents
+            // the set (h-1) | k, which is at index h-1+k+1 = h+k.
+            if MResetX(products[h + k]) == One {
+                // If we measure 1, qubit representing set k needs to be included in targets.
+                CZ(products[h - 1], products[k]);
+            }
+        }
+        // Done with qubit at index h-1. Go to next original qubit.
+        h = h / 2;
+    }
+}
+
+function GetAuxCountForPP(nQubits : Int) : Int {
+    Fact(nQubits >= 0, "Number of qubits for power product construction must be non-negative.");
+    // Number of power products is 2^n - 1 (this excludes the empty product).
+    // Number of original qubits is n.
+    // Aux qubits needed is (2^n - 1) - n = 2^n - n - 1.
+    (1 <<< nQubits) - nQubits - 1
+}
+
+// =============================
+// Tests
+
+internal operation ConstructDestructPowerProducts(qs : Qubit[]) : Unit {
+    // For monomials with more than one variable we need auxilliary qubits.
+    use aux_qubits = Qubit[GetAuxCountForPP(Length(qs))];
+
+    // Construct/destruct should leave qs unchanged.
+    let products = ConstructPowerProducts(qs, aux_qubits);
+    DestructPowerProducts(products);
+}
+
+@Test()
+operation TestCreateDestructPowerProducts() : Unit {
+    // Check that construction and destruction of power products does not affect the register.
+    for i in 0..5 {
+        let success = CheckOperationsAreEqual(
+            i,
+            qs => ConstructDestructPowerProducts(qs),
+            qs => {}
+        );
+        Fact(success, $"Construction/Destruction of power products must be identity for {i} qubits.");
+    }
+}
+
+internal operation CheckPowerProducts(nQubits : Int, address_value : Int) : Unit {
+    // Prepare qubit register.
+    Fact(nQubits >= 0, "Number of qubits must be non-negative.");
+    use qs = Qubit[nQubits];
+    let address_space = 1 <<< nQubits;
+
+    // Prepare random basis state in qs.
+    Fact(address_value >= 0 and address_value < address_space, "Value must fit in the number of qubits.");
+    let state = Std.Convert.IntAsBoolArray(address_value, nQubits);
+    ApplyPauliFromBitString(PauliX, true, state, qs);
+
+    // Construct power products.
+    use aux_qubits = Qubit[GetAuxCountForPP(nQubits)];
+    let products = ConstructPowerProducts(qs, aux_qubits);
+    Fact(Length(products) == address_space - 1, $"Power product length should be {address_space - 1}.");
+
+    // Verify that each product qubit is correct.
+    for index in 0..address_space-2 {
+        // Shift by 1 since empty product is not included.
+        let monomial_index = index + 1;
+        mutable expected_value = true;
+        for bit_position in 0..nQubits-1 {
+            if ((monomial_index &&& (1 <<< bit_position)) != 0) {
+                // This qubit is included in the product.
+                set expected_value = expected_value and state[bit_position];
+            }
+        }
+        within {
+            if (expected_value) {
+                // Invert if expected value is 1 - we'll check for |0⟩ state.
+                X(products[index]);
+            }
+        } apply {
+            Fact(CheckZero(products[index]), $"Power product at index {index} should match expected value {expected_value}.");
+        }
+    }
+
+    // Destruct power products to reset aux qubits.
+    DestructPowerProducts(products);
+
+    // Reset original qubits.
+    ApplyPauliFromBitString(PauliX, true, state, qs);
+
+    // All qubits should be back to |0⟩ state at this point.
+}
+
+@Test()
+operation TestPowerProductsExhaustive() : Unit {
+    // Test power products construction for various numbers of qubits and basis states.
+    for nQubits in 0..5 {
+        let address_space = 1 <<< nQubits;
+        for value in 0..address_space-1 {
+            CheckPowerProducts(nQubits, value);
+        }
+    }
+}
diff --git a/library/table_lookup/src/RecursiveSelect.qs b/library/table_lookup/src/RecursiveSelect.qs
new file mode 100644
index 0000000000..aa21201171
--- /dev/null
+++ b/library/table_lookup/src/RecursiveSelect.qs
@@ -0,0 +1,368 @@
+// Copyright (c) Microsoft Corporation.
+// Licensed under the MIT License.
+
+import Std.Arrays.*;
+import Std.Diagnostics.*;
+import Std.Math.*;
+import Std.Convert.*;
+
+
+/// # Summary
+/// Performs table lookup using a SELECT network respecting longer addresses.
+///
+/// # 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.
+/// The implementation of the SELECT network is based on unary encoding as
+/// presented in [1]. The recursive implementation differs from the one
+/// presented in [2] by allowing addresses beyond the length of `data`.
+///
+/// # Input
+/// ## data
+/// The classical table lookup data which is prepared in `target` with
+/// respect to the state in `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
+///
+/// # References
+/// 1. [arXiv:1805.03662](https://arxiv.org/abs/1805.03662)
+///    "Encoding Electronic Spectra in Quantum Circuits with Linear T
+///    Complexity", Section A.
+/// 2. [arXiv:2211.01133](https://arxiv.org/abs/2211.01133)
+///    "Space-time optimized table lookup", Section 2.
+operation RecursiveLookup(
+    useAnd : Bool,
+    data : Bool[][],
+    address : Qubit[],
+    target : Qubit[]
+) : Unit {
+    let data_length = Length(data);
+    if data_length == 0 {
+        return ();
+    }
+    let address_size = Length(address);
+    if address_size == 0 {
+        ApplyPauliFromBitString(PauliX, true, data[0], target);
+        return ();
+    }
+    let addressable_space = 1 <<< address_size;
+    let data_length = MinI(data_length, addressable_space);
+    let data = data[...data_length - 1];
+    let highest_address_qubit = Tail(address);
+    let lower_address_qubits = Most(address);
+    let parts = Partitioned([MinI(addressable_space / 2, data_length)], data);
+
+    within {
+        X(highest_address_qubit);
+    } apply {
+        ControlledRecursiveSelect(useAnd, highest_address_qubit, parts[0], lower_address_qubits, target);
+    }
+    ControlledRecursiveSelect(useAnd, highest_address_qubit, parts[1], lower_address_qubits, target);
+}
+
+/// # Summary
+/// Performs table lookup using a SELECT network assuming no addresses beyond data length.
+operation RecursiveLookupOpt(
+    useAnd : Bool,
+    data : Bool[][],
+    address : Qubit[],
+    target : Qubit[]
+) : Unit {
+    let data_length = Length(data);
+    if data_length == 0 {
+        // If there's no data, there's nothing to apply.
+        return ();
+    }
+    if data_length == 1 {
+        // Base case: always apply data value if data length is 1.
+        // This version doesn't support address values beyond data length and some value needs to be applied.
+        // Since this is the only data value, it is the one to be applied.
+        ApplyPauliFromBitString(PauliX, true, data[0], target);
+        return ();
+    }
+    let addressable_space = 1 <<< Length(address);
+    if addressable_space == 1 {
+        return ();
+    }
+
+    let data_length = MinI(data_length, addressable_space);
+    let data = data[...data_length - 1];
+    let address_size_needed = BitSizeI(data_length - 1);
+    let (lower_address_qubits, highest_address_qubit) = MostAndTail(address[...address_size_needed - 1]);
+    let parts = Partitioned([2^(address_size_needed - 1)], data);
+
+    within {
+        X(highest_address_qubit);
+    } apply {
+        ControlledRecursiveSelectOpt(useAnd, highest_address_qubit, parts[0], lower_address_qubits, target);
+    }
+    ControlledRecursiveSelectOpt(useAnd, highest_address_qubit, parts[1], lower_address_qubits, target);
+}
+
+// Complete version of recursive select network that ignores address values
+// beyond data length. This is equivalent to padding data with false values
+// to cover the entire addressable space.
+// If data length is 1, single data value is used only if address is zero.
+operation ControlledRecursiveSelect(
+    useAnd : Bool,
+    control : Qubit,
+    data : Bool[][],
+    address : Qubit[],
+    target : Qubit[]
+) : Unit {
+    let data_length = Length(data);
+    if (data_length == 0) {
+        // If there's no data, there's nothing to do.
+        return ();
+    }
+
+    let address_size = Length(address);
+    if address_size == 0 {
+        // Base case. Use CX on qubits where data is true.
+        Fact(data_length == 1, "Data length must be 1 when address size is 0.");
+        Controlled ApplyPauliFromBitString([control], (PauliX, true, data[0], target));
+        return ();
+    }
+
+    let address_space = 1 <<< address_size;
+    Fact(data_length <= address_space, "Data length must not exceed addressable space.");
+
+    let highest_address_qubit = Tail(address);
+    let lower_address_qubits = Most(address);
+    let data_parts = Partitioned([address_space / 2], data);
+
+    use aux = Qubit();
+    within {
+        X(highest_address_qubit);
+    } apply {
+        if useAnd {
+            AND(control, highest_address_qubit, aux);
+        } else {
+            CCNOT(control, highest_address_qubit, aux);
+        }
+    }
+    ControlledRecursiveSelect(useAnd, aux, data_parts[0], lower_address_qubits, target);
+    CNOT(control, aux);
+    ControlledRecursiveSelect(useAnd, aux, data_parts[1], lower_address_qubits, target);
+    if useAnd {
+        Adjoint AND(control, highest_address_qubit, aux);
+    } else {
+        Adjoint CCNOT(control, highest_address_qubit, aux);
+    }
+}
+
+// Optimized version of recursive select network that expects all address values
+// to be within data length. If address value exceeds data length, behavior is undefined.
+// If data length is 1, single data value is always used.
+operation ControlledRecursiveSelectOpt(
+    useAnd : Bool,
+    control : Qubit,
+    data : Bool[][],
+    address : Qubit[],
+    target : Qubit[]
+) : Unit {
+    let data_length = Length(data);
+    Fact(data_length > 0, "ControlledRecursiveSelectOpt: Data cannot be empty.");
+
+    let address_size_needed = BitSizeI(data_length - 1);
+    Fact(Length(address) >= address_size_needed, "ControlledRecursiveSelectOpt: Address register is too short.");
+
+    if data_length == 1 {
+        // Base case: always apply data value if data length is 1.
+        Controlled ApplyPauliFromBitString([control], (PauliX, true, data[0], target));
+    } else {
+        use helper = Qubit();
+
+        // Get just enough address qubits to address all data and split data.
+        let (lower_address_qubits, highest_address_qubit) = MostAndTail(address[...address_size_needed - 1]);
+        let parts = Partitioned([1 <<< (address_size_needed - 1)], data);
+
+        within {
+            X(highest_address_qubit);
+        } apply {
+            if useAnd {
+                AND(control, highest_address_qubit, helper);
+            } else {
+                CCNOT(control, highest_address_qubit, helper);
+            }
+        }
+        ControlledRecursiveSelectOpt(useAnd, helper, parts[0], lower_address_qubits, target);
+        CNOT(control, helper);
+        ControlledRecursiveSelectOpt(useAnd, helper, parts[1], lower_address_qubits, target);
+        if useAnd {
+            Adjoint AND(control, highest_address_qubit, helper);
+        } else {
+            Adjoint CCNOT(control, highest_address_qubit, helper);
+        }
+    }
+}
+
+// =============================
+// Tests
+
+@Test()
+operation CheckRecursiveLookup() : Unit {
+    let n = 3;
+    let data = [[true, false, false], [false, true, false], [false, false, true], [false, false, false], [true, true, false], [false, true, true], [true, false, true], [true, true, true]];
+
+    use addr = Qubit[n];
+    use target = Qubit[3];
+
+    // Check that data at all indices is looked up correctly.
+    for i in IndexRange(data) {
+        ApplyXorInPlace(i, addr);
+        RecursiveLookup(true, data, addr, target);
+
+        ApplyPauliFromBitString(PauliX, true, data[i], target);
+        let zero = CheckAllZero(target);
+        Fact(zero, $"Target should match {data[i]} at index {i}.");
+        ResetAll(addr);
+    }
+}
+
+@Test()
+operation CheckRecursiveLookupOpt() : Unit {
+    let n = 3;
+    let data = [[true, false, false], [false, true, false], [false, false, true], [false, false, false], [true, true, false], [false, true, true], [true, false, true], [true, true, true]];
+
+    use addr = Qubit[n];
+    use target = Qubit[3];
+
+    // Check that data at all indices is looked up correctly.
+    for i in IndexRange(data) {
+        ApplyXorInPlace(i, addr);
+        RecursiveLookupOpt(true, data, addr, target);
+
+        ApplyPauliFromBitString(PauliX, true, data[i], target);
+        let zero = CheckAllZero(target);
+        Fact(zero, $"Target should match {data[i]} at index {i}.");
+        ResetAll(addr);
+    }
+}
+
+@Test()
+operation CheckRecursiveLookupShorterData() : Unit {
+    let n = 3;
+    let width = 3;
+    let data = [[true, false, false], [false, true, false], [false, false, true]];
+
+    use addr = Qubit[n];
+    use target = Qubit[width];
+
+    // Check that shorter data at all indices is looked up correctly.
+    // This works for all addresses even beyond data length.
+    for i in 0..2^n-1 {
+        ApplyXorInPlace(i, addr);
+        RecursiveLookup(true, data, addr, target);
+
+        mutable expected_data = [false, false, false];
+        if i < Length(data) {
+            ApplyPauliFromBitString(PauliX, true, data[i], target);
+            set expected_data = data[i];
+        } else {
+            // For out-of-bounds indices, target should remain |0...0⟩.
+        }
+        let zero = CheckAllZero(target);
+        Fact(zero, $"Target should match {expected_data} at index {i}.");
+        ResetAll(addr);
+    }
+}
+
+@Test()
+operation CheckRecursiveLookupShorterDataOpt() : Unit {
+    let n = 3;
+    let width = 3;
+    let data = [[true, false, false], [false, true, false], [false, false, true]];
+
+    use addr = Qubit[n];
+    use target = Qubit[width];
+
+    // Check that shorter data at all indices is looked up correctly.
+    // This only works up to data length.
+    for i in IndexRange(data) {
+        ApplyXorInPlace(i, addr);
+        RecursiveLookupOpt(true, data, addr, target);
+
+        ApplyPauliFromBitString(PauliX, true, data[i], target);
+        let expected_data = data[i];
+        let zero = CheckAllZero(target);
+        Fact(zero, $"Target should match {expected_data} at index {i}.");
+        ResetAll(addr);
+    }
+}
+
+@Test()
+operation CheckRecursiveLookupLongerData() : Unit {
+    let n = 2;
+    let width = 3;
+    let data = [[true, false, false], [false, true, false], [false, false, true], [false, false, false], [true, true, false], [false, true, true], [true, true, true]];
+
+    use addr = Qubit[n];
+    use target = Qubit[width];
+
+    // Check that longer data at all available indices is looked up correctly.
+    for i in 0..2^n-1 {
+        ApplyXorInPlace(i, addr);
+        RecursiveLookup(true, data, addr, target);
+
+        ApplyPauliFromBitString(PauliX, true, data[i], target);
+        let zero = CheckAllZero(target);
+        Fact(zero, $"Target should match {data[i]} at index {i}.");
+        ResetAll(addr);
+    }
+}
+
+@Test()
+operation CheckRecursiveLookupLongerDataOpt() : Unit {
+    let n = 2;
+    let width = 3;
+    let data = [[true, false, false], [false, true, false], [false, false, true], [false, false, false], [true, true, false], [false, true, true], [true, true, true]];
+
+    use addr = Qubit[n];
+    use target = Qubit[width];
+
+    // Check that longer data at all available indices is looked up correctly.
+    for i in 0..2^n-1 {
+        ApplyXorInPlace(i, addr);
+        RecursiveLookupOpt(true, data, addr, target);
+
+        ApplyPauliFromBitString(PauliX, true, data[i], target);
+        let zero = CheckAllZero(target);
+        Fact(zero, $"Target should match {data[i]} at index {i}.");
+        ResetAll(addr);
+    }
+}
+
+@Test()
+operation TestRecursiveLookupMatchesStd() : Unit {
+    let n = 3;
+    let width = 4;
+    let data = [[true, false, false, false], [false, true, false, false], [false, false, true, false], [false, false, false, false], [true, true, false, false], [false, true, true, false], [true, false, true, true], [true, true, true, true]];
+
+    // Use adjoint Std.TableLookup.Select because this check takes adjoint of that.
+    let equal = CheckOperationsAreEqual(
+        n + width,
+        qs => RecursiveLookup(true, data, qs[0..n-1], qs[n...]),
+        qs => Adjoint Std.TableLookup.Select(data, qs[0..n-1], qs[n...])
+    );
+    Fact(equal, "RecursiveLookup should match Std.TableLookup.Select.");
+}
+
+@Test()
+operation TestRecursiveLookupMatchesStdOpt() : Unit {
+    let n = 3;
+    let width = 4;
+    let data = [[true, false, false, false], [false, true, false, false], [false, false, true, false], [false, false, false, false], [true, true, false, false], [false, true, true, false], [true, false, true, true], [true, true, true, true]];
+
+    // Use adjoint Std.TableLookup.Select because this check takes adjoint of that.
+    let equal = CheckOperationsAreEqual(
+        n + width,
+        qs => RecursiveLookupOpt(true, data, qs[0..n-1], qs[n...]),
+        qs => Adjoint Std.TableLookup.Select(data, qs[0..n-1], qs[n...])
+    );
+    Fact(equal, "RecursiveLookupOpt should match Std.TableLookup.Select.");
+}
diff --git a/library/table_lookup/src/Tests.qs b/library/table_lookup/src/Tests.qs
new file mode 100644
index 0000000000..3b3742d99c
--- /dev/null
+++ b/library/table_lookup/src/Tests.qs
@@ -0,0 +1,84 @@
+// Copyright (c) Microsoft Corporation.
+// Licensed under the MIT License.
+
+import Std.Diagnostics.*;
+
+import Main.*;
+
+internal operation MatchLookupToStd(
+    options : LookupOptions
+) : Unit {
+    let n = 3;
+    let width = 4;
+    let data = [[true, false, false, false], [false, true, false, false], [false, false, true, false], [false, false, false, false], [true, true, false, false], [false, true, true, false], [true, false, true, true], [true, true, true, true]];
+
+    // Use adjoint Std.TableLookup.Select because this check takes adjoint of that.
+    let equal = CheckOperationsAreEqual(
+        n + width,
+        qs => Lookup(options, data, qs[0..n-1], qs[n...]),
+        qs => Adjoint Std.TableLookup.Select(data, qs[0..n-1], qs[n...])
+    );
+    Fact(equal, "Lookup should match Std.TableLookup.Select.");
+}
+
+internal operation MatchControlledLookupToMCX(
+    options : LookupOptions
+) : Unit {
+    let n = 2;
+    let width = 3;
+    let data = [[true, false, false], [false, true, false], [false, false, true], [true, true, true]];
+
+
+    // CheckOperationsAreEqual uses adjoint variant of the reference operation (seond operation).
+    // Select from the standard library uses assumptions that the target is in zero state,
+    // so its adjoint always returns target to zero state. So it won't work for CheckOperationsAreEqual directly.
+    // Instead, we compare controlled Select to controlled LookupViaMCX, which works in all cases.
+    let equal = CheckOperationsAreEqual(
+        1 + n + width,
+        qs => Controlled Lookup(
+            [qs[0]],
+            (options, data, qs[1..n], qs[n + 1...])
+        ),
+        qs => Controlled LookupViaMCX(
+            [qs[0]],
+            (data, qs[1..n], qs[n + 1...])
+        )
+    );
+    Fact(equal, "Controlled Lookup should match controlled LookupViaMCX.");
+}
+
+internal operation TestOnAllAlgorithms(op : LookupOptions => Unit) : Unit {
+    let algorithms = [
+        DoStdLookup(),
+        DoMCXLookup(),
+        DoRecursiveSelectLookup(),
+        DoPPLookup(),
+        DoSplitPPLookup()
+    ];
+    for algorithm in algorithms {
+        let options = new LookupOptions {
+            lookupAlgorithm = algorithm,
+            unlookupAlgorithm = DoUnlookupViaLookup(),
+            failOnLongData = false,
+            failOnShortData = false,
+            respectExcessiveAddress = false,
+            preferMeasurementBasedUncomputation = true,
+        };
+        op(options);
+    }
+}
+
+@Test()
+operation TestDefaultLookupMatchesStd() : Unit {
+    MatchLookupToStd(DefaultLookupOptions());
+}
+
+@Test()
+operation TestLookupMatchesStd() : Unit {
+    TestOnAllAlgorithms(MatchLookupToStd);
+}
+
+@Test()
+operation TestControlledLookupMatchesMCX() : Unit {
+    TestOnAllAlgorithms(MatchControlledLookupToMCX);
+}
diff --git a/library/table_lookup/src/Utils.qs b/library/table_lookup/src/Utils.qs
new file mode 100644
index 0000000000..2d4cb48f24
--- /dev/null
+++ b/library/table_lookup/src/Utils.qs
@@ -0,0 +1,255 @@
+// Copyright (c) Microsoft Corporation.
+// Licensed under the MIT License.
+
+import Std.Diagnostics.*;
+import Std.Math.*;
+import Std.Arrays.*;
+import Std.Convert.*;
+import Std.Logical.Xor;
+
+import Lookup.*;
+
+struct AddressAndData {
+    // Lower part of the address needed to index into data.
+    fitAddress : Qubit[],
+    // Data padded or trimmed to fit needed address space.
+    fitData : Bool[][],
+}
+
+function PrepareAddressAndData(
+    options : LookupOptions,
+    address : Qubit[],
+    data : Bool[][]
+) : AddressAndData {
+    let address_size = Length(address);
+    let address_space = 1 <<< address_size;
+    let data_length = Length(data);
+
+    if (data_length == address_space) {
+        // Data length match address space, nothing to adjust.
+        return new AddressAndData {
+            fitAddress = address,
+            fitData = data,
+        };
+    }
+
+    if (data_length > address_space) {
+        // Truncate longer data if needed.
+        if options.failOnLongData {
+            fail $"Data length {data_length} exceeds address space {address_space}.";
+        }
+        return new AddressAndData {
+            fitAddress = address,
+            fitData = data[...address_space-1]
+        };
+    }
+
+    // Data is shorter than addressable space. Truncate excessive address if needed.
+
+    if (not options.failOnShortData) {
+        fail $"Data length {data_length} is shorter than address space {address_space}.";
+    }
+
+    if (options.respectExcessiveAddress) {
+        return new AddressAndData {
+            fitAddress = address,
+            fitData = data,
+        };
+    }
+
+    if (data_length <= 1) {
+        // No address qubits are needed for data length 0.
+        // Case data_length == 1 is for compatibility with earlier behavior.
+        return new AddressAndData {
+            fitAddress = [],
+            fitData = data,
+        };
+    }
+
+    let address_size_needed = BitSizeI(data_length - 1);
+    Fact(address_size_needed <= address_size, "Internal error: address_size_needed should be at most address_size.");
+
+    let address_space_needed = 1 <<< address_size_needed;
+    Fact(address_space_needed >= data_length, "Internal error: address_space_needed should be at least data_length.");
+
+    return new AddressAndData {
+        // Trim address qubits to needed size.
+        fitAddress = address[...address_size_needed - 1],
+        // Shorter data in this case will be handled later.
+        fitData = data,
+    };
+}
+
+/// # Summary
+/// Computes the Fast Möbius Transform of a boolean array over GF(2).
+/// Also known as the Walsh-Hadamard Transform or subset sum transform.
+///
+/// # Description
+/// This transform converts minterm coefficients to monomial coefficients.
+/// For each position i in the result, it computes the XOR (sum over GF(2)) of all
+/// input elements at positions that are subsets of i (when i is interpreted as a bitmask).
+///
+/// This is equivalent to multiplying the input vector by a triangular matrix
+/// where entry (i,j) is 1 if j is a subset of i (as bitmasks), and 0 otherwise.
+///
+/// # Input
+/// ## qs
+/// Boolean array of minterm coefficients of length 2^n for some integer n ≥ 0.
+///
+/// # Output
+/// Boolean array of the same length as input containing monomial coefficients.
+///
+/// # Remarks
+/// This function is the classical preprocessing step for quantum phase lookup operations,
+/// converting phase data from standard basis coefficients to power product coefficients.
+/// The transformation is its own inverse when applied twice.
+function FastMobiusTransform(coefficients : Bool[]) : Bool[] {
+    let len = Length(coefficients);
+    Fact((len &&& (len-1)) == 0, "Length of a qubit register should be a power of 2");
+    let n = BitSizeI(len)-1;
+
+    mutable result = coefficients;
+    // For each bit position (from least to most significant).
+    for i in 0..n-1 {
+        let step = 2^i;
+        // For each pair of positions that differ only in that bit.
+        for j in 0..(step * 2)..len-1 {
+            for k in 0..step-1 {
+                // XOR the "upper" position with the "lower" position.
+                result[j + k + step] = Xor(result[j + k + step], result[j + k]);
+            }
+        }
+    }
+    return result;
+}
+
+/// # Summary
+/// Measures all qubits in the `target` register in the X basis. Resets them to |0⟩.
+/// Computes and pads resulting phase data to cover the entire address space `2^address_size`.
+operation MeasureAndComputePhaseData(
+    target : Qubit[],
+    data : Bool[][],
+    address_size : Int
+) : Bool[] {
+    // Measure target register in X basis.
+    mutable measurements = [];
+    for qubit in target {
+        set measurements += [MResetX(qubit) == One];
+    }
+
+    // Get phasing data via parity checks.
+    mutable phaseData = [];
+    for x in data {
+        set phaseData += [BinaryInnerProduct(x, measurements)];
+    }
+
+    // Pad phase data at the end to cover the entire address space.
+    Padded(-2^address_size, false, phaseData)
+}
+
+/// # Summary
+/// Computes dot (inner) product of two vectors over GF(2) field.
+/// This isn't a proper inner product as it is not positive-definite.
+///
+/// It is used to see if a phase correction is needed for a bit string `data`
+/// after obtaining a measurement result `measurements`.
+function BinaryInnerProduct(data : Bool[], measurements : Bool[]) : Bool {
+    mutable sum = false;
+    for i in IndexRange(measurements) {
+        set sum = Xor(sum, (data[i] and measurements[i]));
+    }
+    sum
+}
+
+/// # Summary
+/// Combines multiple control qubits into a single control qubit using auxiliary qubits.
+/// Logarithmic depth and linear number of auxiliary qubits are used.
+operation CombineControls(controls : Qubit[], aux : Qubit[]) : Unit is Adj {
+    Fact(Length(controls) >= 1, "CombineControls: controls must not be empty.");
+    Fact(Length(controls) == Length(aux) + 1, "CombineControls: control and aux length mismatch.");
+    let combined = controls + aux;
+    let aux_offset = Length(controls);
+    for i in 0..aux_offset-2 {
+        AND(combined[i * 2], combined[i * 2 + 1], combined[aux_offset + i]);
+    }
+}
+
+/// # Summary
+/// Retrieves the combined control qubit after CombineControls operation.
+function GetCombinedControl(controls : Qubit[], aux : Qubit[]) : Qubit {
+    Fact(Length(controls) >= 1, "GetCombinedControl: controls must not be empty.");
+    Fact(Length(controls) == Length(aux) + 1, "GetCombinedControl: control and aux length mismatch.");
+    if Length(controls) == 1 {
+        return Head(controls);
+    } else {
+        return Tail(aux);
+    }
+}
+
+// =============================
+// Tests
+
+@Test()
+function TestFastMobiusTransform() : Unit {
+    // Test cases for FastMobiusTransform.
+    let testCases = [
+        ([], []),
+        ([false], [false]),
+        ([true], [true]),
+        ([false, false], [false, false]),
+        ([false, true], [false, true]),
+        ([true, false], [true, true]),
+        ([true, true], [true, false]),
+        ([false, false, false, false], [false, false, false, false]),
+        ([false, false, false, true], [false, false, false, true]),
+        ([false, false, true, false], [false, false, true, true]),
+        ([false, false, true, true], [false, false, true, false]),
+        ([true, true, true, true], [true, false, false, false]),
+        ([true, false, false, false], [true, true, true, true]),
+    ];
+    for (input, expected) in testCases {
+        let output = FastMobiusTransform(input);
+        Fact(output == expected, $"FastMobiusTransform({input}) should be {expected}, got {output}");
+        // Test that applying the transform twice returns the original input.
+        let roundTrip = FastMobiusTransform(output);
+        Fact(roundTrip == input, $"FastMobiusTransform(FastMobiusTransform({input})) should be {input}, got {roundTrip}");
+    }
+}
+
+internal operation TestCombineControlsForN(n : Int) : Unit {
+    let all_ones = 2^n - 1;
+    use controls = Qubit[n];
+    use aux = Qubit[n - 1];
+
+    // Test all combinations of control qubits.
+    for i in 0..all_ones {
+        ApplyXorInPlace(i, controls);
+
+        // Combine controls.
+        within {
+            CombineControls(controls, aux);
+        } apply {
+            let combined = GetCombinedControl(controls, aux);
+            // Check that combined control is |1⟩ iff all controls are |1⟩.
+            if i == all_ones {
+                within {
+                    // Ensure combined control is |1⟩.
+                    X(combined);
+                } apply {
+                    Fact(CheckZero(combined), $"Combined control should be |1⟩ when all {n} controls are |1⟩.");
+                }
+            } else {
+                Fact(CheckZero(combined), $"Combined control should be |0⟩ when some of {n} controls are |0⟩.");
+            }
+        }
+        ApplyXorInPlace(i, controls);
+        // Check that all qubits are reset to |0⟩.
+        Fact(CheckAllZero(controls + aux), "All qubits should be reset to |0⟩ after CombineControls adjoint.");
+    }
+}
+
+@Test()
+operation TestCombineControls() : Unit {
+    TestCombineControlsForN(1);
+    TestCombineControlsForN(4);
+}