Add ECAdd() Bloq (#1425)

* Initial commit of ec add waiting on equals to be merged.

* Working on tests for ECAdd

* ECAdd implementation and tests

* remove modmul typo

* Fix mypy errors

* Better bugfix for ModAdd

* Change mod inv classical impl to use monttgomery inv

* Fix pytest error

* Fix some comments

* ECAdd lots of testing

* Add comments about bugs to be fixed

* Reduce complexity by keeping intermediate values mod p

* Stash qmontgomery tests

* Address comments

* Fix montgomery prod/inv calculations + pylint/mypy

---------

Co-authored-by: Noureldin <noureldinyosri@google.com>
Co-authored-by: Matthew Harrigan <mpharrigan@google.com>

Author: 3 people

qualtran/_infra/data_types.py

@@ -772,9 +772,14 @@ class QMontgomeryUInt(QDType):
         bitsize: The number of qubits used to represent the integer.
 
     References:
-        [Montgomery modular multiplication](https://en.wikipedia.org/wiki/Montgomery_modular_multiplication)
+        [Montgomery modular multiplication](https://en.wikipedia.org/wiki/Montgomery_modular_multiplication).
+
+        [Performance Analysis of a Repetition Cat Code Architecture: Computing 256-bit Elliptic Curve Logarithm in 9 Hours with 126133 Cat Qubits](https://arxiv.org/abs/2302.06639).
+        Gouzien et al. 2023.
+        We follow Montgomery form as described in the above paper; namely, r = 2^bitsize.
     """
 
+    # TODO(https://github.com/quantumlib/Qualtran/issues/1471): Add modulus p as a class member.
     bitsize: SymbolicInt
 
     @property
@@ -810,6 +815,43 @@ def assert_valid_classical_val_array(
         if np.any(val_array >= 2**self.bitsize):
             raise ValueError(f"Too-large classical values encountered in {debug_str}")
 
+    def montgomery_inverse(self, xm: int, p: int) -> int:
+        """Returns the modular inverse of an integer in montgomery form.
+
+        Args:
+            xm: An integer in montgomery form.
+            p: The modulus of the finite field.
+        """
+        return ((pow(xm, -1, p)) * pow(2, 2 * self.bitsize, p)) % p
+
+    def montgomery_product(self, xm: int, ym: int, p: int) -> int:
+        """Returns the modular product of two integers in montgomery form.
+
+        Args:
+            xm: The first montgomery form integer for the product.
+            ym: The second montgomery form integer for the product.
+            p: The modulus of the finite field.
+        """
+        return (xm * ym * pow(2, -self.bitsize, p)) % p
+
+    def montgomery_to_uint(self, xm: int, p: int) -> int:
+        """Converts an integer in montgomery form to a normal form integer.
+
+        Args:
+            xm: An integer in montgomery form.
+            p: The modulus of the finite field.
+        """
+        return (xm * pow(2, -self.bitsize, p)) % p
+
+    def uint_to_montgomery(self, x: int, p: int) -> int:
+        """Converts an integer into montgomery form.
+
+        Args:
+            x: An integer.
+            p: The modulus of the finite field.
+        """
+        return (x * pow(2, int(self.bitsize), p)) % p
+
 
 @attrs.frozen
 class QGF(QDType):

qualtran/_infra/data_types_test.py

@@ -135,6 +135,32 @@ def test_qmontgomeryuint():
     assert is_symbolic(QMontgomeryUInt(sympy.Symbol('x')))
 
 
+@pytest.mark.parametrize('p', [13, 17, 29])
+@pytest.mark.parametrize('val', [1, 5, 7, 9])
+def test_qmontgomeryuint_operations(val, p):
+    qmontgomeryuint_8 = QMontgomeryUInt(8)
+    # Convert value to montgomery form and get the modular inverse.
+    val_m = qmontgomeryuint_8.uint_to_montgomery(val, p)
+    mod_inv = qmontgomeryuint_8.montgomery_inverse(val_m, p)
+
+    # Calculate the product in montgomery form and convert back to normal form for assertion.
+    assert (
+        qmontgomeryuint_8.montgomery_to_uint(
+            qmontgomeryuint_8.montgomery_product(val_m, mod_inv, p), p
+        )
+        == 1
+    )
+
+
+@pytest.mark.parametrize('p', [13, 17, 29])
+@pytest.mark.parametrize('val', [1, 5, 7, 9])
+def test_qmontgomeryuint_conversions(val, p):
+    qmontgomeryuint_8 = QMontgomeryUInt(8)
+    assert val == qmontgomeryuint_8.montgomery_to_uint(
+        qmontgomeryuint_8.uint_to_montgomery(val, p), p
+    )
+
+
 def test_qgf():
     qgf_256 = QGF(characteristic=2, degree=8)
     assert str(qgf_256) == 'QGF(2**8)'

qualtran/bloqs/arithmetic/_shims.py

@@ -22,8 +22,11 @@
 
 from attrs import frozen
 
-from qualtran import Bloq, QBit, QUInt, Register, Signature
+from qualtran import Bloq, QBit, QMontgomeryUInt, QUInt, Register, Signature
+from qualtran.bloqs.arithmetic.bitwise import BitwiseNot
+from qualtran.bloqs.arithmetic.controlled_addition import CAdd
 from qualtran.bloqs.basic_gates import Toffoli
+from qualtran.bloqs.basic_gates.swap import TwoBitCSwap
 from qualtran.resource_counting import BloqCountDictT, SympySymbolAllocator
 
 
@@ -39,6 +42,20 @@ def build_call_graph(self, ssa: SympySymbolAllocator) -> BloqCountDictT:
         return {Toffoli(): self.n - 2}
 
 
+@frozen
+class CSub(Bloq):
+    n: int
+
+    @cached_property
+    def signature(self) -> 'Signature':
+        return Signature(
+            [Register('ctrl', QBit()), Register('x', QUInt(self.n)), Register('y', QUInt(self.n))]
+        )
+
+    def build_call_graph(self, ssa: SympySymbolAllocator) -> BloqCountDictT:
+        return {CAdd(QMontgomeryUInt(self.n)): 1, BitwiseNot(QMontgomeryUInt(self.n)): 3}
+
+
 @frozen
 class Lt(Bloq):
     n: int
@@ -62,3 +79,6 @@ class CHalf(Bloq):
     @cached_property
     def signature(self) -> 'Signature':
         return Signature([Register('ctrl', QBit()), Register('x', QUInt(self.n))])
+
+    def build_call_graph(self, ssa: SympySymbolAllocator) -> BloqCountDictT:
+        return {TwoBitCSwap(): self.n}

qualtran/bloqs/factoring/ecc/ec_add.ipynb

@@ -41,16 +41,22 @@
     "This takes elliptic curve points given by (a, b) and (x, y)\n",
     "and outputs the sum (x_r, y_r) in the second pair of registers.\n",
     "\n",
+    "Because the decomposition of this Bloq is complex, we split it into six separate parts\n",
+    "corresponding to the parts described in figure 10 of the Litinski paper cited below. We follow\n",
+    "the signature from figure 5 and break down the further decompositions based on the steps in\n",
+    "figure 10.\n",
+    "\n",
     "#### Parameters\n",
     " - `n`: The bitsize of the two registers storing the elliptic curve point\n",
-    " - `mod`: The modulus of the field in which we do the addition. \n",
+    " - `mod`: The modulus of the field in which we do the addition.\n",
+    " - `window_size`: The number of bits in the ModMult window. \n",
     "\n",
     "#### Registers\n",
-    " - `a`: The x component of the first input elliptic curve point of bitsize `n`.\n",
-    " - `b`: The y component of the first input elliptic curve point of bitsize `n`.\n",
-    " - `x`: The x component of the second input elliptic curve point of bitsize `n`, which will contain the x component of the resultant curve point.\n",
-    " - `y`: The y component of the second input elliptic curve point of bitsize `n`, which will contain the y component of the resultant curve point.\n",
-    " - `lam`: The precomputed lambda slope used in the addition operation. \n",
+    " - `a`: The x component of the first input elliptic curve point of bitsize `n` in montgomery form.\n",
+    " - `b`: The y component of the first input elliptic curve point of bitsize `n` in montgomery form.\n",
+    " - `x`: The x component of the second input elliptic curve point of bitsize `n` in montgomery form, which will contain the x component of the resultant curve point.\n",
+    " - `y`: The y component of the second input elliptic curve point of bitsize `n` in montgomery form, which will contain the y component of the resultant curve point.\n",
+    " - `lam_r`: The precomputed lambda slope used in the addition operation if (a, b) = (x, y) in montgomery form. \n",
     "\n",
     "#### References\n",
     " - [How to compute a 256-bit elliptic curve private key with only 50 million Toffoli gates](https://arxiv.org/abs/2306.08585). Litinski. 2023. Fig 5.\n"
@@ -91,6 +97,18 @@
     "ec_add = ECAdd(n, mod=p)"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "170da165",
+   "metadata": {
+    "cq.autogen": "ECAdd.ec_add_small"
+   },
+   "outputs": [],
+   "source": [
+    "ec_add_small = ECAdd(5, mod=7)"
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "39210af4",
@@ -111,8 +129,8 @@
    "outputs": [],
    "source": [
     "from qualtran.drawing import show_bloqs\n",
-    "show_bloqs([ec_add],\n",
-    "           ['`ec_add`'])"
+    "show_bloqs([ec_add, ec_add_small],\n",
+    "           ['`ec_add`', '`ec_add_small`'])"
    ]
   },
   {
@@ -157,7 +175,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.7"
+   "version": "3.10.13"
   }
  },
  "nbformat": 4,