Merge branch 'main' into comparison

qualtran/bloqs/gf_arithmetic/gf2_inverse.ipynb

@@ -58,15 +58,22 @@
     "    a^{-1} = a^{2^m - 2}\n",
     "$$\n",
     "\n",
-    "Thus, the inverse can be obtained via $m - 1$ squaring and multiplication operations.\n",
+    "The exponential $a^{2^m - 2}$ is computed using $\\mathcal{O}(m)$ squaring and\n",
+    "$\\mathcal{O}(\\log_2(m))$ multiplications via Itoh-Tsujii inversion. The algorithm is described on\n",
+    "page 4 and 5 of Ref[1] and resembles binary exponentiation. The inverse is computed as $B_{n-1}^2$,\n",
+    "where $B_1 = x$ and $B_{i+j} = B_i B_j^{2^i}$.\n",
     "\n",
     "#### Parameters\n",
     " - `bitsize`: The degree $m$ of the galois field $GF(2^m)$. Also corresponds to the number of qubits in the input register whose inverse should be calculated. \n",
     "\n",
     "#### Registers\n",
     " - `x`: Input THRU register of size $m$ that stores elements from $GF(2^m)$.\n",
     " - `result`: Output RIGHT register of size $m$ that stores $x^{-1}$ from $GF(2^m)$.\n",
-    " - `junk`: Output RIGHT register of size $m$ and shape ($m - 2$) that stores results from intermediate multiplications.\n"
+    " - `junk`: Output RIGHT register of size $m$ and shape ($m - 2$) that stores results from intermediate multiplications. \n",
+    "\n",
+    "#### References\n",
+    " - [Efficient quantum circuits for binary elliptic curve arithmetic: reducing T -gate complexity](https://arxiv.org/abs/1209.6348). Section 2.3\n",
+    " - [Structure of parallel multipliers for a class of fields GF(2^m)](https://doi.org/10.1016/0890-5401(89)90045-X)\n"
    ]
   },
   {

qualtran/bloqs/gf_arithmetic/gf2_inverse.py

@@ -30,10 +30,11 @@
 from qualtran.bloqs.gf_arithmetic.gf2_addition import GF2Addition
 from qualtran.bloqs.gf_arithmetic.gf2_multiplication import GF2Multiplication
 from qualtran.bloqs.gf_arithmetic.gf2_square import GF2Square
-from qualtran.symbolics import is_symbolic, SymbolicInt
+from qualtran.resource_counting.generalizers import ignore_alloc_free, ignore_split_join
+from qualtran.symbolics import bit_length, ceil, is_symbolic, log2, SymbolicInt
 
 if TYPE_CHECKING:
-    from qualtran import BloqBuilder, Soquet
+    from qualtran import BloqBuilder, Soquet, SoquetT
     from qualtran.resource_counting import BloqCountDictT, BloqCountT, CostKey, SympySymbolAllocator
     from qualtran.simulation.classical_sim import ClassicalValT
 
@@ -62,7 +63,10 @@ class GF2Inverse(Bloq):
         a^{-1} = a^{2^m - 2}
     $$
 
-    Thus, the inverse can be obtained via $m - 1$ squaring and multiplication operations.
+    The exponential $a^{2^m - 2}$ is computed using $\mathcal{O}(m)$ squaring and
+    $\mathcal{O}(\log_2(m))$ multiplications via Itoh-Tsujii inversion. The algorithm is described on
+    page 4 and 5 of Ref[1] and resembles binary exponentiation. The inverse is computed as $B_{n-1}^2$,
+    where $B_1 = x$ and $B_{i+j} = B_i B_j^{2^i}$.
 
     Args:
         bitsize: The degree $m$ of the galois field $GF(2^m)$. Also corresponds to the number of
@@ -73,15 +77,21 @@ class GF2Inverse(Bloq):
         result: Output RIGHT register of size $m$ that stores $x^{-1}$ from $GF(2^m)$.
         junk: Output RIGHT register of size $m$ and shape ($m - 2$) that stores
             results from intermediate multiplications.
+
+    References:
+        [Efficient quantum circuits for binary elliptic curve arithmetic: reducing T -gate complexity](https://arxiv.org/abs/1209.6348).
+        Section 2.3
+
+        [Structure of parallel multipliers for a class of fields GF(2^m)](https://doi.org/10.1016/0890-5401(89)90045-X)
     """
 
     bitsize: SymbolicInt
 
     @cached_property
     def signature(self) -> 'Signature':
         junk_reg = (
-            [Register('junk', dtype=self.qgf, shape=(self.bitsize - 2,), side=Side.RIGHT)]
-            if is_symbolic(self.bitsize) or self.bitsize > 2
+            [Register('junk', dtype=self.qgf, shape=(self.n_junk_regs,), side=Side.RIGHT)]
+            if is_symbolic(self.bitsize) or self.bitsize > 1
             else []
         )
         return Signature(
@@ -96,60 +106,111 @@ def signature(self) -> 'Signature':
     def qgf(self) -> QGF:
         return QGF(characteristic=2, degree=self.bitsize)
 
+    @cached_property
+    def n_junk_regs(self) -> SymbolicInt:
+        return 2 * bit_length(self.bitsize - 1) + self.bitsize_hamming_weight
+
+    @cached_property
+    def bitsize_hamming_weight(self) -> SymbolicInt:
+        """Hamming weight of self.bitsize - 1"""
+        return (
+            bit_length(self.bitsize - 1)
+            if is_symbolic(self.bitsize)
+            else int(self.bitsize - 1).bit_count()
+        )
+
     def my_static_costs(self, cost_key: 'CostKey'):
+        from qualtran._infra.gate_with_registers import total_bits
         from qualtran.resource_counting import QubitCount
 
         if isinstance(cost_key, QubitCount):
-            return self.signature.n_qubits()
+            return total_bits(self.signature.rights())
 
         return NotImplemented
 
-    def build_composite_bloq(self, bb: 'BloqBuilder', *, x: 'Soquet') -> Dict[str, 'Soquet']:
+    def build_composite_bloq(self, bb: 'BloqBuilder', *, x: 'Soquet') -> Dict[str, 'SoquetT']:
         if is_symbolic(self.bitsize):
             raise DecomposeTypeError(f"Cannot decompose symbolic {self}")
+
         result = bb.allocate(dtype=self.qgf)
         if self.bitsize == 1:
             x, result = bb.add(GF2Addition(self.bitsize), x=x, y=result)
             return {'x': x, 'result': result}
 
-        x = bb.add(GF2Square(self.bitsize), x=x)
-        x, result = bb.add(GF2Addition(self.bitsize), x=x, y=result)
-
         junk = []
-        for i in range(2, self.bitsize):
-            x = bb.add(GF2Square(self.bitsize), x=x)
-            x, result, new_result = bb.add(GF2Multiplication(self.bitsize), x=x, y=result)
-            junk.append(result)
-            result = new_result
-        x = bb.add(GF2Square(self.bitsize), x=x)
-        return {'x': x, 'result': result} | ({'junk': np.array(junk)} if junk else {})
+        beta = bb.allocate(dtype=self.qgf)
+        x, beta = bb.add(GF2Addition(self.bitsize), x=x, y=beta)
+        is_first = True
+        bitsize_minus_one = int(self.bitsize - 1)
+        for i in range(bitsize_minus_one.bit_length()):
+            if (1 << i) & bitsize_minus_one:
+                if is_first:
+                    beta, result = bb.add(GF2Addition(self.bitsize), x=beta, y=result)
+                    is_first = False
+                else:
+                    for j in range(2**i):
+                        result = bb.add(GF2Square(self.bitsize), x=result)
+                    beta, result, new_result = bb.add(
+                        GF2Multiplication(self.bitsize), x=beta, y=result
+                    )
+                    junk.append(result)
+                    result = new_result
+            beta_squared = bb.allocate(dtype=self.qgf)
+            beta, beta_squared = bb.add(GF2Addition(self.bitsize), x=beta, y=beta_squared)
+            for j in range(2**i):
+                beta_squared = bb.add(GF2Square(self.bitsize), x=beta_squared)
+            beta, beta_squared, beta_new = bb.add(
+                GF2Multiplication(self.bitsize), x=beta, y=beta_squared
+            )
+            junk.extend([beta, beta_squared])
+            beta = beta_new
+        junk.append(beta)
+        result = bb.add(GF2Square(self.bitsize), x=result)
+        assert len(junk) == self.n_junk_regs, f'{len(junk)=}, {self.n_junk_regs=}'
+        return {'x': x, 'result': result, 'junk': np.array(junk)}
 
     def build_call_graph(
         self, ssa: 'SympySymbolAllocator'
     ) -> Union['BloqCountDictT', Set['BloqCountT']]:
-        if is_symbolic(self.bitsize) or self.bitsize > 2:
-            return {
-                GF2Addition(self.bitsize): 1,
-                GF2Square(self.bitsize): self.bitsize - 1,
-                GF2Multiplication(self.bitsize): self.bitsize - 2,
-            }
-        return {GF2Addition(self.bitsize): 1} | (
-            {GF2Square(self.bitsize): 1} if self.bitsize == 2 else {}
-        )
+        if not is_symbolic(self.bitsize) and self.bitsize == 1:
+            return {GF2Addition(self.bitsize): 1}
+        square_count = self.bitsize + 2 ** ceil(log2(self.bitsize)) - 1
+        if not is_symbolic(self.bitsize):
+            n = self.bitsize - 1
+            square_count -= n & (-n)
+        return {
+            GF2Addition(self.bitsize): 2 + ceil(log2(self.bitsize)),
+            GF2Square(self.bitsize): square_count,
+            GF2Multiplication(self.bitsize): ceil(log2(self.bitsize))
+            + self.bitsize_hamming_weight
+            - 1,
+        }
 
     def on_classical_vals(self, *, x) -> Dict[str, 'ClassicalValT']:
         assert isinstance(x, self.qgf.gf_type)
-        x_temp = x**2
-        result = x_temp
         junk = []
-        for i in range(2, int(self.bitsize)):
-            junk.append(result)
-            x_temp = x_temp * x_temp
-            result = result * x_temp
+        bitsize_minus_one = int(self.bitsize - 1)
+        beta = x
+        result = self.qgf.gf_type(0)
+        is_first = True
+        for i in range(bitsize_minus_one.bit_length()):
+            if (1 << i) & bitsize_minus_one:
+                if is_first:
+                    is_first = False
+                    result = beta
+                else:
+                    for j in range(2**i):
+                        result = result**2
+                    junk.append(result)
+                    result = result * beta
+            beta_squared = beta ** (2 ** (2**i))
+            junk.extend([beta, beta_squared])
+            beta = beta * beta_squared
+        junk.append(beta)
         return {'x': x, 'result': x ** (-1), 'junk': np.array(junk)}
 
 
-@bloq_example
+@bloq_example(generalizer=[ignore_split_join, ignore_alloc_free])
 def _gf16_inverse() -> GF2Inverse:
     gf16_inverse = GF2Inverse(4)
     return gf16_inverse

qualtran/bloqs/gf_arithmetic/gf2_inverse_test.py

@@ -22,7 +22,9 @@
     GF2Inverse,
 )
 from qualtran.resource_counting import get_cost_value, QECGatesCost, QubitCount
-from qualtran.testing import assert_consistent_classical_action
+from qualtran.resource_counting.generalizers import ignore_alloc_free, ignore_split_join
+from qualtran.symbolics import ceil, log2
+from qualtran.testing import assert_consistent_classical_action, assert_equivalent_bloq_counts
 
 
 def test_gf16_inverse(bloq_autotester):
@@ -36,8 +38,10 @@ def test_gf2_inverse_symbolic(bloq_autotester):
 def test_gf2_inverse_symbolic_toffoli_complexity():
     bloq = _gf2_inverse_symbolic.make()
     m = bloq.bitsize
-    assert get_cost_value(bloq, QECGatesCost()).total_toffoli_only() - m**2 * (m - 2) == 0
-    assert sympy.simplify(get_cost_value(bloq, QubitCount()) - m**2) == 0
+    expected_expr = m**2 * (2 * ceil(log2(m)) - 1)
+    assert get_cost_value(bloq, QECGatesCost()).total_toffoli_only() - expected_expr == 0
+    expected_expr = m * (3 * ceil(log2(m)) + 2)
+    assert sympy.simplify(get_cost_value(bloq, QubitCount()) - expected_expr) == 0
 
 
 def test_gf2_inverse_classical_sim_quick():
@@ -53,3 +57,9 @@ def test_gf2_inverse_classical_sim(m):
     bloq = GF2Inverse(m)
     GFM = GF(2**m)
     assert_consistent_classical_action(bloq, x=GFM.elements[1:])
+
+
+@pytest.mark.parametrize('m', [*range(1, 12)])
+def test_gf2_equivalent_bloq_counts(m):
+    bloq = GF2Inverse(m)
+    assert_equivalent_bloq_counts(bloq, generalizer=[ignore_split_join, ignore_alloc_free])

qualtran/bloqs/gf_arithmetic/gf2_multiplication.ipynb

@@ -56,12 +56,13 @@
     "gates.\n",
     "\n",
     "#### Parameters\n",
-    " - `bitsize`: The degree $m$ of the galois field $GF(2^m)$. Also corresponds to the number of qubits in each of the two input registers $a$ and $b$ that should be multiplied. \n",
+    " - `bitsize`: The degree $m$ of the galois field $GF(2^m)$. Also corresponds to the number of qubits in each of the two input registers $a$ and $b$ that should be multiplied.\n",
+    " - `plus_equal_prod`: If True, implements the `PlusEqualProduct` version that applies the map $|x\\rangle |y\\rangle |z\\rangle \\rightarrow |x\\rangle |y\\rangle |x + z\\rangle$. \n",
     "\n",
     "#### Registers\n",
     " - `x`: Input THRU register of size $m$ that stores elements from $GF(2^m)$.\n",
     " - `y`: Input THRU register of size $m$ that stores elements from $GF(2^m)$.\n",
-    " - `result`: Output RIGHT register of size $m$ that stores the product $x * y$ in $GF(2^m)$.  \n",
+    " - `result`: Register of size $m$ that stores the product $x * y$ in $GF(2^m)$. If plus_equal_prod is True - result is a THRU register and stores $result + x * y$. If plus_equal_prod is False - result is a RIGHT register and stores $x * y$.  \n",
     "\n",
     "#### References\n",
     " - [On the Design and Optimization of a Quantum Polynomial-Time Attack on Elliptic Curve Cryptography](https://arxiv.org/abs/0710.1093). \n",
@@ -99,7 +100,7 @@
    },
    "outputs": [],
    "source": [
-    "gf16_multiplication = GF2Multiplication(4)"
+    "gf16_multiplication = GF2Multiplication(4, plus_equal_prod=True)"
    ]
   },
   {
@@ -114,7 +115,7 @@
     "import sympy\n",
     "\n",
     "m = sympy.Symbol('m')\n",
-    "gf2_multiplication_symbolic = GF2Multiplication(m)"
+    "gf2_multiplication_symbolic = GF2Multiplication(m, plus_equal_prod=False)"
    ]
   },
   {