Add arithmetic library

build.py

@@ -440,7 +440,11 @@ def run_python_tests(cwd, interpreter, pip_env):
 def run_python_integration_tests(cwd, interpreter):
     # don't check to see if pip succeeds. We'll see if the import works later.
     # If it doesn't, we'll skip the tests.
-    command_args = [interpreter, "-m", "pytest"]
+    paths = [
+        cwd,
+        os.path.join(root_dir, "library", "arithmetic", "test"),
+    ]
+    command_args = [interpreter, "-m", "pytest", *paths]
     subprocess.run(command_args, check=True, text=True, cwd=cwd)
 
 

library/arithmetic/README.md

@@ -0,0 +1,38 @@
+# Arithmetic
+
+This library contains advanced quantum arithmetic algorithms.
+
+Unless otherwise noted, all quantum inputs are interpreted as unsigned little-endian
+integers and represented by `Qubit[]`.
+
+The library contains the following algorithms:
+* `Add.qs` - quantum-quantum in-place addition and subtraction (modulo `2^n`).
+* `AddConst.qs` - quantum-classical addition (modulo `2^n`).
+* `AddLookup.qs` - computes `x += table[i]` where `x`, `i` are quantum registers and
+    `table` is a classical table. Supports modular and non-modular addition.
+* `Compare.qs` - compares two unsigned quantum integers, writing the result to an output
+   qubit. Supports inequality and equality.
+* Modular arithmetic:
+  * `ModAdd.qs` - quantum-quantum in-place addition modulo classical constant.
+  * `ModDiv.qs` - division of two quantum numbers modulo classical constant.
+    This is based on the Extended Euclidean Algorithm, so some restrictions apply (in
+    particular, the divisor must be mutually prime with the modulus). It can also be
+    used for modular multiplication and modular inversion.
+  * `ModExp.qs` - modular exponentiation (computes `t:=(t*b^x)%m` where
+    `t`, `x` are quantum and `b`, `m` are classical).
+  * `ModMul.qs` - modular multiplication and square.
+  * `ModNegate.qs` - modular negation.
+
+### Space-optimized and time-optimized variants
+
+Some algorithms (addition, constant addition, comparison) are implemented with two
+different circuit variants: space-optimized (minimizing the number of qubits used) and
+time-optimized (minimizing the number of certain gates). The variants are functionally
+equivalent but have different resource usage.
+
+To select a variant, use the Q# configuration `"optimize"` with the value `"space"` or
+`"time"`. For example, when creating a QDK Context using the Python API:
+`qdk.Context(..., qdk_config={"optimize": "space"})`.
+
+The variant used by default is unspecified (but it is usually the space-optimized one).
+

library/arithmetic/qsharp.json

@@ -0,0 +1,19 @@
+{
+    "license": "MIT",
+    "author": "Microsoft Quantum",
+    "files": [
+        "src/Add.qs",
+        "src/AddConst.qs",
+        "src/AddLookup.qs",
+        "src/ClassicalMath.qs",
+        "src/Compare.qs",
+        "src/ModAdd.qs",
+        "src/ModDiv.qs",
+        "src/ModExp.qs",
+        "src/ModMul.qs",
+        "src/ModNegate.qs",
+        "src/MultiControl.qs",
+        "src/ResourceEstimation.qs",
+        "src/Utils.qs"
+    ]
+}

library/arithmetic/src/Add.qs

@@ -0,0 +1,40 @@
+// Copyright (c) Microsoft Corporation.
+// Licensed under the MIT License.
+
+import Std.Arithmetic.RippleCarryCGIncByLE;
+import Std.Arithmetic.RippleCarryTTKIncByLE;
+
+/// This file re-exports addition algorithms from Std.Arithmetic.
+
+/// Computes y += x (mod 2^n).
+operation Add(x : Qubit[], y : Qubit[]) : Unit is Ctl + Adj {
+    body (...) {
+        let optimize = Std.Core.ConfigValue("optimize", "");
+        if (optimize == "space") {
+            RippleCarryTTKIncByLE(x, y);
+        } elif (optimize == "time") {
+            RippleCarryCGIncByLE(x, y);
+        } else {
+            RippleCarryTTKIncByLE(x, y);
+        }
+    }
+    controlled (controls, ...) {
+        let optimize = Std.Core.ConfigValue("optimize", "");
+        if (Length(controls) == 0) {
+            Add(x, y);
+        } elif (optimize == "space") {
+            Controlled RippleCarryTTKIncByLE(controls, (x, y));
+        } elif (optimize == "time") {
+            Controlled RippleCarryCGIncByLE(controls, (x, y));
+        } else {
+            Controlled RippleCarryTTKIncByLE(controls, (x, y));
+        }
+    }
+}
+
+/// Computes y -= x (mod 2^n).
+operation Subtract(x : Qubit[], y : Qubit[]) : Unit is Ctl + Adj {
+    Adjoint Add(x, y);
+}
+
+export Add, Subtract;

library/arithmetic/src/AddConst.qs

@@ -0,0 +1,121 @@
+// Copyright (c) Microsoft Corporation.
+// Licensed under the MIT License.
+
+import Std.Arithmetic.RippleCarryCGAddLE;
+import Std.Arithmetic.RippleCarryCGIncByLE;
+import Std.Convert.BigIntAsBoolArray;
+import Std.Math.TrailingZeroCountL;
+
+import ClassicalMath.SafeMod;
+
+/// # Summary
+/// Computes `input := (input + constant) % 2^Length(input)`.
+///
+/// # Reference
+/// - [1](https://arxiv.org/pdf/2007.07391) "Compilation of Fault-Tolerant Quantum
+///   Heuristics for Combinatorial Optimization", Sanders et al. (Fig. 18).
+///
+/// # Resources
+/// Uses n-1 auxiliary qubits and 2n-2 Toffoli gates.
+///
+/// # Input
+/// ## constant
+/// Classical constant to add.
+/// ## inp
+/// Target register updated in place.
+operation AddConstantSanders(constant : BigInt, inp : Qubit[]) : Unit is Adj + Ctl {
+    body (...) {
+        Controlled AddConstantSanders([], (constant, inp));
+    }
+    controlled (ctrl, ...) {
+        let n = Length(inp);
+        let constant_bits = BigIntAsBoolArray(constant, n);
+
+        if n == 1 {
+            // Base case: single qubit addition.
+            if (constant_bits[0]) { Controlled X(ctrl, (inp[0])); }
+        } else {
+            use ancillas = Qubit[n - 1];
+
+            if (constant_bits[0]) { CNOT(inp[0], ancillas[0]); }
+
+            for i in 1..n - 2 {
+                let j = i - 1;
+                Controlled CNOT(ctrl, (ancillas[j], inp[i]));
+                within {
+                    if (constant_bits[i]) { X(ancillas[j]); }
+                } apply {
+                    AND(ancillas[j], inp[i], ancillas[i]);
+                }
+                CNOT(ancillas[j], ancillas[i]);
+            }
+
+            Controlled CNOT(ctrl, (ancillas[n - 2], inp[n - 1]));
+
+            for i in n - 2..-1..1 {
+                let j = i - 1;
+                CNOT(ancillas[j], ancillas[i]);
+                within {
+                    if (constant_bits[i]) { X(ancillas[j]); }
+                } apply {
+                    Adjoint AND(ancillas[j], inp[i], ancillas[i]);
+                }
+            }
+
+            if (constant_bits[0]) { CNOT(inp[0], ancillas[0]); }
+
+            for i in 0..n - 1 {
+                if (constant_bits[i]) { Controlled X(ctrl, (inp[i])); }
+            }
+        }
+    }
+}
+
+/// # Summary
+/// Constant adder using the Gidney ripple-carry adder.
+///
+/// # Reference
+/// - [1](https://arxiv.org/abs/1709.06648) "Halving the cost of quantum addition",
+///   Craig Gidney.
+operation AddConstantUsingCGAdd(constant : BigInt, input : Qubit[]) : Unit is Adj + Ctl {
+    body (...) {
+        Controlled AddConstantUsingCGAdd([], (constant, input));
+    }
+    controlled (ctrl, ...) {
+        use anc = Qubit[Length(input)];
+        within {
+            Controlled ApplyXorInPlaceL(ctrl, (constant, anc));
+        } apply {
+            RippleCarryCGIncByLE(anc, input);
+        }
+    }
+}
+
+/// # Summary
+/// Computes `input := (input + constant) % 2^Length(input)`.
+///
+/// # Input
+/// ## constant
+/// Classical constant to add.
+/// ## input
+/// Target register updated in place.
+operation AddConstant(constant : BigInt, input : Qubit[]) : Unit is Adj + Ctl {
+    body (...) {
+        Controlled AddConstant([], (constant, input));
+    }
+    controlled (ctrl, ...) {
+        let n = Length(input);
+        let constant = SafeMod(constant, 1L <<< n);
+        if (constant != 0L) {
+            let tz = TrailingZeroCountL(constant);
+            let optimize = Std.Core.ConfigValue("optimize", "");
+            if (optimize == "time" and Length(ctrl) > 0) {
+                Controlled AddConstantUsingCGAdd(ctrl, (constant >>> tz, input[tz...]));
+            } else {
+                Controlled AddConstantSanders(ctrl, (constant >>> tz, input[tz...]));
+            }
+        }
+    }
+}
+
+export AddConstantSanders, AddConstant;

library/arithmetic/src/AddLookup.qs

@@ -0,0 +1,160 @@
+// Copyright (c) Microsoft Corporation.
+// Licensed under the MIT License.
+
+import Std.Arithmetic.RippleCarryCGIncByLE;
+import Std.Arrays.Mapped;
+import Std.Convert.BigIntAsBoolArray;
+import Std.Diagnostics.Fact;
+import Std.ResourceEstimation.IsResourceEstimating;
+import Std.ResourceEstimation.RepeatEstimates;
+import Std.TableLookup.Select;
+
+import ClassicalMath.SafeMod;
+import ModAdd.ModAdd;
+
+/// References:
+/// - Thomas Haener, Vadym Kliuchnikov, Martin Roetteler, Mathias Soeken,
+///   "Space-time optimized table lookup", 2022.
+///   https://arxiv.org/abs/2211.01133
+/// - Dominic W. Berry, Craig Gidney, Mario Motta, Jarrod R. McClean, Ryan Babbush,
+///   "Qubitization of Arbitrary Basis Quantum Chemistry Leveraging Sparsity and
+///   Low Rank Factorization" (Appendix C), 2019.
+///   https://arxiv.org/abs/1902.02134
+
+function _CombineTables(tables : BigInt[][], num_bits : Int, modulus : BigInt) : BigInt[] {
+    mutable combined = Mapped(x -> SafeMod(x, modulus), tables[0]);
+    for i in 1..Length(tables) - 1 {
+        let shift = i * num_bits;
+        for j in 0..Length(combined) - 1 {
+            let shifted_value = SafeMod(tables[i][j], modulus) <<< shift;
+            set combined w/= j <- combined[j] ||| shifted_value;
+        }
+    }
+    return combined;
+}
+
+function DataRowsAsBits(data : BigInt[], num_bits : Int) : Bool[][] {
+    let wrap = 1L <<< num_bits;
+    return Mapped(x -> BigIntAsBoolArray(SafeMod(x, wrap), num_bits), data);
+}
+
+// Requires 0 <= data[i] < 2^Length(address).
+operation _Lookup(data : BigInt[], address : Qubit[], target : Qubit[]) : Unit is Adj {
+    let address_size = Length(address);
+    let can_use_formula = (address_size >= 3) and (Length(data) == 2^address_size);
+
+    if (can_use_formula and IsResourceEstimating()) {
+        let num_ancilla = address_size - 1;
+        use q_anc = Qubit[num_ancilla];
+        within {
+            RepeatEstimates(2^address_size - 2);
+        } apply {
+            AND(address[0], address[1], q_anc[0]);
+            Adjoint AND(address[0], address[1], q_anc[0]);
+        }
+    } else {
+        let data_bits = DataRowsAsBits(data, Length(target));
+        Select(data_bits, address, target);
+    }
+}
+
+/// # Summary
+/// Computes `q_result := (q_result + data[q_address]) % 2^n`.
+/// Requires `0 <= data[i] < 2^n`.
+///
+/// # Input
+/// ## q_address
+/// Register encoding the table index `q_address`.
+/// ## q_result
+/// Target register updated in place.
+/// ## data
+/// Lookup table of values to add.
+operation AddLookup(
+    q_address : Qubit[],
+    q_result : Qubit[],
+    data : BigInt[]
+) : Unit {
+    use q_select_output = Qubit[Length(q_result)];
+
+    within {
+        _Lookup(data, q_address, q_select_output);
+    } apply {
+        RippleCarryCGIncByLE(q_select_output, q_result);
+    }
+}
+
+/// # Summary
+/// Computes `q_result := (q_result + data[q_address]) % modulus`.
+///
+/// # Input
+/// ## q_address
+/// Register encoding the table index `q_address`.
+/// ## q_result
+/// Target register updated in place.
+/// ## data
+/// Lookup table of values to add.
+/// ## modulus
+/// Modulus used for modular addition.
+operation ModAddLookup(
+    q_address : Qubit[],
+    q_result : Qubit[],
+    data : BigInt[],
+    modulus : BigInt
+) : Unit is Adj {
+    use q_select_output = Qubit[Length(q_result)];
+    let data_modded = Mapped(x -> SafeMod(x, modulus), data);
+    within {
+        _Lookup(data_modded, q_address, q_select_output);
+    } apply {
+        ModAdd(q_select_output, q_result, modulus);
+    }
+}
+
+/// # Summary
+/// Computes `q_result[i] := (q_result[i] + tables[i][q_address]) % modulus` for each `i`.
+///
+/// # Input
+/// ## q_address
+/// Register encoding the table index `q_address`.
+/// ## q_result
+/// Array of target registers, each updated in place.
+/// ## tables
+/// Array of lookup tables matched one-to-one with `q_result`.
+/// ## modulus
+/// Modulus used for each modular addition.
+operation ParallelModAddLookup(
+    q_address : Qubit[],
+    q_result : Qubit[][],
+    tables : BigInt[][],
+    modulus : BigInt
+) : Unit {
+    let num_tables = Length(tables);
+    Fact(num_tables == Length(q_result), "Size mismatch.");
+    Fact(num_tables > 0, "Must provide at least one table.");
+    let result_size = Length(q_result[0]);
+    let table_length = Length(tables[0]);
+    for i in 0..num_tables - 1 {
+        Fact(Length(q_result[i]) == result_size, "All target registers must have equal size.");
+        Fact(Length(tables[i]) == table_length, "All tables must have the same length.");
+    }
+
+    let optimize = Std.Core.ConfigValue("optimize", "");
+    if (optimize == "space") {
+        for i in 0..num_tables - 1 {
+            ModAddLookup(q_address, q_result[i], tables[i], modulus);
+        }
+    } else {
+        use q_select_output = Qubit[num_tables * result_size];
+        let tables_combined = _CombineTables(tables, result_size, modulus);
+        within {
+            _Lookup(tables_combined, q_address, q_select_output);
+        } apply {
+            for i in 0..num_tables - 1 {
+                let q_temp = q_select_output[i * result_size..(i + 1) * result_size - 1];
+                ModAdd(q_temp, q_result[i], modulus);
+            }
+        }
+    }
+}
+
+export AddLookup, ModAddLookup, ParallelModAddLookup;