Add more numerics functionality (#102)

* Adding high-level integer functionality and polynomial evaluation in fixed-point.

* Restructuring Numerics.

* Adapt to naming convention.

* Add Numerics to pack step.

* Tabs to spaces.

* Primitive -> Intrinsic

* Fix naming of 2's complement inversion.

* Improved controlled poly. evaluation.

* Update Numerics/src/FixedPoint/Types.qs

* Applying first round of suggestions.

* Fixed a typo in one test.

* Last couple changes from review and discussion.

Build/step-build-libs.yml

@@ -12,6 +12,7 @@ steps:
     projects: |
       $(LibrariesRootFolder)/Standard.sln
       $(LibrariesRootFolder)/Chemistry.sln
+      $(LibrariesRootFolder)/Numerics.sln
     arguments: >
         -c $(BuildConfiguration)
         -v n
@@ -32,6 +33,7 @@ steps:
      $(LibrariesRootFolder)/Chemistry/tests/ChemistryTests/ChemistryTests.csproj
      $(LibrariesRootFolder)/Chemistry/tests/SystemTests/SystemTests.csproj
      $(LibrariesRootFolder)/Chemistry/tests/DataModelTests/DataModelTests.csproj
+     $(LibrariesRootFolder)/Numerics/tests/NumericsTests.csproj
     arguments: >
         -c $(BuildConfiguration)
         -v n
@@ -50,10 +52,10 @@ steps:
     projects: |
      $(LibrariesRootFolder)/Standard/src/Standard.csproj
      $(LibrariesRootFolder)/Chemistry/src/DataModel/DataModel.csproj
+     $(LibrariesRootFolder)/Numerics/src/Numerics.csproj
     arguments: >
       --no-build
       -c $(BuildConfiguration)
       -v n
       -o $(System.DefaultWorkingDirectory)
       /p:PackageVersion=$(Nuget.Version)
-

Numerics.sln

@@ -0,0 +1,37 @@
+
+Microsoft Visual Studio Solution File, Format Version 12.00
+# Visual Studio 15
+VisualStudioVersion = 15.0.28307.539
+MinimumVisualStudioVersion = 10.0.40219.1
+Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "Numerics", "Numerics\src\Numerics.csproj", "{C7E6F1F9-5DB7-41D6-ADA4-3E8F5BAEF731}"
+EndProject
+Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "NumericsTests", "Numerics\tests\NumericsTests.csproj", "{B8E99F28-CF37-4AA1-BDEC-16D3573CAEFD}"
+EndProject
+Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "Standard", "Standard\src\Standard.csproj", "{4EE69CC7-1158-49A2-A28F-9443B06D4A32}"
+EndProject
+Global
+	GlobalSection(SolutionConfigurationPlatforms) = preSolution
+		Debug|Any CPU = Debug|Any CPU
+		Release|Any CPU = Release|Any CPU
+	EndGlobalSection
+	GlobalSection(ProjectConfigurationPlatforms) = postSolution
+		{C7E6F1F9-5DB7-41D6-ADA4-3E8F5BAEF731}.Debug|Any CPU.ActiveCfg = Debug|Any CPU
+		{C7E6F1F9-5DB7-41D6-ADA4-3E8F5BAEF731}.Debug|Any CPU.Build.0 = Debug|Any CPU
+		{C7E6F1F9-5DB7-41D6-ADA4-3E8F5BAEF731}.Release|Any CPU.ActiveCfg = Release|Any CPU
+		{C7E6F1F9-5DB7-41D6-ADA4-3E8F5BAEF731}.Release|Any CPU.Build.0 = Release|Any CPU
+		{B8E99F28-CF37-4AA1-BDEC-16D3573CAEFD}.Debug|Any CPU.ActiveCfg = Debug|Any CPU
+		{B8E99F28-CF37-4AA1-BDEC-16D3573CAEFD}.Debug|Any CPU.Build.0 = Debug|Any CPU
+		{B8E99F28-CF37-4AA1-BDEC-16D3573CAEFD}.Release|Any CPU.ActiveCfg = Release|Any CPU
+		{B8E99F28-CF37-4AA1-BDEC-16D3573CAEFD}.Release|Any CPU.Build.0 = Release|Any CPU
+		{4EE69CC7-1158-49A2-A28F-9443B06D4A32}.Debug|Any CPU.ActiveCfg = Debug|Any CPU
+		{4EE69CC7-1158-49A2-A28F-9443B06D4A32}.Debug|Any CPU.Build.0 = Debug|Any CPU
+		{4EE69CC7-1158-49A2-A28F-9443B06D4A32}.Release|Any CPU.ActiveCfg = Release|Any CPU
+		{4EE69CC7-1158-49A2-A28F-9443B06D4A32}.Release|Any CPU.Build.0 = Release|Any CPU
+	EndGlobalSection
+	GlobalSection(SolutionProperties) = preSolution
+		HideSolutionNode = FALSE
+	EndGlobalSection
+	GlobalSection(ExtensibilityGlobals) = postSolution
+		SolutionGuid = {62FF308E-50FB-430C-B3EE-BE8E142835C9}
+	EndGlobalSection
+EndGlobal

Numerics/src/FixedPoint/Addition.qs

@@ -0,0 +1,51 @@
+// Copyright (c) Microsoft Corporation. All rights reserved.
+// Licensed under the MIT License.
+
+namespace Microsoft.Quantum.Arithmetic {
+    open Microsoft.Quantum.Canon;
+
+    /// # Summary
+    /// Adds a classical constant to a quantum fixed-point number.
+    ///
+    /// # Input
+    /// ## constant
+    /// Constant to add to the quantum fixed-point number.
+    /// ## fp
+    /// Fixed-point number to which the constant will
+    /// be added.
+    operation AddConstantFxP(constant : Double, fp : FixedPoint) : Unit is Adj + Ctl {
+        let (px, xs) = fp!;
+        let n = Length(xs);
+        using (ys = Qubit[n]){
+            let tmpFp = FixedPoint(px, ys);
+            ApplyWithCA(PrepareFxP(constant, _), AddFxP(_, fp), tmpFp);
+        }
+    }
+
+    /// # Summary
+    /// Adds two fixed-point numbers stored in quantum registers.
+    ///
+    /// # Description
+    /// Given two fixed-point registers respectively in states $\ket{f_1}$ and $\ket{f_2}$,
+    /// performs the operation $\ket{f_1} \ket{f_2} \mapsto \ket{f_1} \ket{f_1 + f_2}$.
+    ///
+    /// # Input
+    /// ## fp1
+    /// First fixed-point number
+    /// ## fp2
+    /// Second fixed-point number, will be updated to contain the sum of the
+    /// two inputs.
+    ///
+    /// # Remarks
+    /// The current implementation requires the two fixed-point numbers
+    /// to have the same point position counting from the least-significant
+    /// bit, i.e., $n_i$ and $p_i$ must be equal.
+    operation AddFxP(fp1 : FixedPoint, fp2 : FixedPoint) : Unit is Adj + Ctl {
+        let (px, xs) = fp1!;
+        let (py, ys) = fp2!;
+
+        IdenticalPointPosFactFxP([fp1, fp2]);
+
+        AddI(LittleEndian(xs), LittleEndian(ys));
+    }
+}

Numerics/src/FixedPoint/Comparison.qs

@@ -0,0 +1,31 @@
+// Copyright (c) Microsoft Corporation. All rights reserved.
+// Licensed under the MIT License.
+
+namespace Microsoft.Quantum.Arithmetic {
+
+    /// # Summary
+    /// Compares two fixed-point numbers stored in quantum registers, and
+    /// controls a flip on the result.
+    ///
+    /// # Input
+    /// ## fp1
+    /// First fixed-point number to be compared.
+    /// ## fp2
+    /// Second fixed-point number to be compared.
+    /// ## result
+    /// Result of the comparison. Will be flipped if `fp1 > fp2`.
+    ///
+    /// # Remarks
+    /// The current implementation requires the two fixed-point numbers
+    /// to have the same point position and the same number of qubits.
+    operation CompareGreaterThanFxP(fp1 : FixedPoint, fp2 : FixedPoint,
+                                    result : Qubit) : Unit is Adj + Ctl {
+        let (px, xs) = fp1!;
+        let (py, ys) = fp2!;
+
+        IdenticalFormatFactFxP([fp1, fp2]);
+        CompareGTSI(SignedLittleEndian(LittleEndian(xs)),
+                    SignedLittleEndian(LittleEndian(ys)),
+                    result);
+    }
+}

Numerics/src/FixedPoint/Facts.qs

@@ -0,0 +1,63 @@
+// Copyright (c) Microsoft Corporation. All rights reserved.
+// Licensed under the MIT License.
+
+namespace Microsoft.Quantum.Arithmetic {
+    open Microsoft.Quantum.Diagnostics;
+    open Microsoft.Quantum.Arrays;
+
+    /// # Summary
+    /// Helper function to assert that a quantum fixed-point number is
+    /// initialized to zero, i.e., all qubits are in state $\ket{0}$.
+    operation AssertAllZeroFxP(fp : FixedPoint) : Unit is Adj + Ctl {
+        let (p, xs) = fp!;
+        AssertAllZero(xs);
+    }
+
+    /// # Summary
+    /// Assert that all fixed-point numbers in the provided array
+    /// have identical point positions and qubit numbers.
+    ///
+    /// # Input
+    /// ## fixedPoints
+    /// Array of quantum fixed-point numbers that will be checked for
+    /// compatibility (using assertions).
+    function IdenticalFormatFactFxP(fixedPoints : FixedPoint[]) : Unit {
+        if (Length(fixedPoints) == 0) {
+            return ();
+        }
+        let (position, register) = fixedPoints[0]!;
+        Fact(position > 0, "Point position must be greater than zero.");
+        let n = Length(register);
+        for (fp in Most(fixedPoints)) {
+            let (pos, reg) = fp!;
+            EqualityFactI(pos, position,
+                "FixedPoint numbers must have identical binary point position.");
+            EqualityFactI(Length(reg), n,
+                "FixedPoint numbers must have identical number of qubits.");
+        }
+    }
+
+    /// # Summary
+    /// Assert that all fixed-point numbers in the provided array
+    /// have identical point positions when counting from the least-
+    /// significant bit. I.e., number of bits minus point position must
+    /// be constant for all fixed-point numbers in the array.
+    ///
+    /// # Input
+    /// ## fixedPoints
+    /// Array of quantum fixed-point numbers that will be checked for
+    /// compatibility (using assertions).
+    function IdenticalPointPosFactFxP(fixedPoints : FixedPoint[]) : Unit {
+        if (Length(fixedPoints) == 0) {
+            return ();
+        }
+        let (position, register) = fixedPoints[0]!;
+        Fact(position > 0, "Point position must be greater than zero.");
+        let n = Length(register);
+        for (fp in Most(fixedPoints)) {
+            let (pos, reg) = fp!;
+            EqualityFactI(Length(reg) - pos, n - position,
+                "FixedPoint numbers must have identical point alignment.");
+        }
+    }
+}

Numerics/src/FixedPoint/Init.qs

@@ -0,0 +1,41 @@
+// Copyright (c) Microsoft Corporation. All rights reserved.
+// Licensed under the MIT License.
+
+namespace Microsoft.Quantum.Arithmetic {
+    open Microsoft.Quantum.Convert;
+    open Microsoft.Quantum.Intrinsic;
+    open Microsoft.Quantum.Math;
+
+    /// # Summary
+    /// Initialize a quantum fixed-point number to a classical constant.
+    ///
+    /// # Input
+    /// ## constant
+    /// Constant to which to initialize the quantum fixed-point number.
+    /// ## fp
+    /// Fixed-point number (of type FixedPoint) to initialize.
+    operation PrepareFxP(constant : Double, fp : FixedPoint) : Unit{
+        body (...) {
+            let (p, q) = fp!;
+            let n = Length(q);
+            let sign = constant < 0.;
+            mutable rescaledConstant = PowD(2., IntAsDouble(n-p)) * AbsD(constant) + 0.5;
+            mutable keepAdding = sign;
+            for (i in 0..n-1) {
+                let intConstant = Floor(rescaledConstant);
+                set rescaledConstant = 0.5 * rescaledConstant;
+                mutable currentBit = (intConstant &&& 1) == (sign ? 0 | 1);
+                if (keepAdding) {
+                    set keepAdding = currentBit;
+                    set currentBit = not currentBit;
+                }
+                if (currentBit) {
+                    X(q[i]);
+                }
+            }
+        }
+        controlled auto;
+        adjoint self;
+        adjoint controlled auto;
+    }
+}

Numerics/src/FixedPoint/Measurement.qs

@@ -0,0 +1,32 @@
+// Copyright (c) Microsoft Corporation. All rights reserved.
+// Licensed under the MIT License.
+
+namespace Microsoft.Quantum.Arithmetic {
+    open Microsoft.Quantum.Convert;
+    open Microsoft.Quantum.Measurement;
+    open Microsoft.Quantum.Math;
+
+    /// # Summary
+    /// Measure a fixed-point number, returns its value as Double, and resets
+    /// all the register to zero.
+    ///
+    /// # Input
+    /// ## fp
+    /// Fixed-point number to measure.
+    operation MeasureFxP(fp : FixedPoint) : Double {
+        let (p, xs) = fp!;
+        let n = Length(xs);
+        let sign = MResetZ(xs[n-1]) == One;
+        mutable keepAdding = sign;
+        mutable fpAsDouble = 0.;
+        for (i in 0..n-2) {
+            mutable currentRes = MResetZ(xs[i]) == (sign ? Zero | One);
+            if (keepAdding) {
+                set keepAdding = currentRes;
+                set currentRes = not currentRes;
+            }
+            set fpAsDouble = fpAsDouble * 0.5 + (currentRes == true ? 1. | 0.);
+        }
+        return (sign ? -1.0 | 1.0) * fpAsDouble * PowD(2.0, IntAsDouble(p-2));
+    }
+}

Numerics/src/FixedPoint/Multiplication.qs

@@ -0,0 +1,85 @@
+// Copyright (c) Microsoft Corporation. All rights reserved.
+// Licensed under the MIT License.
+
+namespace Microsoft.Quantum.Arithmetic {
+    open Microsoft.Quantum.Arrays;
+    open Microsoft.Quantum.Intrinsic;
+    open Microsoft.Quantum.Canon;
+
+    /// # Summary
+    /// Multiplies two fixed-point numbers in quantum registers.
+    ///
+    /// # Input
+    /// ## fp1
+    /// First fixed-point number.
+    /// ## fp2
+    /// Second fixed-point number.
+    /// ## result
+    /// Result fixed-point number, must be in state $\ket{0}$ initially.
+    ///
+    /// # Remarks
+    /// The current implementation requires the three fixed-point numbers
+    /// to have the same point position and the same number of qubits.
+    operation MultiplyFxP(fp1 : FixedPoint, fp2 : FixedPoint,
+                                       result : FixedPoint) : Unit is Adj {
+
+        body(...) {
+            (Controlled MultiplyFxP) (new Qubit[0],
+                                                   (fp1, fp2, result));
+        }
+        controlled (controls, ...){
+            IdenticalFormatFactFxP([fp1, fp2, result]);
+            AssertAllZeroFxP(result);
+            let (px, xs) = fp1!;
+            let (py, ys) = fp2!;
+            let (pz, zs) = result!;
+            let n = Length(xs);
+
+            using (tmpResult = Qubit[2*n]){
+                let xsInt = SignedLittleEndian(LittleEndian(xs));
+                let ysInt = SignedLittleEndian(LittleEndian(ys));
+                let tmpResultInt = SignedLittleEndian(
+                    LittleEndian(tmpResult));
+                MultiplySI(xsInt, ysInt, tmpResultInt);
+                (Controlled ApplyToEachCA)(controls,
+                                           (CNOT,
+                                            Zip(tmpResult[n-px..2*n-px-1], zs)));
+                (Adjoint MultiplySI)(xsInt, ysInt, tmpResultInt);
+            }
+        }
+    }
+
+    /// # Summary
+    /// Squares a fixed-point number.
+    ///
+    /// # Input
+    /// ## fp
+    /// Fixed-point number.
+    /// ## result
+    /// Result fixed-point number,
+    /// must be in state $\ket{0}$ initially.
+    operation SquareFxP(fp : FixedPoint, result : FixedPoint) : Unit is Adj {
+        body(...) {
+            (Controlled SquareFxP) (new Qubit[0],
+                                           (fp, result));
+        }
+        controlled (controls, ...){
+            IdenticalFormatFactFxP([fp, result]);
+            AssertAllZeroFxP(result);
+            let (px, xs) = fp!;
+            let (py, ys) = result!;
+            let n = Length(xs);
+
+            using (tmpResult = Qubit[2*n]){
+                let xsInt = SignedLittleEndian(LittleEndian(xs));
+                let tmpResultInt = SignedLittleEndian(
+                    LittleEndian(tmpResult));
+                SquareSI(xsInt, tmpResultInt);
+                (Controlled ApplyToEachCA)(controls,
+                                           (CNOT,
+                                            Zip(tmpResult[n-px..2*n-px-1], ys)));
+                (Adjoint SquareSI)(xsInt, tmpResultInt);
+            }
+        }
+    }
+}

Numerics/src/FixedPoint/Polynomial.qs

@@ -0,0 +1,148 @@
+// Copyright (c) Microsoft Corporation. All rights reserved.
+// Licensed under the MIT License.
+
+namespace Microsoft.Quantum.Arithmetic {
+    open Microsoft.Quantum.Canon;
+
+    /// # Summary
+    /// Evaluates a polynomial in a fixed-point representation.
+    ///
+    /// # Input
+    /// ## coefficients
+    /// Coefficients of the polynomial as a double array, i.e., the array
+    /// $[a_0, a_1, ..., a_d]$ for the polynomial
+    /// $P(x) = a_0 + a_1 x + \cdots + a_d x^d$.
+    /// ## fpx
+    /// Input fixed-point number for which to evaluate the polynomial.
+    /// ## result
+    /// Output fixed-point number which will hold $P(x)$. Must be in state
+    /// $\ket{0}$ initially.
+    operation EvaluatePolynomialFxP(coefficients : Double[], fpx : FixedPoint,
+                                   result : FixedPoint) : Unit is Adj {
+        body (...) {
+            (Controlled EvaluatePolynomialFxP) (new Qubit[0],
+                (coefficients, fpx, result));
+        }
+        controlled (controls, ...) {
+            IdenticalFormatFactFxP([fpx, result]);
+            AssertAllZeroFxP(result);
+            let degree = Length(coefficients) - 1;
+            let (p, q) = fpx!;
+            let n = Length(q);
+            if (degree == 0){
+                (Controlled PrepareFxP)(controls,
+                    (coefficients[0], result));
+            }
+            elif (degree > 0) {
+                // initialize ancillary register to a_d
+                using (qubits = Qubit[n * degree]){
+                    let firstIterate = FixedPoint(p,
+                        qubits[(degree-1)*n..degree*n-1]);
+                    PrepareFxP(coefficients[degree], firstIterate);
+                    for (d in degree..(-1)..2) {
+                        let currentIterate = FixedPoint(p, qubits[(d-1)*n..d*n-1]);
+                        let nextIterate = FixedPoint(p, qubits[(d-2)*n..(d-1)*n-1]);
+                        // multiply by x and then add current coefficient
+                        MultiplyFxP(currentIterate, fpx, nextIterate);
+                        AddConstantFxP(coefficients[d-1], nextIterate);
+                    }
+                    let finalIterate = FixedPoint(p, qubits[0..n-1]);
+                    // final multiplication into the result register
+                    (Controlled MultiplyFxP)(controls, (finalIterate, fpx, result));
+                    // add a_0 to complete polynomial evaluation and
+                    (Controlled AddConstantFxP)(controls,
+                        (coefficients[0], result));
+                    // uncompute intermediate results
+                    for (d in 2..degree) {
+                        let currentIterate = FixedPoint(p, qubits[(d-1)*n..d*n-1]);
+                        let nextIterate = FixedPoint(p, qubits[(d-2)*n..(d-1)*n-1]);
+                        (Adjoint AddConstantFxP)(coefficients[d-1], nextIterate);
+                        (Adjoint MultiplyFxP)(currentIterate, fpx,
+                                              nextIterate);
+                    }
+                    PrepareFxP(coefficients[degree], firstIterate);
+                }
+            }
+        }
+    }
+
+    /// # Summary
+    /// Evaluates an even polynomial in a fixed-point representation.
+    ///
+    /// # Input
+    /// ## coefficients
+    /// Coefficients of the even polynomial as a double array, i.e., the array
+    /// $[a_0, a_1, ..., a_k]$ for the even polynomial
+    /// $P(x) = a_0 + a_1 x^2 + \cdots + a_k x^{2k}$.
+    /// ## fpx
+    /// Input fixed-point number for which to evaluate the polynomial.
+    /// ## result
+    /// Output fixed-point number which will hold $P(x)$. Must be in state
+    /// $\ket{0}$ initially.
+    operation EvaluateEvenPolynomialFxP(coefficients : Double[], fpx : FixedPoint,
+                                       result : FixedPoint) : Unit is Adj {
+        body (...) {
+            (Controlled EvaluateEvenPolynomialFxP) (new Qubit[0],
+                (coefficients, fpx, result));
+        }
+        controlled (controls, ...) {
+            IdenticalFormatFactFxP([fpx, result]);
+            AssertAllZeroFxP(result);
+            let halfDegree = Length(coefficients) - 1;
+            let (p, q) = fpx!;
+            let n = Length(q);
+
+            if (halfDegree == 0){
+                (Controlled PrepareFxP)(controls,
+                    (coefficients[0], result));
+            }
+            elif (halfDegree > 0) {
+                // initialize ancillary register to a_d
+                using (xsSquared = Qubit[n]){
+                    let fpxSquared = FixedPoint(p, xsSquared);
+                    ApplyWithCA(SquareFxP(fpx, _),
+                        (Controlled EvaluatePolynomialFxP)(controls,
+                            (coefficients, _, result)),
+                        fpxSquared);
+                }
+            }
+        }
+    }
+
+    /// # Summary
+    /// Evaluates an odd polynomial in a fixed-point representation.
+    ///
+    /// # Input
+    /// ## coefficients
+    /// Coefficients of the odd polynomial as a double array, i.e., the array
+    /// $[a_0, a_1, ..., a_k]$ for the odd polynomial
+    /// $P(x) = a_0 x + a_1 x^3 + \cdots + a_k x^{2k+1}$.
+    /// ## fpx
+    /// Input fixed-point number for which to evaluate the polynomial.
+    /// ## result
+    /// Output fixed-point number which will hold P(x). Must be in state
+    /// $\ket{0}$ initially.
+    operation EvaluateOddPolynomialFxP(coefficients : Double[], fpx : FixedPoint,
+                                      result : FixedPoint) : Unit is Adj {
+        body (...) {
+            (Controlled EvaluateOddPolynomialFxP) (new Qubit[0],
+                (coefficients, fpx, result));
+        }
+        controlled (controls, ...) {
+            IdenticalFormatFactFxP([fpx, result]);
+            AssertAllZeroFxP(result);
+            let halfDegree = Length(coefficients) - 1;
+            let (p, q) = fpx!;
+            let n = Length(q);
+            if (halfDegree >= 0) {
+                using (tmpResult = Qubit[n]) {
+                    let tmpResultFp = FixedPoint(p, tmpResult);
+                    ApplyWithCA(EvaluateEvenPolynomialFxP(coefficients, _, _),
+                           (Controlled MultiplyFxP)(controls,
+                                                    (_, _, result)),
+                           (fpx, tmpResultFp));
+                }
+            }
+        }
+    }
+}

Numerics/src/FixedPoint/Reciprocal.qs

@@ -0,0 +1,44 @@
+// Copyright (c) Microsoft Corporation. All rights reserved.
+// Licensed under the MIT License.
+
+namespace Microsoft.Quantum.Arithmetic {
+    open Microsoft.Quantum.Math;
+    open Microsoft.Quantum.Canon;
+    open Microsoft.Quantum.Arrays;
+    open Microsoft.Quantum.Intrinsic;
+    open Microsoft.Quantum.Diagnostics;
+
+    /// # Summary
+    /// Computes $1/x$ for a fixed-point number $x$.
+    ///
+    /// # Input
+    /// ## x
+    /// Fixed-point number to be inverted.
+    /// ## result
+    /// Fixed-point number that will hold the result. Must be initialized to $\ket{0.0}$.
+    operation ComputeReciprocalFxP(x : FixedPoint, result : FixedPoint) : Unit is Adj {
+        body (...) {
+            (Controlled ComputeReciprocalFxP) (new Qubit[0], (x, result));
+        }
+        controlled (controls, ...) {
+            let (p, xs) = x!;
+            let (pRes, rs) = result!;
+            let n = Length(xs);
+            AssertAllZero(rs);
+            Fact(p + pRes - 1 + n >= Length(rs), "Output register is too wide.");
+            using ((sign, tmpRes) = (Qubit(), Qubit[2*n])) {
+                CNOT(Tail(xs), sign);
+                (Controlled Invert2sSI)
+                    ([sign], SignedLittleEndian(LittleEndian(xs)));
+                ComputeReciprocalI(LittleEndian(xs), LittleEndian(tmpRes));
+                (Controlled ApplyToEachCA)(controls,
+                    (CNOT, Zip(tmpRes[p+pRes-1+n-Length(rs)..Min([n+p+pRes, 2*n-1])], rs)));
+                (Controlled Invert2sSI)([sign], SignedLittleEndian(LittleEndian(rs)));
+                (Adjoint ComputeReciprocalI)(LittleEndian(xs), LittleEndian(tmpRes));
+                (Controlled Adjoint Invert2sSI)
+                    ([sign], SignedLittleEndian(LittleEndian(xs)));
+                CNOT(Tail(xs), sign);
+            }
+        }
+    }
+}

Numerics/src/FixedPoint/Types.qs

@@ -0,0 +1,10 @@
+// Copyright (c) Microsoft Corporation. All rights reserved.
+// Licensed under the MIT License.
+
+namespace Microsoft.Quantum.Arithmetic {
+    /// # Summary
+    /// Represents a register of qubits encoding a fixed-point number. Consists of an integer that is equal to the number of
+    /// qubits to the left of the binary point, i.e., qubits of weight greater
+    /// than or equal to 1, and a quantum register.
+    newtype FixedPoint = (Int, Qubit[]);
+}

Numerics/src/Integer/ArithmeticWrappers.qs

@@ -0,0 +1,71 @@
+// Copyright (c) Microsoft Corporation. All rights reserved.
+// Licensed under the MIT License.
+
+namespace Microsoft.Quantum.Arithmetic {
+    open Microsoft.Quantum.Intrinsic;
+    open Microsoft.Quantum.Arrays;
+
+    /// # Summary
+    /// Automatically chooses between addition with
+    /// carry and without, depending on the register size of `ys`.
+    ///
+    /// # Input
+    /// ## xs
+    /// $n$-bit addend.
+    /// ## ys
+    /// Addend with at least $n$ qubits. Will hold the result.
+    operation AddI (xs: LittleEndian, ys: LittleEndian) : Unit is Adj + Ctl {
+        if (Length(xs!) == Length(ys!)) {
+            RippleCarryAdderNoCarryTTK(xs, ys);
+        }
+        elif (Length(ys!) > Length(xs!)) {
+            using (qs = Qubit[Length(ys!) - Length(xs!) - 1]){
+                RippleCarryAdderTTK(LittleEndian(xs! + qs),
+                                    LittleEndian(Most(ys!)), Tail(ys!));
+            }
+        }
+        else {
+            fail "xs must not contain more qubits than ys!";
+        }
+    }
+
+    /// # Summary
+    /// Wrapper for integer comparison: `result = x > y`.
+    ///
+    /// # Input
+    /// ## xs
+    /// First $n$-bit number
+    /// ## ys
+    /// Second $n$-bit number
+    /// ## result
+    /// Will be flipped if $x > y$
+    operation CompareGTI (xs: LittleEndian, ys: LittleEndian,
+                            result: Qubit) : Unit is Adj + Ctl {
+        GreaterThan(xs, ys, result);
+    }
+
+    /// # Summary
+    /// Wrapper for signed integer comparison: `result = xs > ys`.
+    ///
+    /// # Input
+    /// ## xs
+    /// First $n$-bit number
+    /// ## ys
+    /// Second $n$-bit number
+    /// ## result
+    /// Will be flipped if $xs > ys$
+    operation CompareGTSI (xs: SignedLittleEndian,
+                           ys: SignedLittleEndian,
+                           result: Qubit) : Unit is Adj + Ctl {
+        using (tmp = Qubit()) {
+            CNOT(Tail(xs!!), tmp);
+            CNOT(Tail(ys!!), tmp);
+            X(tmp);
+            (Controlled CompareGTI)([tmp], (xs!, ys!, result));
+            X(tmp);
+            CCNOT(tmp, Tail(ys!!), result);
+            CNOT(Tail(xs!!), tmp);
+            CNOT(Tail(ys!!), tmp);
+        }
+    }
+}

Numerics/src/Integer/Division.qs

@@ -0,0 +1,56 @@
+// Copyright (c) Microsoft Corporation. All rights reserved.
+// Licensed under the MIT License.
+
+namespace Microsoft.Quantum.Arithmetic {
+    open Microsoft.Quantum.Intrinsic;
+    open Microsoft.Quantum.Diagnostics;
+
+    /// # Summary
+    /// Divides two quantum integers.
+    ///
+    /// # Description
+    /// `xs` will hold the
+    /// remainder `xs - floor(xs/ys) * ys` and `result` will hold
+    /// `floor(xs/ys)`.
+    ///
+    /// # Input
+    /// ## xs
+    /// $n$-bit dividend, will be replaced by the remainder.
+    /// ## ys
+    /// $n$-bit divisor
+    /// ## result
+    /// $n$-bit result, must be in state $\ket{0}$ initially
+    /// and will be replaced by the result of the integer division.
+    ///
+    /// # Remarks
+    /// Uses a standard shift-and-subtract approach to implement the division.
+    /// The controlled version is specialized such the subtraction does not
+    /// require additional controls.
+    operation DivideI (xs: LittleEndian, ys: LittleEndian,
+                               result: LittleEndian) : Unit {
+        body (...) {
+            (Controlled DivideI) (new Qubit[0], (xs, ys, result));
+        }
+        controlled (controls, ...) {
+            let n = Length(result!);
+
+            EqualityFactI(n, Length(ys!), "Integer division requires
+                           equally-sized registers ys and result.");
+            EqualityFactI(n, Length(xs!), "Integer division
+                            requires an n-bit dividend registers.");
+            AssertAllZero(result!);
+
+            let xpadded = LittleEndian(xs! + result!);
+
+            for (i in (n-1)..(-1)..0) {
+                let xtrunc = LittleEndian(xpadded![i..i+n-1]);
+                (Controlled CompareGTI) (controls, (ys, xtrunc, result![i]));
+                // if ys > xtrunc, we don't subtract:
+                (Controlled X) (controls, result![i]);
+                (Controlled Adjoint AddI) ([result![i]], (ys, xtrunc));
+            }
+        }
+        adjoint auto;
+        adjoint controlled auto;
+    }
+}

Numerics/src/Integer/Inversion.qs

@@ -0,0 +1,75 @@
+// Copyright (c) Microsoft Corporation. All rights reserved.
+// Licensed under the MIT License.
+
+namespace Microsoft.Quantum.Arithmetic {
+    open Microsoft.Quantum.Arrays;
+    open Microsoft.Quantum.Diagnostics;
+    open Microsoft.Quantum.Canon;
+    open Microsoft.Quantum.Intrinsic;
+
+    /// # Summary
+    /// Inverts a given integer modulo 2's complement.
+    ///
+    /// # Input
+    /// ## xs
+    /// n-bit signed integer (SignedLittleEndian), will be inverted modulo
+    /// 2's complement.
+    operation Invert2sSI (xs: SignedLittleEndian) : Unit {
+        body (...) {
+            (Controlled Invert2sSI) (new Qubit[0], xs);
+        }
+        controlled (controls, ...) {
+            ApplyToEachCA((Controlled X)(controls, _), xs!!);
+
+            using (ancillas = Qubit[Length(xs!!)]) {
+                (Controlled X)(controls, ancillas[0]);
+                AddI(LittleEndian(ancillas), xs!);
+                (Controlled X)(controls, ancillas[0]);
+            }
+        }
+        adjoint auto;
+        adjoint controlled auto;
+    }
+
+    /// # Summary
+    /// Computes the reciprocal 1/x for an unsigned integer x
+    /// using integer division. The result, interpreted as an integer,
+    /// will be `floor(2^(2*n-1) / x)`.
+    ///
+    /// # Input
+    /// ## xs
+    /// n-bit unsigned integer
+    /// ## result
+    /// 2n-bit output, must be in $\ket{0}$ initially.
+    ///
+    /// # Remarks
+    /// For the input x=0, the output will be all-ones.
+    operation ComputeReciprocalI (xs: LittleEndian,
+                                  result: LittleEndian) : Unit {
+        body (...) {
+            (Controlled ComputeReciprocalI) (new Qubit[0], (xs, result));
+        }
+        controlled (controls, ...) {
+            let n = Length(xs!);
+            AssertIntEqual(Length(result!), 2*n,
+                           "Result register must contain 2n qubits.");
+            AssertAllZero(result!);
+            using ((lhs, padding) = (Qubit[2*n], Qubit[n])) {
+                let paddedxs = LittleEndian(xs! + padding);
+                X(Tail(lhs)); // initialize left-hand side to 2^{2n-1}
+                // ... and divide:
+                (Controlled DivideI) (controls,
+                    (LittleEndian(lhs), paddedxs, result));
+                // uncompute lhs
+                for (i in 0..2*n-1) {
+                    (Controlled AddI) ([result![i]],
+                        (LittleEndian(paddedxs![0..2*n-1-i]),
+                         LittleEndian(lhs[i..2*n-1])));
+                }
+                X(Tail(lhs));
+            }
+        }
+        adjoint auto;
+        adjoint controlled auto;
+    }
+}

Numerics/src/Integer/Multiplication.qs

@@ -0,0 +1,177 @@
+// Copyright (c) Microsoft Corporation. All rights reserved.
+// Licensed under the MIT License.
+
+namespace Microsoft.Quantum.Arithmetic {
+    open Microsoft.Quantum.Arrays;
+    open Microsoft.Quantum.Intrinsic;
+    open Microsoft.Quantum.Diagnostics;
+
+    /// # Summary
+    /// Multiply integer `xs` by integer `ys` and store the result in `result`,
+    /// which must be zero initially.
+    ///
+    /// # Input
+    /// ## xs
+    /// $n$-bit multiplicand (LittleEndian)
+    /// ## ys
+    /// $n$-bit multiplier (LittleEndian)
+    /// ## result
+    /// $2n$-bit result (LittleEndian), must be in state $\ket{0}$ initially.
+    ///
+    /// # Remarks
+    /// Uses a standard shift-and-add approach to implement the multiplication.
+    /// The controlled version was improved by copying out $x_i$ to an ancilla
+    /// qubit conditioned on the control qubits, and then controlling the
+    /// addition on the ancilla qubit.
+    operation MultiplyI (xs: LittleEndian, ys: LittleEndian,
+                         result: LittleEndian) : Unit {
+        body (...) {
+            let n = Length(xs!);
+
+            EqualityFactI(n, Length(ys!), "Integer multiplication requires
+                           equally-sized registers xs and ys.");
+            EqualityFactI(2 * n, Length(result!), "Integer multiplication
+                            requires a 2n-bit result registers.");
+            AssertAllZero(result!);
+
+            for (i in 0..n-1) {
+                (Controlled AddI) ([xs![i]], (ys, LittleEndian(result![i..i+n])));
+            }
+        }
+        controlled (controls, ...) {
+            let n = Length(xs!);
+
+            EqualityFactI(n, Length(ys!), "Integer multiplication requires
+                           equally-sized registers xs and ys.");
+            EqualityFactI(2 * n, Length(result!), "Integer multiplication
+                            requires a 2n-bit result registers.");
+            AssertAllZero(result!);
+
+            using (anc = Qubit()) {
+                for (i in 0..n-1) {
+                    (Controlled CNOT) (controls, (xs![i], anc));
+                    (Controlled AddI) ([anc], (ys, LittleEndian(result![i..i+n])));
+                    (Controlled CNOT) (controls, (xs![i], anc));
+                }
+            }
+        }
+        adjoint auto;
+        adjoint controlled auto;
+    }
+
+    /// # Summary
+    /// Computes the square of the integer `xs` into `result`,
+    /// which must be zero initially.
+    ///
+    /// # Input
+    /// ## xs
+    /// $n$-bit number to square (LittleEndian)
+    /// ## result
+    /// $2n$-bit result (LittleEndian), must be in state $\ket{0}$ initially.
+    ///
+    /// # Remarks
+    /// Uses a standard shift-and-add approach to compute the square. Saves
+    /// $n-1$ qubits compared to the straight-forward solution which first
+    /// copies out xs before applying a regular multiplier and then undoing
+    /// the copy operation.
+    operation SquareI (xs: LittleEndian, result: LittleEndian) : Unit {
+        body (...) {
+            (Controlled SquareI) (new Qubit[0], (xs, result));
+        }
+        controlled (controls, ...) {
+            let n = Length(xs!);
+
+            EqualityFactI(2 * n, Length(result!), "Integer multiplication
+                            requires a 2n-bit result registers.");
+            AssertAllZero(result!);
+
+            using (anc = Qubit()) {
+                for (i in 0..n-1) {
+                    (Controlled CNOT) (controls, (xs![i], anc));
+                    (Controlled AddI) ([anc], (xs,
+                        LittleEndian(result![i..i+n])));
+                    (Controlled CNOT) (controls, (xs![i], anc));
+                }
+            }
+        }
+        adjoint auto;
+        adjoint controlled auto;
+    }
+
+    /// # Summary
+    /// Multiply signed integer `xs` by signed integer `ys` and store
+    /// the result in `result`, which must be zero initially.
+    ///
+    /// # Input
+    /// ## xs
+    /// n-bit multiplicand (SignedLittleEndian)
+    /// ## ys
+    /// n-bit multiplier (SignedLittleEndian)
+    /// ## result
+    /// 2n-bit result (SignedLittleEndian), must be in state $\ket{0}$
+    /// initially.
+    operation MultiplySI (xs: SignedLittleEndian,
+                          ys: SignedLittleEndian,
+                          result: SignedLittleEndian): Unit {
+        body (...) {
+            (Controlled MultiplySI) (new Qubit[0], (xs, ys, result));
+        }
+        controlled (controls, ...) {
+            let n = Length(xs!!);
+            using ((signx, signy) = (Qubit(), Qubit())) {
+                CNOT(Tail(xs!!), signx);
+                CNOT(Tail(ys!!), signy);
+                (Controlled Invert2sSI)([signx], xs);
+                (Controlled Invert2sSI)([signy], ys);
+
+                (Controlled MultiplyI) (controls, (xs!, ys!, result!));
+                CNOT(signx, signy);
+                // No controls required since `result` will still be zero
+                // if we did not perform the multiplication above.
+                (Controlled Invert2sSI)([signy], result);
+                CNOT(signx, signy);
+
+                (Controlled Adjoint Invert2sSI)([signx], xs);
+                (Controlled Adjoint Invert2sSI)([signy], ys);
+                CNOT(Tail(xs!!), signx);
+                CNOT(Tail(ys!!), signy);
+            }
+        }
+        adjoint auto;
+        adjoint controlled auto;
+    }
+
+    /// # Summary
+    /// Square signed integer `xs` and store
+    /// the result in `result`, which must be zero initially.
+    ///
+    /// # Input
+    /// ## xs
+    /// n-bit integer to square (SignedLittleEndian)
+    /// ## result
+    /// 2n-bit result (SignedLittleEndian), must be in state $\ket{0}$
+    /// initially.
+    ///
+    /// # Remarks
+    /// The implementation relies on IntegerSquare.
+    operation SquareSI (xs: SignedLittleEndian,
+                        result: SignedLittleEndian): Unit {
+        body (...) {
+            (Controlled SquareSI) (new Qubit[0], (xs, result));
+        }
+        controlled (controls, ...) {
+            let n = Length(xs!!);
+            using ((signx, signy) = (Qubit(), Qubit())) {
+                CNOT(Tail(xs!!), signx);
+                (Controlled Invert2sSI)([signx], xs);
+
+                (Controlled SquareI) (controls, (xs!, result!));
+
+                (Controlled Adjoint Invert2sSI)([signx], xs);
+                CNOT(Tail(xs!!), signx);
+            }
+        }
+        adjoint auto;
+        adjoint controlled auto;
+    }
+}

Numerics/src/Integer/Types.qs

@@ -0,0 +1,8 @@
+// Copyright (c) Microsoft Corporation. All rights reserved.
+// Licensed under the MIT License.
+
+namespace Microsoft.Quantum.Arithmetic {
+    /// # Summary
+    /// Type of a signed integer stored in little endian (see LittleEndian).
+    newtype SignedLittleEndian = LittleEndian;
+}

Numerics/src/Numerics.csproj

@@ -0,0 +1,28 @@
+<Project Sdk="Microsoft.NET.Sdk">
+  <PropertyGroup>
+    <TargetFramework>netstandard2.0</TargetFramework>
+    <AssemblyName>Microsoft.Quantum.Numerics</AssemblyName>
+    <PlatformTarget>x64</PlatformTarget>
+  </PropertyGroup>
+
+  <PropertyGroup>
+    <NoWarn>0162</NoWarn>
+    <RunQDocGen>True</RunQDocGen>
+    <Authors>Microsoft</Authors>
+    <Description>Microsoft's Quantum numerics libraries.</Description>
+    <Copyright>© Microsoft Corporation. All rights reserved.</Copyright>
+    <PackageReleaseNotes>See: https://docs.microsoft.com/en-us/quantum/relnotes/</PackageReleaseNotes>
+    <PackageLicenseUrl>https://github.com/Microsoft/QuantumLibraries/raw/master/LICENSE.txt</PackageLicenseUrl>
+    <PackageProjectUrl>https://github.com/Microsoft/QuantumLibraries</PackageProjectUrl>
+    <PackageIconUrl>https://secure.gravatar.com/avatar/bd1f02955b2853ba0a3b1cdc2434e8ec.png</PackageIconUrl>
+    <PackageTags>Quantum Q# Qsharp</PackageTags>
+  </PropertyGroup>
+
+  <ItemGroup>
+    <ProjectReference Include="..\..\Standard\src\Standard.csproj" />
+  </ItemGroup>
+
+  <ItemGroup>
+    <PackageReference Include="Microsoft.Quantum.Development.Kit" Version="0.6.1905.301" />
+  </ItemGroup>
+</Project>

Numerics/tests/FixedPointTests.qs

@@ -0,0 +1,376 @@
+// Copyright (c) Microsoft Corporation. All rights reserved.
+// Licensed under the MIT License.
+
+namespace Microsoft.Quantum.Numerics.ToffoliTests {
+    open Microsoft.Quantum.Arrays;
+    open Microsoft.Quantum.Convert;
+    open Microsoft.Quantum.Math;
+    open Microsoft.Quantum.Intrinsic;
+    open Microsoft.Quantum.Arithmetic;
+    open Microsoft.Quantum.Diagnostics;
+
+    operation PrepareFxPTest() : Unit {
+        for (a in [1.2, 3.9, 3.14159, -0.6, -4.5, -3.1931, 0.0]){
+            using (xs = Qubit[10]) {
+                let fp = FixedPoint(4, xs);
+                PrepareFxP(a, fp);
+                let measured = MeasureFxP(fp);
+                EqualityFactB(AbsD(measured - a) <= 1./IntAsDouble(2^7), true,
+                    $"FixedPoint initialized to {a} but measured {measured}.");
+                ResetAll(xs);
+            }
+        }
+    }
+
+    operation CompareGreaterThanFxPTest() : Unit {
+        for (a in [1.2, 3.9, 3.14159, -0.6, -4.5, -3.1931, 0.0]){
+            for (b in [1.1, 3.95, 3.14259, -0.4, -4.6, -3.931, 0.1]) {
+                using ((xs, ys, res) = (Qubit[10], Qubit[10], Qubit())) {
+                    let fp1 = FixedPoint(4, xs);
+                    let fp2 = FixedPoint(4, ys);
+                    PrepareFxP(a, fp1);
+                    PrepareFxP(b, fp2);
+                    CompareGreaterThanFxP(fp1, fp2, res);
+                    let measured = M(res);
+                    EqualityFactB(a > b, measured == One,
+                        $"FixedPoint comparison: {a} > {b} != {measured}.");
+                    ResetAll(xs + ys + [res]);
+                }
+            }
+        }
+    }
+
+    operation AddConstantFxPTest() : Unit {
+        for (a in [1.2, 3.9, 3.14159, -0.6, -4.5, -3.1931, 0.0]){
+            for (b in [1.2, 3.9, 3.14159, -0.6, -4.5, -3.1931, 0.0]){
+                using (xs = Qubit[11]) {
+                    let fp = FixedPoint(5, xs);
+                    PrepareFxP(a, fp);
+                    AddConstantFxP(b, fp);
+                    let measured = MeasureFxP(fp);
+                    EqualityWithinToleranceFact(measured, (a+b), 1. / IntAsDouble(2^6));
+                }
+            }
+        }
+    }
+
+    operation AddFxPTest() : Unit {
+        for (a in [1.2, 3.9, 3.14159, -0.6, -4.5, -3.1931, 0.0]){
+            for (b in [1.2, 3.9, 3.14159, -0.6, -4.5, -3.1931, 0.0]){
+                using ((xs, ys) = (Qubit[11], Qubit[11])) {
+                    let fp1 = FixedPoint(5, xs);
+                    let fp2 = FixedPoint(5, ys);
+                    PrepareFxP(a, fp1);
+                    PrepareFxP(b, fp2);
+                    AddFxP(fp1, fp2);
+                    let measured = MeasureFxP(fp2);
+                    EqualityWithinToleranceFact(measured, (a+b), 1. / IntAsDouble(2^6));
+                    ResetAll(xs + ys);
+                }
+            }
+        }
+    }
+
+    operation MultiplyFxPTest() : Unit {
+        for (pos in 5..8) {
+            for (a in [1.2, 3.9, 3.14159, -0.6, -3.5, -3.1931, 0.0]){
+                for (b in [1.2, 3.9, 3.14159, -0.6, -3.5, -3.1931, 0.0]){
+                    using ((xs, ys, zs) = (Qubit[13], Qubit[13], Qubit[13])) {
+                        let fp1 = FixedPoint(pos, xs);
+                        let fp2 = FixedPoint(pos, ys);
+                        let fp3 = FixedPoint(pos, zs);
+                        PrepareFxP(a, fp1);
+                        PrepareFxP(b, fp2);
+                        MultiplyFxP(fp1, fp2, fp3);
+                        let measured = MeasureFxP(fp3);
+                        let eps = 1./IntAsDouble(2^(13-pos));
+                        let epsTotal = AbsD(a) * eps + AbsD(b) * eps + eps * eps;
+                        EqualityWithinToleranceFact(measured, a * b, epsTotal);
+                        ResetAll(xs + ys + zs);
+                    }
+                }
+            }
+        }
+    }
+
+    operation SquareFxPTest() : Unit {
+        for (pos in 5..8) {
+            for (a in [1.2, 3.9, 3.14159, -0.6, -3.5, -3.1931, 0.0]){
+                using ((xs, ys) = (Qubit[13], Qubit[13])) {
+                    let fp1 = FixedPoint(pos, xs);
+                    let fp2 = FixedPoint(pos, ys);
+                    PrepareFxP(a, fp1);
+                    SquareFxP(fp1, fp2);
+                    let measured = MeasureFxP(fp2);
+                    let eps = 1./IntAsDouble(2^(13-pos));
+                    let epsTotal = 2. * AbsD(a) * eps + eps * eps;
+                    EqualityWithinToleranceFact(measured, a * a, epsTotal);
+                    ResetAll(xs + ys);
+                }
+            }
+        }
+    }
+
+    function _computeReciprocal(a : Double, n : Int, pos : Int, pos2 : Int) : Double {
+        let p = pos;
+        let intA = a >= 0. ? Floor(AbsD(a) * IntAsDouble(2^(n-p)) + 0.5)
+                           | Ceiling(AbsD(a) * IntAsDouble(2^(n-p)) - 0.5);
+        let intDiv = 2^(2*n-1) / intA;
+        let aReciprUnsigned = IntAsDouble((intDiv >>> (p+pos2-1)) &&& (2^n-1)) / IntAsDouble(2^(n-pos2));
+        return (a >= 0. ? 1. | -1.) * aReciprUnsigned;
+    }
+
+    operation ComputeReciprocalFxPTest() : Unit {
+        for (pos in 5..8) {
+            for (pos2 in pos-1..pos+3) {
+                for (a in [1.2, 3.9, -0.314159, -0.6, -3.5, -3.1931, 0.127]){
+                    let n = 20;
+                    using ((xs, ys) = (Qubit[n], Qubit[n])) {
+                        let fp1 = FixedPoint(pos, xs);
+                        let fp2 = FixedPoint(pos, ys);
+                        PrepareFxP(a, fp1);
+                        ComputeReciprocalFxP(fp1, fp2);
+                        let measured = MeasureFxP(fp2);
+                        let eps = 1./IntAsDouble(2^(n-pos));
+                        let eps2 = 1./IntAsDouble(2^(n-pos2));
+                        let aEpsLarger = a + (a>=0. ? eps | -eps);
+                        let aEpsSmaller = a - (a>=0. ? eps | -eps);
+                        let res1 = _computeReciprocal(a+eps,n,pos,pos2);
+                        let res2 = _computeReciprocal(a-eps,n,pos,pos2);
+                        let minRes = MinD(res1, res2) - eps2;
+                        let maxRes = MaxD(res1, res2) + eps2;
+                        let isWithinTol = minRes <= measured and
+                                          maxRes >= measured;
+                        EqualityFactB(isWithinTol,
+                            true,
+                            $"FixedPoint reciprocal 1/{a}: {measured} is not within [{minRes},{maxRes}].");
+                        ResetAll(xs + ys);
+                    }
+                }
+            }
+        }
+    }
+
+    operation SquareFxPCtrlTest() : Unit {
+        for (ctrl in 0..3) {
+            for (pos in 5..8) {
+                for (a in [1.2, 3.9, 3.14159, -0.6, -3.5, -3.1931, 0.0]){
+                    using ((xs, ys, cs) = (Qubit[13], Qubit[13], Qubit[2])) {
+                        ApplyXorInPlace(ctrl, LittleEndian(cs));
+                        let fp1 = FixedPoint(pos, xs);
+                        let fp2 = FixedPoint(pos, ys);
+                        PrepareFxP(a, fp1);
+                        (Controlled SquareFxP)(cs, (fp1, fp2));
+                        let measured = MeasureFxP(fp2);
+                        let eps = 1./IntAsDouble(2^(13-pos));
+                        let epsTotal = 2. * AbsD(a) * eps + eps * eps;
+                        if (ctrl == 3) {
+                            EqualityWithinToleranceFact(measured, a * a, epsTotal);
+                        }
+                        else {
+                            let measuredI = MeasureInteger(LittleEndian(ys));
+                            EqualityFactI(measuredI, 0,
+                                "Controlled FixedPoint square changed the result register!");
+                        }
+                        ResetAll(xs + ys + cs);
+                    }
+                }
+            }
+        }
+    }
+
+    operation MultiplyFxPCtrlTest() : Unit {
+        for (ctrl in 0..3) {
+            for (pos in 5..8) {
+                for (a in [1.2, 3.9, 3.14159, -0.6, -3.5, -3.1931, 0.0]){
+                    for (b in [1.2, 3.9, 3.14159, -0.6, -3.5, -3.1931, 0.0]){
+                        using ((xs, ys, zs, cs) = (Qubit[13], Qubit[13], Qubit[13], Qubit[2])) {
+                            ApplyXorInPlace(ctrl, LittleEndian(cs));
+                            let fp1 = FixedPoint(pos, xs);
+                            let fp2 = FixedPoint(pos, ys);
+                            let fp3 = FixedPoint(pos, zs);
+                            PrepareFxP(a, fp1);
+                            PrepareFxP(b, fp2);
+                            (Controlled MultiplyFxP)(cs, (fp1, fp2, fp3));
+                            let measured = MeasureFxP(fp3);
+                            let eps = 1./IntAsDouble(2^(13-pos));
+                            let epsTotal = AbsD(a) * eps + AbsD(b) * eps + eps * eps;
+                            if (ctrl == 3) {
+                                EqualityWithinToleranceFact(measured, a * b, epsTotal);
+                            }
+                            else {
+                                let measuredI = MeasureInteger(LittleEndian(zs));
+                                EqualityFactI(measuredI, 0,
+                                    "Controlled FixedPoint multiplication changed the output register!");
+                            }
+                            ResetAll(xs + ys + zs + cs);
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    operation EvaluatePolynomialFxPTest() : Unit {
+        for (pos in 4..5) {
+            for (coeffs in [[1.3, -2.4, 1.9],
+                            [-0.3, -0.2],
+                            [0.1, 1.1, -0.1, 0.2],
+                            [0.1, -0.1, 0.2, 0.2, -0.1, 0.3],
+                            [0.2]]){
+                for (a in [0.0, 0.1, -0.1, 0.2, -0.2, 1.3, -1.3]){
+                    let n = 20;
+                    using ((xs, ys) = (Qubit[n], Qubit[n])) {
+                        let fp1 = FixedPoint(pos, xs);
+                        let fp2 = FixedPoint(pos, ys);
+                        PrepareFxP(a, fp1);
+                        EvaluatePolynomialFxP(coeffs, fp1, fp2);
+                        let measured = MeasureFxP(fp2);
+                        let eps = 1./IntAsDouble(2^(n-pos));
+                        mutable epsTotal = 0.;
+                        mutable errX = eps;
+                        mutable result = Tail(coeffs);
+                        set epsTotal = epsTotal + eps;
+                        for (coeff in coeffs[Length(coeffs)-2..(-1)..0]) {
+                            set epsTotal = epsTotal + AbsD(result) * eps
+                                           + AbsD(a) * epsTotal + eps * epsTotal;
+                            set result = result * a + coeff;
+                            set epsTotal = epsTotal + eps;
+                        }
+                        EqualityWithinToleranceFact(measured, result, epsTotal);
+                        ResetAll(xs + ys);
+                    }
+                }
+            }
+        }
+    }
+
+    operation EvaluatePolynomialFxPCtrlTest() : Unit {
+        for (ctrl in 0..3) {
+            for (pos in 4..5) {
+                for (coeffs in [[1.3, -2.4, 1.9],
+                                [-0.3, -0.2],
+                                [0.2]]){
+                    for (a in [0.0, 0.1, -0.2, 1.3, -1.3]){
+                        let n = 20;
+                        using ((xs, ys, ctrls) = (Qubit[n], Qubit[n], Qubit[2])) {
+                            let fp1 = FixedPoint(pos, xs);
+                            let fp2 = FixedPoint(pos, ys);
+                            ApplyXorInPlace(ctrl, LittleEndian(ctrls));
+                            PrepareFxP(a, fp1);
+                            (Controlled EvaluatePolynomialFxP)(ctrls,
+                                (coeffs, fp1, fp2));
+                            let measured = MeasureFxP(fp2);
+                            let eps = 1./IntAsDouble(2^(n-pos));
+                            mutable epsTotal = 0.;
+                            mutable errX = eps;
+                            mutable result = Tail(coeffs);
+                            set epsTotal = epsTotal + eps;
+                            for (coeff in coeffs[Length(coeffs)-2..(-1)..0]) {
+                                set epsTotal = epsTotal + AbsD(result) * eps
+                                               + AbsD(a) * epsTotal + eps * epsTotal;
+                                set result = result * a + coeff;
+                                set epsTotal = epsTotal + eps;
+                            }
+                            if (ctrl == 3) {
+                                EqualityWithinToleranceFact(measured, result, epsTotal);
+                            }
+                            else{
+                                let measuredI = MeasureInteger(LittleEndian(ys));
+                                EqualityFactI(measuredI,
+                                    0,
+                                    $"Controlled FixedPoint polynomial evaluation changed the output register!");
+                            }
+                            ResetAll(xs + ys + ctrls);
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    operation EvaluateOddPolynomialFxPTest() : Unit {
+        for (pos in 4..5) {
+            for (coeffs in [[1.3, -2.4, 1.9],
+                            [-0.3],
+                            [0.1, -0.1, 0.2, 0.2]]){
+                for (a in [0.0, 0.1, -0.1, 0.2, -0.2, 1.3, -1.3]){
+                    let n = 20;
+                    using ((xs, ys) = (Qubit[n], Qubit[n])) {
+                        let fp1 = FixedPoint(pos, xs);
+                        let fp2 = FixedPoint(pos, ys);
+                        PrepareFxP(a, fp1);
+                        EvaluateOddPolynomialFxP(coeffs, fp1, fp2);
+                        let measured = MeasureFxP(fp2);
+                        let eps = 1./IntAsDouble(2^(n-pos));
+                        mutable epsTotal = 0.;
+                        mutable errX = eps;
+                        mutable result = Tail(coeffs);
+                        set epsTotal = epsTotal + eps;
+                        let aSquare = a * a;
+                        for (coeff in coeffs[Length(coeffs)-2..(-1)..0]) {
+                            set epsTotal = epsTotal + AbsD(result) * eps
+                                           + AbsD(aSquare) * epsTotal + eps * epsTotal;
+                            set result = result * aSquare + coeff;
+                            set epsTotal = epsTotal + eps;
+                        }
+                        set epsTotal = epsTotal + AbsD(result) * eps + AbsD(a) * epsTotal
+                                       + eps * epsTotal;
+                        set result = result * a;
+
+                        EqualityWithinToleranceFact(measured, result, epsTotal);
+                        ResetAll(xs + ys);
+                    }
+                }
+            }
+        }
+    }
+
+    operation EvaluateOddPolynomialFxPCtrlTest() : Unit {
+        for (ctrl in 0..3) {
+            for (pos in 4..5) {
+                for (coeffs in [[1.3, -2.4, 1.9],
+                                [-0.3, -0.2],
+                                [0.1, 1.1, -0.1, 0.2],
+                                [0.1, -0.1, 0.2, 0.2, -0.1, 0.3],
+                                [0.2]]){
+                    for (a in [0.0, 0.1, -0.1, 0.2, -0.2, 1.3, -1.3]){
+                        let n = 20;
+                        using ((xs, ys, ctrls) = (Qubit[n], Qubit[n], Qubit[2])) {
+                            let fp1 = FixedPoint(pos, xs);
+                            let fp2 = FixedPoint(pos, ys);
+                            PrepareFxP(a, fp1);
+                            EvaluateOddPolynomialFxP(coeffs, fp1, fp2);
+                            let measured = MeasureFxP(fp2);
+                            let eps = 1./IntAsDouble(2^(n-pos));
+                            mutable epsTotal = 0.;
+                            mutable errX = eps;
+                            mutable result = Tail(coeffs);
+                            set epsTotal = epsTotal + eps;
+                            let aSquare = a * a;
+                            for (coeff in coeffs[Length(coeffs)-2..(-1)..0]) {
+                                set epsTotal = epsTotal + AbsD(result) * eps
+                                               + AbsD(aSquare) * epsTotal + eps * epsTotal;
+                                set result = result * aSquare + coeff;
+                                set epsTotal = epsTotal + eps;
+                            }
+                            set epsTotal = epsTotal + AbsD(result) * eps + AbsD(a) * epsTotal
+                                           + eps * epsTotal;
+                            set result = result * a;
+                            if (ctrl == 3) {
+                                EqualityWithinToleranceFact(measured, result, epsTotal);
+                            }
+                            else{
+                                let measuredI = MeasureInteger(LittleEndian(ys));
+                                EqualityFactI(measuredI,
+                                    0,
+                                    $"Controlled FixedPoint polynomial evaluation changed the output register!");
+                            }
+                            ResetAll(xs + ys + ctrls);
+                        }
+                    }
+                }
+            }
+        }
+    }
+}

Numerics/tests/IntegerHighLevelTests.qs

@@ -0,0 +1,321 @@
+// Copyright (c) Microsoft Corporation. All rights reserved.
+// Licensed under the MIT License.
+
+namespace Microsoft.Quantum.Numerics.ToffoliTests {
+    open Microsoft.Quantum.Diagnostics;
+    open Microsoft.Quantum.Intrinsic;
+    open Microsoft.Quantum.Arithmetic;
+    open Microsoft.Quantum.Extensions.Math;
+
+    operation MultiplyIExhaustiveTest () : Unit {
+        ExhaustiveTestHelper2Args(IntegerMultiplicationRun(false, _, _, _, _));
+    }
+
+    operation SquareIExhaustiveTest () : Unit {
+        ExhaustiveTestHelper1Arg(IntegerSquareRun(false, _, _, _));
+    }
+
+    operation DivideIExhaustiveTest () : Unit {
+        ExhaustiveTestHelper2Args(IntegerDivisionRun);
+    }
+
+    operation SquareSIExhaustiveTest () : Unit {
+        ExhaustiveTestHelper1Arg(IntegerSquareRun(true, _, _, _));
+    }
+
+    operation CompareGTSIExhaustiveTest () : Unit {
+        ExhaustiveTestHelper2Args(IntegerGreaterThanRun(true, _, _, _, _));
+    }
+
+    operation MultiplySIExhaustiveTest () : Unit {
+        ExhaustiveTestHelper2Args(IntegerMultiplicationRun(true, _, _, _, _));
+    }
+
+    operation ComputeReciprocalIExhaustiveTest () : Unit {
+        ExhaustiveTestHelper1Arg(IntegerReciprocalRun(false, _, _, _));
+    }
+
+    operation IntegerGreaterThanRun(signed: Bool, a: Int, b: Int,
+                                    n: Int, numCtrl: Int) : Unit {
+        using ((aqs, bqs, result, ctrlqs) = (Qubit[n], Qubit[n],
+                                             Qubit(), Qubit[numCtrl])) {
+            ApplyXorInPlace(a, LittleEndian(aqs));
+            ApplyXorInPlace(b, LittleEndian(bqs));
+            if (signed) {
+                CompareGTSI(
+                    SignedLittleEndian(LittleEndian(aqs)),
+                    SignedLittleEndian(LittleEndian(bqs)),
+                    result);
+            }
+            else {
+                CompareGTI(LittleEndian(aqs),
+                    LittleEndian(bqs),
+                    result);
+            }
+            mutable asigned = a;
+            mutable bsigned = b;
+            if (signed and a >= 2^(n-1)) {
+                set asigned = -2^n+a;
+            }
+            if (signed and b >= 2^(n-1)) {
+                set bsigned = -2^n+b;
+            }
+            mutable res = asigned > bsigned;
+            mutable resMeasured = M(result);
+            EqualityFactB(res, resMeasured == One,
+                $"Integer comparison failed:
+                  {asigned} > {bsigned} = {res} != {resMeasured} [n={n}]");
+            ResetAll(aqs + bqs + [result]);
+            for (ctrlState in 0..2^numCtrl-1) {
+                ApplyXorInPlace(ctrlState, LittleEndian(ctrlqs));
+                ApplyXorInPlace(a, LittleEndian(aqs));
+                ApplyXorInPlace(b, LittleEndian(bqs));
+                if (signed) {
+                    (Controlled CompareGTSI) (ctrlqs,
+                        (SignedLittleEndian(LittleEndian(aqs)),
+                         SignedLittleEndian(LittleEndian(bqs)),
+                         result));
+                }
+                else {
+                    (Controlled CompareGTI) (ctrlqs,
+                        (LittleEndian(aqs),
+                         LittleEndian(bqs),
+                         result));
+                }
+                set res = asigned > bsigned;
+                if (ctrlState < 2^numCtrl-1) {
+                    set res = false;
+                }
+                set resMeasured = M(result);
+                EqualityFactB(res, resMeasured == One,
+                    $"Controlled integer comparison failed.");
+                ResetAll(aqs + bqs + [result] + ctrlqs);
+            }
+        }
+    }
+
+    operation IntegerMultiplicationRun(signed: Bool, a: Int, b: Int,
+                                             n: Int, numCtrl: Int) : Unit {
+        using ((aqs, bqs, cqs, ctrlqs) = (Qubit[n], Qubit[n],
+                                          Qubit[2*n], Qubit[numCtrl])) {
+            ApplyXorInPlace(a, LittleEndian(aqs));
+            ApplyXorInPlace(b, LittleEndian(bqs));
+            if (signed) {
+                MultiplySI(
+                        SignedLittleEndian(LittleEndian(aqs)),
+                        SignedLittleEndian(LittleEndian(bqs)),
+                        SignedLittleEndian(LittleEndian(cqs)));
+            }
+            else {
+                MultiplyI (LittleEndian(aqs),
+                        LittleEndian(bqs),
+                        LittleEndian(cqs));
+            }
+            mutable asigned = a;
+            mutable bsigned = b;
+            if (signed and a >= 2^(n-1)) {
+                set asigned = -2^n+a;
+            }
+            if (signed and b >= 2^(n-1)) {
+                set bsigned = -2^n+b;
+            }
+            mutable c = asigned * bsigned;
+            mutable cMeasured = MeasureInteger(LittleEndian(cqs));
+            if (signed and cMeasured >= 2^(2*n-1)){
+                set cMeasured = -2^(2*n) + cMeasured;
+            }
+            EqualityFactI(c, cMeasured,
+                $"Multiplication did not yield the correct result:
+                  {asigned} * {bsigned} = {c} != {cMeasured} [n={n}]");
+            ResetAll(aqs + bqs + cqs);
+            for (ctrlState in 0..2^numCtrl-1) {
+                ApplyXorInPlace(ctrlState, LittleEndian(ctrlqs));
+                ApplyXorInPlace(a, LittleEndian(aqs));
+                ApplyXorInPlace(b, LittleEndian(bqs));
+                if (signed) {
+                    (Controlled MultiplySI) (ctrlqs,
+                        (SignedLittleEndian(LittleEndian(aqs)),
+                         SignedLittleEndian(LittleEndian(bqs)),
+                         SignedLittleEndian(LittleEndian(cqs))));
+                }
+                else {
+                    (Controlled MultiplyI) (ctrlqs,
+                        (LittleEndian(aqs),
+                         LittleEndian(bqs),
+                         LittleEndian(cqs)));
+                }
+                set c = asigned * bsigned;
+                if (ctrlState != 2^numCtrl-1) {
+                    set c = 0;
+                }
+                set cMeasured = MeasureInteger(LittleEndian(cqs));
+                if (signed and cMeasured >= 2^(2*n-1)){
+                    set cMeasured = -2^(2*n) + cMeasured;
+                }
+                EqualityFactI(c, cMeasured,
+                    "Controlled multiplication did not yield the correct result.");
+                ResetAll(aqs + bqs + cqs + ctrlqs);
+            }
+        }
+    }
+
+    operation IntegerSquareRun(signed: Bool, a: Int,
+                               n: Int, numCtrl: Int) : Unit {
+        using ((aqs, cqs, ctrlqs) = (Qubit[n], Qubit[2*n],
+                                     Qubit[numCtrl])) {
+            ApplyXorInPlace(a, LittleEndian(aqs));
+            if (signed) {
+                SquareSI(SignedLittleEndian(LittleEndian(aqs)),
+                         SignedLittleEndian(LittleEndian(cqs)));
+            }
+            else {
+                SquareI(LittleEndian(aqs), LittleEndian(cqs));
+            }
+            mutable signeda = a;
+            if (signed and a >= 2^(n-1)) {
+                set signeda = -2^n + a;
+            }
+            mutable c = signeda * signeda;
+            mutable cMeasured = MeasureInteger(LittleEndian(cqs));
+            if (signed and cMeasured >= 2^(2*n-1)){
+                set cMeasured = -2^(2*n) + cMeasured;
+            }
+            EqualityFactI(c, cMeasured,
+                "Square did not yield the correct result.");
+            ResetAll(aqs + cqs);
+            for (ctrlState in 0..2^numCtrl-1) {
+                ApplyXorInPlace(ctrlState, LittleEndian(ctrlqs));
+                ApplyXorInPlace(a, LittleEndian(aqs));
+                if (signed) {
+                    (Controlled SquareSI) (ctrlqs,
+                        (SignedLittleEndian(LittleEndian(aqs)),
+                         SignedLittleEndian(LittleEndian(cqs))));
+                }
+                else {
+                    (Controlled SquareI) (ctrlqs,
+                        (LittleEndian(aqs), LittleEndian(cqs)));
+                }
+                set c = signeda * signeda;
+                if (ctrlState != 2^numCtrl-1) {
+                    set c = 0;
+                }
+                set cMeasured = MeasureInteger(LittleEndian(cqs));
+                if (signed and cMeasured >= 2^(2*n-1)){
+                    set cMeasured = -2^(2*n) + cMeasured;
+                }
+                EqualityFactI(c, cMeasured,
+                    "Controlled square did not yield the correct result.");
+                ResetAll(aqs + cqs + ctrlqs);
+            }
+        }
+    }
+
+    operation IntegerReciprocalRun(signed: Bool, a: Int,
+                                   n: Int, numCtrl: Int) : Unit {
+        using ((aqs, cqs, ctrlqs) = (Qubit[n], Qubit[2*n],
+                                     Qubit[numCtrl])) {
+            ApplyXorInPlace(a, LittleEndian(aqs));
+            ComputeReciprocalI(LittleEndian(aqs), LittleEndian(cqs));
+            mutable c = a > 0 ? 2^(2*n-1) / a | (2^(2*n)-1);
+            mutable cMeasured = MeasureInteger(LittleEndian(cqs));
+            mutable aMeasured = MeasureInteger(LittleEndian(aqs));
+            EqualityFactI(a, aMeasured,
+                "Reciprocal modified the input.");
+            EqualityFactI(c, cMeasured,
+                $"Reciprocal did not yield the correct result:
+                  1/{a} = {c} != {cMeasured}.");
+            ResetAll(aqs + cqs);
+            for (ctrlState in 0..2^numCtrl-1) {
+                ApplyXorInPlace(ctrlState, LittleEndian(ctrlqs));
+                ApplyXorInPlace(a, LittleEndian(aqs));
+                (Controlled ComputeReciprocalI) (ctrlqs,
+                    (LittleEndian(aqs), LittleEndian(cqs)));
+                set c = a > 0 ? 2^(2*n-1) / a | (2^(2*n)-1);
+                if (ctrlState != 2^numCtrl-1) {
+                    set c = 0;
+                }
+                set cMeasured = MeasureInteger(LittleEndian(cqs));
+                EqualityFactI(c, cMeasured,
+                    $"Controlled reciprocal did not yield the correct result:
+                     1/{a} = {c} != {cMeasured}.");
+                set aMeasured = MeasureInteger(LittleEndian(aqs));
+                EqualityFactI(a, aMeasured,
+                    "Reciprocal modified the input.");
+                ResetAll(aqs + cqs + ctrlqs);
+            }
+        }
+    }
+
+    operation IntegerDivisionRun(a: Int, b: Int, n: Int, numCtrl: Int): Unit {
+        using ((aqs, bqs, cqs, ctrlqs) = (Qubit[n], Qubit[n],
+                                          Qubit[n], Qubit[numCtrl])) {
+            mutable c = 0;
+            if (b > 0) {
+                set c = a / b;
+            }
+            else{
+                set c = 2^n - 1;
+            }
+            let rem = a - c * b;
+            ApplyXorInPlace(a, LittleEndian(aqs));
+            ApplyXorInPlace(b, LittleEndian(bqs));
+            DivideI (LittleEndian(aqs), LittleEndian(bqs),
+                             LittleEndian(cqs));
+            mutable resMeasured = MeasureInteger(LittleEndian(cqs));
+            mutable remMeasured = MeasureInteger(LittleEndian(aqs));
+            EqualityFactI(c, resMeasured,
+                $"Controlled division did not yield the correct result:
+                    {a} / {b} = {c} != {resMeasured}");
+            EqualityFactI(remMeasured, rem,
+                "Controlled division did not yield the correct remainder:
+                    rem = {rem} != {remMeasured} for {a} / {b}");
+            ResetAll(aqs + bqs + cqs + ctrlqs);
+            for (ctrlState in 0..2^numCtrl-1) {
+                ApplyXorInPlace(ctrlState, LittleEndian(ctrlqs));
+                ApplyXorInPlace(a, LittleEndian(aqs));
+                ApplyXorInPlace(b, LittleEndian(bqs));
+                (Controlled DivideI) (ctrlqs,
+                    (LittleEndian(aqs), LittleEndian(bqs), LittleEndian(cqs)));
+                set resMeasured = MeasureInteger(LittleEndian(cqs));
+                set remMeasured = MeasureInteger(LittleEndian(aqs));
+                if (ctrlState == 2^numCtrl-1) {
+                    EqualityFactI(c, resMeasured,
+                        $"Controlled division did not yield the correct result:
+                            {a} / {b} = {c} != {resMeasured}");
+                    EqualityFactI(remMeasured, rem,
+                        "Controlled division did not yield the correct remainder:
+                            rem = {rem} != {remMeasured} for {a} / {b}");
+                }
+                else {
+                    EqualityFactI(0, resMeasured,
+                        "Controlled division was not trivial.");
+                    EqualityFactI(a, remMeasured,
+                        "Controlled division was not trivial.");
+                }
+                ResetAll(aqs + bqs + cqs + ctrlqs);
+            }
+        }
+    }
+
+    operation ExhaustiveTestHelper1Arg (TestFunction: ((Int, Int, Int) => Unit)) : Unit {
+        for (numCtrlQubits in 0..2) {
+            for (numQubits in 1..5) {
+                for (a in 0..2^numQubits-1) {
+                    TestFunction(a, numQubits, numCtrlQubits);
+                }
+            }
+        }
+    }
+
+    operation ExhaustiveTestHelper2Args (TestFunction: ((Int, Int, Int, Int) => Unit)) : Unit {
+        for (numCtrlQubits in 0..2) {
+            for (numQubits in 1..5) {
+                for (a in 0..2^numQubits-1) {
+                    for (b in 0..2^numQubits-1) {
+                        TestFunction(a, b, numQubits, numCtrlQubits);
+                    }
+                }
+            }
+        }
+    }
+}

Numerics/tests/NumericsTests.cs

@@ -0,0 +1,52 @@
+// Copyright (c) Microsoft Corporation. All rights reserved.
+// Licensed under the MIT License.
+
+using Microsoft.Quantum.Simulation.XUnit;
+using Microsoft.Quantum.Simulation.Simulators;
+using Xunit.Abstractions;
+using System.Diagnostics;
+
+namespace Microsoft.Quantum.Numerics.Tests
+{
+    public class NumericsTests
+    {
+        private readonly ITestOutputHelper output;
+
+        public NumericsTests(ITestOutputHelper output)
+        {
+            this.output = output;
+        }
+
+        /// <summary>
+        /// This driver will run all Q# tests (operations named "...Test")
+        /// that are located inside Microsoft.Quantum.Numerics.Tests using the quantum
+        /// simulator.
+        /// </summary>
+        [OperationDriver(TestNamespace = "Microsoft.Quantum.Numerics.Tests",
+                         TestCasePrefix = "QSim:")]
+        public void QSimTests(TestOperation op)
+        {
+            var sim = new QuantumSimulator();
+            // OnLog defines action(s) performed when Q# test calls function Message
+            sim.OnLog += (msg) => { output.WriteLine(msg); };
+            sim.OnLog += (msg) => { Debug.WriteLine(msg); };
+            op.TestOperationRunner(sim);
+        }
+
+        /// <summary>
+        /// This driver will run all Q# tests (operations named "...Test")
+        /// that are located inside Microsoft.Quantum.Numerics.ToffoliTests using the
+        /// Toffoli simulator.
+        /// </summary>
+        [OperationDriver(TestNamespace = "Microsoft.Quantum.Numerics.ToffoliTests",
+                         TestCasePrefix = "ToffSim:")]
+        public void ToffoliSimTests(TestOperation op)
+        {
+            var sim = new ToffoliSimulator();
+            // OnLog defines action(s) performed when Q# test calls function Message
+            sim.OnLog += (msg) => { output.WriteLine(msg); };
+            sim.OnLog += (msg) => { Debug.WriteLine(msg); };
+            op.TestOperationRunner(sim);
+        }
+    }
+}

Numerics/tests/NumericsTests.csproj

@@ -0,0 +1,26 @@
+<Project Sdk="Microsoft.NET.Sdk">
+  <PropertyGroup>
+    <TargetFramework>netcoreapp2.0</TargetFramework>
+    <PlatformTarget>x64</PlatformTarget>
+    <IsPackable>false</IsPackable>
+    <LangVersion>latest</LangVersion>
+  </PropertyGroup>
+
+  <PropertyGroup>
+    <NoWarn>0162</NoWarn>
+  </PropertyGroup>
+
+  <ItemGroup>
+    <ProjectReference Include="..\src\Numerics.csproj" />
+    <ProjectReference Include="..\..\Standard\src\Standard.csproj" />
+  </ItemGroup>
+
+  <ItemGroup>
+    <PackageReference Include="Microsoft.Quantum.Development.Kit" Version="0.6.1905.301" />
+    <PackageReference Include="Microsoft.Quantum.Xunit" Version="0.6.1905.301" />
+    <PackageReference Include="Microsoft.NET.Test.Sdk" Version="15.8.0" />
+    <PackageReference Include="xunit" Version="2.4.0" />
+    <PackageReference Include="xunit.runner.visualstudio" Version="2.4.0" />
+    <DotNetCliToolReference Include="dotnet-xunit" Version="2.3.1" />
+  </ItemGroup>
+</Project>