From 8ff6b5f12cccb33fba3f4bc03bc0ece2e885f5f8 Mon Sep 17 00:00:00 2001
From: Anurudh Peduri <anurudhp97@gmail.com>
Date: Sat, 10 Feb 2024 21:43:49 +0100
Subject: [PATCH 01/34] wip: negative powers

---
 qualtran/bloqs/generalized_qsp.py      | 81 ++++++++++++++++++--------
 qualtran/bloqs/generalized_qsp_test.py |  8 ++-
 2 files changed, 63 insertions(+), 26 deletions(-)

diff --git a/qualtran/bloqs/generalized_qsp.py b/qualtran/bloqs/generalized_qsp.py
index b430a3233..69e7e9fd5 100644
--- a/qualtran/bloqs/generalized_qsp.py
+++ b/qualtran/bloqs/generalized_qsp.py
@@ -12,11 +12,11 @@
 #  See the License for the specific language governing permissions and
 #  limitations under the License.
 from functools import cached_property
-from typing import Sequence, Tuple
+from typing import Optional, Sequence, Tuple
 
 import cirq
 import numpy as np
-from attrs import frozen
+from attrs import field, frozen
 from numpy.polynomial import Polynomial
 from numpy.typing import NDArray
 
@@ -123,6 +123,7 @@ def qsp_complementary_polynomial(
     smaller_roots: list[complex] = []  # roots r s.t. \abs{r} < 1
 
     for r in roots:
+        assert not (np.abs(r) <= 1e-5), "zero root!"
         if np.allclose(np.abs(r), 1):
             units.append(r)
         elif np.abs(r) > 1:
@@ -130,22 +131,6 @@ def qsp_complementary_polynomial(
         else:
             smaller_roots.append(r)
 
-    if verify:
-        # verify that the non-unit roots indeed occur in conjugate pairs.
-        def is_permutation(A, B):
-            assert len(A) == len(B)
-            A = list(A)
-            for z in B:
-                for w in A:
-                    if np.allclose(z, w):
-                        A.remove(w)
-                        break
-                else:
-                    return False
-            return True
-
-        assert is_permutation(smaller_roots, 1 / np.array(larger_roots).conj())
-
     # pair up roots in `units`, claimed in Eq. 40 and the explanation preceding it.
     # all unit roots must have even multiplicity.
     paired_units: list[complex] = []
@@ -157,14 +142,45 @@ def is_permutation(A, B):
                 matched_z = w
                 break
 
-        if matched_z:
+        if matched_z is not None:
             paired_units.append(z)
             unpaired_units.remove(matched_z)
         else:
             unpaired_units.append(z)
 
+    unpaired_conj_units: list[complex] = []
+    for z in unpaired_units:
+        matched_z_conj = None
+        for w in unpaired_conj_units:
+            if np.allclose(z.conjugate(), w):
+                matched_z_conj = w
+                break
+
+        if matched_z_conj is not None:
+            smaller_roots.append(z)
+            larger_roots.append(matched_z_conj)
+            unpaired_conj_units.remove(matched_z_conj)
+        else:
+            unpaired_conj_units.append(z)
+
     if verify:
-        assert len(unpaired_units) == 0
+        assert len(unpaired_conj_units) == 0
+
+        # verify that the non-unit roots indeed occur in conjugate pairs.
+        def assert_is_permutation(A, B):
+            assert len(A) == len(B)
+            A = list(A)
+            unmatched = []
+            for z in B:
+                for w in A:
+                    if np.allclose(z, w, rtol=1e-5, atol=1e-5):
+                        A.remove(w)
+                        break
+                else:
+                    unmatched.append(z)
+            assert len(unmatched) == 0
+
+        assert_is_permutation(smaller_roots, 1 / np.array(larger_roots).conj())
 
     # Q = G \hat{G}, where
     # - \hat{G}^2 is the monomials which are unit roots of R, which occur in pairs.
@@ -239,9 +255,11 @@ def safe_angle(x):
 class GeneralizedQSP(GateWithRegisters):
     r"""Applies a QSP polynomial $P$ to a unitary $U$ to obtain a block-encoding of $P(U)$.
 
+    Can optionally provide a negative power offset $k$ (defaults to 0),
+    to obtain $U^{-k} P(U)$. (Theorem 6)
     This gate represents the following unitary:
 
-        $$ \begin{bmatrix} P(U) & \cdot \\ \cdot & \cdot \end{bmatrix} $$
+        $$ \begin{bmatrix} U^{-k} P(U) & \cdot \\ \cdot & \cdot \end{bmatrix} $$
 
     The polynomial $P$ must satisfy:
     $\abs{P(e^{i \theta})}^2 \le 1$ for every $\theta \in \mathbb{R}$.
@@ -251,14 +269,17 @@ class GeneralizedQSP(GateWithRegisters):
     Args:
         U: Unitary operation.
         P: Co-efficients of a complex polynomial.
+        negative_power: value of $k$, which effectively applies $z^{-k} P(z)$. defaults to 0.
 
     References:
         [Generalized Quantum Signal Processing](https://arxiv.org/abs/2308.01501)
-            Motlagh and Wiebe. (2023). Theorem 3; Figure 2.
+            Motlagh and Wiebe. (2023). Theorem 3; Figure 2; Theorem 6.
     """
 
     U: GateWithRegisters
     P: Sequence[complex]
+    _precomputed_Q: Optional[Sequence[complex]] = None
+    negative_power: int = field(default=0, kw_only=True)
 
     @cached_property
     def signature(self) -> Signature:
@@ -266,6 +287,8 @@ def signature(self) -> Signature:
 
     @cached_property
     def Q(self):
+        if self._precomputed_Q is not None:
+            return self._precomputed_Q
         return qsp_complementary_polynomial(self.P)
 
     @cached_property
@@ -290,7 +313,19 @@ def decompose_from_registers(
         assert len(signal) == 1
         signal_qubit = signal[0]
 
+        num_inverse_applications = self.negative_power
+
         yield SU2RotationGate(self._theta[0], self._phi[0], self._lambda).on(signal_qubit)
         for theta, phi in zip(self._theta[1:], self._phi[1:]):
-            yield self.U.on_registers(**quregs).controlled_by(signal_qubit, control_values=[0])
+            if num_inverse_applications > 0:
+                # apply C-U^\dagger
+                yield self.U.adjoint().on_registers(**quregs).controlled_by(signal_qubit)
+                num_inverse_applications -= 1
+            else:
+                # apply C[0]-U
+                yield self.U.on_registers(**quregs).controlled_by(signal_qubit, control_values=[0])
             yield SU2RotationGate(theta, phi, 0).on(signal_qubit)
+
+        while num_inverse_applications > 0:
+            yield self.U.adjoint().on_registers(**quregs)
+            num_inverse_applications -= 1
diff --git a/qualtran/bloqs/generalized_qsp_test.py b/qualtran/bloqs/generalized_qsp_test.py
index 7462618fe..346261676 100644
--- a/qualtran/bloqs/generalized_qsp_test.py
+++ b/qualtran/bloqs/generalized_qsp_test.py
@@ -12,7 +12,7 @@
 #  See the License for the specific language governing permissions and
 #  limitations under the License.
 from functools import cached_property
-from typing import Sequence, Union
+from typing import Optional, Sequence, Union
 
 import cirq
 import numpy as np
@@ -146,10 +146,12 @@ def assert_matrices_almost_equal(A: NDArray, B: NDArray):
     assert np.linalg.norm(A - B) <= 1e-5
 
 
-def verify_generalized_qsp(U: GateWithRegisters, P: Sequence[complex]):
+def verify_generalized_qsp(
+    U: GateWithRegisters, P: Sequence[complex], Q: Optional[Sequence[complex]] = None
+):
     input_unitary = cirq.unitary(U)
     N = input_unitary.shape[0]
-    gqsp_U = GeneralizedQSP(U, P)
+    gqsp_U = GeneralizedQSP(U, P, Q)
     result_unitary = cirq.unitary(gqsp_U)
 
     expected_top_left = evaluate_polynomial_of_matrix(P, input_unitary)

From c56d7084b06aa7464dbe38e7335103b021fd30c3 Mon Sep 17 00:00:00 2001
From: Anurudh Peduri <anurudhp97@gmail.com>
Date: Thu, 15 Feb 2024 16:40:02 +0100
Subject: [PATCH 02/34] notebook: Hamiltonian simulation example.

---
 qualtran/bloqs/generalized_qsp.ipynb | 103 +++++++++++++++++++++++----
 qualtran/bloqs/generalized_qsp.py    |  21 +++---
 2 files changed, 101 insertions(+), 23 deletions(-)

diff --git a/qualtran/bloqs/generalized_qsp.ipynb b/qualtran/bloqs/generalized_qsp.ipynb
index 4ca71bdef..ff3a6d875 100644
--- a/qualtran/bloqs/generalized_qsp.ipynb
+++ b/qualtran/bloqs/generalized_qsp.ipynb
@@ -15,12 +15,33 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from qualtran.bloqs.qubitization_walk_operator_test import get_walk_operator_for_1d_Ising_model\n",
+    "import numpy as np\n",
+    "import scipy\n",
     "from qualtran.drawing import show_bloq\n",
     "\n",
+    "from qualtran.bloqs.qubitization_walk_operator_test import get_walk_operator_for_1d_Ising_model\n",
+    "from qualtran.bloqs.qubitization_walk_operator import QubitizationWalkOperator\n",
+    "\n",
     "from qualtran.bloqs.generalized_qsp import GeneralizedQSP"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "6963c30f339d42de",
+   "metadata": {},
+   "source": [
+    "`GeneralizedQSP` implements the Quantum Eigenvalue Transform on a unitary $U$ using QSP. Given a complex GQSP polynomial $P$ (and its complement $Q$), it implements the unitary:\n",
+    "$$U' = \\begin{bmatrix} P(U) & \\cdot \\\\ Q(U) & \\cdot \\end{bmatrix}$$\n",
+    "\n",
+    "Here, the polynomials $P, Q$ must satisfy the following constraint:\n",
+    "\n",
+    "$$\\left| P(e^{i\\theta}) \\right|^2 + \\left| Q(e^{i\\theta}) \\right|^2 = 1 ~~\\text{for every}~ \\theta \\in [0, 2\\pi]$$\n",
+    "\n",
+    "A polynomial $P$ is said to be a GQSP polynomial iff it satisfies $\\left| P(e^{i\\theta}) \\right|^2 \\le 1$ for every $\\theta \\in [0, 2\\pi]$. \n",
+    "\n",
+    "Reference: https://doi.org/10.48550/arXiv.2308.01501"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -32,42 +53,94 @@
     "show_bloq(U.decompose_bloq())"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "dc7100fd94d58c6b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pU = GeneralizedQSP(U, (0.5, 0.5))\n",
+    "show_bloq(pU.decompose_bloq())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "78cd3857297f092b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pU = GeneralizedQSP(U, (0.5, 0, 0.5))\n",
+    "show_bloq(pU.decompose_bloq())"
+   ]
+  },
   {
    "cell_type": "markdown",
-   "id": "6963c30f339d42de",
+   "id": "c84e5443-483c-4294-8a79-1d2bed70e1f5",
    "metadata": {},
    "source": [
-    "`GeneralizedQSP` implements the Quantum Eigenvalue Transform on a unitary $U$ using QSP. Given a complex QSP polynomial $P$ (and its complement $Q$), it implements the unitary:\n",
-    "$$U' = \\begin{bmatrix} P(U) & \\cdot \\\\ Q(U) & \\cdot \\end{bmatrix}$$\n",
+    "## Hamiltonian Simulation by GQSP\n",
     "\n",
-    "Here, the polynomials $P, Q$ must satisfy the following constraint:\n",
+    "Given the Szegedy Quantum Walk Operator for a Hamiltonian $H$, one can construct the walk operator for $e^{-iHt}$ using GQSP (Corollary 8).\n",
+    "\n",
+    "### Recap:\n",
     "\n",
-    "$$\\left\\mid P(e^{i\\theta}) \\right\\mid^2 + \\left\\mid Q(e^{i\\theta}) \\right\\mid^2 = 1 ~~\\text{for every}~ \\theta \\in [0, 2\\pi]$$\n",
+    "For a Hamiltonian $H = \\sum_i \\alpha_i U_i$, given the SELECT and PREPARE oracles\n",
+    "$$ \\text{SELECT} = \\sum_j \\ketbra{j}{j} \\otimes U_j $$\n",
+    "$$ \\text{PREPARE} \\ket{0} = \\sum_j \\frac{\\sqrt{\\alpha_j}}{\\|\\alpha\\|_1} \\ket{j} $$\n",
+    "we can implement the [QubitizationWalkOperator](qubitization_walk_operator.ipynb) that encodes the spectrum of $H$ in the eigenphases of the walk operator $W$.\n",
     "\n",
+    "### Approximating $\\cos$\n",
+    "We can use the Jacobi-Anger expansion to obtain low-degree polynomial approximations for the $\\cos$ function:\n",
     "\n",
-    "Reference: https://arxiv.org/abs/2308.01501"
+    "$$e^{it\\cos\\theta} = \\sum_{n = -\\infty}^{\\infty} i^n J_n(t) (e^{i\\theta})^n$$\n",
+    "\n",
+    "We can cutoff at $d = O(t + \\log(1/\\epsilon) / \\log\\log(1/\\epsilon))$ to get an $\\epsilon$-approximation (Theorem 7):\n",
+    "\n",
+    "$$P[t](z) = \\sum_{n = -d}^d i^n J_n(t) z^n$$\n",
+    "\n",
+    "### Obtaining $e^{iHt}$\n",
+    "\n",
+    "As the eigenphases of the walk operator are $e^{-i\\arccos(E_k / \\|\\alpha\\|_1)}$, we can use the GQSP polynomial with $P = P[\\|\\alpha\\|_1 t]$ (and complementary $Q = 0$) to obtain $P(U) = e^{-iHt}$.\n",
+    "The obtained GQSP operator $W'$ can then be used with the PREPARE oracle to simulate the hamiltonian:\n",
+    "$$(\\text{PREPARE}^\\dagger \\otimes I) W' (\\text{PREPARE} \\otimes I) |0\\rangle|\\psi\\rangle = |0\\rangle e^{-iHt}|\\psi\\rangle$$"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "dc7100fd94d58c6b",
-   "metadata": {},
+   "id": "721cd78c6d6e60c4",
+   "metadata": {
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
    "outputs": [],
    "source": [
-    "pU = GeneralizedQSP(U, (0.5, 0.5))\n",
-    "show_bloq(pU.decompose_bloq())"
+    "def hamiltonian_simulation_by_gqsp(\n",
+    "        W: QubitizationWalkOperator, t: float, *, precision=1e-5, max_degree=None\n",
+    ") -> GeneralizedQSP:\n",
+    "    degree = t + 3 * np.log(1/precision) / np.log(np.log(1/precision))\n",
+    "    if max_degree is not None:\n",
+    "        degree = min(max_degree, degree)\n",
+    "    degree = int(np.round(degree))\n",
+    "\n",
+    "    coeff_indices = np.arange(-degree, degree + 1)\n",
+    "    approx_cos = 1j**coeff_indices * scipy.special.jv(coeff_indices, t)\n",
+    "\n",
+    "    return GeneralizedQSP(W, approx_cos, np.zeros(2*degree + 1), negative_power=degree)"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "78cd3857297f092b",
+   "id": "c4ccd9c8-92ce-4b53-92bc-722f263a32ad",
    "metadata": {},
    "outputs": [],
    "source": [
-    "pU = GeneralizedQSP(U, (0.5, 0, 0.5))\n",
-    "show_bloq(pU.decompose_bloq())"
+    "W_e_iHt = hamiltonian_simulation_by_gqsp(U, 5)\n",
+    "show_bloq(W_e_iHt.decompose_bloq())"
    ]
   }
  ],
@@ -87,7 +160,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.9"
+   "version": "3.11.6"
   }
  },
  "nbformat": 4,
diff --git a/qualtran/bloqs/generalized_qsp.py b/qualtran/bloqs/generalized_qsp.py
index 69e7e9fd5..3daf70f1a 100644
--- a/qualtran/bloqs/generalized_qsp.py
+++ b/qualtran/bloqs/generalized_qsp.py
@@ -20,7 +20,7 @@
 from numpy.polynomial import Polynomial
 from numpy.typing import NDArray
 
-from qualtran import GateWithRegisters, QBit, Register, Signature
+from qualtran import GateWithRegisters, QBit, Register, Signature, Adjoint
 
 
 @frozen
@@ -317,15 +317,20 @@ def decompose_from_registers(
 
         yield SU2RotationGate(self._theta[0], self._phi[0], self._lambda).on(signal_qubit)
         for theta, phi in zip(self._theta[1:], self._phi[1:]):
-            if num_inverse_applications > 0:
-                # apply C-U^\dagger
-                yield self.U.adjoint().on_registers(**quregs).controlled_by(signal_qubit)
-                num_inverse_applications -= 1
-            else:
-                # apply C[0]-U
-                yield self.U.on_registers(**quregs).controlled_by(signal_qubit, control_values=[0])
+            # TODO debug controlled adjoint bloq serialization
+            # if num_inverse_applications > 0:
+            #     # apply C-U^\dagger
+            #     yield self.U.adjoint().on_registers(**quregs).controlled_by(signal_qubit)
+            #     num_inverse_applications -= 1
+            # else:
+            # apply C[0]-U
+            yield self.U.on_registers(**quregs).controlled_by(signal_qubit, control_values=[0])
             yield SU2RotationGate(theta, phi, 0).on(signal_qubit)
 
         while num_inverse_applications > 0:
             yield self.U.adjoint().on_registers(**quregs)
             num_inverse_applications -= 1
+
+    def __hash__(self):
+        # TODO hash properly
+        return hash((self.U, *list(self.P), *list(self.Q), self.negative_power))

From 845f03a24a884c5d04663900743d139b2a269e8b Mon Sep 17 00:00:00 2001
From: Anurudh Peduri <anurudhp97@gmail.com>
Date: Thu, 15 Feb 2024 16:58:23 +0100
Subject: [PATCH 03/34] notebook

---
 qualtran/bloqs/generalized_qsp.ipynb | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/qualtran/bloqs/generalized_qsp.ipynb b/qualtran/bloqs/generalized_qsp.ipynb
index ff3a6d875..7c6b35c0b 100644
--- a/qualtran/bloqs/generalized_qsp.ipynb
+++ b/qualtran/bloqs/generalized_qsp.ipynb
@@ -104,7 +104,7 @@
     "\n",
     "As the eigenphases of the walk operator are $e^{-i\\arccos(E_k / \\|\\alpha\\|_1)}$, we can use the GQSP polynomial with $P = P[\\|\\alpha\\|_1 t]$ (and complementary $Q = 0$) to obtain $P(U) = e^{-iHt}$.\n",
     "The obtained GQSP operator $W'$ can then be used with the PREPARE oracle to simulate the hamiltonian:\n",
-    "$$(\\text{PREPARE}^\\dagger \\otimes I) W' (\\text{PREPARE} \\otimes I) |0\\rangle|\\psi\\rangle = |0\\rangle e^{-iHt}|\\psi\\rangle$$"
+    "$$(\\langle0| \\otimes \\text{PREPARE}^\\dagger \\otimes I) W' (|0\\rangle \\otimes \\text{PREPARE} \\otimes I) |0\\rangle|\\psi\\rangle = |0\\rangle e^{-iHt}|\\psi\\rangle$$"
    ]
   },
   {

From d7e35c703c468a5e58169bd8e6e161082b55279e Mon Sep 17 00:00:00 2001
From: Anurudh Peduri <anurudhp97@gmail.com>
Date: Thu, 15 Feb 2024 17:36:43 +0100
Subject: [PATCH 04/34] lint + TODO comment

---
 qualtran/bloqs/generalized_qsp.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/qualtran/bloqs/generalized_qsp.py b/qualtran/bloqs/generalized_qsp.py
index 3daf70f1a..cab51fb19 100644
--- a/qualtran/bloqs/generalized_qsp.py
+++ b/qualtran/bloqs/generalized_qsp.py
@@ -22,6 +22,8 @@
 
 from qualtran import GateWithRegisters, QBit, Register, Signature, Adjoint
 
+# TODO refactor using Bloq instead of cirq-ft infra
+
 
 @frozen
 class SU2RotationGate(GateWithRegisters):

From 814cd51f179c732c5209ab71d9bb783a606811fa Mon Sep 17 00:00:00 2001
From: Anurudh Peduri <anurudhp97@gmail.com>
Date: Thu, 15 Feb 2024 18:14:19 +0100
Subject: [PATCH 05/34] nb: explain negative power terms

---
 qualtran/bloqs/generalized_qsp.ipynb | 22 ++++++++++++++++++++++
 1 file changed, 22 insertions(+)

diff --git a/qualtran/bloqs/generalized_qsp.ipynb b/qualtran/bloqs/generalized_qsp.ipynb
index 7c6b35c0b..7acfb310c 100644
--- a/qualtran/bloqs/generalized_qsp.ipynb
+++ b/qualtran/bloqs/generalized_qsp.ipynb
@@ -75,6 +75,28 @@
     "show_bloq(pU.decompose_bloq())"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "a58f06ba-9287-435d-92a3-256f747024c2",
+   "metadata": {},
+   "source": [
+    "### Negative degree terms\n",
+    "\n",
+    "To apply GQSP for a polynomial $P'(z) = z^{-k} P(z)$, we can just pass the polynomial $P$ along with negative power $k$.\n",
+    "The QSP angle sequence is the same for both, and $P'$ can be achieved by running $(U^\\dagger)^k$ at any point in the circuit."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ee60e95f-979f-420d-91b1-b6955b9c5a3b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pU = GeneralizedQSP(U, (0.5, 0, 0.5), negative_power=1)\n",
+    "show_bloq(pU.decompose_bloq())"
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "c84e5443-483c-4294-8a79-1d2bed70e1f5",

From 96046135dba9796b7f0eac19ec00a9a6bfb58e3a Mon Sep 17 00:00:00 2001
From: Anurudh Peduri <anurudhp97@gmail.com>
Date: Fri, 16 Feb 2024 20:01:19 +0100
Subject: [PATCH 06/34] call graph for GQSP

---
 qualtran/bloqs/generalized_qsp.ipynb   |  58 ++++++++------
 qualtran/bloqs/generalized_qsp.py      | 104 ++++++++++++++++++-------
 qualtran/bloqs/generalized_qsp_test.py |   8 +-
 3 files changed, 115 insertions(+), 55 deletions(-)

diff --git a/qualtran/bloqs/generalized_qsp.ipynb b/qualtran/bloqs/generalized_qsp.ipynb
index 7acfb310c..70bcd2839 100644
--- a/qualtran/bloqs/generalized_qsp.ipynb
+++ b/qualtran/bloqs/generalized_qsp.ipynb
@@ -19,10 +19,13 @@
     "import scipy\n",
     "from qualtran.drawing import show_bloq\n",
     "\n",
+    "import qualtran\n",
     "from qualtran.bloqs.qubitization_walk_operator_test import get_walk_operator_for_1d_Ising_model\n",
     "from qualtran.bloqs.qubitization_walk_operator import QubitizationWalkOperator\n",
+    "from qualtran.cirq_interop import CirqGateAsBloq\n",
     "\n",
-    "from qualtran.bloqs.generalized_qsp import GeneralizedQSP"
+    "from qualtran.bloqs.generalized_qsp import GeneralizedQSP, HamiltonianSimulationByGQSP\n",
+    "from qualtran.bloqs.generalized_qsp_test import RandomGate"
    ]
   },
   {
@@ -60,10 +63,18 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "pU = GeneralizedQSP(U, (0.5, 0.5))\n",
+    "pU = GeneralizedQSP(U, (0.5, 0.5), (-0.5, 0.5))\n",
     "show_bloq(pU.decompose_bloq())"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "935a03f7-5843-4b11-abe6-5eb9048c0ab5",
+   "metadata": {},
+   "source": [
+    "There is also a method that directly computes $Q$ from $P$:"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -71,7 +82,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "pU = GeneralizedQSP(U, (0.5, 0, 0.5))\n",
+    "pU = GeneralizedQSP.from_qsp_polynomial(U, (0.5, 0, 0.5))\n",
     "show_bloq(pU.decompose_bloq())"
    ]
   },
@@ -93,7 +104,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "pU = GeneralizedQSP(U, (0.5, 0, 0.5), negative_power=1)\n",
+    "pU = GeneralizedQSP.from_qsp_polynomial(U, (0.5, 0, 0.5), negative_power=1)\n",
     "show_bloq(pU.decompose_bloq())"
    ]
   },
@@ -118,6 +129,8 @@
     "\n",
     "$$e^{it\\cos\\theta} = \\sum_{n = -\\infty}^{\\infty} i^n J_n(t) (e^{i\\theta})^n$$\n",
     "\n",
+    "where $J_n$ is the $n$-th [Bessel function of the first kind](https://en.wikipedia.org/wiki/Bessel_function#Bessel_functions_of_the_first_kind), which is provided by `scipy.special.jv`.\n",
+    "\n",
     "We can cutoff at $d = O(t + \\log(1/\\epsilon) / \\log\\log(1/\\epsilon))$ to get an $\\epsilon$-approximation (Theorem 7):\n",
     "\n",
     "$$P[t](z) = \\sum_{n = -d}^d i^n J_n(t) z^n$$\n",
@@ -132,37 +145,32 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "721cd78c6d6e60c4",
-   "metadata": {
-    "jupyter": {
-     "outputs_hidden": false
-    }
-   },
+   "id": "c4ccd9c8-92ce-4b53-92bc-722f263a32ad",
+   "metadata": {},
    "outputs": [],
    "source": [
-    "def hamiltonian_simulation_by_gqsp(\n",
-    "        W: QubitizationWalkOperator, t: float, *, precision=1e-5, max_degree=None\n",
-    ") -> GeneralizedQSP:\n",
-    "    degree = t + 3 * np.log(1/precision) / np.log(np.log(1/precision))\n",
-    "    if max_degree is not None:\n",
-    "        degree = min(max_degree, degree)\n",
-    "    degree = int(np.round(degree))\n",
-    "\n",
-    "    coeff_indices = np.arange(-degree, degree + 1)\n",
-    "    approx_cos = 1j**coeff_indices * scipy.special.jv(coeff_indices, t)\n",
-    "\n",
-    "    return GeneralizedQSP(W, approx_cos, np.zeros(2*degree + 1), negative_power=degree)"
+    "W_e_iHt = HamiltonianSimulationByGQSP(U, t=5, alpha=1, precision=1e-3)\n",
+    "\n",
+    "def decomp_pred(inst):\n",
+    "    from qualtran.bloqs.select_and_prepare import PrepareOracle, SelectOracle\n",
+    "\n",
+    "    return any(isinstance(inst.bloq, cls) for cls in [HamiltonianSimulationByGQSP, GeneralizedQSP, SelectOracle, PrepareOracle])\n",
+    "    \n",
+    "show_bloq(\n",
+    "    W_e_iHt\n",
+    "    .decompose_bloq()\n",
+    "    .flatten(pred=decomp_pred)\n",
+    ")"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "c4ccd9c8-92ce-4b53-92bc-722f263a32ad",
+   "id": "c6627cb6-446d-4f76-b727-8ed2f67bc382",
    "metadata": {},
    "outputs": [],
    "source": [
-    "W_e_iHt = hamiltonian_simulation_by_gqsp(U, 5)\n",
-    "show_bloq(W_e_iHt.decompose_bloq())"
+    "W_e_iHt.call_graph()"
    ]
   }
  ],
diff --git a/qualtran/bloqs/generalized_qsp.py b/qualtran/bloqs/generalized_qsp.py
index cab51fb19..da2886cae 100644
--- a/qualtran/bloqs/generalized_qsp.py
+++ b/qualtran/bloqs/generalized_qsp.py
@@ -12,21 +12,26 @@
 #  See the License for the specific language governing permissions and
 #  limitations under the License.
 from functools import cached_property
-from typing import Optional, Sequence, Tuple
+from typing import Dict, Sequence, Set, Tuple, TYPE_CHECKING
 
 import cirq
 import numpy as np
+import scipy
 from attrs import field, frozen
 from numpy.polynomial import Polynomial
 from numpy.typing import NDArray
 
-from qualtran import GateWithRegisters, QBit, Register, Signature, Adjoint
+from qualtran import Bloq, GateWithRegisters, QBit, Register, Signature
+from qualtran.bloqs.basic_gates import Ry, ZPowGate
+from qualtran.bloqs.qubitization_walk_operator import QubitizationWalkOperator
 
-# TODO refactor using Bloq instead of cirq-ft infra
+if TYPE_CHECKING:
+    from qualtran import BloqBuilder, SoquetT
+    from qualtran.resource_counting import BloqCountT, SympySymbolAllocator
 
 
 @frozen
-class SU2RotationGate(GateWithRegisters):
+class SU2RotationGate(Bloq):
     theta: float
     phi: float
     lambd: float
@@ -65,15 +70,12 @@ def rotation_matrix(self):
             ]
         )
 
-    def decompose_from_registers(
-        self, *, context: cirq.DecompositionContext, q: NDArray[cirq.Qid]
-    ) -> cirq.OP_TREE:
-        qubit = q[0]
-
-        yield cirq.Rz(rads=np.pi - self.lambd).on(qubit)
-        yield cirq.Ry(rads=2 * self.theta).on(qubit)
-        yield cirq.Rz(rads=-self.phi).on(qubit)
-        yield cirq.GlobalPhaseGate(np.exp(1j * (np.pi + self.lambd + self.phi) / 2)).on()
+    def build_composite_bloq(self, bb: 'BloqBuilder', q: 'SoquetT') -> Dict[str, 'SoquetT']:
+        q = bb.add(ZPowGate(exponent=2, global_shift=0.5), q=q)
+        q = bb.add(ZPowGate(exponent=1 - self.lambd / np.pi, global_shift=-1), q=q)
+        q = bb.add(Ry(angle=2 * self.theta), q=q)
+        q = bb.add(ZPowGate(exponent=-self.phi / np.pi, global_shift=-1), q=q)
+        return {'q': q}
 
 
 def qsp_complementary_polynomial(
@@ -261,7 +263,7 @@ class GeneralizedQSP(GateWithRegisters):
     to obtain $U^{-k} P(U)$. (Theorem 6)
     This gate represents the following unitary:
 
-        $$ \begin{bmatrix} U^{-k} P(U) & \cdot \\ \cdot & \cdot \end{bmatrix} $$
+        $$ \begin{bmatrix} U^{-k} P(U) & \cdot \\ Q(U) & \cdot \end{bmatrix} $$
 
     The polynomial $P$ must satisfy:
     $\abs{P(e^{i \theta})}^2 \le 1$ for every $\theta \in \mathbb{R}$.
@@ -280,18 +282,19 @@ class GeneralizedQSP(GateWithRegisters):
 
     U: GateWithRegisters
     P: Sequence[complex]
-    _precomputed_Q: Optional[Sequence[complex]] = None
+    Q: Sequence[complex]
     negative_power: int = field(default=0, kw_only=True)
 
     @cached_property
     def signature(self) -> Signature:
         return Signature([Register('signal', QBit()), *self.U.signature])
 
-    @cached_property
-    def Q(self):
-        if self._precomputed_Q is not None:
-            return self._precomputed_Q
-        return qsp_complementary_polynomial(self.P)
+    @staticmethod
+    def from_qsp_polynomial(
+        U: GateWithRegisters, P: Sequence[complex], *, negative_power: int = 0
+    ) -> 'GeneralizedQSP':
+        Q = qsp_complementary_polynomial(P)
+        return GeneralizedQSP(U, P, Q, negative_power=negative_power)
 
     @cached_property
     def _qsp_phases(self) -> Tuple[Sequence[float], Sequence[float], float]:
@@ -309,16 +312,24 @@ def _phi(self) -> Sequence[float]:
     def _lambda(self) -> float:
         return self._qsp_phases[2]
 
+    @cached_property
+    def signal_rotations(self) -> NDArray[SU2RotationGate]:
+        return np.array(
+            [
+                SU2RotationGate(theta, phi, self._lambda if i == 0 else 0)
+                for i, (theta, phi) in enumerate(zip(self._theta, self._phi))
+            ]
+        )
+
     def decompose_from_registers(
         self, *, context: cirq.DecompositionContext, signal, **quregs: NDArray[cirq.Qid]
     ) -> cirq.OP_TREE:
-        assert len(signal) == 1
         signal_qubit = signal[0]
 
         num_inverse_applications = self.negative_power
 
-        yield SU2RotationGate(self._theta[0], self._phi[0], self._lambda).on(signal_qubit)
-        for theta, phi in zip(self._theta[1:], self._phi[1:]):
+        yield self.signal_rotations[0].on(signal_qubit)
+        for signal_rotation in self.signal_rotations[1:]:
             # TODO debug controlled adjoint bloq serialization
             # if num_inverse_applications > 0:
             #     # apply C-U^\dagger
@@ -327,12 +338,49 @@ def decompose_from_registers(
             # else:
             # apply C[0]-U
             yield self.U.on_registers(**quregs).controlled_by(signal_qubit, control_values=[0])
-            yield SU2RotationGate(theta, phi, 0).on(signal_qubit)
+            yield signal_rotation.on(signal_qubit)
 
-        while num_inverse_applications > 0:
+        for _ in range(num_inverse_applications):
             yield self.U.adjoint().on_registers(**quregs)
-            num_inverse_applications -= 1
 
     def __hash__(self):
-        # TODO hash properly
-        return hash((self.U, *list(self.P), *list(self.Q), self.negative_power))
+        return hash((self.U, *self.signal_rotations, self.negative_power))
+
+    def build_call_graph(self, ssa: 'SympySymbolAllocator') -> Set['BloqCountT']:
+        return {(self.U, max(len(self.P), self.negative_power))} | {
+            (rotation, 1) for rotation in self.signal_rotations
+        }
+
+
+@frozen
+class HamiltonianSimulationByGQSP(GateWithRegisters):
+    WalkOperator: QubitizationWalkOperator
+    t: float = field(kw_only=True)
+    alpha: float = field(kw_only=True)
+    precision: float = field(kw_only=True)
+
+    @cached_property
+    def degree(self) -> int:
+        log_precision = np.log(1 / self.precision)
+        degree = self.t * self.alpha + log_precision / np.log(log_precision)
+        return int(np.ceil(degree))
+
+    @cached_property
+    def gqsp(self) -> GeneralizedQSP:
+        coeff_indices = np.arange(-self.degree, self.degree + 1)
+        approx_cos = 1j**coeff_indices * scipy.special.jv(coeff_indices, self.t * self.alpha)
+
+        return GeneralizedQSP(
+            self.WalkOperator, approx_cos, np.zeros(2 * self.degree + 1), negative_power=self.degree
+        )
+
+    @cached_property
+    def signature(self) -> 'Signature':
+        return self.gqsp.signature
+
+    def build_composite_bloq(self, bb: 'BloqBuilder', **soqs: 'SoquetT') -> Dict[str, 'SoquetT']:
+        # TODO call prepare oracle
+        # soqs |= bb.add_d(self.WalkOperator.prepare, selection=soqs['selection'])
+        soqs = bb.add_d(self.gqsp, **soqs)
+        # soqs |= bb.add_d(self.WalkOperator.prepare.adjoint(), selection=soqs['selection'])
+        return soqs
diff --git a/qualtran/bloqs/generalized_qsp_test.py b/qualtran/bloqs/generalized_qsp_test.py
index 346261676..bac9508aa 100644
--- a/qualtran/bloqs/generalized_qsp_test.py
+++ b/qualtran/bloqs/generalized_qsp_test.py
@@ -30,6 +30,7 @@
     qsp_phase_factors,
     SU2RotationGate,
 )
+from qualtran.cirq_interop import BloqAsCirqGate
 
 
 def assert_angles_almost_equal(
@@ -50,7 +51,7 @@ def test_cirq_decompose_SU2_to_single_qubit_pauli_gates():
         gate = SU2RotationGate(theta, phi, lambd)
 
         expected = gate.rotation_matrix
-        actual = cirq.unitary(gate)
+        actual = cirq.unitary(BloqAsCirqGate(gate))
         np.testing.assert_allclose(actual, expected)
 
 
@@ -151,7 +152,10 @@ def verify_generalized_qsp(
 ):
     input_unitary = cirq.unitary(U)
     N = input_unitary.shape[0]
-    gqsp_U = GeneralizedQSP(U, P, Q)
+    if Q is None:
+        gqsp_U = GeneralizedQSP.from_qsp_polynomial(U, P)
+    else:
+        gqsp_U = GeneralizedQSP(U, P, Q)
     result_unitary = cirq.unitary(gqsp_U)
 
     expected_top_left = evaluate_polynomial_of_matrix(P, input_unitary)

From 47e0dfb32e2754658b59991c024e609dcb14eef7 Mon Sep 17 00:00:00 2001
From: Anurudh Peduri <anurudhp97@gmail.com>
Date: Sun, 18 Feb 2024 21:48:24 +0100
Subject: [PATCH 07/34] test negative power poly

---
 qualtran/bloqs/generalized_qsp_test.py | 23 ++++++++++++++++-------
 1 file changed, 16 insertions(+), 7 deletions(-)

diff --git a/qualtran/bloqs/generalized_qsp_test.py b/qualtran/bloqs/generalized_qsp_test.py
index bac9508aa..847af852d 100644
--- a/qualtran/bloqs/generalized_qsp_test.py
+++ b/qualtran/bloqs/generalized_qsp_test.py
@@ -129,10 +129,12 @@ def _unitary_(self):
         return self.matrix
 
 
-def evaluate_polynomial_of_matrix(P: Sequence[complex], U: NDArray) -> NDArray:
+def evaluate_polynomial_of_matrix(
+    P: Sequence[complex], U: NDArray, *, negative_power: int = 0
+) -> NDArray:
     assert U.ndim == 2 and U.shape[0] == U.shape[1]
 
-    pow_U = np.identity(U.shape[0], dtype=U.dtype)
+    pow_U = np.linalg.matrix_power(U.conjugate().T, negative_power)
     result = np.zeros(U.shape, dtype=U.dtype)
 
     for c in P:
@@ -148,14 +150,18 @@ def assert_matrices_almost_equal(A: NDArray, B: NDArray):
 
 
 def verify_generalized_qsp(
-    U: GateWithRegisters, P: Sequence[complex], Q: Optional[Sequence[complex]] = None
+    U: GateWithRegisters,
+    P: Sequence[complex],
+    Q: Optional[Sequence[complex]] = None,
+    *,
+    negative_power: int = 0,
 ):
     input_unitary = cirq.unitary(U)
     N = input_unitary.shape[0]
     if Q is None:
-        gqsp_U = GeneralizedQSP.from_qsp_polynomial(U, P)
+        gqsp_U = GeneralizedQSP.from_qsp_polynomial(U, P, negative_power=negative_power)
     else:
-        gqsp_U = GeneralizedQSP(U, P, Q)
+        gqsp_U = GeneralizedQSP(U, P, Q, negative_power=negative_power)
     result_unitary = cirq.unitary(gqsp_U)
 
     expected_top_left = evaluate_polynomial_of_matrix(P, input_unitary)
@@ -182,13 +188,16 @@ def test_generalized_qsp_with_real_poly_on_random_unitaries(bitsize: int, degree
 @pytest.mark.slow
 @pytest.mark.parametrize("bitsize", [1, 2, 3])
 @pytest.mark.parametrize("degree", [2, 3, 4, 5, 50, 100, 120])
-def test_generalized_qsp_with_complex_poly_on_random_unitaries(bitsize: int, degree: int):
+@pytest.mark.parametrize("negative_power", [0, 1, 2])
+def test_generalized_qsp_with_complex_poly_on_random_unitaries(
+    bitsize: int, degree: int, negative_power: int
+):
     random_state = np.random.RandomState(42)
 
     for _ in range(10):
         U = RandomGate.create(bitsize, random_state=random_state)
         P = random_qsp_polynomial(degree, random_state=random_state)
-        verify_generalized_qsp(U, P)
+        verify_generalized_qsp(U, P, negative_power=negative_power)
 
 
 @define(slots=False)

From 38bfd2f243a45b15b44ef8dfc6a3870242978bb3 Mon Sep 17 00:00:00 2001
From: Anurudh Peduri <anurudhp97@gmail.com>
Date: Sun, 18 Feb 2024 21:48:58 +0100
Subject: [PATCH 08/34] call prepare oracle correctly

---
 qualtran/bloqs/generalized_qsp.py | 29 ++++++++++++++++++++++++++---
 1 file changed, 26 insertions(+), 3 deletions(-)

diff --git a/qualtran/bloqs/generalized_qsp.py b/qualtran/bloqs/generalized_qsp.py
index da2886cae..1fd7a74f3 100644
--- a/qualtran/bloqs/generalized_qsp.py
+++ b/qualtran/bloqs/generalized_qsp.py
@@ -379,8 +379,31 @@ def signature(self) -> 'Signature':
         return self.gqsp.signature
 
     def build_composite_bloq(self, bb: 'BloqBuilder', **soqs: 'SoquetT') -> Dict[str, 'SoquetT']:
-        # TODO call prepare oracle
-        # soqs |= bb.add_d(self.WalkOperator.prepare, selection=soqs['selection'])
+        prepare = self.WalkOperator.prepare
+        state_prep_ancilla = {
+            reg.name: bb.allocate(reg.total_bits()) for reg in prepare.junk_registers
+        }
+
+        # PREPARE
+        prepare_soqs = bb.add_d(
+            self.WalkOperator.prepare, selection=soqs['selection'], **state_prep_ancilla
+        )
+        soqs['selection'] = prepare_soqs['selection']
+        del prepare_soqs['selection']
+        state_prep_ancilla = prepare_soqs
+
+        # GQSP
         soqs = bb.add_d(self.gqsp, **soqs)
-        # soqs |= bb.add_d(self.WalkOperator.prepare.adjoint(), selection=soqs['selection'])
+
+        # PREPARE†
+        prepare_soqs = bb.add_d(
+            self.WalkOperator.prepare.adjoint(), selection=soqs['selection'], **state_prep_ancilla
+        )
+        soqs['selection'] = prepare_soqs['selection']
+        del prepare_soqs['selection']
+        state_prep_ancilla = prepare_soqs
+
+        for soq in state_prep_ancilla.values():
+            bb.free(soq)
+
         return soqs

From bc6352104e750af992a1e0709b18355605533bac Mon Sep 17 00:00:00 2001
From: Anurudh Peduri <anurudhp97@gmail.com>
Date: Sun, 18 Feb 2024 21:59:58 +0100
Subject: [PATCH 09/34] add adjoint for random gate

---
 qualtran/bloqs/generalized_qsp_test.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/qualtran/bloqs/generalized_qsp_test.py b/qualtran/bloqs/generalized_qsp_test.py
index 847af852d..19a671542 100644
--- a/qualtran/bloqs/generalized_qsp_test.py
+++ b/qualtran/bloqs/generalized_qsp_test.py
@@ -128,6 +128,9 @@ def signature(self) -> Signature:
     def _unitary_(self):
         return self.matrix
 
+    def adjoint(self) -> GateWithRegisters:
+        return RandomGate(self.bitsize, self.matrix.conj().T)
+
 
 def evaluate_polynomial_of_matrix(
     P: Sequence[complex], U: NDArray, *, negative_power: int = 0

From 8ccb722924b4fc4aa360710cbb8b1a86784ff2b2 Mon Sep 17 00:00:00 2001
From: Anurudh Peduri <anurudhp97@gmail.com>
Date: Sun, 18 Feb 2024 22:18:57 +0100
Subject: [PATCH 10/34] fix negative power test

---
 qualtran/bloqs/generalized_qsp_test.py | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/qualtran/bloqs/generalized_qsp_test.py b/qualtran/bloqs/generalized_qsp_test.py
index 19a671542..a2cbb0696 100644
--- a/qualtran/bloqs/generalized_qsp_test.py
+++ b/qualtran/bloqs/generalized_qsp_test.py
@@ -137,7 +137,7 @@ def evaluate_polynomial_of_matrix(
 ) -> NDArray:
     assert U.ndim == 2 and U.shape[0] == U.shape[1]
 
-    pow_U = np.linalg.matrix_power(U.conjugate().T, negative_power)
+    pow_U = np.linalg.matrix_power(U.conj().T, negative_power)
     result = np.zeros(U.shape, dtype=U.dtype)
 
     for c in P:
@@ -167,11 +167,15 @@ def verify_generalized_qsp(
         gqsp_U = GeneralizedQSP(U, P, Q, negative_power=negative_power)
     result_unitary = cirq.unitary(gqsp_U)
 
-    expected_top_left = evaluate_polynomial_of_matrix(P, input_unitary)
+    expected_top_left = evaluate_polynomial_of_matrix(
+        P, input_unitary, negative_power=negative_power
+    )
     actual_top_left = result_unitary[:N, :N]
     assert_matrices_almost_equal(expected_top_left, actual_top_left)
 
-    expected_bottom_left = evaluate_polynomial_of_matrix(gqsp_U.Q, input_unitary)
+    expected_bottom_left = evaluate_polynomial_of_matrix(
+        gqsp_U.Q, input_unitary, negative_power=negative_power
+    )
     actual_bottom_left = result_unitary[N:, :N]
     assert_matrices_almost_equal(expected_bottom_left, actual_bottom_left)
 

From 962b7003c181d17180b2527f888f31eebf1ba105 Mon Sep 17 00:00:00 2001
From: Anurudh Peduri <anurudhp97@gmail.com>
Date: Sun, 18 Feb 2024 22:35:19 +0100
Subject: [PATCH 11/34] gqsp - fewer uses of U

---
 qualtran/bloqs/generalized_qsp.py | 26 ++++++++++++++------------
 1 file changed, 14 insertions(+), 12 deletions(-)

diff --git a/qualtran/bloqs/generalized_qsp.py b/qualtran/bloqs/generalized_qsp.py
index 1fd7a74f3..93a32784f 100644
--- a/qualtran/bloqs/generalized_qsp.py
+++ b/qualtran/bloqs/generalized_qsp.py
@@ -21,7 +21,7 @@
 from numpy.polynomial import Polynomial
 from numpy.typing import NDArray
 
-from qualtran import Bloq, GateWithRegisters, QBit, Register, Signature
+from qualtran import Bloq, GateWithRegisters, QBit, Register, Signature, CtrlSpec
 from qualtran.bloqs.basic_gates import Ry, ZPowGate
 from qualtran.bloqs.qubitization_walk_operator import QubitizationWalkOperator
 
@@ -330,14 +330,13 @@ def decompose_from_registers(
 
         yield self.signal_rotations[0].on(signal_qubit)
         for signal_rotation in self.signal_rotations[1:]:
-            # TODO debug controlled adjoint bloq serialization
-            # if num_inverse_applications > 0:
-            #     # apply C-U^\dagger
-            #     yield self.U.adjoint().on_registers(**quregs).controlled_by(signal_qubit)
-            #     num_inverse_applications -= 1
-            # else:
-            # apply C[0]-U
-            yield self.U.on_registers(**quregs).controlled_by(signal_qubit, control_values=[0])
+            if num_inverse_applications > 0:
+                # apply C-U^\dagger
+                yield self.U.adjoint().on_registers(**quregs).controlled_by(signal_qubit)
+                num_inverse_applications -= 1
+            else:
+                # apply C[0]-U
+                yield self.U.on_registers(**quregs).controlled_by(signal_qubit, control_values=[0])
             yield signal_rotation.on(signal_qubit)
 
         for _ in range(num_inverse_applications):
@@ -347,9 +346,12 @@ def __hash__(self):
         return hash((self.U, *self.signal_rotations, self.negative_power))
 
     def build_call_graph(self, ssa: 'SympySymbolAllocator') -> Set['BloqCountT']:
-        return {(self.U, max(len(self.P), self.negative_power))} | {
-            (rotation, 1) for rotation in self.signal_rotations
-        }
+        degree = len(self.P)
+        return {
+            (self.U.adjoint().controlled(), min(degree, self.negative_power)),
+            (self.U.controlled(control_values=[0]), max(0, degree - self.negative_power)),
+            (self.U.adjoint(), max(0, self.negative_power - degree)),
+        } | {(rotation, 1) for rotation in self.signal_rotations}
 
 
 @frozen

From fa9e4768f2126c4d25f51424fc5033e3146ee94c Mon Sep 17 00:00:00 2001
From: Anurudh Peduri <anurudhp97@gmail.com>
Date: Sun, 18 Feb 2024 22:46:00 +0100
Subject: [PATCH 12/34] hash random gate to support call graph

---
 qualtran/bloqs/generalized_qsp_test.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/qualtran/bloqs/generalized_qsp_test.py b/qualtran/bloqs/generalized_qsp_test.py
index a2cbb0696..a597da8c9 100644
--- a/qualtran/bloqs/generalized_qsp_test.py
+++ b/qualtran/bloqs/generalized_qsp_test.py
@@ -131,6 +131,9 @@ def _unitary_(self):
     def adjoint(self) -> GateWithRegisters:
         return RandomGate(self.bitsize, self.matrix.conj().T)
 
+    def __hash__(self):
+        return hash(tuple(np.ravel(self.matrix)))
+
 
 def evaluate_polynomial_of_matrix(
     P: Sequence[complex], U: NDArray, *, negative_power: int = 0

From 95dd66dd609af5d343c65829df4776ae026371ac Mon Sep 17 00:00:00 2001
From: Anurudh Peduri <anurudhp97@gmail.com>
Date: Sun, 18 Feb 2024 22:48:58 +0100
Subject: [PATCH 13/34] format

---
 qualtran/bloqs/generalized_qsp.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/qualtran/bloqs/generalized_qsp.py b/qualtran/bloqs/generalized_qsp.py
index 93a32784f..b73fe8805 100644
--- a/qualtran/bloqs/generalized_qsp.py
+++ b/qualtran/bloqs/generalized_qsp.py
@@ -21,7 +21,7 @@
 from numpy.polynomial import Polynomial
 from numpy.typing import NDArray
 
-from qualtran import Bloq, GateWithRegisters, QBit, Register, Signature, CtrlSpec
+from qualtran import Bloq, GateWithRegisters, QBit, Register, Signature
 from qualtran.bloqs.basic_gates import Ry, ZPowGate
 from qualtran.bloqs.qubitization_walk_operator import QubitizationWalkOperator
 

From a0ae7b09cc170729a7d0baa04c7148f12b168a93 Mon Sep 17 00:00:00 2001
From: Anurudh Peduri <anurudhp97@gmail.com>
Date: Tue, 20 Feb 2024 20:53:29 +0100
Subject: [PATCH 14/34] rename

---
 qualtran/bloqs/generalized_qsp.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/qualtran/bloqs/generalized_qsp.py b/qualtran/bloqs/generalized_qsp.py
index b73fe8805..5b0a21e52 100644
--- a/qualtran/bloqs/generalized_qsp.py
+++ b/qualtran/bloqs/generalized_qsp.py
@@ -356,7 +356,7 @@ def build_call_graph(self, ssa: 'SympySymbolAllocator') -> Set['BloqCountT']:
 
 @frozen
 class HamiltonianSimulationByGQSP(GateWithRegisters):
-    WalkOperator: QubitizationWalkOperator
+    walk_operator: QubitizationWalkOperator
     t: float = field(kw_only=True)
     alpha: float = field(kw_only=True)
     precision: float = field(kw_only=True)
@@ -373,7 +373,7 @@ def gqsp(self) -> GeneralizedQSP:
         approx_cos = 1j**coeff_indices * scipy.special.jv(coeff_indices, self.t * self.alpha)
 
         return GeneralizedQSP(
-            self.WalkOperator, approx_cos, np.zeros(2 * self.degree + 1), negative_power=self.degree
+            self.walk_operator, approx_cos, np.zeros(2 * self.degree + 1), negative_power=self.degree
         )
 
     @cached_property
@@ -381,14 +381,14 @@ def signature(self) -> 'Signature':
         return self.gqsp.signature
 
     def build_composite_bloq(self, bb: 'BloqBuilder', **soqs: 'SoquetT') -> Dict[str, 'SoquetT']:
-        prepare = self.WalkOperator.prepare
+        prepare = self.walk_operator.prepare
         state_prep_ancilla = {
             reg.name: bb.allocate(reg.total_bits()) for reg in prepare.junk_registers
         }
 
         # PREPARE
         prepare_soqs = bb.add_d(
-            self.WalkOperator.prepare, selection=soqs['selection'], **state_prep_ancilla
+            self.walk_operator.prepare, selection=soqs['selection'], **state_prep_ancilla
         )
         soqs['selection'] = prepare_soqs['selection']
         del prepare_soqs['selection']
@@ -399,7 +399,7 @@ def build_composite_bloq(self, bb: 'BloqBuilder', **soqs: 'SoquetT') -> Dict[str
 
         # PREPARE†
         prepare_soqs = bb.add_d(
-            self.WalkOperator.prepare.adjoint(), selection=soqs['selection'], **state_prep_ancilla
+            self.walk_operator.prepare.adjoint(), selection=soqs['selection'], **state_prep_ancilla
         )
         soqs['selection'] = prepare_soqs['selection']
         del prepare_soqs['selection']

From 595577b075673c11302eeb01c4bc59f8024bcf99 Mon Sep 17 00:00:00 2001
From: Anurudh Peduri <anurudhp97@gmail.com>
Date: Fri, 23 Feb 2024 21:30:42 +0100
Subject: [PATCH 15/34] wip call graph for ham. sim. by gqsp

---
 qualtran/bloqs/generalized_qsp.py      | 19 ++++++++++++++++---
 qualtran/bloqs/generalized_qsp_test.py |  5 ++++-
 2 files changed, 20 insertions(+), 4 deletions(-)

diff --git a/qualtran/bloqs/generalized_qsp.py b/qualtran/bloqs/generalized_qsp.py
index 5b0a21e52..780ce8a9d 100644
--- a/qualtran/bloqs/generalized_qsp.py
+++ b/qualtran/bloqs/generalized_qsp.py
@@ -17,11 +17,12 @@
 import cirq
 import numpy as np
 import scipy
+import sympy
 from attrs import field, frozen
 from numpy.polynomial import Polynomial
 from numpy.typing import NDArray
 
-from qualtran import Bloq, GateWithRegisters, QBit, Register, Signature
+from qualtran import Bloq, CtrlSpec, GateWithRegisters, QBit, Register, Signature
 from qualtran.bloqs.basic_gates import Ry, ZPowGate
 from qualtran.bloqs.qubitization_walk_operator import QubitizationWalkOperator
 
@@ -349,7 +350,7 @@ def build_call_graph(self, ssa: 'SympySymbolAllocator') -> Set['BloqCountT']:
         degree = len(self.P)
         return {
             (self.U.adjoint().controlled(), min(degree, self.negative_power)),
-            (self.U.controlled(control_values=[0]), max(0, degree - self.negative_power)),
+            (self.U.controlled(ctrl_spec=CtrlSpec(cvs=0)), max(0, degree - self.negative_power)),
             (self.U.adjoint(), max(0, self.negative_power - degree)),
         } | {(rotation, 1) for rotation in self.signal_rotations}
 
@@ -373,7 +374,10 @@ def gqsp(self) -> GeneralizedQSP:
         approx_cos = 1j**coeff_indices * scipy.special.jv(coeff_indices, self.t * self.alpha)
 
         return GeneralizedQSP(
-            self.walk_operator, approx_cos, np.zeros(2 * self.degree + 1), negative_power=self.degree
+            self.walk_operator,
+            approx_cos,
+            np.zeros(2 * self.degree + 1),
+            negative_power=self.degree,
         )
 
     @cached_property
@@ -409,3 +413,12 @@ def build_composite_bloq(self, bb: 'BloqBuilder', **soqs: 'SoquetT') -> Dict[str
             bb.free(soq)
 
         return soqs
+
+    def build_call_graph(self, ssa: 'SympySymbolAllocator') -> Set['BloqCountT']:
+        t = ssa.new_symbol('t')
+        alpha = ssa.new_symbol('alpha')
+        precision = ssa.new_symbol('precision')
+        d = t * alpha + sympy.log(precision) / sympy.log(sympy.log(precision))
+
+        # TODO account for SU2 rotation gates
+        return {(self.walk_operator, d)}
diff --git a/qualtran/bloqs/generalized_qsp_test.py b/qualtran/bloqs/generalized_qsp_test.py
index a597da8c9..0dee7c5fb 100644
--- a/qualtran/bloqs/generalized_qsp_test.py
+++ b/qualtran/bloqs/generalized_qsp_test.py
@@ -23,7 +23,7 @@
 from numpy.polynomial import Polynomial
 from numpy.typing import NDArray
 
-from qualtran import GateWithRegisters, Signature
+from qualtran import Bloq, GateWithRegisters, Signature
 from qualtran.bloqs.generalized_qsp import (
     GeneralizedQSP,
     qsp_complementary_polynomial,
@@ -131,6 +131,9 @@ def _unitary_(self):
     def adjoint(self) -> GateWithRegisters:
         return RandomGate(self.bitsize, self.matrix.conj().T)
 
+    def controlled(self, *args, **kwargs):
+        return Bloq.controlled(self, *args, **kwargs)
+
     def __hash__(self):
         return hash(tuple(np.ravel(self.matrix)))
 

From 07a784e09b71ac9832709b1937afcb80425a4f45 Mon Sep 17 00:00:00 2001
From: Anurudh Peduri <anurudhp97@gmail.com>
Date: Sun, 25 Feb 2024 23:05:25 +0100
Subject: [PATCH 16/34] add tests

---
 qualtran/_infra/adjoint.py             |   2 +-
 qualtran/bloqs/generalized_qsp_test.py | 158 ++++++++++++++++++++++++-
 2 files changed, 157 insertions(+), 3 deletions(-)

diff --git a/qualtran/_infra/adjoint.py b/qualtran/_infra/adjoint.py
index ca790c005..825b31478 100644
--- a/qualtran/_infra/adjoint.py
+++ b/qualtran/_infra/adjoint.py
@@ -147,7 +147,7 @@ def decompose_from_registers(
         self, *, context: cirq.DecompositionContext, **quregs: NDArray[cirq.Qid]
     ) -> cirq.OP_TREE:
         if isinstance(self.subbloq, GateWithRegisters):
-            return cirq.inverse(self.subbloq.decompose_from_registers(context=context, **quregs))
+            return cirq.inverse(*self.subbloq.decompose_from_registers(context=context, **quregs))
         return super().decompose_from_registers(context=context, **quregs)
 
     def supports_decompose_bloq(self) -> bool:
diff --git a/qualtran/bloqs/generalized_qsp_test.py b/qualtran/bloqs/generalized_qsp_test.py
index 0dee7c5fb..9a98abf0f 100644
--- a/qualtran/bloqs/generalized_qsp_test.py
+++ b/qualtran/bloqs/generalized_qsp_test.py
@@ -12,24 +12,38 @@
 #  See the License for the specific language governing permissions and
 #  limitations under the License.
 from functools import cached_property
-from typing import Optional, Sequence, Union
+from typing import Optional, Sequence, Tuple, Union
 
+import attrs
 import cirq
 import numpy as np
 import pytest
+import scipy
 import sympy
 from attrs import define, frozen
 from cirq.testing import random_unitary
 from numpy.polynomial import Polynomial
 from numpy.typing import NDArray
 
-from qualtran import Bloq, GateWithRegisters, Signature
+from qualtran import (
+    Bloq,
+    BloqBuilder,
+    BoundedQUInt,
+    GateWithRegisters,
+    QBit,
+    Register,
+    Signature,
+    SoquetT,
+)
 from qualtran.bloqs.generalized_qsp import (
     GeneralizedQSP,
+    HamiltonianSimulationByGQSP,
     qsp_complementary_polynomial,
     qsp_phase_factors,
     SU2RotationGate,
 )
+from qualtran.bloqs.qubitization_walk_operator import QubitizationWalkOperator
+from qualtran.bloqs.select_and_prepare import PrepareOracle, SelectOracle
 from qualtran.cirq_interop import BloqAsCirqGate
 
 
@@ -279,3 +293,143 @@ def test_generalized_real_qsp_with_symbolic_signal_matrix(degree: int):
     for _ in range(10):
         P = random_qsp_polynomial(degree, random_state=random_state)
         SymbolicGQSP(P).verify()
+
+
+@frozen
+class RandomPrepareOracle(PrepareOracle):
+    U: RandomGate
+
+    @property
+    def selection_registers(self) -> tuple[Register, ...]:
+        return (Register('selection', BoundedQUInt(bitsize=self.U.bitsize)),)
+
+    @staticmethod
+    def create(bitsize: int, *, random_state: np.random.RandomState):
+        matrix = random_unitary(2**bitsize, random_state=random_state)
+        alpha = matrix[:, 0]
+        matrix = matrix * (alpha.conj() / np.abs(alpha))[:, None]
+        return RandomPrepareOracle(RandomGate(bitsize, matrix))
+
+    def build_composite_bloq(self, bb: BloqBuilder, selection: SoquetT) -> dict[str, SoquetT]:
+        selection = bb.add(self.U, q=selection)
+        return {'selection': selection}
+
+    def __pow__(self, power):
+        if power == -1:
+            return self.U.adjoint()
+        return NotImplemented
+
+    @cached_property
+    def alphas(self):
+        np.testing.assert_almost_equal(np.imag(self.U.matrix[:, 0]), 0)
+        return np.sqrt(self.U.matrix[:, 0])
+
+
+@frozen
+class PauliSelectOracle(SelectOracle):
+    select_bitsize: int
+    target_bitsize: int
+    select_unitaries: tuple[cirq.DensePauliString, ...]
+    control_val: Optional[int] = None
+
+    @property
+    def control_registers(self) -> Tuple[Register, ...]:
+        return () if self.control_val is None else (Register('control', QBit()),)
+
+    @property
+    def selection_registers(self) -> Tuple[Register, ...]:
+        return (Register('selection', BoundedQUInt(bitsize=self.select_bitsize)),)
+
+    @property
+    def target_registers(self) -> Tuple[Register, ...]:
+        return (Register('target', BoundedQUInt(bitsize=self.target_bitsize)),)
+
+    def adjoint(self):
+        return self
+
+    def __pow__(self, power):
+        if abs(power) == 1:
+            return self
+        return NotImplemented
+
+    def controlled(
+        self,
+        num_controls: Optional[int] = None,
+        control_values=None,
+        control_qid_shape: Optional[Tuple[int, ...]] = None,
+    ) -> 'cirq.Gate':
+        if num_controls is None:
+            num_controls = 1
+        if control_values is None:
+            control_values = [1] * num_controls
+        if (
+            isinstance(control_values, Sequence)
+            and isinstance(control_values[0], int)
+            and len(control_values) == 1
+            and self.control_val is None
+        ):
+            return attrs.evolve(self, control_val=control_values[0])
+        raise NotImplementedError()
+
+    def decompose_from_registers(
+        self,
+        *,
+        context: cirq.DecompositionContext,
+        selection: NDArray[cirq.Qid],
+        target: NDArray[cirq.Qid],
+        **quregs: NDArray[cirq.Qid],
+    ) -> cirq.OP_TREE:
+        if self.control_val is not None:
+            selection = np.concatenate([selection, quregs['control']])
+
+        for cv, U in enumerate(self.select_unitaries):
+            bits = tuple(map(int, bin(cv)[2:].zfill(self.select_bitsize)))
+            if self.control_val is not None:
+                bits = (*bits, self.control_val)
+            yield U.on(*target).controlled_by(*selection, control_values=bits)
+
+
+def random_qubitization_walk_operator(
+    select_bitsize: int, target_bitsize: int, *, random_state: np.random.RandomState
+) -> tuple[QubitizationWalkOperator, cirq.PauliSum]:
+    prepare = RandomPrepareOracle.create(select_bitsize, random_state=random_state)
+
+    dps = tuple(
+        cirq.DensePauliString(random_state.random_integers(0, 3, size=target_bitsize))
+        for _ in range(2**select_bitsize)
+    )
+    select = PauliSelectOracle(select_bitsize, target_bitsize, dps)
+
+    ham = cirq.PauliSum.from_pauli_strings(
+        [
+            dp.on(*cirq.LineQubit.range(target_bitsize)) * alpha
+            for dp, alpha in zip(dps, prepare.alphas)
+        ]
+    )
+
+    return QubitizationWalkOperator(prepare=prepare, select=select), ham
+
+
+def verify_hamiltonian_simulation_by_gqsp(
+    W: QubitizationWalkOperator, H: NDArray[np.complex_], *, t: float, alpha: float = 1
+):
+    N = H.shape[0]
+    W_e_iHt = HamiltonianSimulationByGQSP(W, t=t, alpha=alpha, precision=1e-10)
+    result_unitary = cirq.unitary(W_e_iHt)
+
+    expected_top_left = scipy.linalg.expm(1j * H * t)
+    actual_top_left = result_unitary[:N, :N]
+    assert_matrices_almost_equal(expected_top_left, actual_top_left)
+
+
+@pytest.mark.slow
+@pytest.mark.parametrize("select_bitsize", [1, 2, 3])
+@pytest.mark.parametrize("target_bitsize", [1, 2, 3])
+def test_hamiltonian_simulation_by_gqsp(select_bitsize: int, target_bitsize: int):
+    random_state = np.random.RandomState(42)
+
+    for _ in range(5):
+        W, H = random_qubitization_walk_operator(
+            select_bitsize, target_bitsize, random_state=random_state
+        )
+        verify_hamiltonian_simulation_by_gqsp(W, H.matrix(), t=5, alpha=1)

From 57416ce0dc3f989aad4898c54009407210850160 Mon Sep 17 00:00:00 2001
From: Anurudh Peduri <anurudhp97@gmail.com>
Date: Mon, 26 Feb 2024 22:47:41 +0100
Subject: [PATCH 17/34] fix bug in using adjoint - use pow(-1) for GWR

---
 qualtran/bloqs/generalized_qsp.py      |  2 +-
 qualtran/bloqs/generalized_qsp_test.py | 11 +++++++----
 2 files changed, 8 insertions(+), 5 deletions(-)

diff --git a/qualtran/bloqs/generalized_qsp.py b/qualtran/bloqs/generalized_qsp.py
index 780ce8a9d..24269611f 100644
--- a/qualtran/bloqs/generalized_qsp.py
+++ b/qualtran/bloqs/generalized_qsp.py
@@ -333,7 +333,7 @@ def decompose_from_registers(
         for signal_rotation in self.signal_rotations[1:]:
             if num_inverse_applications > 0:
                 # apply C-U^\dagger
-                yield self.U.adjoint().on_registers(**quregs).controlled_by(signal_qubit)
+                yield (self.U**-1).on_registers(**quregs).controlled_by(signal_qubit)
                 num_inverse_applications -= 1
             else:
                 # apply C[0]-U
diff --git a/qualtran/bloqs/generalized_qsp_test.py b/qualtran/bloqs/generalized_qsp_test.py
index 9a98abf0f..437d9374d 100644
--- a/qualtran/bloqs/generalized_qsp_test.py
+++ b/qualtran/bloqs/generalized_qsp_test.py
@@ -26,7 +26,6 @@
 from numpy.typing import NDArray
 
 from qualtran import (
-    Bloq,
     BloqBuilder,
     BoundedQUInt,
     GateWithRegisters,
@@ -145,9 +144,6 @@ def _unitary_(self):
     def adjoint(self) -> GateWithRegisters:
         return RandomGate(self.bitsize, self.matrix.conj().T)
 
-    def controlled(self, *args, **kwargs):
-        return Bloq.controlled(self, *args, **kwargs)
-
     def __hash__(self):
         return hash(tuple(np.ravel(self.matrix)))
 
@@ -306,8 +302,15 @@ def selection_registers(self) -> tuple[Register, ...]:
     @staticmethod
     def create(bitsize: int, *, random_state: np.random.RandomState):
         matrix = random_unitary(2**bitsize, random_state=random_state)
+
+        # make the first column (weights alpha_i) all reals
         alpha = matrix[:, 0]
         matrix = matrix * (alpha.conj() / np.abs(alpha))[:, None]
+
+        # verify that it is still unitary
+        np.testing.assert_allclose(matrix @ matrix.conj().T, np.eye(2**bitsize), atol=1e-10)
+        np.testing.assert_allclose(matrix.conj().T @ matrix, np.eye(2**bitsize), atol=1e-10)
+
         return RandomPrepareOracle(RandomGate(bitsize, matrix))
 
     def build_composite_bloq(self, bb: BloqBuilder, selection: SoquetT) -> dict[str, SoquetT]:

From b3e9cef95f419f138b38399ea363968c4121d403 Mon Sep 17 00:00:00 2001
From: Anurudh Peduri <anurudhp97@gmail.com>
Date: Mon, 26 Feb 2024 22:49:39 +0100
Subject: [PATCH 18/34] cleanup notebook

---
 qualtran/bloqs/generalized_qsp.ipynb | 14 ++------------
 1 file changed, 2 insertions(+), 12 deletions(-)

diff --git a/qualtran/bloqs/generalized_qsp.ipynb b/qualtran/bloqs/generalized_qsp.ipynb
index 70bcd2839..d051cba3b 100644
--- a/qualtran/bloqs/generalized_qsp.ipynb
+++ b/qualtran/bloqs/generalized_qsp.ipynb
@@ -150,17 +150,7 @@
    "outputs": [],
    "source": [
     "W_e_iHt = HamiltonianSimulationByGQSP(U, t=5, alpha=1, precision=1e-3)\n",
-    "\n",
-    "def decomp_pred(inst):\n",
-    "    from qualtran.bloqs.select_and_prepare import PrepareOracle, SelectOracle\n",
-    "\n",
-    "    return any(isinstance(inst.bloq, cls) for cls in [HamiltonianSimulationByGQSP, GeneralizedQSP, SelectOracle, PrepareOracle])\n",
-    "    \n",
-    "show_bloq(\n",
-    "    W_e_iHt\n",
-    "    .decompose_bloq()\n",
-    "    .flatten(pred=decomp_pred)\n",
-    ")"
+    "show_bloq(W_e_iHt.decompose_bloq())"
    ]
   },
   {
@@ -190,7 +180,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.6"
+   "version": "3.11.7"
   }
  },
  "nbformat": 4,

From 94dac1779686576e70c6258d231a97cd4de45d9f Mon Sep 17 00:00:00 2001
From: Anurudh Peduri <anurudhp97@gmail.com>
Date: Mon, 26 Feb 2024 23:06:25 +0100
Subject: [PATCH 19/34] fix expression for alpha

---
 qualtran/bloqs/generalized_qsp_test.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/qualtran/bloqs/generalized_qsp_test.py b/qualtran/bloqs/generalized_qsp_test.py
index 437d9374d..265005149 100644
--- a/qualtran/bloqs/generalized_qsp_test.py
+++ b/qualtran/bloqs/generalized_qsp_test.py
@@ -325,7 +325,7 @@ def __pow__(self, power):
     @cached_property
     def alphas(self):
         np.testing.assert_almost_equal(np.imag(self.U.matrix[:, 0]), 0)
-        return np.sqrt(self.U.matrix[:, 0])
+        return self.U.matrix[:, 0] ** 2
 
 
 @frozen
@@ -403,6 +403,8 @@ def random_qubitization_walk_operator(
     )
     select = PauliSelectOracle(select_bitsize, target_bitsize, dps)
 
+    np.testing.assert_allclose(np.linalg.norm(prepare.alphas, 1), 1)
+
     ham = cirq.PauliSum.from_pauli_strings(
         [
             dp.on(*cirq.LineQubit.range(target_bitsize)) * alpha

From 671fe13ae6740b3b89730e7a9873649f1854260a Mon Sep 17 00:00:00 2001
From: Anurudh Peduri <anurudhp97@gmail.com>
Date: Mon, 26 Feb 2024 23:06:34 +0100
Subject: [PATCH 20/34] cleanup bloq construction

---
 qualtran/bloqs/generalized_qsp.py | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/qualtran/bloqs/generalized_qsp.py b/qualtran/bloqs/generalized_qsp.py
index 24269611f..87f326ec0 100644
--- a/qualtran/bloqs/generalized_qsp.py
+++ b/qualtran/bloqs/generalized_qsp.py
@@ -394,8 +394,7 @@ def build_composite_bloq(self, bb: 'BloqBuilder', **soqs: 'SoquetT') -> Dict[str
         prepare_soqs = bb.add_d(
             self.walk_operator.prepare, selection=soqs['selection'], **state_prep_ancilla
         )
-        soqs['selection'] = prepare_soqs['selection']
-        del prepare_soqs['selection']
+        soqs['selection'] = prepare_soqs.pop('selection')
         state_prep_ancilla = prepare_soqs
 
         # GQSP
@@ -405,8 +404,7 @@ def build_composite_bloq(self, bb: 'BloqBuilder', **soqs: 'SoquetT') -> Dict[str
         prepare_soqs = bb.add_d(
             self.walk_operator.prepare.adjoint(), selection=soqs['selection'], **state_prep_ancilla
         )
-        soqs['selection'] = prepare_soqs['selection']
-        del prepare_soqs['selection']
+        soqs['selection'] = prepare_soqs.pop('selection')
         state_prep_ancilla = prepare_soqs
 
         for soq in state_prep_ancilla.values():

From f576a5df5e3868047e44b3cbcba8cc63d5d4fe76 Mon Sep 17 00:00:00 2001
From: Anurudh Peduri <anurudhp97@gmail.com>
Date: Mon, 26 Feb 2024 23:17:34 +0100
Subject: [PATCH 21/34] implement `**-1` for RandomGate

---
 qualtran/bloqs/generalized_qsp_test.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/qualtran/bloqs/generalized_qsp_test.py b/qualtran/bloqs/generalized_qsp_test.py
index 265005149..98a1a0034 100644
--- a/qualtran/bloqs/generalized_qsp_test.py
+++ b/qualtran/bloqs/generalized_qsp_test.py
@@ -144,6 +144,11 @@ def _unitary_(self):
     def adjoint(self) -> GateWithRegisters:
         return RandomGate(self.bitsize, self.matrix.conj().T)
 
+    def __pow__(self, power):
+        if power == -1:
+            return self.adjoint()
+        return NotImplemented
+
     def __hash__(self):
         return hash(tuple(np.ravel(self.matrix)))
 

From 53bd323068c8e37e8eee5ea00267783c5881e63e Mon Sep 17 00:00:00 2001
From: Anurudh Peduri <anurudhp97@gmail.com>
Date: Tue, 27 Feb 2024 11:20:20 +0100
Subject: [PATCH 22/34] cleanup unused fix

---
 qualtran/_infra/adjoint.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/qualtran/_infra/adjoint.py b/qualtran/_infra/adjoint.py
index 825b31478..ca790c005 100644
--- a/qualtran/_infra/adjoint.py
+++ b/qualtran/_infra/adjoint.py
@@ -147,7 +147,7 @@ def decompose_from_registers(
         self, *, context: cirq.DecompositionContext, **quregs: NDArray[cirq.Qid]
     ) -> cirq.OP_TREE:
         if isinstance(self.subbloq, GateWithRegisters):
-            return cirq.inverse(*self.subbloq.decompose_from_registers(context=context, **quregs))
+            return cirq.inverse(self.subbloq.decompose_from_registers(context=context, **quregs))
         return super().decompose_from_registers(context=context, **quregs)
 
     def supports_decompose_bloq(self) -> bool:

From 52a980b5a10c244ca3bf1b3c9dbd74012e7609fe Mon Sep 17 00:00:00 2001
From: Anurudh Peduri <anurudhp97@gmail.com>
Date: Tue, 27 Feb 2024 11:20:47 +0100
Subject: [PATCH 23/34] add test for `cos` approximation

---
 qualtran/bloqs/generalized_qsp.py      | 13 ++++++----
 qualtran/bloqs/generalized_qsp_test.py | 34 +++++++++++++++++++++++---
 2 files changed, 38 insertions(+), 9 deletions(-)

diff --git a/qualtran/bloqs/generalized_qsp.py b/qualtran/bloqs/generalized_qsp.py
index 87f326ec0..b2240bc56 100644
--- a/qualtran/bloqs/generalized_qsp.py
+++ b/qualtran/bloqs/generalized_qsp.py
@@ -290,9 +290,9 @@ class GeneralizedQSP(GateWithRegisters):
     def signature(self) -> Signature:
         return Signature([Register('signal', QBit()), *self.U.signature])
 
-    @staticmethod
+    @classmethod
     def from_qsp_polynomial(
-        U: GateWithRegisters, P: Sequence[complex], *, negative_power: int = 0
+        cls, U: GateWithRegisters, P: Sequence[complex], *, negative_power: int = 0
     ) -> 'GeneralizedQSP':
         Q = qsp_complementary_polynomial(P)
         return GeneralizedQSP(U, P, Q, negative_power=negative_power)
@@ -365,17 +365,20 @@ class HamiltonianSimulationByGQSP(GateWithRegisters):
     @cached_property
     def degree(self) -> int:
         log_precision = np.log(1 / self.precision)
-        degree = self.t * self.alpha + log_precision / np.log(log_precision)
+        degree = self.t * self.alpha + 3 * log_precision / np.log(log_precision)
         return int(np.ceil(degree))
 
     @cached_property
-    def gqsp(self) -> GeneralizedQSP:
+    def _approx_cos(self) -> NDArray[np.complex_]:
         coeff_indices = np.arange(-self.degree, self.degree + 1)
         approx_cos = 1j**coeff_indices * scipy.special.jv(coeff_indices, self.t * self.alpha)
+        return approx_cos
 
+    @cached_property
+    def gqsp(self) -> GeneralizedQSP:
         return GeneralizedQSP(
             self.walk_operator,
-            approx_cos,
+            self._approx_cos,
             np.zeros(2 * self.degree + 1),
             negative_power=self.degree,
         )
diff --git a/qualtran/bloqs/generalized_qsp_test.py b/qualtran/bloqs/generalized_qsp_test.py
index 98a1a0034..0c9f9813a 100644
--- a/qualtran/bloqs/generalized_qsp_test.py
+++ b/qualtran/bloqs/generalized_qsp_test.py
@@ -73,6 +73,7 @@ def check_polynomial_pair_on_random_points_on_unit_circle(
     Q: Union[Sequence[complex], Polynomial],
     *,
     random_state: np.random.RandomState,
+    rtol: float = 1e-7,
     n_points: int = 1000,
 ):
     P = Polynomial(P)
@@ -80,7 +81,7 @@ def check_polynomial_pair_on_random_points_on_unit_circle(
 
     for _ in range(n_points):
         z = np.exp(random_state.random() * np.pi * 2j)
-        np.testing.assert_allclose(np.abs(P(z)) ** 2 + np.abs(Q(z)) ** 2, 1)
+        np.testing.assert_allclose(np.abs(P(z)) ** 2 + np.abs(Q(z)) ** 2, 1, rtol=rtol)
 
 
 def random_qsp_polynomial(
@@ -391,7 +392,7 @@ def decompose_from_registers(
             selection = np.concatenate([selection, quregs['control']])
 
         for cv, U in enumerate(self.select_unitaries):
-            bits = tuple(map(int, bin(cv)[2:].zfill(self.select_bitsize)))
+            bits = tuple(map(int, bin(cv)[2:].zfill(self.select_bitsize)))[::-1]
             if self.control_val is not None:
                 bits = (*bits, self.control_val)
             yield U.on(*target).controlled_by(*selection, control_values=bits)
@@ -432,9 +433,34 @@ def verify_hamiltonian_simulation_by_gqsp(
     assert_matrices_almost_equal(expected_top_left, actual_top_left)
 
 
+@pytest.mark.parametrize("precision", [1e-5, 1e-10])
+def test_cos_approximation_fast(precision: float):
+    random_state = np.random.RandomState(42)
+
+    for t in [1, 2, 3]:
+        for alpha in [0.5, 1]:
+            bloq = HamiltonianSimulationByGQSP(None, t=t, alpha=alpha, precision=precision)
+            check_polynomial_pair_on_random_points_on_unit_circle(
+                bloq._approx_cos, [0], random_state=random_state, rtol=precision
+            )
+
+
 @pytest.mark.slow
-@pytest.mark.parametrize("select_bitsize", [1, 2, 3])
-@pytest.mark.parametrize("target_bitsize", [1, 2, 3])
+@pytest.mark.parametrize("precision", [1e-5, 1e-7, 1e-10])
+def test_cos_approximation(precision):
+    random_state = np.random.RandomState(42)
+
+    for t in random_state.random(10):
+        for alpha in random_state.random(5):
+            bloq = HamiltonianSimulationByGQSP(None, t=t * 10, alpha=alpha, precision=precision)
+            check_polynomial_pair_on_random_points_on_unit_circle(
+                bloq._approx_cos, [0], random_state=random_state, rtol=precision
+            )
+
+
+# @pytest.mark.slow
+@pytest.mark.parametrize("select_bitsize", [1])
+@pytest.mark.parametrize("target_bitsize", [1])
 def test_hamiltonian_simulation_by_gqsp(select_bitsize: int, target_bitsize: int):
     random_state = np.random.RandomState(42)
 

From b2e5a6e2acd94e26c3808b12b4db5c0e136e0956 Mon Sep 17 00:00:00 2001
From: Anurudh Peduri <anurudhp97@gmail.com>
Date: Tue, 27 Feb 2024 11:28:38 +0100
Subject: [PATCH 24/34] docstring

---
 qualtran/bloqs/generalized_qsp.py | 22 +++++++++++++++++++++-
 1 file changed, 21 insertions(+), 1 deletion(-)

diff --git a/qualtran/bloqs/generalized_qsp.py b/qualtran/bloqs/generalized_qsp.py
index b2240bc56..fbe9e4a86 100644
--- a/qualtran/bloqs/generalized_qsp.py
+++ b/qualtran/bloqs/generalized_qsp.py
@@ -357,6 +357,24 @@ def build_call_graph(self, ssa: 'SympySymbolAllocator') -> Set['BloqCountT']:
 
 @frozen
 class HamiltonianSimulationByGQSP(GateWithRegisters):
+    r"""Hamiltonian simulation using Generalized QSP given a qubitized quantum walk operator.
+
+    Implements Hamiltonian simulation given a walk operator from SELECT and PREPARE oracles.
+
+    We can use the Jacobi-Anger expansion to obtain low-degree polynomial approximations for the $\cos$ function:
+
+        $$ e^{it\cos\theta} = \sum_{n = -\infty}^{\infty} i^n J_n(t) (e^{i\theta})^n $$
+
+    where $J_n$ is the $n$-th [Bessel function of the first kind](https://en.wikipedia.org/wiki/Bessel_function#Bessel_functions_of_the_first_kind), which is provided by `scipy.special.jv`.
+    We can cutoff at $d = O(t + \log(1/\epsilon) / \log\log(1/\epsilon))$ to get an $\epsilon$-approximation (Theorem 7):
+
+        $$ P[t](z) = \sum_{n = -d}^d i^n J_n(t) z^n $$
+
+    References:
+        [Generalized Quantum Signal Processing](https://arxiv.org/abs/2308.01501)
+            Motlagh and Wiebe. (2023). Theorem 7.
+    """
+
     walk_operator: QubitizationWalkOperator
     t: float = field(kw_only=True)
     alpha: float = field(kw_only=True)
@@ -364,12 +382,14 @@ class HamiltonianSimulationByGQSP(GateWithRegisters):
 
     @cached_property
     def degree(self) -> int:
+        r"""degree of the polynomial to approximate the function e^{it\cos(\theta)}"""
         log_precision = np.log(1 / self.precision)
         degree = self.t * self.alpha + 3 * log_precision / np.log(log_precision)
         return int(np.ceil(degree))
 
     @cached_property
     def _approx_cos(self) -> NDArray[np.complex_]:
+        r"""polynomial approximation for e^{it\cos(\theta)} in terms of monomials (e^{i\theta})^n"""
         coeff_indices = np.arange(-self.degree, self.degree + 1)
         approx_cos = 1j**coeff_indices * scipy.special.jv(coeff_indices, self.t * self.alpha)
         return approx_cos
@@ -419,7 +439,7 @@ def build_call_graph(self, ssa: 'SympySymbolAllocator') -> Set['BloqCountT']:
         t = ssa.new_symbol('t')
         alpha = ssa.new_symbol('alpha')
         precision = ssa.new_symbol('precision')
-        d = t * alpha + sympy.log(precision) / sympy.log(sympy.log(precision))
+        d = t * alpha + 3 * sympy.log(precision) / sympy.log(sympy.log(precision))
 
         # TODO account for SU2 rotation gates
         return {(self.walk_operator, d)}

From a1afb1fadafc256913bcc4314cfa4c4e8d116728 Mon Sep 17 00:00:00 2001
From: Anurudh Peduri <anurudhp97@gmail.com>
Date: Tue, 27 Feb 2024 12:06:18 +0100
Subject: [PATCH 25/34] use adjoint

---
 qualtran/bloqs/generalized_qsp.ipynb | 2 +-
 qualtran/bloqs/generalized_qsp.py    | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/qualtran/bloqs/generalized_qsp.ipynb b/qualtran/bloqs/generalized_qsp.ipynb
index d051cba3b..3402c1287 100644
--- a/qualtran/bloqs/generalized_qsp.ipynb
+++ b/qualtran/bloqs/generalized_qsp.ipynb
@@ -105,7 +105,7 @@
    "outputs": [],
    "source": [
     "pU = GeneralizedQSP.from_qsp_polynomial(U, (0.5, 0, 0.5), negative_power=1)\n",
-    "show_bloq(pU.decompose_bloq())"
+    "show_bloq(pU)"
    ]
   },
   {
diff --git a/qualtran/bloqs/generalized_qsp.py b/qualtran/bloqs/generalized_qsp.py
index fbe9e4a86..d551d0f95 100644
--- a/qualtran/bloqs/generalized_qsp.py
+++ b/qualtran/bloqs/generalized_qsp.py
@@ -333,7 +333,7 @@ def decompose_from_registers(
         for signal_rotation in self.signal_rotations[1:]:
             if num_inverse_applications > 0:
                 # apply C-U^\dagger
-                yield (self.U**-1).on_registers(**quregs).controlled_by(signal_qubit)
+                yield self.U.adjoint().on_registers(**quregs).controlled_by(signal_qubit)
                 num_inverse_applications -= 1
             else:
                 # apply C[0]-U

From f9e425d338ab0f073de6c5d03a736871301e49f6 Mon Sep 17 00:00:00 2001
From: Anurudh Peduri <anurudhp97@gmail.com>
Date: Tue, 27 Feb 2024 17:04:38 +0100
Subject: [PATCH 26/34] reorganize code

---
 qualtran/bloqs/generalized_qsp_test.py | 50 +++++++++++++-------------
 1 file changed, 25 insertions(+), 25 deletions(-)

diff --git a/qualtran/bloqs/generalized_qsp_test.py b/qualtran/bloqs/generalized_qsp_test.py
index 0c9f9813a..5038a6095 100644
--- a/qualtran/bloqs/generalized_qsp_test.py
+++ b/qualtran/bloqs/generalized_qsp_test.py
@@ -297,6 +297,31 @@ def test_generalized_real_qsp_with_symbolic_signal_matrix(degree: int):
         SymbolicGQSP(P).verify()
 
 
+@pytest.mark.parametrize("precision", [1e-5, 1e-10])
+def test_cos_approximation_fast(precision: float):
+    random_state = np.random.RandomState(42)
+
+    for t in [1, 2, 3]:
+        for alpha in [0.5, 1]:
+            bloq = HamiltonianSimulationByGQSP(None, t=t, alpha=alpha, precision=precision)
+            check_polynomial_pair_on_random_points_on_unit_circle(
+                bloq._approx_cos, [0], random_state=random_state, rtol=precision
+            )
+
+
+@pytest.mark.slow
+@pytest.mark.parametrize("precision", [1e-5, 1e-7, 1e-10])
+def test_cos_approximation(precision):
+    random_state = np.random.RandomState(42)
+
+    for t in random_state.random(10):
+        for alpha in random_state.random(5):
+            bloq = HamiltonianSimulationByGQSP(None, t=t * 10, alpha=alpha, precision=precision)
+            check_polynomial_pair_on_random_points_on_unit_circle(
+                bloq._approx_cos, [0], random_state=random_state, rtol=precision
+            )
+
+
 @frozen
 class RandomPrepareOracle(PrepareOracle):
     U: RandomGate
@@ -433,31 +458,6 @@ def verify_hamiltonian_simulation_by_gqsp(
     assert_matrices_almost_equal(expected_top_left, actual_top_left)
 
 
-@pytest.mark.parametrize("precision", [1e-5, 1e-10])
-def test_cos_approximation_fast(precision: float):
-    random_state = np.random.RandomState(42)
-
-    for t in [1, 2, 3]:
-        for alpha in [0.5, 1]:
-            bloq = HamiltonianSimulationByGQSP(None, t=t, alpha=alpha, precision=precision)
-            check_polynomial_pair_on_random_points_on_unit_circle(
-                bloq._approx_cos, [0], random_state=random_state, rtol=precision
-            )
-
-
-@pytest.mark.slow
-@pytest.mark.parametrize("precision", [1e-5, 1e-7, 1e-10])
-def test_cos_approximation(precision):
-    random_state = np.random.RandomState(42)
-
-    for t in random_state.random(10):
-        for alpha in random_state.random(5):
-            bloq = HamiltonianSimulationByGQSP(None, t=t * 10, alpha=alpha, precision=precision)
-            check_polynomial_pair_on_random_points_on_unit_circle(
-                bloq._approx_cos, [0], random_state=random_state, rtol=precision
-            )
-
-
 # @pytest.mark.slow
 @pytest.mark.parametrize("select_bitsize", [1])
 @pytest.mark.parametrize("target_bitsize", [1])

From f4b1ddf13c2e5e535e25b7f4b94537bc2e1ce0fe Mon Sep 17 00:00:00 2001
From: Anurudh Peduri <anurudhp97@gmail.com>
Date: Tue, 27 Feb 2024 17:21:02 +0100
Subject: [PATCH 27/34] fix warning

---
 qualtran/bloqs/generalized_qsp_test.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/qualtran/bloqs/generalized_qsp_test.py b/qualtran/bloqs/generalized_qsp_test.py
index 5038a6095..ce851d332 100644
--- a/qualtran/bloqs/generalized_qsp_test.py
+++ b/qualtran/bloqs/generalized_qsp_test.py
@@ -429,7 +429,7 @@ def random_qubitization_walk_operator(
     prepare = RandomPrepareOracle.create(select_bitsize, random_state=random_state)
 
     dps = tuple(
-        cirq.DensePauliString(random_state.random_integers(0, 3, size=target_bitsize))
+        cirq.DensePauliString(random_state.randint(0, 4, size=target_bitsize))
         for _ in range(2**select_bitsize)
     )
     select = PauliSelectOracle(select_bitsize, target_bitsize, dps)

From 830a2fd2d96d5d29f306b4c14a2a23caa7a0f990 Mon Sep 17 00:00:00 2001
From: Anurudh Peduri <anurudhp97@gmail.com>
Date: Mon, 4 Mar 2024 12:44:25 +0100
Subject: [PATCH 28/34] compute exp-cos approx accurately

---
 qualtran/bloqs/generalized_qsp.py      | 35 ++++++++++++------
 qualtran/bloqs/generalized_qsp_test.py | 49 +++++++++++++++++---------
 2 files changed, 56 insertions(+), 28 deletions(-)

diff --git a/qualtran/bloqs/generalized_qsp.py b/qualtran/bloqs/generalized_qsp.py
index d551d0f95..aac4dd8c1 100644
--- a/qualtran/bloqs/generalized_qsp.py
+++ b/qualtran/bloqs/generalized_qsp.py
@@ -237,7 +237,7 @@ def qsp_phase_factors(
     lambd = 0
 
     def safe_angle(x):
-        return 0 if np.isclose(x, 0) else np.angle(x)
+        return 0 if np.isclose(x, 0, atol=1e-10) else np.angle(x)
 
     for d in reversed(range(n)):
         assert S.shape == (2, d + 1)
@@ -245,7 +245,7 @@ def safe_angle(x):
         a, b = S[:, d]
         theta[d] = np.arctan2(np.abs(b), np.abs(a))
         # \phi_d = arg(a / b)
-        phi[d] = 0 if np.isclose(np.abs(b), 0) else safe_angle(a * np.conj(b))
+        phi[d] = 0 if np.isclose(np.abs(b), 0, atol=1e-10) else safe_angle(a * np.conj(b))
 
         if d == 0:
             lambd = safe_angle(b)
@@ -383,23 +383,31 @@ class HamiltonianSimulationByGQSP(GateWithRegisters):
     @cached_property
     def degree(self) -> int:
         r"""degree of the polynomial to approximate the function e^{it\cos(\theta)}"""
-        log_precision = np.log(1 / self.precision)
-        degree = self.t * self.alpha + 3 * log_precision / np.log(log_precision)
-        return int(np.ceil(degree))
+        return len(self.approx_cos) // 2
 
     @cached_property
-    def _approx_cos(self) -> NDArray[np.complex_]:
-        r"""polynomial approximation for e^{it\cos(\theta)} in terms of monomials (e^{i\theta})^n"""
-        coeff_indices = np.arange(-self.degree, self.degree + 1)
+    def approx_cos(self) -> NDArray[np.complex_]:
+        r"""polynomial approximation for $$e^{i\theta} \mapsto e^{it\cos(\theta)}$$"""
+        d = 0
+        while True:
+            term = scipy.special.jv(d + 1, self.t * self.alpha)
+            if np.isclose(term, 0, atol=self.precision / 2):
+                break
+            d += 1
+
+        coeff_indices = np.arange(-d, d + 1)
         approx_cos = 1j**coeff_indices * scipy.special.jv(coeff_indices, self.t * self.alpha)
         return approx_cos
 
     @cached_property
     def gqsp(self) -> GeneralizedQSP:
+        # return GeneralizedQSP.from_qsp_polynomial(
+        #     self.walk_operator, self.approx_cos / np.sqrt(2), negative_power=self.degree
+        # )
         return GeneralizedQSP(
             self.walk_operator,
-            self._approx_cos,
-            np.zeros(2 * self.degree + 1),
+            self.approx_cos / np.sqrt(2),
+            self.approx_cos / np.sqrt(2),
             negative_power=self.degree,
         )
 
@@ -439,7 +447,12 @@ def build_call_graph(self, ssa: 'SympySymbolAllocator') -> Set['BloqCountT']:
         t = ssa.new_symbol('t')
         alpha = ssa.new_symbol('alpha')
         precision = ssa.new_symbol('precision')
-        d = t * alpha + 3 * sympy.log(precision) / sympy.log(sympy.log(precision))
+        d = sympy.O(
+            t * alpha + sympy.log(precision) / sympy.log(sympy.log(precision)),
+            (t, sympy.oo),
+            (alpha, sympy.oo),
+            (precision, 0),
+        )
 
         # TODO account for SU2 rotation gates
         return {(self.walk_operator, d)}
diff --git a/qualtran/bloqs/generalized_qsp_test.py b/qualtran/bloqs/generalized_qsp_test.py
index ce851d332..e98aa12b0 100644
--- a/qualtran/bloqs/generalized_qsp_test.py
+++ b/qualtran/bloqs/generalized_qsp_test.py
@@ -297,29 +297,43 @@ def test_generalized_real_qsp_with_symbolic_signal_matrix(degree: int):
         SymbolicGQSP(P).verify()
 
 
-@pytest.mark.parametrize("precision", [1e-5, 1e-10])
-def test_cos_approximation_fast(precision: float):
+@pytest.mark.parametrize("precision", [1e-5, 1e-7, 1e-10])
+def test_cos_approximation(precision: float):
     random_state = np.random.RandomState(42)
 
-    for t in [1, 2, 3]:
-        for alpha in [0.5, 1]:
+    for t in [2, 3, 5, 10]:
+        for alpha in [0.5, 1, 2]:
             bloq = HamiltonianSimulationByGQSP(None, t=t, alpha=alpha, precision=precision)
-            check_polynomial_pair_on_random_points_on_unit_circle(
-                bloq._approx_cos, [0], random_state=random_state, rtol=precision
+
+            P = Polynomial(bloq.approx_cos)
+            theta = 2 * np.pi * random_state.random(1000)
+            e_itheta = np.exp(1j * theta)
+            np.testing.assert_allclose(
+                P(e_itheta) * e_itheta ** (-bloq.degree),
+                np.exp(1j * t * alpha * np.cos(theta)),
+                rtol=precision * 2,
             )
 
 
 @pytest.mark.slow
-@pytest.mark.parametrize("precision", [1e-5, 1e-7, 1e-10])
-def test_cos_approximation(precision):
+@pytest.mark.parametrize("bitsize", [1, 2])
+@pytest.mark.parametrize("t", [2, 3, 5])
+@pytest.mark.parametrize("alpha", [1, 2, 3])
+@pytest.mark.parametrize("precision", [1e-7, 1e-10])
+def test_generalized_qsp_with_exp_cos_approx_on_random_unitaries(
+    bitsize: int, t: float, alpha: float, precision: float
+):
     random_state = np.random.RandomState(42)
 
-    for t in random_state.random(10):
-        for alpha in random_state.random(5):
-            bloq = HamiltonianSimulationByGQSP(None, t=t * 10, alpha=alpha, precision=precision)
-            check_polynomial_pair_on_random_points_on_unit_circle(
-                bloq._approx_cos, [0], random_state=random_state, rtol=precision
-            )
+    for _ in range(5):
+        U = RandomGate.create(bitsize, random_state=random_state)
+        P = HamiltonianSimulationByGQSP(
+            None, t=t, alpha=alpha, precision=precision
+        ).approx_cos / np.sqrt(2)
+        check_polynomial_pair_on_random_points_on_unit_circle(
+            P, P, random_state=random_state, rtol=2 * precision
+        )
+        verify_generalized_qsp(U, P, P, negative_power=len(P) // 2)
 
 
 @frozen
@@ -450,15 +464,16 @@ def verify_hamiltonian_simulation_by_gqsp(
     W: QubitizationWalkOperator, H: NDArray[np.complex_], *, t: float, alpha: float = 1
 ):
     N = H.shape[0]
-    W_e_iHt = HamiltonianSimulationByGQSP(W, t=t, alpha=alpha, precision=1e-10)
+
+    W_e_iHt = HamiltonianSimulationByGQSP(W, t=t, alpha=alpha, precision=1e-7)
     result_unitary = cirq.unitary(W_e_iHt)
 
-    expected_top_left = scipy.linalg.expm(1j * H * t)
+    expected_top_left = scipy.linalg.expm(-1j * H * t)
     actual_top_left = result_unitary[:N, :N]
     assert_matrices_almost_equal(expected_top_left, actual_top_left)
 
 
-# @pytest.mark.slow
+@pytest.mark.slow
 @pytest.mark.parametrize("select_bitsize", [1])
 @pytest.mark.parametrize("target_bitsize", [1])
 def test_hamiltonian_simulation_by_gqsp(select_bitsize: int, target_bitsize: int):

From 916b6746be7fc311157aca394c0d22f9bb77a4f8 Mon Sep 17 00:00:00 2001
From: Anurudh Peduri <anurudhp97@gmail.com>
Date: Mon, 4 Mar 2024 14:20:11 +0100
Subject: [PATCH 29/34] fix sympy: all vars need to have same limit

---
 qualtran/bloqs/generalized_qsp.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/qualtran/bloqs/generalized_qsp.py b/qualtran/bloqs/generalized_qsp.py
index aac4dd8c1..91ce3205c 100644
--- a/qualtran/bloqs/generalized_qsp.py
+++ b/qualtran/bloqs/generalized_qsp.py
@@ -446,12 +446,12 @@ def build_composite_bloq(self, bb: 'BloqBuilder', **soqs: 'SoquetT') -> Dict[str
     def build_call_graph(self, ssa: 'SympySymbolAllocator') -> Set['BloqCountT']:
         t = ssa.new_symbol('t')
         alpha = ssa.new_symbol('alpha')
-        precision = ssa.new_symbol('precision')
+        inv_precision = ssa.new_symbol('1/precision')
         d = sympy.O(
-            t * alpha + sympy.log(precision) / sympy.log(sympy.log(precision)),
+            t * alpha + sympy.log(1 / inv_precision) / sympy.log(sympy.log(1 / inv_precision)),
             (t, sympy.oo),
             (alpha, sympy.oo),
-            (precision, 0),
+            (inv_precision, sympy.oo),
         )
 
         # TODO account for SU2 rotation gates

From dc0f2865452e967c1fae9858c9ddd4d6e8ee921e Mon Sep 17 00:00:00 2001
From: Anurudh Peduri <anurudhp97@gmail.com>
Date: Mon, 4 Mar 2024 14:26:51 +0100
Subject: [PATCH 30/34] mark slow test

---
 qualtran/bloqs/generalized_qsp_test.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/qualtran/bloqs/generalized_qsp_test.py b/qualtran/bloqs/generalized_qsp_test.py
index e98aa12b0..3feee5e1f 100644
--- a/qualtran/bloqs/generalized_qsp_test.py
+++ b/qualtran/bloqs/generalized_qsp_test.py
@@ -114,6 +114,7 @@ def test_complementary_polynomial(degree: int):
         check_polynomial_pair_on_random_points_on_unit_circle(P, Q, random_state=random_state)
 
 
+@pytest.mark.slow
 @pytest.mark.parametrize("degree", [3, 4, 5, 10, 20, 30, 100])
 def test_real_polynomial_has_real_complementary_polynomial(degree: int):
     random_state = np.random.RandomState(42)

From 6f8d34b1b351fe463467085825e6e10bb51065e4 Mon Sep 17 00:00:00 2001
From: Anurudh Peduri <anurudhp97@gmail.com>
Date: Tue, 5 Mar 2024 15:03:05 +0100
Subject: [PATCH 31/34] update ham-sim poly test

---
 qualtran/bloqs/generalized_qsp_test.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/qualtran/bloqs/generalized_qsp_test.py b/qualtran/bloqs/generalized_qsp_test.py
index 3feee5e1f..c34b4ed04 100644
--- a/qualtran/bloqs/generalized_qsp_test.py
+++ b/qualtran/bloqs/generalized_qsp_test.py
@@ -328,13 +328,13 @@ def test_generalized_qsp_with_exp_cos_approx_on_random_unitaries(
 
     for _ in range(5):
         U = RandomGate.create(bitsize, random_state=random_state)
-        P = HamiltonianSimulationByGQSP(
-            None, t=t, alpha=alpha, precision=precision
-        ).approx_cos / np.sqrt(2)
+        gqsp = HamiltonianSimulationByGQSP(None, t=t, alpha=alpha, precision=precision).gqsp
+        P, Q = gqsp.P, gqsp.Q
+
         check_polynomial_pair_on_random_points_on_unit_circle(
-            P, P, random_state=random_state, rtol=2 * precision
+            P, Q, random_state=random_state, rtol=2 * precision
         )
-        verify_generalized_qsp(U, P, P, negative_power=len(P) // 2)
+        verify_generalized_qsp(U, P, Q, negative_power=len(P) // 2)
 
 
 @frozen

From 8b332149ebd2fd712adc14ee51d7bbd78f7e8eba Mon Sep 17 00:00:00 2001
From: Anurudh Peduri <anurudhp97@gmail.com>
Date: Tue, 5 Mar 2024 15:08:18 +0100
Subject: [PATCH 32/34] use exact poly as in paper

---
 qualtran/bloqs/generalized_qsp.py | 10 ++--------
 1 file changed, 2 insertions(+), 8 deletions(-)

diff --git a/qualtran/bloqs/generalized_qsp.py b/qualtran/bloqs/generalized_qsp.py
index 91ce3205c..8af535ede 100644
--- a/qualtran/bloqs/generalized_qsp.py
+++ b/qualtran/bloqs/generalized_qsp.py
@@ -401,14 +401,8 @@ def approx_cos(self) -> NDArray[np.complex_]:
 
     @cached_property
     def gqsp(self) -> GeneralizedQSP:
-        # return GeneralizedQSP.from_qsp_polynomial(
-        #     self.walk_operator, self.approx_cos / np.sqrt(2), negative_power=self.degree
-        # )
-        return GeneralizedQSP(
-            self.walk_operator,
-            self.approx_cos / np.sqrt(2),
-            self.approx_cos / np.sqrt(2),
-            negative_power=self.degree,
+        return GeneralizedQSP.from_qsp_polynomial(
+            self.walk_operator, self.approx_cos, negative_power=self.degree
         )
 
     @cached_property

From 026df56cc5c803cc9276e18627f2367a96656cd0 Mon Sep 17 00:00:00 2001
From: Anurudh Peduri <anurudhp97@gmail.com>
Date: Wed, 6 Mar 2024 17:07:17 +0100
Subject: [PATCH 33/34] cleanup code

---
 qualtran/bloqs/generalized_qsp.py | 11 +++++------
 1 file changed, 5 insertions(+), 6 deletions(-)

diff --git a/qualtran/bloqs/generalized_qsp.py b/qualtran/bloqs/generalized_qsp.py
index 8af535ede..95a8305a6 100644
--- a/qualtran/bloqs/generalized_qsp.py
+++ b/qualtran/bloqs/generalized_qsp.py
@@ -383,19 +383,18 @@ class HamiltonianSimulationByGQSP(GateWithRegisters):
     @cached_property
     def degree(self) -> int:
         r"""degree of the polynomial to approximate the function e^{it\cos(\theta)}"""
-        return len(self.approx_cos) // 2
-
-    @cached_property
-    def approx_cos(self) -> NDArray[np.complex_]:
-        r"""polynomial approximation for $$e^{i\theta} \mapsto e^{it\cos(\theta)}$$"""
         d = 0
         while True:
             term = scipy.special.jv(d + 1, self.t * self.alpha)
             if np.isclose(term, 0, atol=self.precision / 2):
                 break
             d += 1
+        return d
 
-        coeff_indices = np.arange(-d, d + 1)
+    @cached_property
+    def approx_cos(self) -> NDArray[np.complex_]:
+        r"""polynomial approximation for $$e^{i\theta} \mapsto e^{it\cos(\theta)}$$"""
+        coeff_indices = np.arange(-self.degree, self.degree + 1)
         approx_cos = 1j**coeff_indices * scipy.special.jv(coeff_indices, self.t * self.alpha)
         return approx_cos
 

From fc01bc914fd276bcc4b433c11de383867c2482ea Mon Sep 17 00:00:00 2001
From: Anurudh Peduri <anurudhp97@gmail.com>
Date: Fri, 15 Mar 2024 09:39:01 +0100
Subject: [PATCH 34/34] remove custom hash

---
 qualtran/bloqs/generalized_qsp.py | 7 ++-----
 1 file changed, 2 insertions(+), 5 deletions(-)

diff --git a/qualtran/bloqs/generalized_qsp.py b/qualtran/bloqs/generalized_qsp.py
index 95a8305a6..0fe5f2a96 100644
--- a/qualtran/bloqs/generalized_qsp.py
+++ b/qualtran/bloqs/generalized_qsp.py
@@ -282,8 +282,8 @@ class GeneralizedQSP(GateWithRegisters):
     """
 
     U: GateWithRegisters
-    P: Sequence[complex]
-    Q: Sequence[complex]
+    P: Tuple[complex, ...] = field(converter=tuple)
+    Q: Tuple[complex, ...] = field(converter=tuple)
     negative_power: int = field(default=0, kw_only=True)
 
     @cached_property
@@ -343,9 +343,6 @@ def decompose_from_registers(
         for _ in range(num_inverse_applications):
             yield self.U.adjoint().on_registers(**quregs)
 
-    def __hash__(self):
-        return hash((self.U, *self.signal_rotations, self.negative_power))
-
     def build_call_graph(self, ssa: 'SympySymbolAllocator') -> Set['BloqCountT']:
         degree = len(self.P)
         return {