Implement square root.

qsharp.json

@@ -16,6 +16,7 @@
       "src/QuantumArithmetic/LYY2021.qs",
       "src/QuantumArithmetic/LYY2021Test.qs",
       "src/QuantumArithmetic/MCT2017.qs",
+      "src/QuantumArithmetic/MCT2017_SquareRoot.qs",
       "src/QuantumArithmetic/NZLS2023.qs",
       "src/QuantumArithmetic/NZLS2023Test.qs",
       "src/QuantumArithmetic/OFOSG2023.qs",

src/QuantumArithmetic/MCT2017_SquareRoot.qs

@@ -0,0 +1,108 @@
+import QuantumArithmetic.GKDKH2021.QuantumFullAdder;
+import QuantumArithmetic.Utils.SWAPViaRelabel;
+/// Square Root algorithm, presented in the paper:
+///   T-count and Qubit Optimized Quantum Circuit Design of the Non-Restoring Square Root Algorithm
+///   Edgard Muñoz-Coreas, Himanshu Thapliyal, 2017.
+///   https://arxiv.org/abs/1712.08254
+/// All numbers are unsigned integers, little-endian.
+import Std.Diagnostics.Fact;
+import QuantumArithmetic.Utils.DivCeil;
+
+// Computes ys-=xs if ctrl=1, and ys+=xs if ctrl=0.
+operation AddSub(ctrl : Qubit, xs : Qubit[], ys : Qubit[]) : Unit is Adj + Ctl {
+    let config = new QuantumArithmetic.TMVH2019.Config { Adder = Std.Arithmetic.RippleCarryCGIncByLE };
+    QuantumArithmetic.TMVH2019.AddSub(ctrl, xs, ys, config);
+}
+
+// Computes ys+=xs if ctrl=1, does nothing if ctrl=0.
+operation CtrlAdd(ctrl : Qubit, xs : Qubit[], ys : Qubit[]) : Unit is Adj + Ctl {
+    Controlled Std.Arithmetic.RippleCarryCGIncByLE([ctrl], (xs, ys));
+}
+
+/// Computes R;Ans = R-Sqrt(R)^2;Sqrt(R).
+/// R and Ans must be of the same size.
+/// This is the implementation from the paper, but it is incorrect when the
+/// highest bit of R is 1.
+operation SquareRootInternal(R : Qubit[], Ans : Qubit[]) : Unit is Adj + Ctl {
+    let n = Length(R);
+    Fact(n % 2 == 0, "n must be even");
+    Fact(n >= 4, "n is too small");
+    Fact(Length(Ans) == n, "Size mismatch");
+    let m = n / 2;
+    use z = Qubit();
+    let F = Ans[n-2..n-1] + Ans[0..n-3];
+
+    X(F[0]); // Set F:=1.
+
+    // Part 1: Initial Substraction.
+    X(R[n-2]); // Step 1.
+    CNOT(R[n-2], R[n-1]);
+    CNOT(R[n-1], F[1]);
+    X(R[n-1]);
+    CNOT(R[n-1], z);
+    CNOT(R[n-1], F[2]);
+    X(R[n-1]);
+    AddSub(z, F[0..3], R[n-4..n-1]);
+
+    // Part 2: Conditional Addition or Subtraction.
+    for i in 2..m-1 {
+        X(z);
+        CNOT(z, F[1]);
+        X(z);
+        CNOT(F[2], z);
+        CNOT(R[n-1], F[1]);
+        X(R[n-1]);
+        CNOT(R[n-1], z);
+        CNOT(R[n-1], F[i + 1]);
+        X(R[n-1]);
+        for j in i + 1..-1..3 {
+            SWAP(F[j], F[j-1]);
+        }
+        AddSub(z, F[0..2 * i + 1], R[n-2 * i-2..n-1]);
+    }
+
+    // Part 3: Remainder Restoration.
+    X(z);
+    CNOT(z, F[1]);
+    X(z);
+    CNOT(F[2], z);
+    X(R[n-1]);
+    CNOT(R[n-1], z);
+    CNOT(R[n-1], F[m + 1]);
+    X(R[n-1]);
+    X(z);
+    CtrlAdd(z, F, R);
+    X(z);
+    for j in m + 1..-1..3 {
+        SWAP(F[j], F[j-1]);
+    }
+    CNOT(F[2], z);
+
+
+    X(F[0]);
+}
+
+/// Computes R;Ans = R-Sqrt(R)^2;Sqrt(R).
+/// R can be of any size.
+/// Must be Length(Ans)>=⌈Length(R)/2⌉.
+/// Ans must be prepared in zero state.
+operation SquareRoot(R : Qubit[], Ans : Qubit[]) : Unit {
+    let n = Length(R);
+    Fact(Length(Ans) >= DivCeil(n, 2), "Ans is to small.");
+    if (n == 1) {
+        SWAP(R[0], Ans[0]);
+    } else {
+        // Add minimal necessary padding, so R has even length and its highest
+        // bit is zero.
+        let pad_R_size = 2 -(n % 2);
+        use pad_R = Qubit[pad_R_size];
+        if (Length(Ans) > n + pad_R_size) {
+            SquareRootInternal(R + pad_R, Ans[0..n + pad_R_size-1]);
+        } else {
+            use pad_Ans = Qubit[n + pad_R_size-Length(Ans)];
+            SquareRootInternal(R + pad_R, Ans + pad_Ans);
+        }
+    }
+}
+
+export SquareRoot;

src/TestUtils.qs

@@ -65,6 +65,25 @@ operation BinaryOpInPlace(n : Int, x_val : BigInt, y_val : BigInt, op : (Qubit[]
     return ans;
 }
 
+// Tests arbitrary operation acting on 2 registers.
+// Applies `op` on two registers of sizes nx, ny, populated with x,y.
+// Returns new values of the registers.
+operation BinaryOpArb(
+    nx : Int,
+    ny : Int,
+    x : BigInt,
+    y : BigInt,
+    op : (Qubit[], Qubit[]) => Unit
+) : (BigInt, BigInt) {
+    use X = Qubit[nx];
+    use Y = Qubit[ny];
+    ApplyBigInt(x, X);
+    ApplyBigInt(y, Y);
+    op(X, Y);
+    return (MeasureBigInt(X), MeasureBigInt(Y));
+}
+
+
 // Tests arbitrary operation acting on 3 registers.
 // Applies `op` on three registers of sizes nx, ny, nz, populated with x,y,z.
 // Returns new values of the registers.

test/MCT2017_SquareRoot_test.py

@@ -0,0 +1,32 @@
+import pytest
+from qsharp import init, eval
+import random
+import math
+
+
+@pytest.fixture(scope="session", autouse=True)
+def setup():
+    init(project_root='.')
+
+
+@pytest.mark.parametrize("n", [1, 2, 3, 4, 5, 6, 7, 8])
+def test_SquareRoot_Exhaustive(n: int):
+    op = "QuantumArithmetic.MCT2017_SquareRoot.SquareRoot"
+    for x in range(2**n):
+        ans = eval(f"TestUtils.BinaryOpArb({n},{n},{x}L,0L,{op})")
+        true_root = math.floor(math.sqrt(x))
+        expected = (x - true_root**2), true_root
+        assert ans == expected
+
+
+@pytest.mark.parametrize("n1,n2", [
+    (2, 1), (2, 3), (5, 3), (5, 10), (10, 5), (10, 6), (10, 11), (16, 8),
+    (32, 16), (32, 32), (64, 32), (100, 50)])
+def test_SquareRoot(n1: int, n2: int):
+    op = "QuantumArithmetic.MCT2017_SquareRoot.SquareRoot"
+    for _ in range(5):
+        x = random.randint(0, 2**n1-1)
+        ans = eval(f"TestUtils.BinaryOpArb({n1},{n2},{x}L,0L,{op})")
+        true_root = math.isqrt(x)
+        expected = (x - true_root**2), true_root
+        assert ans == expected