diff --git a/build.py b/build.py
index b05cf59fb50..4768b6000f2 100755
--- a/build.py
+++ b/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)
 
 
diff --git a/library/arithmetic/README.md b/library/arithmetic/README.md
new file mode 100644
index 00000000000..060bc9d1aad
--- /dev/null
+++ b/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).
+
diff --git a/library/arithmetic/qsharp.json b/library/arithmetic/qsharp.json
new file mode 100644
index 00000000000..44387d53555
--- /dev/null
+++ b/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"
+    ]
+}
diff --git a/library/arithmetic/src/Add.qs b/library/arithmetic/src/Add.qs
new file mode 100644
index 00000000000..dfd078727ce
--- /dev/null
+++ b/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;
diff --git a/library/arithmetic/src/AddConst.qs b/library/arithmetic/src/AddConst.qs
new file mode 100644
index 00000000000..8ed01bfe4bd
--- /dev/null
+++ b/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;
diff --git a/library/arithmetic/src/AddLookup.qs b/library/arithmetic/src/AddLookup.qs
new file mode 100644
index 00000000000..add4172b6f6
--- /dev/null
+++ b/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;
diff --git a/library/arithmetic/src/ClassicalMath.qs b/library/arithmetic/src/ClassicalMath.qs
new file mode 100644
index 00000000000..42cbb8a295d
--- /dev/null
+++ b/library/arithmetic/src/ClassicalMath.qs
@@ -0,0 +1,34 @@
+// Copyright (c) Microsoft Corporation.
+// Licensed under the MIT License.
+
+/// Classical mathematical functions used in arithmetic algorithms.
+
+import Std.Math.Log;
+import Std.Math.LogOf2;
+
+/// Computes x%mod with result in [0, mod).
+function SafeMod(x : BigInt, mod : BigInt) : BigInt {
+    let remainder = x % mod;
+    return remainder < 0L ? remainder + mod | remainder;
+}
+
+/// Computes ⌈a/b⌉.
+function DivCeil(a : Int, b : Int) : Int {
+    return ((a + b - 1) / b);
+}
+
+/// Computes base-2 logarithm of x.
+function Log2(x : Double) : Double {
+    return Log(x) / LogOf2();
+}
+
+/// Computes [(a^(2^i))%N for i in 0..n-1].
+function ComputeSequentialSquares(a : BigInt, N : BigInt, n : Int) : BigInt[] {
+    mutable ans : BigInt[] = [((a % N) + N) % N];
+    for i in 1..n-1 {
+        set ans += [(ans[i-1] * ans[i-1]) % N];
+    }
+    return ans;
+}
+
+export SafeMod, DivCeil, Log2, ComputeSequentialSquares;
diff --git a/library/arithmetic/src/Compare.qs b/library/arithmetic/src/Compare.qs
new file mode 100644
index 00000000000..93dfacfd517
--- /dev/null
+++ b/library/arithmetic/src/Compare.qs
@@ -0,0 +1,113 @@
+// Copyright (c) Microsoft Corporation.
+// Licensed under the MIT License.
+
+import Std.Arithmetic.MAJ;
+import Std.Arithmetic.ApplyIfGreaterLE;
+import Std.Diagnostics.Fact;
+
+import Utils.ParallelCNOT;
+import Utils.CheckIfAllZero;
+
+/// # Summary
+/// Flips `q_answer` iff `q_input > q_out` using a Cuccaro-style ripple comparator.
+/// The registers `q_input` and `q_out` are restored to their original values.
+/// In controlled form, this operation supports at most one control qubit.
+///
+/// # Reference
+/// - [1](https://arxiv.org/abs/1202.6614) "Constant-Optimized Quantum Circuits for
+///   Modular Multiplication and Exponentiation", Igor L. Markov, Mehdi Saeedi.
+///   (Fig. 6).
+///
+/// # Input
+/// ## q_input
+/// First input register.
+/// ## q_out
+/// Second input register with the same length as `q_input`.
+/// ## q_answer
+/// Target qubit flipped when `q_input > q_out`.
+operation CompareCuccaro(
+    q_input : Qubit[],
+    q_out : Qubit[],
+    q_answer : Qubit
+) : Unit is Adj + Ctl {
+    body (...) {
+        Controlled CompareCuccaro([], (q_input, q_out, q_answer));
+    }
+    controlled (controls, ...) {
+        let n = Length(q_input);
+        Fact(n >= 1, "Input register must contain at least one qubit.");
+        Fact(Length(q_out) == n, "Input size mismatch.");
+        Fact(Length(controls) <= 1, "Operation takes at most 1 control.");
+
+        use q_anc = Qubit();
+
+        // Negate q_out so q_input + q_out becomes q_input - q_out.
+        ApplyToEach(X, q_out);
+
+        MAJ(q_anc, q_input[0], q_out[0]);
+        for i in 1..n - 1 {
+            MAJ(q_out[i - 1], q_input[i], q_out[i]);
+        }
+
+        Controlled CNOT(controls, (q_out[n - 1], q_answer));
+
+        for i in n - 1..-1..1 {
+            Adjoint MAJ(q_out[i - 1], q_input[i], q_out[i]);
+        }
+        Adjoint MAJ(q_anc, q_input[0], q_out[0]);
+
+        ApplyToEach(X, q_out);
+    }
+    adjoint self;
+}
+
+operation CompareRippleCarry(x : Qubit[], y : Qubit[], result : Qubit) : Unit is Adj + Ctl {
+    body (...) {
+        ApplyIfGreaterLE(X, x, y, result);
+    }
+    controlled (controls, ...) {
+        ApplyIfGreaterLE(Controlled X(controls, _), x, y, result);
+    }
+    adjoint self;
+}
+
+/// Flips `result` iff x > y.
+operation CompareGT(x : Qubit[], y : Qubit[], result : Qubit) : Unit is Adj + Ctl {
+    let optimize = Std.Core.ConfigValue("optimize", "");
+    if (optimize == "space") {
+        CompareCuccaro(x, y, result);
+    } elif (optimize == "time") {
+        CompareRippleCarry(x, y, result);
+    } else {
+        CompareCuccaro(x, y, result);
+    }
+}
+
+/// Flips `result` iff x < y.
+operation CompareLT(x : Qubit[], y : Qubit[], result : Qubit) : Unit is Adj + Ctl {
+    CompareGT(y, x, result);
+}
+
+/// Flips `result` iff x <= y.
+operation CompareLE(x : Qubit[], y : Qubit[], result : Qubit) : Unit is Adj + Ctl {
+    CompareGT(x, y, result);
+    X(result);
+}
+
+/// Flips `result` iff x >= y.
+operation CompareGE(x : Qubit[], y : Qubit[], result : Qubit) : Unit is Adj + Ctl {
+    CompareGT(y, x, result);
+    X(result);
+}
+
+/// Flips `result` iff x == y.
+operation CompareEQ(x : Qubit[], y : Qubit[], result : Qubit) : Unit is Adj + Ctl {
+    within {
+        ParallelCNOT(x, y);
+    } apply {
+        CheckIfAllZero(y, result);
+    }
+}
+
+
+export CompareCuccaro, CompareRippleCarry, CompareGT, CompareLT, CompareLE, CompareGE, CompareEQ;
diff --git a/library/arithmetic/src/ModAdd.qs b/library/arithmetic/src/ModAdd.qs
new file mode 100644
index 00000000000..b8584e34333
--- /dev/null
+++ b/library/arithmetic/src/ModAdd.qs
@@ -0,0 +1,104 @@
+// Copyright (c) Microsoft Corporation.
+// Licensed under the MIT License.
+
+import Std.Diagnostics.Fact;
+
+import Add.Add;
+import Add.Subtract;
+import AddConst.AddConstant;
+import Compare.CompareGT;
+
+/// Modular addition: y = (x+y) % modulus.
+/// Uses measurement-based uncomputation, which makes it non-adjointable.
+operation _ModAddMBU(x : Qubit[], y : Qubit[], modulus : BigInt) : Unit {
+    body (...) {
+        Controlled _ModAddMBU([], (x, y, modulus));
+    }
+    controlled (controls, ...) {
+        let n = Length(x);
+        Fact(Length(y) == n, "Input size mismatch.");
+        Fact(modulus >= 2L, "Modulus must be at least 2.");
+        Fact(modulus < (1L <<< n), "Modulus is too large for number of bits.");
+
+        use carry = Qubit();
+
+        // Step 1. Add.
+        Controlled Add(controls, (x, y + [carry]));
+
+        // Step 2. "-p".
+        AddConstant((1L <<< n) - modulus, y + [carry]);
+        X(carry);
+
+        // Step 3. Controlled "+p".
+        Controlled AddConstant([carry], (modulus, y));
+
+        // Step 4. Uncompute carry with comparator (measurement-based).
+        H(carry);
+        let qflag = M(carry);
+        if (qflag == One) {
+            H(carry);
+            Controlled CompareGT(controls, (x, y, carry));
+            H(carry);
+            X(carry);
+        }
+    }
+}
+
+/// Computes (x, y) := x, ((x+y)%p).
+/// See https://arxiv.org/pdf/1706.06752, Fig. 3.
+operation _ModAddNoMBU(x : Qubit[], y : Qubit[], modulus : BigInt) : Unit is Ctl + Adj {
+    body (...) {
+        Controlled _ModAddNoMBU([], (x, y, modulus));
+    }
+    controlled (controls, ...) {
+        let n = Length(x);
+        Fact(Length(y) == n, "Registers must have the same length");
+        use carry = Qubit();
+
+        // Step 1. Add.
+        Controlled Add(controls, (x, y + [carry]));
+
+        // Step 2. "-p".
+        AddConstant((1L <<< n) - modulus, y + [carry]);
+        X(carry);
+
+        // Step 3. Controlled "+p".
+        Controlled AddConstant([carry], (modulus, y));
+
+        // Step 4. Uncompute carry with comparator.
+        Controlled CompareGT(controls, (x, y, carry));
+        X(carry);
+    }
+}
+
+/// # Summary
+/// In-place modular addition of two equal-length little-endian registers.
+///
+/// Given registers `x` and `y` encoding integers with `0 <= x, y < modulus`,
+/// this operation computes `(x, y) -> (x, (x + y) mod modulus)`.
+///
+/// # References
+/// - [1](https://arxiv.org/abs/1706.06752) "Quantum resource estimates for computing
+///   elliptic curve discrete logarithms", Martin Roetteler, Michael Naehrig, Krysta M.
+///   Svore, Kristin Lauter. (Fig. 3).
+/// - [2](https://arxiv.org/abs/2407.20167) "Measurement-based uncomputation of quantum
+///   circuits for modular arithmetic", Alessandro Luongo, Antonio Michele Miti,
+///   Varun Narasimhachar, Adithya Sireesh.
+///
+/// # Input
+/// ## x
+/// Addend register. This register is preserved.
+/// ## y
+/// Accumulator register. This register is updated in place.
+/// ## modulus
+/// Classical modulus. Must satisfy `2 <= modulus < 2^Length(x)`.
+operation ModAdd(x : Qubit[], y : Qubit[], modulus : BigInt) : Unit is Adj + Ctl {
+    body (...) {
+        _ModAddMBU(x, y, modulus);
+    }
+    adjoint (...) {
+        Adjoint _ModAddNoMBU(x, y, modulus);
+    }
+}
+
+export ModAdd;
diff --git a/library/arithmetic/src/ModDiv.qs b/library/arithmetic/src/ModDiv.qs
new file mode 100644
index 00000000000..8547444febe
--- /dev/null
+++ b/library/arithmetic/src/ModDiv.qs
@@ -0,0 +1,322 @@
+// Copyright (c) Microsoft Corporation.
+// Licensed under the MIT License.
+
+import Std.Diagnostics.Fact;
+import Std.Convert.IntAsDouble;
+import Std.Arrays.Reversed;
+import Std.Math;
+
+import Add.Subtract;
+import ClassicalMath.DivCeil;
+import ClassicalMath.Log2;
+import Compare.CompareGT;
+import ModAdd.ModAdd;
+import ModMul.ModDouble;
+import ResourceEstimation.LoopA;
+import Utils.CheckIfAllZero;
+import Utils.ParallelSWAP;
+import Utils.RotateRight;
+
+/// This file contains an implementation of modular multiplication, division, and inversion
+/// based on the Extended Euclidean Algorithm, as described in the paper:
+/// Andre Schrottenloher, "Optimized Point Addition Circuits for Elliptic Curve Discrete
+/// Logarithms", 2026, https://arxiv.org/abs/2606.02235.
+
+/// Compresses six garbage qubits into five qubits.
+///
+/// Under the assumption that pairs `(q[2*i], q[2*i+1])` are never `(Zero, One)`
+/// for `i in 0..2`, this operation maps states so that `q[5]` is always `Zero`.
+/// This matches Fig. 1.
+operation CompressGarbage(qs : Qubit[]) : Unit is Adj + Ctl {
+    Fact(Length(qs) == 6, "CompressGarbage requires exactly 6 qubits.");
+
+    CNOT(qs[1], qs[0]);
+    CNOT(qs[3], qs[2]);
+    CNOT(qs[5], qs[4]);
+    CNOT(qs[0], qs[2]);
+    CNOT(qs[5], qs[3]);
+    X(qs[4]);
+    CCNOT(qs[1], qs[3], qs[5]);
+    CNOT(qs[1], qs[4]);
+    X(qs[2]);
+    CCNOT(qs[3], qs[4], qs[5]);
+    CCNOT(qs[5], qs[4], qs[1]);
+    CCNOT(qs[5], qs[2], qs[0]);
+    CCNOT(qs[0], qs[1], qs[5]);
+}
+
+// Size to which we can truncate registers on i-th iteration of GCD.
+function ActiveRegisterSize(n : Int, iterationIdx : Int) : Int {
+    let cPad = 2.3;
+    let regSizeApprox = IntAsDouble(n)
+        - 0.5 * Log2(8.0 / 3.0) * IntAsDouble(iterationIdx + 1) + cPad * Math.Sqrt(IntAsDouble(n));
+    return Math.MinI(n, Math.Ceiling(regSizeApprox));
+}
+
+/// Number of GCD iterations from Schrottenloher (2026), Section 3.1.
+///
+/// The paper shows that the required number of iterations follows a normal
+/// distribution with mean 1.413 and standard deviation 0.6*sqrt(n).
+/// Choosing c_iter=2.4 corresponds to 4 standard deviations, so the probability of
+/// error (not doing enough iterations) is 3e-5.
+function NumGcdIterations(n : Int) : Int {
+    let cIter = 2.4;
+    let estimate = 1.413 * IntAsDouble(n) + cIter * Math.Sqrt(IntAsDouble(n));
+    return Math.MinI(2 * n, Math.Ceiling(estimate));
+}
+
+/// Size of the garbage vector with compression.
+function GarbageVectorSize(n : Int) : Int {
+    let numIterations = DivCeil(NumGcdIterations(n), 3);
+    return 5 * numIterations + 1;
+}
+
+/// Auxiliary registers used by algorithms.
+/// qubits = g ∪ t.
+/// g and t may overlap if using qubit sharing.
+struct _AuxQubits {
+    // All auxiliary qubits.
+    qubits : Qubit[],
+    // Garbage register.
+    g : Qubit[],
+    // Register for the first argument to _Gcd.
+    t : Qubit[],
+}
+
+
+// Allocates auxiliary registers with qubit sharing.
+// Qubit sharing: share qubits between t and garbage vector in such a way that on every
+// iteration of GCD used qubits don't overlap.
+// We can do this because as we need to use more qubits of g, we need to use less
+// qubits of t (as number in t becomes smaller).
+// We share a suffix of g and a suffix of t.
+operation _AllocateAuxQubits(n : Int) : _AuxQubits {
+    let garbageSize = GarbageVectorSize(n);
+    let numIterations = DivCeil(NumGcdIterations(n), 3);
+    let minActiveSize = ActiveRegisterSize(n, 3 * numIterations - 1);
+    let numSharedQubits = n - minActiveSize;
+
+    let qubits = QIR.Runtime.AllocateQubitArray(garbageSize + minActiveSize);
+    let g = qubits[0..garbageSize-1];
+    let tExtra = qubits[garbageSize..garbageSize + minActiveSize - 1];
+    let tShared = Reversed(g)[0..numSharedQubits - 1];
+    let t = tExtra + tShared;
+    return _AuxQubits(qubits, g, t);
+}
+
+operation _ReleaseAuxQubits(aux : _AuxQubits) : Unit {
+    QIR.Runtime.ReleaseQubitArray(aux.qubits);
+}
+
+/// Single iteration of GCD algorithm from Schrottenloher (2026), Alg. 2.
+///
+/// This operation updates only the active prefixes of `u` and `v`.
+operation _GcdIteration(iteration_idx : Int, u : Qubit[], v : Qubit[], b0 : Qubit, b01 : Qubit) : Unit is Adj {
+    let n = Length(u);
+    let regSize = ActiveRegisterSize(n, iteration_idx);
+    let uReg = u[0..regSize - 1];
+    let vReg = v[0..regSize - 1];
+
+    // b0 := v % 2.
+    CNOT(vReg[0], b0);
+    // b01 := b0 & (u > v).
+    (Controlled CompareGT)([b0], (uReg, vReg, b01));
+    // if (b0 & b1) swap(u, v).
+    (Controlled ParallelSWAP)([b01], (uReg, vReg));
+    // if (b0) v -= u.
+    (Controlled Subtract)([b0], (uReg, vReg));
+    // v := v / 2.
+    RotateRight(vReg);
+}
+
+/// GCD algorithm from Schrottenloher (2026), Alg. 2, with garbage compression enabled.
+///
+/// Preconditions: `u` is odd, `g` is initialized to `0`, and `gcd(u, v) = 1`.
+/// Postcondition: `u = 1`, `v = 0`, and `g` stores compressed garbage bits.
+operation _Gcd(u : Qubit[], v : Qubit[], g : Qubit[]) : Unit is Adj {
+    let n = Length(u);
+    Fact(n > 0, "Registers must be non-empty.");
+    Fact(Length(v) == n, "Registers u and v must have the same length.");
+
+    let numIterations = DivCeil(NumGcdIterations(n), 3);
+    Fact(
+        Length(g) == GarbageVectorSize(n),
+        "Garbage register size mismatch for compressed GCD."
+    );
+
+    LoopA(numIterations, idx => {
+        let page = g[5 * idx..5 * idx + 5];
+        for i in 0..2 {
+            _GcdIteration(3 * idx + i, u, v, page[2 * i], page[2 * i + 1]);
+        }
+        CompressGarbage(page);
+    });
+}
+
+/// Single reconstruction iteration, corresponds to a single GCD iteration.
+operation _ReconstructionIteration(
+    r : Qubit[],
+    s : Qubit[],
+    b0 : Qubit,
+    b01 : Qubit,
+    modulus : BigInt
+) : Unit is Adj {
+    // s := (2*s) % modulus.
+    ModDouble(s, modulus);
+    // if (b0) s += r (mod modulus).
+    Controlled ModAdd([b0], (r, s, modulus));
+    // if (b0 & b1) swap(r, s).
+    (Controlled ParallelSWAP)([b01], (r, s));
+}
+
+/// Restores Bezout coefficients from compressed garbage vector produced by GCD.
+///
+/// Algorithm 3 from Schrottenloher (2026), using compressed pages.
+operation _BezoutReconstruction(modulus : BigInt, r : Qubit[], s : Qubit[], g : Qubit[]) : Unit is Adj {
+    let n = Length(r);
+    Fact(Length(s) == n, "Registers r and s must have the same length.");
+
+    let numIterations = DivCeil(NumGcdIterations(n), 3);
+    Fact(
+        Length(g) == GarbageVectorSize(n),
+        "Garbage register size mismatch for compressed reconstruction."
+    );
+
+    LoopA(numIterations, idx => {
+        let reverseIdx = numIterations - 1 - idx;
+        let page = g[5 * reverseIdx..5 * reverseIdx + 4] + [g[Length(g) - 1]];
+
+        Adjoint CompressGarbage(page);
+        _ReconstructionIteration(r, s, page[4], page[5], modulus);
+        _ReconstructionIteration(r, s, page[2], page[3], modulus);
+        _ReconstructionIteration(r, s, page[0], page[1], modulus);
+        CompressGarbage(page);
+    });
+}
+
+/// Computes `(x, y) -> (x, (x * y) % modulus)`.
+/// Takes all ancillas as inputs.
+/// Returns ancillas (g, t) in the zero state if x != 0.
+/// If x == 0, they will not be returned in the zero state and must be reset.
+operation _ModMul(x : Qubit[], y : Qubit[], modulus : BigInt, aux : _AuxQubits) : Unit is Adj {
+    let n = Length(x);
+    Fact(Length(y) == n, "Registers x and y must have the same length.");
+    Fact(3L <= modulus and modulus < (1L <<< n), "Modulus must satisfy 3 <= modulus < 2^n.");
+    Fact(modulus % 2L == 1L, "Modulus must be odd.");
+    Fact(Length(aux.g) == GarbageVectorSize(n), "Garbage size mismatch");
+
+    // Prepare t = modulus, run GCD on (t, x).
+    within {
+        ApplyXorInPlaceL(modulus, aux.t);
+        _Gcd(aux.t, x, aux.g);
+        ApplyXorInPlaceL(1L, aux.t);
+    } apply {
+        _BezoutReconstruction(modulus, y, x, aux.g);
+        ParallelSWAP(x, y);
+    }
+}
+
+
+/// # Summary
+/// Computes `(x, y) -> (x, (x * y) % modulus)`.
+/// Requires `0 < x < modulus` and `gcd(x, modulus) = 1`.
+/// If `x == 0`, this operation calls `Reset` on non-zero qubits.
+///
+/// # Input
+/// ## x
+/// Register storing `x`.
+/// ## y
+/// Register storing `y` and updated in place.
+/// ## modulus
+/// Odd modulus satisfying `3 <= modulus < 2^Length(x)`.
+operation ModMul(x : Qubit[], y : Qubit[], modulus : BigInt) : Unit {
+    let aux = _AllocateAuxQubits(Length(x));
+    _ModMul(x, y, modulus, aux);
+    ResetAll(aux.qubits);
+    _ReleaseAuxQubits(aux);
+}
+
+/// # Summary
+/// Computes `(x, y) -> (x, (y * x^-1) % modulus)`.
+/// Requires `0 < x < modulus` and `gcd(x, modulus) = 1`.
+/// If `x == 0`, this operation calls `Reset` on non-zero qubits.
+///
+/// # Input
+/// ## x
+/// Register storing `x`.
+/// ## y
+/// Register storing `y` and updated in place.
+/// ## modulus
+/// Odd modulus satisfying `3 <= modulus < 2^Length(x)`.
+operation ModDiv(x : Qubit[], y : Qubit[], modulus : BigInt) : Unit {
+    let aux = _AllocateAuxQubits(Length(x));
+    Adjoint _ModMul(x, y, modulus, aux);
+    ResetAll(aux.qubits);
+    _ReleaseAuxQubits(aux);
+}
+
+/// # Summary
+/// Computes `(x, y) -> (x, (x * y) % modulus)` if `x != 0`.
+/// Computes `(0, y) -> (0, y)` if `x == 0`.
+/// Requires `0 <= x < modulus` and `gcd(x, modulus) = 1`.
+///
+/// # Input
+/// ## x
+/// Register storing `x`.
+/// ## y
+/// Register storing `y` and updated in place.
+/// ## modulus
+/// Odd modulus satisfying `3 <= modulus < 2^Length(x)`.
+operation SafeModMul(x : Qubit[], y : Qubit[], modulus : BigInt) : Unit {
+    use isXZero = Qubit();
+    within {
+        CheckIfAllZero(x, isXZero);
+        CNOT(isXZero, x[0]);
+    } apply {
+        let aux = _AllocateAuxQubits(Length(x));
+        _ModMul(x, y, modulus, aux);
+        _ReleaseAuxQubits(aux);
+    }
+}
+
+/// # Summary
+/// Computes `(x, y) -> (x, (y * x^-1) % modulus)` if `x != 0`.
+/// Computes `(0, y) -> (0, y)` if `x == 0`.
+/// Requires `0 <= x < modulus` and `gcd(x, modulus) = 1`.
+///
+/// # Input
+/// ## x
+/// Register storing `x`.
+/// ## y
+/// Register storing `y` and updated in place.
+/// ## modulus
+/// Odd modulus satisfying `3 <= modulus < 2^Length(x)`.
+operation SafeModDiv(x : Qubit[], y : Qubit[], modulus : BigInt) : Unit {
+    use isXZero = Qubit();
+    within {
+        CheckIfAllZero(x, isXZero);
+        CNOT(isXZero, x[0]);
+    } apply {
+        let aux = _AllocateAuxQubits(Length(x));
+        Adjoint _ModMul(x, y, modulus, aux);
+        _ReleaseAuxQubits(aux);
+    }
+}
+
+/// # Summary
+/// Computes `(x, 0) -> (x, x^-1 (mod modulus))`.
+/// Requires `x != 0` and `gcd(x, modulus) = 1`.
+///
+/// # Input
+/// ## x
+/// Register storing `x`.
+/// ## ans
+/// Zero-initialized register that is updated to `x^-1 (mod modulus)`.
+/// ## modulus
+/// Odd modulus satisfying `3 <= modulus < 2^Length(x)`.
+operation ModInv(x : Qubit[], ans : Qubit[], modulus : BigInt) : Unit {
+    X(ans[0]); // ans:=1.
+    ModDiv(x, ans, modulus);
+}
+
+export ModMul, ModDiv, SafeModMul, SafeModDiv, ModInv;
diff --git a/library/arithmetic/src/ModExp.qs b/library/arithmetic/src/ModExp.qs
new file mode 100644
index 00000000000..000c4b2871e
--- /dev/null
+++ b/library/arithmetic/src/ModExp.qs
@@ -0,0 +1,167 @@
+// Copyright (c) Microsoft Corporation.
+// Licensed under the MIT License.
+
+// Q# operations for windowed modular exponentiation.
+
+import Std.Arrays.Mapped;
+import Std.Convert.IntAsBigInt;
+import Std.Diagnostics.Fact;
+import Std.Math;
+import Std.Math.*;
+
+import AddLookup.ModAddLookup;
+import ClassicalMath.*;
+import ResourceEstimation.LoopA;
+import Utils.ParallelSWAP;
+
+/// Computes powers of base: base_pows[i] = base^i % modulus.
+function ComputeBasePowers(base : BigInt, modulus : BigInt, num_powers : Int) : BigInt[] {
+    mutable base_pows = Repeated(0L, num_powers);
+    base_pows[0] = 1L;
+    let base_mod = base % modulus;
+    for a in 1..num_powers - 1 {
+        base_pows[a] = (base_pows[a - 1] * base_mod) % modulus;
+    }
+    return base_pows;
+}
+
+/// Generates a modular-exponentiation lookup table.
+/// For each a in [0,num_a), b in [0, num_b):
+/// data[(b * num_a) + a] = (factor * b * base^a * sign) % modulus.
+function ModExpLookupTable(
+    factor : BigInt,
+    exp_length : Int,
+    mul_length : Int,
+    base : BigInt,
+    modulus : BigInt,
+    sign : BigInt
+) : BigInt[] {
+    let num_a = 1 <<< exp_length;
+    let num_b = 1 <<< mul_length;
+    let base_pows = ComputeBasePowers(base, modulus, num_a);
+
+    mutable data = Repeated(0L, num_a * num_b);
+    for b in 0..num_b - 1 {
+        let fb : BigInt = SafeMod(sign * factor * IntAsBigInt(b), modulus);
+        for a in 0..num_a - 1 {
+            data[(b * num_a) + a] = (fb * base_pows[a]) % modulus;
+        }
+    }
+    return data;
+}
+
+/// Internal helper for windowed modular multiplication.
+operation WindowModularMultiply(
+    q_exponent_window : Qubit[],
+    q_source : Qubit[],
+    q_target : Qubit[],
+    multiply_window_size : Int,
+    adjusted_base : BigInt,
+    modulus : BigInt,
+    sign : BigInt
+) : Unit is Adj {
+    let num_multiply_windows = (Length(q_source) + multiply_window_size - 1) / multiply_window_size;
+
+    LoopA(num_multiply_windows, j => {
+        let window_start = j * multiply_window_size;
+        let window_end = Math.MinI(window_start + multiply_window_size, Length(q_source));
+        let q_multiply_window = q_source[window_start..window_end - 1];
+
+        let data = ModExpLookupTable(
+            1L <<< window_start,
+            Length(q_exponent_window),
+            Length(q_multiply_window),
+            adjusted_base,
+            modulus,
+            sign
+        );
+
+        let q_address = q_exponent_window + q_multiply_window;
+        ModAddLookup(q_address, q_target, data, modulus);
+    });
+}
+
+/// # Summary
+/// Windowed modular exponentiation over a little-endian exponent register.
+///
+/// Computes `q_target := (q_target * base^exponent)%modulus`
+/// using alternating forward and inverse windowed modular-multiply updates.
+///
+/// Reference:
+/// - [1](https://arxiv.org/abs/1905.07682) "Windowed quantum arithmetic", Craig Gidney.
+///   Section 3.5.
+///
+/// # Input
+/// ## q_target
+/// Target register that is multiplied in place by `base^exponent (mod modulus)`.
+/// ## q_exponent
+/// Little-endian exponent register controlling the modular exponentiation.
+/// ## base
+/// Classical base of the exponentiation.
+/// ## modulus
+/// Classical modulus used for all modular arithmetic.
+/// ## multiply_window_size
+/// Number of source-register bits used per multiplication lookup window.
+/// ## exponent_window_size
+/// Number of exponent bits used per exponentiation window.
+operation WindowModularExp(
+    q_target : Qubit[],
+    q_exponent : Qubit[],
+    base : BigInt,
+    modulus : BigInt,
+    multiply_window_size : Int,
+    exponent_window_size : Int
+) : Unit {
+    let result_size = Length(q_target);
+    let exponent_size = Length(q_exponent);
+    Fact(multiply_window_size <= result_size, "multiply_window_size too large");
+    Fact(exponent_window_size <= exponent_size, "exponent_window_size too large");
+
+    let num_exponent_windows = DivCeil(exponent_size, exponent_window_size);
+    let base_sqs = ComputeSequentialSquares(base, modulus, exponent_size);
+    use q_minus_reg = Qubit[result_size];
+
+    LoopA(num_exponent_windows, window_idx => {
+        let i = window_idx * exponent_window_size;
+        let window_end = Math.MinI(i + exponent_window_size, exponent_size);
+        let q_exponent_window = q_exponent[i..window_end - 1];
+        let adjusted_base = base_sqs[i];  // base^(2^i) % modulus.
+
+        // Determine q_plus and q_minus based on window parity.
+        let (q_plus, q_minus) = if (window_idx % 2 == 0) {
+            (q_target, q_minus_reg)
+        } else {
+            (q_minus_reg, q_target)
+        };
+
+        // Forward pass: q_minus += (q_plus * adjusted_base) % modulus.
+        WindowModularMultiply(
+            q_exponent_window,
+            q_plus,
+            q_minus,
+            multiply_window_size,
+            adjusted_base,
+            modulus,
+            1L
+        );
+
+        // Inverse pass: q_plus -= q_minus * inv(adjusted_base) (mod modulus).
+        let adjusted_base_inv = Math.InverseModL(adjusted_base, modulus);
+        WindowModularMultiply(
+            q_exponent_window,
+            q_minus,
+            q_plus,
+            multiply_window_size,
+            adjusted_base_inv,
+            modulus,
+            -1L
+        );
+    });
+
+    // If num_exponent_windows is odd, swap q_minus_reg back to q_target.
+    if (num_exponent_windows % 2 == 1) {
+        ParallelSWAP(q_minus_reg, q_target);
+    }
+}
+
+export ComputeBasePowers, WindowModularMultiply, WindowModularExp;
diff --git a/library/arithmetic/src/ModMul.qs b/library/arithmetic/src/ModMul.qs
new file mode 100644
index 00000000000..99de1326b5f
--- /dev/null
+++ b/library/arithmetic/src/ModMul.qs
@@ -0,0 +1,115 @@
+// Copyright (c) Microsoft Corporation.
+// Licensed under the MIT License.
+
+import Std.Diagnostics.Fact;
+
+import AddConst.AddConstant;
+import ModAdd.ModAdd;
+import ResourceEstimation.LoopA;
+import Utils.RotateLeft;
+
+/// # Summary
+/// Computes `x := (2 * x) % modulus` in place.
+/// Requires an odd modulus satisfying `0 < modulus < 2^Length(x)`.
+///
+/// # Reference
+/// - [1](https://arxiv.org/abs/1706.06752) "Quantum resource estimates for computing
+///   elliptic curve discrete logarithms", Martin Roetteler, Michael Naehrig, Krysta M.
+///   Svore, Kristin Lauter. (Fig. 4).
+///
+/// # Input
+/// ## x
+/// Register storing `x`, updated in place.
+/// ## modulus
+/// Odd classical modulus satisfying `0 < modulus < 2^Length(x)`.
+operation ModDouble(x : Qubit[], modulus : BigInt) : Unit is Adj {
+    let n = Length(x);
+    Fact(modulus % 2L == 1L, "Modulus must be odd.");
+    Fact(modulus > 0L and modulus < (1L <<< n), "Modulus must satisfy 0 < modulus < 2^Length(x).");
+    use carry = Qubit();
+
+    // Step 1. Multiply by 2 by bit shift.
+    RotateLeft(x + [carry]);
+
+    // Step 2. "-p".
+    AddConstant((1L <<< n) - modulus, x + [carry]);
+    X(carry);
+
+    // Step 3. Controlled "+p".
+    Controlled AddConstant([carry], (modulus, x));
+
+    // Step 4. Uncompute carry (with CNOT).
+    // The result 2x mod p is odd iff p was subtracted (since p is odd).
+    CNOT(x[0], carry);
+    X(carry);
+}
+
+/// # Summary
+/// Computes `(x, y, ans) -> (x, y, ((ans << (n - 1)) + x * y) % modulus)`.
+/// Assumes `ans` is initialized to `0` for multiplication semantics.
+/// Requires an odd modulus satisfying `0 < modulus < 2^Length(x)`.
+///
+/// # Reference
+/// - [1](https://arxiv.org/abs/1706.06752) "Quantum resource estimates for computing
+///   elliptic curve discrete logarithms", Martin Roetteler, Michael Naehrig, Krysta M.
+///   Svore, Kristin Lauter. (Fig. 5).
+///
+/// # Input
+/// ## x
+/// First multiplicand register.
+/// ## y
+/// Second multiplicand register.
+/// ## ans
+/// Accumulator register updated in place.
+/// ## modulus
+/// Odd classical modulus satisfying `0 < modulus < 2^Length(x)`.
+operation ModMul(x : Qubit[], y : Qubit[], ans : Qubit[], modulus : BigInt) : Unit is Adj {
+    let n = Length(x);
+    Fact(Length(y) == n, "Registers must have the same length.");
+    Fact(Length(ans) == n, "Registers must have the same length.");
+    Fact(modulus % 2L == 1L, "Modulus must be odd.");
+    Fact(modulus > 0L and modulus < (1L <<< n), "Modulus must satisfy 0 < modulus < 2^Length(x).");
+
+    LoopA(n, idx => {
+        Controlled ModAdd([x[n - 1 - idx]], (y, ans, modulus));
+        if idx != n - 1 {
+            ModDouble(ans, modulus);
+        }
+    });
+}
+
+/// # Summary
+/// Computes `(x, ans) -> (x, ((ans << (n - 1)) + x * x) % modulus)`.
+/// Requires an odd modulus satisfying `0 < modulus < 2^Length(x)`.
+///
+/// # Reference
+/// - [1](https://arxiv.org/abs/1706.06752) "Quantum resource estimates for computing
+///   elliptic curve discrete logarithms", Martin Roetteler, Michael Naehrig, Krysta M.
+///   Svore, Kristin Lauter. (Fig. 6).
+///
+/// # Input
+/// ## x
+/// Register storing the value to be squared.
+/// ## ans
+/// Accumulator register updated in place.
+/// ## modulus
+/// Odd classical modulus satisfying `0 < modulus < 2^Length(x)`.
+operation ModSquare(x : Qubit[], ans : Qubit[], modulus : BigInt) : Unit is Adj {
+    let n = Length(x);
+    Fact(Length(ans) == n, "Registers must have the same length.");
+    Fact(modulus % 2L == 1L, "Modulus must be odd.");
+    Fact(modulus > 0L and modulus < (1L <<< n), "Modulus must satisfy 0 < modulus < 2^Length(x).");
+    use x_copy = Qubit();
+
+    LoopA(n, idx => {
+        let bit_idx = n - 1 - idx;
+        CNOT(x[bit_idx], x_copy);
+        Controlled ModAdd([x_copy], (x, ans, modulus));
+        CNOT(x[bit_idx], x_copy);
+        if bit_idx != 0 {
+            ModDouble(ans, modulus);
+        }
+    });
+}
+
+export ModDouble, ModMul, ModSquare;
diff --git a/library/arithmetic/src/ModNegate.qs b/library/arithmetic/src/ModNegate.qs
new file mode 100644
index 00000000000..07b7f6130a6
--- /dev/null
+++ b/library/arithmetic/src/ModNegate.qs
@@ -0,0 +1,38 @@
+// Copyright (c) Microsoft Corporation.
+// Licensed under the MIT License.
+
+import AddConst.AddConstant;
+import Utils.CheckIfAllOnes;
+import Utils.CheckIfAllZero;
+
+/// # Summary
+/// Reversible, in-place modular negation of an integer modulo a constant
+/// integer modulus. Given an $n$-bit integer $x$ encoded in a little-endian
+/// register `xs` and a constant integer `modulus`, this operation computes
+/// $-x \bmod m$. The result is held in register `xs`.
+///
+/// # Input
+/// ## modulus
+/// Constant integer modulus.
+/// ## xs
+/// Qubit register encoding the integer `x`; replaced with `-x mod modulus`.
+///
+/// # Reference
+/// - [1](https://arxiv.org/pdf/2306.08585) "How to compute a 256-bit elliptic curve
+///   private key with only 50 million Toffoli gates",  Daniel Litinski. (Fig. 4).
+operation ModNegate(xs : Qubit[], modulus : BigInt) : Unit is Adj + Ctl {
+    body (...) {
+        (Controlled ModNegate)([], (xs, modulus));
+    }
+    controlled (controls, ...) {
+        let negModulus = (1L <<< Length(xs)) - modulus - 1L;
+        use isAllZeros = Qubit();
+        (Controlled CheckIfAllZero)(controls, (xs, isAllZeros));
+        (Controlled ApplyXorInPlaceL)([isAllZeros], (modulus, xs));
+        (Controlled AddConstant)(controls, (negModulus, xs));
+        CheckIfAllOnes(controls + xs, isAllZeros);
+        (Controlled ApplyToEachCA)(controls, (X, xs));
+    }
+}
+
+export ModNegate;
diff --git a/library/arithmetic/src/MultiControl.qs b/library/arithmetic/src/MultiControl.qs
new file mode 100644
index 00000000000..056a310da65
--- /dev/null
+++ b/library/arithmetic/src/MultiControl.qs
@@ -0,0 +1,50 @@
+// Copyright (c) Microsoft Corporation.
+// Licensed under the MIT License.
+
+import Std.Diagnostics.Fact;
+import Std.Arithmetic.RippleCarryCGAddLE;
+import Std.Arithmetic.RippleCarryTTKIncByLE;
+
+/// # Summary
+/// Flips the output qubit if and only if all the input qubits are 1.
+///
+/// # Resources
+/// For n>=2, uses n-2 auxiliary qubits, n-2 AND/IAND pairs and 1 CCNOT gate.
+///
+/// # Input
+/// ## controlQubits
+/// Input qubits to check.
+/// ## output
+/// Output qubit to flip.
+///
+/// # Reference
+/// - Michael A. Nielsen and Isaac L. Chuang,
+///   "Quantum Computation and Quantum Information" (Section 4.3, Fig. 4.10).
+operation MultiControl(controlQubits : Qubit[], output : Qubit) : Unit is Adj + Ctl {
+    body (...) {
+        let nQubits = Length(controlQubits);
+        if (nQubits == 0) {
+            // Do nothing.
+        } elif (nQubits == 1) {
+            CNOT(controlQubits[0], output);
+        } elif (nQubits == 2) {
+            CCNOT(controlQubits[0], controlQubits[1], output);
+        } else {
+            use anc = Qubit[nQubits-2];
+            within {
+                AND(controlQubits[0], controlQubits[1], anc[0]);
+                for i in 0..nQubits-4 {
+                    AND(anc[i], controlQubits[i + 2], anc[i + 1]);
+                }
+            } apply {
+                CCNOT(anc[nQubits-3], controlQubits[nQubits-1], output);
+            }
+        }
+    }
+    controlled (controls, ...) {
+        MultiControl(controls + controlQubits, output);
+    }
+    adjoint self;
+}
+
+export MultiControl;
diff --git a/library/arithmetic/src/ResourceEstimation.qs b/library/arithmetic/src/ResourceEstimation.qs
new file mode 100644
index 00000000000..f9623575667
--- /dev/null
+++ b/library/arithmetic/src/ResourceEstimation.qs
@@ -0,0 +1,56 @@
+// Copyright (c) Microsoft Corporation.
+// Licensed under the MIT License.
+
+import Std.ResourceEstimation.BeginEstimateCaching;
+import Std.ResourceEstimation.EndEstimateCaching;
+import Std.ResourceEstimation.IsResourceEstimating;
+import Std.ResourceEstimation.RepeatEstimates;
+
+/// # Summary
+/// Repeats an operation `num_iterations` times with resource-estimation-friendly behavior.
+///
+/// When running under resource estimation, this operation uses `RepeatEstimates`
+/// and executes `iteration(0)` once to model the loop body cost without
+/// classically iterating through all indices. During simulation/execution, it
+/// executes `iteration(i)` for each `i` in `0 .. num_iterations - 1`.
+///
+/// # Input
+/// ## num_iterations
+/// Number of loop iterations.
+/// ## iteration
+/// Operation that implements one loop iteration and receives the iteration index.
+operation Loop(num_iterations : Int, iteration : ((Int) => Unit)) : Unit {
+    if (num_iterations == 0) {
+        // Do nothing.
+    } elif (IsResourceEstimating()) {
+        within {
+            RepeatEstimates(num_iterations);
+        } apply {
+            iteration(0);
+        }
+    } else {
+        for i in 0..num_iterations-1 {
+            iteration(i);
+        }
+    }
+}
+
+/// # Summary
+/// Adjointable variant of `Loop` for adjointable iteration operations.
+operation LoopA(num_iterations : Int, iteration : ((Int) => Unit is Adj)) : Unit is Adj {
+    if (num_iterations == 0) {
+        // Do nothing.
+    } elif (IsResourceEstimating()) {
+        within {
+            RepeatEstimates(num_iterations);
+        } apply {
+            iteration(0);
+        }
+    } else {
+        for i in 0..num_iterations-1 {
+            iteration(i);
+        }
+    }
+}
+
+export Loop, LoopA;
diff --git a/library/arithmetic/src/Utils.qs b/library/arithmetic/src/Utils.qs
new file mode 100644
index 00000000000..f9b6566e2a8
--- /dev/null
+++ b/library/arithmetic/src/Utils.qs
@@ -0,0 +1,58 @@
+// Copyright (c) Microsoft Corporation.
+// Licensed under the MIT License.
+
+/// Simple helpers used in various arithmetic circuits.
+
+import Std.Diagnostics.Fact;
+
+import MultiControl.MultiControl;
+
+/// Swaps two same-length qubit arrays element-wise.
+operation ParallelSWAP(xs : Qubit[], ys : Qubit[]) : Unit is Adj + Ctl {
+    Fact(Length(xs) == Length(ys), "ParallelSWAP: registers must have equal length");
+    for idx in 0..Length(xs) - 1 {
+        SWAP(xs[idx], ys[idx]);
+    }
+}
+
+/// Applies CNOT between two registers.
+operation ParallelCNOT(controls : Qubit[], targets : Qubit[]) : Unit is Ctl + Adj {
+    let n : Int = Length(controls);
+    Fact(Length(targets) == n, "Size mismatch.");
+    for i in 0..n-1 {
+        CNOT(controls[i], targets[i]);
+    }
+}
+
+/// Cyclically rotates qubits left with SWAPs.
+/// If the high qubit is 0, this is multiplication by 2.
+operation RotateLeft(x : Qubit[]) : Unit is Adj + Ctl {
+    for i in Length(x) - 1.. -1..1 {
+        SWAP(x[i], x[i - 1]);
+    }
+}
+
+/// Cyclically rotates qubits right with SWAPs.
+/// If x stores an even number, this is division by 2.
+operation RotateRight(x : Qubit[]) : Unit is Adj + Ctl {
+    Adjoint RotateLeft(x);
+}
+
+/// Flips `output` if all qubits in `xs` are in |1> state.
+operation CheckIfAllOnes(xs : Qubit[], output : Qubit) : Unit is Ctl + Adj {
+    MultiControl(xs, output);
+}
+
+/// Flips `output` if all qubits in `xs` are in |0> state.
+operation CheckIfAllZero(xs : Qubit[], output : Qubit) : Unit is Adj + Ctl {
+    body (...) {
+        (Controlled CheckIfAllZero)([], (xs, output));
+    }
+    controlled (controls, ...) {
+        ApplyToEachCA(X, xs);
+        MultiControl(controls + xs, output);
+        ApplyToEachCA(X, xs);
+    }
+}
+
+export ParallelSWAP, ParallelCNOT, RotateLeft, RotateRight, CheckIfAllOnes, CheckIfAllZero;
diff --git a/library/arithmetic/test/__init__.py b/library/arithmetic/test/__init__.py
new file mode 100644
index 00000000000..e69de29bb2d
diff --git a/library/arithmetic/test/test_add.py b/library/arithmetic/test/test_add.py
new file mode 100644
index 00000000000..089df0ed580
--- /dev/null
+++ b/library/arithmetic/test/test_add.py
@@ -0,0 +1,35 @@
+"""Tests for Add.qs."""
+
+import random
+
+import pytest
+from qdk.test_utils import ArithmeticOpTester
+
+from test.test_utils import get_qdk_context
+
+
+@pytest.mark.parametrize("optimize", ["space", "time"])
+@pytest.mark.parametrize("n", [8, 16])
+def test_add(n: int, optimize: str):
+    """Tests AddConst.AddConstant."""
+    context = get_qdk_context(optimize=optimize)
+    modulus = 2**n
+    tester = ArithmeticOpTester("Add.Add", [n, n], context)
+    for _ in range(10):
+        a = random.randint(0, modulus - 1)
+        b = random.randint(0, modulus - 1)
+        ans = tester.run([a, b])
+        assert ans == [a, (a + b) % modulus]
+
+
+@pytest.mark.parametrize("optimize", ["space", "time"])
+@pytest.mark.parametrize("n", [8, 16])
+def test_add_resource_estimates(n: int, optimize: str):
+    context = get_qdk_context(optimize=optimize)
+    tester = ArithmeticOpTester("Add.Add", [n, n], context)
+    counts = context.logical_counts(tester.test_callable, [0, 0])
+    num_qubits, num_ccz = counts._data["numQubits"], counts._data["cczCount"]
+    if optimize == "space":
+        assert (num_qubits, num_ccz) == (2 * n, 2 * n - 2)
+    elif optimize == "time":
+        assert (num_qubits, num_ccz) == (3 * n, n - 1)
diff --git a/library/arithmetic/test/test_add_const.py b/library/arithmetic/test/test_add_const.py
new file mode 100644
index 00000000000..3addc777610
--- /dev/null
+++ b/library/arithmetic/test/test_add_const.py
@@ -0,0 +1,40 @@
+"""Tests for AddConst.qs."""
+
+import random
+
+import pytest
+from qdk.test_utils import ArithmeticOpTester
+
+from test.test_utils import get_qdk_context
+
+
+@pytest.mark.parametrize("optimize", ["space", "time"])
+@pytest.mark.parametrize("n", [1, 2, 3, 8, 16])
+def test_add_constant(n: int, optimize: str):
+    """Tests AddConst.AddConstant."""
+    context = get_qdk_context(optimize=optimize)
+    modulus = 2**n
+    for _ in range(10):
+        a = random.randint(0, modulus - 1)
+        b = random.randint(0, modulus - 1)
+        op = f"AddConst.AddConstant({a}L,_)"
+        ans = ArithmeticOpTester.run_unary_op(op, n, b, context)
+        assert ans == (a + b) % modulus
+
+
+@pytest.mark.parametrize("optimize", ["space", "time"])
+@pytest.mark.parametrize("n", [1, 2, 3, 8, 16])
+def test_add_constant_controlled(n: int, optimize: str):
+    """Tests controlled AddConst.AddConstant."""
+    context = get_qdk_context(optimize=optimize)
+    modulus = 2**n
+    for _ in range(10):
+        ctrl = random.randint(0, 1)
+        a = random.randint(0, modulus - 1)
+        b = random.randint(0, modulus - 1)
+        op = f"((ctrl, x) => (Controlled AddConst.AddConstant)(ctrl, ({a}L, x)))"
+        ans = ArithmeticOpTester.run_op(op, [1, n], [ctrl, b], context)
+        expected = (a + b) % modulus if ctrl == 1 else b
+        assert ans == [ctrl, expected]
+
+
diff --git a/library/arithmetic/test/test_add_lookup.py b/library/arithmetic/test/test_add_lookup.py
new file mode 100644
index 00000000000..fd37bc3e26d
--- /dev/null
+++ b/library/arithmetic/test/test_add_lookup.py
@@ -0,0 +1,75 @@
+"""Tests for AddLookup.qs."""
+
+import random
+
+from qdk.test_utils import ArithmeticOpTester
+
+from test.test_utils import get_qdk_context
+
+
+def _list_to_qs(data: list[int]) -> str:
+    return "[" + ",".join(f"{x}L" for x in data) + "]"
+
+
+def test_add_lookup_simple():
+    """Simple deterministic smoke test for AddLookup.AddLookup."""
+    context = get_qdk_context()
+    address_size = 2
+    n = 4
+    data = [1, 3, 5, 7]
+    op = f"AddLookup.AddLookup(_,_,{_list_to_qs(data)})"
+    op_tester = ArithmeticOpTester(op, arg_sizes=[address_size, n], context=context)
+
+    for address, y in [(0, 0), (1, 2), (2, 14), (3, 15)]:
+        expected = [address, (y + data[address]) % (2**n)]
+        assert op_tester.run([address, y]) == expected
+
+
+def test_add_lookup_non_mod():
+    """Tests AddLookup.AddLookup for addition modulo 2^n."""
+    context = get_qdk_context()
+    address_size = 3
+    n = 8
+    modulus = 2**n
+    table_length = 2**address_size
+
+    data = [random.randint(0, modulus - 1) for _ in range(table_length)]
+    op = f"AddLookup.AddLookup(_,_,{_list_to_qs(data)})"
+    op_tester = ArithmeticOpTester(op, arg_sizes=[address_size, n], context=context)
+
+    for _ in range(20):
+        address = random.randint(0, table_length - 1)
+        y = random.randint(0, modulus - 1)
+        assert op_tester.run([address, y]) == [address, (y + data[address]) % modulus]
+
+
+def test_parallel_mod_add_lookup():
+    """Tests AddLookup.ParallelModAddLookup for modular table adds."""
+    context = get_qdk_context()
+    address_size = 3
+    n = 8
+    modulus = 211
+    table_length = 2**address_size
+
+    table0 = [random.randint(0, 4 * modulus) for _ in range(table_length)]
+    table1 = [random.randint(0, 4 * modulus) for _ in range(table_length)]
+
+    table0_qs = _list_to_qs(table0)
+    table1_qs = _list_to_qs(table1)
+    op = (
+        "((address, target0, target1) => AddLookup.ParallelModAddLookup("
+        f"address, [target0, target1], [{table0_qs}, {table1_qs}], {modulus}L))"
+    )
+    op_tester = ArithmeticOpTester(op, arg_sizes=[address_size, n, n], context=context)
+
+    for _ in range(20):
+        address = random.randint(0, table_length - 1)
+        y0 = random.randint(0, modulus - 1)
+        y1 = random.randint(0, modulus - 1)
+
+        expected = [
+            address,
+            (y0 + (table0[address] % modulus)) % modulus,
+            (y1 + (table1[address] % modulus)) % modulus,
+        ]
+        assert op_tester.run([address, y0, y1]) == expected
diff --git a/library/arithmetic/test/test_compare.py b/library/arithmetic/test/test_compare.py
new file mode 100644
index 00000000000..3a6eed095d5
--- /dev/null
+++ b/library/arithmetic/test/test_compare.py
@@ -0,0 +1,82 @@
+"""Tests for Compare.qs."""
+
+import random
+from collections.abc import Callable
+
+import pytest
+from qdk.test_utils import ArithmeticOpTester
+
+from test.test_utils import get_qdk_context
+
+COMPARE_CASES = [
+    ("CompareGT", lambda x, y: x > y),
+    ("CompareLT", lambda x, y: x < y),
+    ("CompareGE", lambda x, y: x >= y),
+    ("CompareLE", lambda x, y: x <= y),
+    ("CompareEQ", lambda x, y: x == y),
+]
+
+
+@pytest.mark.parametrize("optimize", ["space", "time"])
+@pytest.mark.parametrize("n", [1, 2, 5, 20])
+@pytest.mark.parametrize(("op_name", "predicate"), COMPARE_CASES)
+def test_compare(
+    n: int, op_name: str, predicate: Callable, optimize: str
+):
+    """Test compare operations."""
+    context = get_qdk_context(optimize=optimize)
+    op = f"((x, y, result) => Compare.{op_name}(x, y, result[0]))"
+    tester = ArithmeticOpTester(op, [n, n, 1], context)
+
+    # Explicitly cover x == y corner case.
+    ans_bit_eq = random.randint(0, 1)
+    x_eq = random.randint(0, 2**n - 1)
+    result_eq = tester.run([x_eq, x_eq, ans_bit_eq])
+    assert result_eq == [
+        x_eq,
+        x_eq,
+        1 - ans_bit_eq if predicate(x_eq, x_eq) else ans_bit_eq,
+    ]
+
+    for _ in range(5):
+        ans_bit = random.randint(0, 1)
+        x = random.randint(0, 2**n - 1)
+        y = random.randint(0, 2**n - 1)
+
+        result = tester.run([x, y, ans_bit])
+        assert result == [x, y, 1 - ans_bit if predicate(x, y) else ans_bit]
+
+
+@pytest.mark.parametrize("optimize", ["space", "time"])
+@pytest.mark.parametrize("n", [1, 2, 5, 20])
+@pytest.mark.parametrize(("op_name", "predicate"), COMPARE_CASES)
+def test_compare_controlled(
+    n: int, op_name: str, predicate: Callable, optimize: str
+):
+    """Test controlled compare operations."""
+    context = get_qdk_context(optimize=optimize)
+    op = f"((ctrl, x, y, result) => (Controlled Compare.{op_name})(ctrl, (x, y, result[0])))"
+    tester = ArithmeticOpTester(op, [1, n, n, 1], context)
+
+    # Explicitly cover x == y corner case.
+    ans_bit_eq = random.randint(0, 1)
+    x_eq = random.randint(0, 2**n - 1)
+    result_eq_0 = tester.run([0, x_eq, x_eq, ans_bit_eq])
+    assert result_eq_0 == [0, x_eq, x_eq, ans_bit_eq]
+    result_eq_1 = tester.run([1, x_eq, x_eq, ans_bit_eq])
+    assert result_eq_1 == [
+        1,
+        x_eq,
+        x_eq,
+        1 - ans_bit_eq if predicate(x_eq, x_eq) else ans_bit_eq,
+    ]
+
+    for _ in range(5):
+        ans_bit = random.randint(0, 1)
+        x = random.randint(0, 2**n - 1)
+        y = random.randint(0, 2**n - 1)
+
+        result_0 = tester.run([0, x, y, ans_bit])
+        assert result_0 == [0, x, y, ans_bit]
+        result_1 = tester.run([1, x, y, ans_bit])
+        assert result_1 == [1, x, y, 1 - ans_bit if predicate(x, y) else ans_bit]
diff --git a/library/arithmetic/test/test_mod_add.py b/library/arithmetic/test/test_mod_add.py
new file mode 100644
index 00000000000..9587a2dcb9b
--- /dev/null
+++ b/library/arithmetic/test/test_mod_add.py
@@ -0,0 +1,74 @@
+import random
+
+import pytest
+from qdk.test_utils import ArithmeticOpTester
+
+from test.test_utils import get_qdk_context
+
+
+@pytest.mark.parametrize("optimize", ["space", "time"])
+@pytest.mark.parametrize("n", [2, 5, 20])
+def test_mod_add(n: int, optimize: str):
+    context = get_qdk_context(optimize=optimize)
+    for _ in range(10):
+        modulus = random.randint(2, 2**n - 1)
+        x = random.randint(0, modulus - 1)
+        y = random.randint(0, modulus - 1)
+        op = f"ModAdd.ModAdd(_,_,{modulus}L)"
+        result = ArithmeticOpTester.run_op(op, [n, n], [x, y], context)
+        assert result == [x, (x + y) % modulus]
+
+
+@pytest.mark.parametrize("optimize", ["space", "time"])
+@pytest.mark.parametrize("n", [2, 5, 20])
+def test_mod_add_controlled(n: int, optimize: str):
+    context = get_qdk_context(optimize=optimize)
+    for _ in range(10):
+        modulus = random.randint(2, 2**n - 1)
+        x = random.randint(0, modulus - 1)
+        y = random.randint(0, modulus - 1)
+
+        # Test with control = 1 (should apply ModAdd)
+        op = f"((c,x,y)=>Controlled ModAdd.ModAdd(c,(x,y,{modulus}L)))"
+        tester = ArithmeticOpTester(op, [1, n, n], context)
+        result = tester.run([1, x, y])
+        assert result == [1, x, (x + y) % modulus]
+
+        # Test with control = 0 (should not apply)
+        result = tester.run([0, x, y])
+        assert result == [0, x, y]
+
+
+@pytest.mark.parametrize("optimize", ["space", "time"])
+@pytest.mark.parametrize("n", [2, 5, 20])
+def test_mod_add_adjoint(n: int, optimize: str):
+    context = get_qdk_context(optimize=optimize)
+    for _ in range(10):
+        modulus = random.randint(2, 2**n - 1)
+        x = random.randint(0, modulus - 1)
+        y = random.randint(0, modulus - 1)
+
+        # Adjoint ModAdd should compute (x, (y - x) % modulus)
+        op = f"(Adjoint ModAdd.ModAdd(_,_,{modulus}L))"
+        result = ArithmeticOpTester.run_op(op, [n, n], [x, y], context)
+        assert result == [x, (y - x) % modulus]
+
+
+@pytest.mark.parametrize("optimize", ["space", "time"])
+@pytest.mark.parametrize("n", [2, 5, 20])
+def test_mod_add_controlled_adjoint(n: int, optimize: str):
+    context = get_qdk_context(optimize=optimize)
+    for _ in range(10):
+        modulus = random.randint(2, 2**n - 1)
+        x = random.randint(0, modulus - 1)
+        y = random.randint(0, modulus - 1)
+
+        # Test with control = 1 (should apply Adjoint ModAdd)
+        op = f"((c,x,y)=>Adjoint Controlled ModAdd.ModAdd(c,(x,y,{modulus}L)))"
+        tester = ArithmeticOpTester(op, [1, n, n], context)
+        result = tester.run([1, x, y])
+        assert result == [1, x, (y - x) % modulus]
+
+        # Test with control = 0 (should not apply)
+        result = tester.run([0, x, y])
+        assert result == [0, x, y]
diff --git a/library/arithmetic/test/test_mod_div.py b/library/arithmetic/test/test_mod_div.py
new file mode 100644
index 00000000000..95c419a4975
--- /dev/null
+++ b/library/arithmetic/test/test_mod_div.py
@@ -0,0 +1,100 @@
+import random
+
+import pytest
+from qdk.test_utils import ArithmeticOpTester
+
+from test.test_utils import get_qdk_context
+
+
+def test_compress_garbage() -> None:
+    """Tests that the compression circuit always leaves last qubit in 0 state."""
+    context = get_qdk_context()
+
+    def _is_valid_input(x: int) -> bool:
+        for i in range(3):
+            b0 = (x >> (2 * i)) & 1
+            b1 = (x >> (2 * i + 1)) & 1
+            if b0 == 0 and b1 == 1:
+                return False
+        return True
+
+    tester = ArithmeticOpTester("ModDiv.CompressGarbage", [6], context)
+    for x in range(2**6):
+        if _is_valid_input(x):
+            result = tester.run([x])[0]
+            assert 0 <= result < 32
+
+
+@pytest.mark.parametrize("num_bits,modulus", [(3, 5), (4, 13), (5, 31), (6, 61)])
+def test_mod_mul(num_bits: int, modulus: int) -> None:
+    context = get_qdk_context()
+    op = f"ModDiv.ModMul(_,_,{modulus}L)"
+    tester = ArithmeticOpTester(op, [num_bits, num_bits], context)
+    for _ in range(5):
+        x = random.randint(1, modulus - 1)
+        y = random.randint(0, modulus - 1)
+        result = tester.run([x, y])
+        assert result == [x, (y * x) % modulus]
+
+
+@pytest.mark.parametrize("num_bits,modulus", [(3, 5), (4, 13), (5, 31), (6, 61)])
+def test_safe_mod_mul(num_bits: int, modulus: int) -> None:
+    context = get_qdk_context()
+    op = f"ModDiv.SafeModMul(_,_,{modulus}L)"
+    tester = ArithmeticOpTester(op, [num_bits, num_bits], context)
+
+    # Special branch: x == 0 should leave y unchanged.
+    for _ in range(5):
+        y = random.randint(0, modulus - 1)
+        result = tester.run([0, y])
+        assert result == [0, y]
+
+    # Regular branch: x > 0 should behave as modular multiplication.
+    for _ in range(5):
+        x = random.randint(1, modulus - 1)
+        y = random.randint(0, modulus - 1)
+        result = tester.run([x, y])
+        assert result == [x, (y * x) % modulus]
+
+
+@pytest.mark.parametrize("num_bits,modulus", [(3, 5), (4, 13), (5, 31), (6, 61)])
+def test_mod_div(num_bits: int, modulus: int) -> None:
+    context = get_qdk_context()
+    op = f"ModDiv.ModDiv(_,_,{modulus}L)"
+    tester = ArithmeticOpTester(op, [num_bits, num_bits], context)
+    for _ in range(5):
+        x = random.randint(1, modulus - 1)
+        y = random.randint(0, modulus - 1)
+        result = tester.run([x, y])
+        assert result == [x, (y * pow(x, -1, modulus)) % modulus]
+
+
+@pytest.mark.parametrize("num_bits,modulus", [(3, 5), (4, 13), (5, 31), (6, 61)])
+def test_safe_mod_div(num_bits: int, modulus: int) -> None:
+    context = get_qdk_context()
+    op = f"ModDiv.SafeModDiv(_,_,{modulus}L)"
+    tester = ArithmeticOpTester(op, [num_bits, num_bits], context)
+
+    # Special branch: x == 0 should leave y unchanged.
+    for _ in range(5):
+        y = random.randint(0, modulus - 1)
+        result = tester.run([0, y])
+        assert result == [0, y]
+
+    # Regular branch: x > 0 should behave as modular division.
+    for _ in range(5):
+        x = random.randint(1, modulus - 1)
+        y = random.randint(0, modulus - 1)
+        result = tester.run([x, y])
+        assert result == [x, (y * pow(x, -1, modulus)) % modulus]
+
+
+@pytest.mark.parametrize("num_bits,modulus", [(3, 5), (4, 13), (5, 31), (6, 61)])
+def test_mod_inv(num_bits: int, modulus: int) -> None:
+    context = get_qdk_context()
+    op = f"ModDiv.ModInv(_,_,{modulus}L)"
+    tester = ArithmeticOpTester(op, [num_bits, num_bits], context)
+    for _ in range(5):
+        x = random.randint(1, modulus - 1)
+        result = tester.run([x, 0])
+        assert result == [x, pow(x, -1, modulus)]
diff --git a/library/arithmetic/test/test_mod_exp.py b/library/arithmetic/test/test_mod_exp.py
new file mode 100644
index 00000000000..c7bfbdd1c7d
--- /dev/null
+++ b/library/arithmetic/test/test_mod_exp.py
@@ -0,0 +1,36 @@
+import math
+import random
+
+import pytest
+from qdk.test_utils import ArithmeticOpTester
+
+from test.test_utils import get_qdk_context
+
+
+def _random_coprime(modulus: int, attempts: int = 100) -> int:
+    for _ in range(attempts):
+        ans = random.randint(2, modulus - 1)
+        if math.gcd(ans, modulus) == 1:
+            return ans
+    raise RuntimeError(f"No coprime is found for {modulus} in {attempts} attempts.")
+
+
+@pytest.mark.parametrize(
+    ("target_size", "exp_size"),
+    [(2, 2), (5, 3), (6, 3), (7, 4), (10, 10)],
+)
+def test_window_modular_exp(target_size: int, exp_size: int):
+    """Tests for WindowModularExp."""
+    context = get_qdk_context()
+    mul_w = 2
+    exp_w = 2
+
+    for _ in range(5):
+        mod = 1 + 2 * random.randint(1, 2 ** (target_size - 1) - 1)
+        x = random.randint(0, mod - 1)
+        y = random.randint(0, 2**exp_size - 1)
+        a = _random_coprime(mod, attempts=100)
+        op_name = "ModExp.WindowModularExp"
+        op = f"{op_name}(_,_,{a}L,{mod}L,{mul_w},{exp_w})"
+        result = ArithmeticOpTester.run_op(op, [target_size, exp_size], [x, y], context)
+        assert result == [(x * (a**y)) % mod, y]
diff --git a/library/arithmetic/test/test_mod_mul.py b/library/arithmetic/test/test_mod_mul.py
new file mode 100644
index 00000000000..f8c6d0fef00
--- /dev/null
+++ b/library/arithmetic/test/test_mod_mul.py
@@ -0,0 +1,46 @@
+import random
+
+from qdk.test_utils import ArithmeticOpTester
+
+from test.test_utils import get_qdk_context
+
+
+def test_mod_double():
+    """Tests for ModDouble."""
+    context = get_qdk_context()
+    n = 20
+    for _ in range(10):
+        # The modulus is an odd number in the range [3, 2^n-1].
+        modulus = random.randint(0, 2 ** (n - 1) - 2) * 2 + 3
+        x = random.randint(0, modulus - 1)
+        op = f"ModMul.ModDouble(_,{modulus}L)"
+        assert ArithmeticOpTester.run_unary_op(op, n, x, context) == (2 * x) % modulus
+
+
+def test_mod_mul():
+    """Tests for ModMul."""
+    context = get_qdk_context()
+    n = 8
+    for _ in range(10):
+        # The modulus is an odd number in the range [3, 2^n-1].
+        modulus = random.randint(0, 2 ** (n - 1) - 2) * 2 + 3
+        x = random.randint(0, modulus - 1)
+        y = random.randint(0, modulus - 1)
+        z = random.randint(0, modulus - 1)
+        op = f"ModMul.ModMul(_,_,_,{modulus}L)"
+        result = ArithmeticOpTester.run_op(op, [n, n, n], [x, y, z], context)
+        assert result == [x, y, ((z << (n - 1)) + x * y) % modulus]
+
+
+def test_mod_square():
+    """Tests for ModSquare."""
+    context = get_qdk_context()
+    n = 8
+    for _ in range(10):
+        # The modulus is an odd number in the range [3, 2^n-1].
+        modulus = random.randint(0, 2 ** (n - 1) - 2) * 2 + 3
+        x = random.randint(0, modulus - 1)
+        y = random.randint(0, modulus - 1)
+        op = f"ModMul.ModSquare(_,_,{modulus}L)"
+        result = ArithmeticOpTester.run_op(op, [n, n], [x, y], context)
+        assert result == [x, ((y << (n - 1)) + (x * x)) % modulus]
diff --git a/library/arithmetic/test/test_mod_negate.py b/library/arithmetic/test/test_mod_negate.py
new file mode 100644
index 00000000000..df182b57b06
--- /dev/null
+++ b/library/arithmetic/test/test_mod_negate.py
@@ -0,0 +1,26 @@
+import random
+
+from qdk.test_utils import ArithmeticOpTester
+
+from test.test_utils import get_qdk_context
+
+
+def test_mod_negate():
+    """Tests for ModNegate."""
+    context = get_qdk_context()
+    n = 16
+    for _ in range(10):
+        modulus = random.randint(3, 2**n - 1)
+        x = random.randint(0, modulus - 1)
+        op = f"ModNegate.ModNegate(_,{modulus}L)"
+        result = ArithmeticOpTester.run_unary_op(op, n, x, context)
+        assert result == (-x) % modulus
+
+
+def test_mod_negate_zero():
+    """Tests for ModNegate when input is 0."""
+    context = get_qdk_context()
+    n = 16
+    modulus = random.randint(3, 2**n - 1)
+    op = f"ModNegate.ModNegate(_,{modulus}L)"
+    assert ArithmeticOpTester.run_unary_op(op, n, 0, context) == 0
diff --git a/library/arithmetic/test/test_multi_control.py b/library/arithmetic/test/test_multi_control.py
new file mode 100644
index 00000000000..c3918fe8981
--- /dev/null
+++ b/library/arithmetic/test/test_multi_control.py
@@ -0,0 +1,27 @@
+"""Tests for MultiControl.qs."""
+
+import random
+
+import pytest
+from qdk.test_utils import ArithmeticOpTester
+
+from test.test_utils import get_qdk_context
+
+
+@pytest.mark.parametrize("n", [1, 2, 3, 4, 5, 6, 7, 8, 10, 20])
+def test_multi_control(n: int):
+    """Correctness tests for MultiControl."""
+    context = get_qdk_context()
+    op = "MultiControl.MultiControl"
+    op = f"((ctrl, target) => {op}(ctrl, target[0]))"
+    values = []
+    if n <= 6:
+        values = list(range(2**n))
+    else:
+        values = [2**n - 1] + [random.randint(0, 2**n - 2) for _ in range(10)]
+    tester = ArithmeticOpTester(op, [n, 1], context)
+    for ctl_val in values:
+        target_val = random.randint(0, 1)
+        result = tester.run([ctl_val, target_val])
+        expected = target_val ^ 1 if (ctl_val == 2**n - 1) else target_val
+        assert result == [ctl_val, expected]
diff --git a/library/arithmetic/test/test_utils.py b/library/arithmetic/test/test_utils.py
new file mode 100644
index 00000000000..77108c641a3
--- /dev/null
+++ b/library/arithmetic/test/test_utils.py
@@ -0,0 +1,13 @@
+import functools
+from pathlib import Path
+
+from qdk import Context
+
+
+@functools.cache
+def get_qdk_context(*, optimize: str = "") -> Context:
+    # Context is cached to avoid re-compiling the library for each test.
+    path = str(Path(__file__).resolve().parents[1])
+    return Context(
+        project_root=path, qdk_config={"optimize": optimize}
+    )
diff --git a/source/vscode/src/registry.json b/source/vscode/src/registry.json
index b59c8c38dd9..986fddf65be 100644
--- a/source/vscode/src/registry.json
+++ b/source/vscode/src/registry.json
@@ -90,6 +90,18 @@
         }
       }
     },
+    {
+      "name": "Arithmetic",
+      "description": "Quantum arithmetic algorithms",
+      "dependency": {
+        "github": {
+          "owner": "microsoft",
+          "repo": "qdk",
+          "refs": [{ "ref": "main", "notes": "nightly, unstable" }],
+          "path": "library/arithmetic"
+        }
+      }
+    },
     {
       "name": "QuantumArithmetic",
       "description": "Quantum arithmetic algorithms",