From 9c7bacc70526a1ed4d3e0a35058fd49551a68f3d Mon Sep 17 00:00:00 2001
From: Charles Yuan <charlesyuan@google.com>
Date: Tue, 2 Jul 2024 15:06:09 -0700
Subject: [PATCH] Add `Product` block encoding

---
 dev_tools/autogenerate-bloqs-notebooks-v2.py  |   1 +
 qualtran/bloqs/block_encoding/__init__.py     |   1 +
 .../bloqs/block_encoding/block_encoding.ipynb | 159 +++++++++++
 qualtran/bloqs/block_encoding/product.py      | 264 ++++++++++++++++++
 qualtran/bloqs/block_encoding/product_test.py | 143 ++++++++++
 qualtran/conftest.py                          |   3 +
 qualtran/serialization/resolver_dict.py       |   2 +
 qualtran/symbolics/math_funcs.py              |  24 +-
 8 files changed, 596 insertions(+), 1 deletion(-)
 create mode 100644 qualtran/bloqs/block_encoding/product.py
 create mode 100644 qualtran/bloqs/block_encoding/product_test.py

diff --git a/dev_tools/autogenerate-bloqs-notebooks-v2.py b/dev_tools/autogenerate-bloqs-notebooks-v2.py
index 4358844cc5..3e1d4af185 100644
--- a/dev_tools/autogenerate-bloqs-notebooks-v2.py
+++ b/dev_tools/autogenerate-bloqs-notebooks-v2.py
@@ -572,6 +572,7 @@
             qualtran.bloqs.block_encoding.chebyshev_polynomial._CHEBYSHEV_BLOQ_DOC,
             qualtran.bloqs.block_encoding.unitary._UNITARY_DOC,
             qualtran.bloqs.block_encoding.tensor_product._TENSOR_PRODUCT_DOC,
+            qualtran.bloqs.block_encoding.product._PRODUCT_DOC,
         ],
         directory=f'{SOURCE_DIR}/bloqs/block_encoding/',
     ),
diff --git a/qualtran/bloqs/block_encoding/__init__.py b/qualtran/bloqs/block_encoding/__init__.py
index 7b05237906..ae05c6c9c3 100644
--- a/qualtran/bloqs/block_encoding/__init__.py
+++ b/qualtran/bloqs/block_encoding/__init__.py
@@ -21,5 +21,6 @@
     LCUBlockEncoding,
     LCUBlockEncodingZeroState,
 )
+from qualtran.bloqs.block_encoding.product import Product
 from qualtran.bloqs.block_encoding.tensor_product import TensorProduct
 from qualtran.bloqs.block_encoding.unitary import Unitary
diff --git a/qualtran/bloqs/block_encoding/block_encoding.ipynb b/qualtran/bloqs/block_encoding/block_encoding.ipynb
index d5400fcbc5..a3228c1c91 100644
--- a/qualtran/bloqs/block_encoding/block_encoding.ipynb
+++ b/qualtran/bloqs/block_encoding/block_encoding.ipynb
@@ -890,6 +890,165 @@
     "show_call_graph(tensor_product_block_encoding_g)\n",
     "show_counts_sigma(tensor_product_block_encoding_sigma)"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "1127108b",
+   "metadata": {
+    "cq.autogen": "Product.bloq_doc.md"
+   },
+   "source": [
+    "## `Product`\n",
+    "Product of a sequence of block encodings.\n",
+    "\n",
+    "Builds the block encoding $B[U_1 * U_2 * \\cdots * U_n]$ given block encodings $B[U_1], \\ldots, B[U_n]$.\n",
+    "\n",
+    "When each $B[U_i]$ is a $(\\alpha_i, a_i, \\epsilon_i)$-block encoding of $U_i$, we have that\n",
+    "$B[U_1 * \\cdots * U_n]$ is a $(\\prod_i \\alpha_i, n - 1 + \\max_i a_i, \\sum_i \\alpha_i \\epsilon_i)$-block\n",
+    "encoding of $U_1 * \\cdots * U_n$.\n",
+    "\n",
+    "#### Parameters\n",
+    " - `block_encodings`: A sequence of block encodings. \n",
+    "\n",
+    "#### Registers\n",
+    " - `system`: The system register.\n",
+    " - `ancilla`: The ancilla register (present only if bitsize > 0).\n",
+    " - `resource`: The resource register (present only if bitsize > 0). \n",
+    "\n",
+    "#### References\n",
+    " - [Quantum algorithms: A survey of applications and end-to-end complexities](https://arxiv.org/abs/2310.03011). Dalzell et al. (2023). Ch. 10.2.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3643f471",
+   "metadata": {
+    "cq.autogen": "Product.bloq_doc.py"
+   },
+   "outputs": [],
+   "source": [
+    "from qualtran.bloqs.block_encoding import Product"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d37c3d81",
+   "metadata": {
+    "cq.autogen": "Product.example_instances.md"
+   },
+   "source": [
+    "### Example Instances"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "96a8d277",
+   "metadata": {
+    "cq.autogen": "Product.product_block_encoding"
+   },
+   "outputs": [],
+   "source": [
+    "from qualtran.bloqs.basic_gates import Hadamard, TGate\n",
+    "from qualtran.bloqs.block_encoding.unitary import Unitary\n",
+    "\n",
+    "product_block_encoding = Product((Unitary(TGate()), Unitary(Hadamard())))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "83b0468a",
+   "metadata": {
+    "cq.autogen": "Product.product_block_encoding_override"
+   },
+   "outputs": [],
+   "source": [
+    "from qualtran.bloqs.basic_gates import Hadamard, TGate\n",
+    "from qualtran.bloqs.block_encoding.unitary import Unitary\n",
+    "\n",
+    "u1 = evolve(Unitary(TGate()), alpha=0.5, ancilla_bitsize=2, resource_bitsize=1, epsilon=0.01)\n",
+    "u2 = evolve(Unitary(Hadamard()), alpha=0.5, ancilla_bitsize=1, resource_bitsize=1, epsilon=0.1)\n",
+    "product_block_encoding_override = Product((u1, u2))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "238b45e4",
+   "metadata": {
+    "cq.autogen": "Product.product_block_encoding_symb"
+   },
+   "outputs": [],
+   "source": [
+    "import sympy\n",
+    "\n",
+    "from qualtran.bloqs.basic_gates import Hadamard, TGate\n",
+    "from qualtran.bloqs.block_encoding.unitary import Unitary\n",
+    "\n",
+    "alpha1 = sympy.Symbol('alpha1')\n",
+    "a1 = sympy.Symbol('a1')\n",
+    "eps1 = sympy.Symbol('eps1')\n",
+    "alpha2 = sympy.Symbol('alpha2')\n",
+    "a2 = sympy.Symbol('a2')\n",
+    "eps2 = sympy.Symbol('eps2')\n",
+    "product_block_encoding_symb = Product(\n",
+    "    (\n",
+    "        Unitary(TGate(), alpha=alpha1, ancilla_bitsize=a1, epsilon=eps1),\n",
+    "        Unitary(Hadamard(), alpha=alpha2, ancilla_bitsize=a2, epsilon=eps2),\n",
+    "    )\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "598fe397",
+   "metadata": {
+    "cq.autogen": "Product.graphical_signature.md"
+   },
+   "source": [
+    "#### Graphical Signature"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5e28152f",
+   "metadata": {
+    "cq.autogen": "Product.graphical_signature.py"
+   },
+   "outputs": [],
+   "source": [
+    "from qualtran.drawing import show_bloqs\n",
+    "show_bloqs([product_block_encoding, product_block_encoding_override, product_block_encoding_symb],\n",
+    "           ['`product_block_encoding`', '`product_block_encoding_override`', '`product_block_encoding_symb`'])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d3175473",
+   "metadata": {
+    "cq.autogen": "Product.call_graph.md"
+   },
+   "source": [
+    "### Call Graph"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a5436f3e",
+   "metadata": {
+    "cq.autogen": "Product.call_graph.py"
+   },
+   "outputs": [],
+   "source": [
+    "from qualtran.resource_counting.generalizers import ignore_split_join\n",
+    "product_block_encoding_g, product_block_encoding_sigma = product_block_encoding.call_graph(max_depth=1, generalizer=ignore_split_join)\n",
+    "show_call_graph(product_block_encoding_g)\n",
+    "show_counts_sigma(product_block_encoding_sigma)"
+   ]
   }
  ],
  "metadata": {
diff --git a/qualtran/bloqs/block_encoding/product.py b/qualtran/bloqs/block_encoding/product.py
new file mode 100644
index 0000000000..ab2b5d3f96
--- /dev/null
+++ b/qualtran/bloqs/block_encoding/product.py
@@ -0,0 +1,264 @@
+#  Copyright 2023 Google LLC
+#
+#  Licensed under the Apache License, Version 2.0 (the "License");
+#  you may not use this file except in compliance with the License.
+#  You may obtain a copy of the License at
+#
+#      https://www.apache.org/licenses/LICENSE-2.0
+#
+#  Unless required by applicable law or agreed to in writing, software
+#  distributed under the License is distributed on an "AS IS" BASIS,
+#  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#  See the License for the specific language governing permissions and
+#  limitations under the License.
+
+from functools import cached_property
+from typing import cast, Dict, Tuple
+
+from attrs import evolve, field, frozen, validators
+from numpy.typing import NDArray
+
+from qualtran import (
+    bloq_example,
+    BloqBuilder,
+    BloqDocSpec,
+    CtrlSpec,
+    DecomposeTypeError,
+    QAny,
+    QBit,
+    Register,
+    Signature,
+    Soquet,
+    SoquetT,
+)
+from qualtran.bloqs.basic_gates.x_basis import XGate
+from qualtran.bloqs.block_encoding import BlockEncoding
+from qualtran.bloqs.block_encoding.lcu_select_and_prepare import PrepareOracle
+from qualtran.bloqs.bookkeeping.partition import Partition
+from qualtran.symbolics import is_symbolic, prod, smax, ssum, SymbolicFloat, SymbolicInt
+
+
+@frozen
+class Product(BlockEncoding):
+    r"""Product of a sequence of block encodings.
+
+    Builds the block encoding $B[U_1 * U_2 * \cdots * U_n]$ given block encodings $B[U_1], \ldots, B[U_n]$.
+
+    When each $B[U_i]$ is a $(\alpha_i, a_i, \epsilon_i)$-block encoding of $U_i$, we have that
+    $B[U_1 * \cdots * U_n]$ is a $(\prod_i \alpha_i, n - 1 + \max_i a_i, \sum_i \alpha_i \epsilon_i)$-block
+    encoding of $U_1 * \cdots * U_n$.
+
+    Args:
+        block_encodings: A sequence of block encodings.
+
+    Registers:
+        system: The system register.
+        ancilla: The ancilla register (present only if bitsize > 0).
+        resource: The resource register (present only if bitsize > 0).
+
+    References:
+        [Quantum algorithms: A survey of applications and end-to-end complexities](https://arxiv.org/abs/2310.03011). Dalzell et al. (2023). Ch. 10.2.
+    """
+
+    block_encodings: Tuple[BlockEncoding, ...] = field(
+        converter=lambda x: x if isinstance(x, tuple) else tuple(x), validator=validators.min_len(1)
+    )
+
+    def __attrs_post_init__(self):
+        if not all(u.system_bitsize == self.system_bitsize for u in self.block_encodings):
+            raise ValueError("All block encodings must have the same system size.")
+
+    @cached_property
+    def signature(self) -> Signature:
+        return Signature.build_from_dtypes(
+            system=QAny(self.system_bitsize),
+            ancilla=QAny(self.ancilla_bitsize),  # if ancilla_bitsize is 0, not present
+            resource=QAny(self.resource_bitsize),  # if resource_bitsize is 0, not present
+        )
+
+    @cached_property
+    def system_bitsize(self) -> SymbolicInt:
+        return self.block_encodings[0].system_bitsize
+
+    def pretty_name(self) -> str:
+        return f"B[{'*'.join(u.pretty_name()[2:-1] for u in self.block_encodings)}]"
+
+    @cached_property
+    def alpha(self) -> SymbolicFloat:
+        return prod(u.alpha for u in self.block_encodings)
+
+    @cached_property
+    def ancilla_bitsize(self) -> SymbolicInt:
+        return smax(u.ancilla_bitsize for u in self.block_encodings) + len(self.block_encodings) - 1
+
+    @cached_property
+    def resource_bitsize(self) -> SymbolicInt:
+        return ssum(u.resource_bitsize for u in self.block_encodings)
+
+    @cached_property
+    def epsilon(self) -> SymbolicFloat:
+        return ssum(u.alpha * u.epsilon for u in self.block_encodings)
+
+    @property
+    def target_registers(self) -> Tuple[Register, ...]:
+        return (self.signature.get_right("system"),)
+
+    @property
+    def junk_registers(self) -> Tuple[Register, ...]:
+        return (self.signature.get_right("resource"),)
+
+    @property
+    def selection_registers(self) -> Tuple[Register, ...]:
+        return (self.signature.get_right("ancilla"),)
+
+    @property
+    def signal_state(self) -> PrepareOracle:
+        """This method will be implemented in the future after PrepareOracle is updated for the BlockEncoding interface."""
+        raise NotImplementedError
+
+    def build_composite_bloq(
+        self, bb: BloqBuilder, system: SoquetT, **soqs: SoquetT
+    ) -> Dict[str, SoquetT]:
+        if (
+            is_symbolic(self.system_bitsize)
+            or is_symbolic(self.ancilla_bitsize)
+            or is_symbolic(self.resource_bitsize)
+        ):
+            raise DecomposeTypeError(f"Cannot decompose symbolic {self=}")
+        assert (
+            isinstance(self.system_bitsize, int)
+            and isinstance(self.ancilla_bitsize, int)
+            and isinstance(self.resource_bitsize, int)
+        )
+
+        n = len(self.block_encodings)
+        res_bits_used = 0
+        for i, u in enumerate(reversed(self.block_encodings)):
+            u_soqs = {"system": system}
+
+            # split ancilla register if necessary
+            if self.ancilla_bitsize > 0:
+                anc_bits = cast(int, u.ancilla_bitsize)
+                flag_bits = Register("flag_bits", dtype=QBit(), shape=(n - 1,))  # type: ignore
+                anc_used = Register("anc_used", dtype=QAny(anc_bits))
+                anc_unused_bits = self.ancilla_bitsize - (n - 1) - anc_bits
+                anc_unused = Register("anc_unused", dtype=QAny(anc_unused_bits))
+                anc_regs = [flag_bits]
+                if anc_bits > 0:
+                    anc_regs.append(anc_used)
+                if anc_unused_bits > 0:
+                    anc_regs.append(anc_unused)
+                anc_part = Partition(self.ancilla_bitsize, tuple(anc_regs))
+                anc_part_soqs = bb.add_d(anc_part, x=soqs["ancilla"])
+                if anc_bits > 0:
+                    u_soqs["ancilla"] = anc_part_soqs["anc_used"]
+
+            # split resource register if necessary
+            res_bits = cast(int, u.resource_bitsize)
+            if res_bits > 0:
+                res_before = Register("res_before", dtype=QAny(res_bits_used))
+                res = Register("res", dtype=QAny(res_bits))
+                res_bits_left = self.resource_bitsize - res_bits_used - res_bits
+                res_after = Register("res_after", dtype=QAny(res_bits_left))
+                res_regs = []
+                if res_bits_used > 0:
+                    res_regs.append(res_before)
+                res_regs.append(res)
+                res_bits_used += res_bits
+                if res_bits_left > 0:
+                    res_regs.append(res_after)
+                res_part = Partition(self.resource_bitsize, tuple(res_regs))
+                res_part_soqs = bb.add_d(res_part, x=soqs["resource"])
+                u_soqs["resource"] = res_part_soqs["res"]
+
+            # connect the constituent bloq
+            u_out_soqs = bb.add_d(u, **u_soqs)
+            system = u_out_soqs["system"]
+
+            # un-partition the resource register
+            if res_bits > 0:
+                res_part_soqs["res"] = u_out_soqs["resource"]
+                soqs["resource"] = cast(
+                    Soquet, bb.add(evolve(res_part, partition=False), **res_part_soqs)
+                )
+
+            # un-partition the ancilla register
+            if self.ancilla_bitsize > 0:
+                flag_bits_soq = cast(NDArray, anc_part_soqs["flag_bits"])
+                if anc_bits > 0:
+                    anc_used_soq = cast(Soquet, u_out_soqs["ancilla"])
+                    if i == n - 1:
+                        anc_part_soqs["anc_used"] = anc_used_soq
+                    else:
+                        # set corresponding flag if ancillas are all zero
+                        ctrl, flag_bits_soq[i] = bb.add_t(
+                            XGate().controlled(CtrlSpec(qdtypes=(anc_used_soq.reg.dtype,), cvs=0)),
+                            ctrl=anc_used_soq,
+                            q=flag_bits_soq[i],
+                        )
+                        flag_bits_soq[i] = bb.add(XGate(), q=flag_bits_soq[i])
+                        anc_part_soqs["anc_used"] = cast(Soquet, ctrl)
+                anc_part_soqs["flag_bits"] = flag_bits_soq
+                soqs["ancilla"] = cast(
+                    Soquet, bb.add(evolve(anc_part, partition=False), **anc_part_soqs)
+                )
+
+        out = {"system": system}
+        if self.ancilla_bitsize > 0:
+            out["ancilla"] = soqs["ancilla"]
+        if self.resource_bitsize > 0:
+            out["resource"] = soqs["resource"]
+        return out
+
+
+@bloq_example
+def _product_block_encoding() -> Product:
+    from qualtran.bloqs.basic_gates import Hadamard, TGate
+    from qualtran.bloqs.block_encoding.unitary import Unitary
+
+    product_block_encoding = Product((Unitary(TGate()), Unitary(Hadamard())))
+    return product_block_encoding
+
+
+@bloq_example
+def _product_block_encoding_override() -> Product:
+    from qualtran.bloqs.basic_gates import Hadamard, TGate
+    from qualtran.bloqs.block_encoding.unitary import Unitary
+
+    u1 = evolve(Unitary(TGate()), alpha=0.5, ancilla_bitsize=2, resource_bitsize=1, epsilon=0.01)
+    u2 = evolve(Unitary(Hadamard()), alpha=0.5, ancilla_bitsize=1, resource_bitsize=1, epsilon=0.1)
+    product_block_encoding_override = Product((u1, u2))
+    return product_block_encoding_override
+
+
+@bloq_example
+def _product_block_encoding_symb() -> Product:
+    import sympy
+
+    from qualtran.bloqs.basic_gates import Hadamard, TGate
+    from qualtran.bloqs.block_encoding.unitary import Unitary
+
+    alpha1 = sympy.Symbol('alpha1')
+    a1 = sympy.Symbol('a1')
+    eps1 = sympy.Symbol('eps1')
+    alpha2 = sympy.Symbol('alpha2')
+    a2 = sympy.Symbol('a2')
+    eps2 = sympy.Symbol('eps2')
+    product_block_encoding_symb = Product(
+        (
+            Unitary(TGate(), alpha=alpha1, ancilla_bitsize=a1, epsilon=eps1),
+            Unitary(Hadamard(), alpha=alpha2, ancilla_bitsize=a2, epsilon=eps2),
+        )
+    )
+    return product_block_encoding_symb
+
+
+_PRODUCT_DOC = BloqDocSpec(
+    bloq_cls=Product,
+    import_line="from qualtran.bloqs.block_encoding import Product",
+    examples=[
+        _product_block_encoding,
+        _product_block_encoding_override,
+        _product_block_encoding_symb,
+    ],
+)
diff --git a/qualtran/bloqs/block_encoding/product_test.py b/qualtran/bloqs/block_encoding/product_test.py
new file mode 100644
index 0000000000..9fe36754d2
--- /dev/null
+++ b/qualtran/bloqs/block_encoding/product_test.py
@@ -0,0 +1,143 @@
+#  Copyright 2023 Google LLC
+#
+#  Licensed under the Apache License, Version 2.0 (the "License");
+#  you may not use this file except in compliance with the License.
+#  You may obtain a copy of the License at
+#
+#      https://www.apache.org/licenses/LICENSE-2.0
+#
+#  Unless required by applicable law or agreed to in writing, software
+#  distributed under the License is distributed on an "AS IS" BASIS,
+#  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#  See the License for the specific language governing permissions and
+#  limitations under the License.
+
+from typing import cast
+
+import cirq
+import numpy as np
+import pytest
+import sympy
+
+from qualtran import BloqBuilder, QAny, Register, Signature, Soquet
+from qualtran.bloqs.basic_gates import CNOT, Hadamard, TGate, XGate, ZeroEffect, ZeroState
+from qualtran.bloqs.block_encoding.product import (
+    _product_block_encoding,
+    _product_block_encoding_override,
+    _product_block_encoding_symb,
+    Product,
+)
+from qualtran.bloqs.block_encoding.unitary import Unitary
+from qualtran.cirq_interop.testing import assert_circuit_inp_out_cirqsim
+
+
+def test_product(bloq_autotester):
+    bloq_autotester(_product_block_encoding)
+
+
+def test_product_signature():
+    assert _product_block_encoding().signature == Signature(
+        [Register("system", QAny(1)), Register("ancilla", QAny(1))]
+    )
+    assert _product_block_encoding_override().signature == Signature(
+        [Register("system", QAny(1)), Register("ancilla", QAny(3)), Register("resource", QAny(2))]
+    )
+    assert _product_block_encoding_symb().signature == Signature(
+        [
+            Register("system", QAny(1)),
+            Register("ancilla", QAny(sympy.Max(sympy.Symbol('a1'), sympy.Symbol('a2')) + 1)),
+        ]
+    )
+
+
+def test_product_checks():
+    with pytest.raises(ValueError):
+        _ = Product(())
+    with pytest.raises(ValueError):
+        _ = Product((Unitary(TGate()), Unitary(CNOT())))
+
+
+def test_product_params():
+    bloq = _product_block_encoding()
+    assert bloq.system_bitsize == 1
+    assert bloq.alpha == 1
+    assert bloq.epsilon == 0
+    assert bloq.ancilla_bitsize == 1
+    assert bloq.resource_bitsize == 0
+
+    bloq = _product_block_encoding_override()
+    assert bloq.system_bitsize == 1
+    assert bloq.alpha == 0.5 * 0.5
+    assert bloq.epsilon == 0.5 * 0.01 + 0.5 * 0.1
+    assert bloq.ancilla_bitsize == max(2, 1) + 1
+    assert bloq.resource_bitsize == 1 + 1
+
+    bloq = _product_block_encoding_symb()
+    assert bloq.system_bitsize == 1
+    assert bloq.alpha == sympy.Symbol('alpha1') * sympy.Symbol('alpha2')
+    assert bloq.epsilon == sympy.Symbol('alpha1') * sympy.Symbol('eps1') + sympy.Symbol(
+        'alpha2'
+    ) * sympy.Symbol('eps2')
+    assert bloq.ancilla_bitsize == sympy.Max(sympy.Symbol('a1'), sympy.Symbol('a2')) + 2 - 1
+    assert bloq.resource_bitsize == 0
+
+
+def test_product_tensors():
+    bb = BloqBuilder()
+    system = bb.add_register("system", 1)
+    ancilla = cast(Soquet, bb.add(ZeroState()))
+    system, ancilla = bb.add_t(_product_block_encoding(), system=system, ancilla=ancilla)
+    bb.add(ZeroEffect(), q=ancilla)
+    bloq = bb.finalize(system=system)
+
+    from_gate = np.matmul(TGate().tensor_contract(), Hadamard().tensor_contract())
+    from_tensors = bloq.tensor_contract()
+    np.testing.assert_allclose(from_gate, from_tensors)
+
+
+def test_product_single_tensors():
+    bb = BloqBuilder()
+    system = bb.add_register("system", 1)
+    system = cast(Soquet, bb.add(Product((Unitary(TGate()),)), system=system))
+    bloq = bb.finalize(system=system)
+
+    from_gate = TGate().tensor_contract()
+    from_tensors = bloq.tensor_contract()
+    np.testing.assert_allclose(from_gate, from_tensors)
+
+
+def test_product_override_tensors():
+    bb = BloqBuilder()
+    system = bb.add_register("system", 1)
+    ancilla = bb.join(np.array([bb.add(ZeroState()), bb.add(ZeroState()), bb.add(ZeroState())]))
+    resource = bb.join(np.array([bb.add(ZeroState()), bb.add(ZeroState())]))
+    system, ancilla, resource = bb.add_t(
+        _product_block_encoding_override(), system=system, ancilla=ancilla, resource=resource
+    )
+    for q in bb.split(cast(Soquet, ancilla)):
+        bb.add(ZeroEffect(), q=q)
+    for q in bb.split(cast(Soquet, resource)):
+        bb.add(ZeroEffect(), q=q)
+    bloq = bb.finalize(system=system)
+
+    from_gate = np.matmul(TGate().tensor_contract(), Hadamard().tensor_contract())
+    from_tensors = bloq.tensor_contract()
+    np.testing.assert_allclose(from_gate, from_tensors)
+
+
+def test_product_cirq():
+    qubits = cirq.LineQubit.range(3)
+    op = Product((Product((Unitary(XGate()), Unitary(XGate()))), Unitary(XGate()))).on_registers(
+        system=qubits[:1], ancilla=qubits[1:]
+    )
+    circuit = cirq.Circuit(
+        cirq.decompose_once(
+            op,
+            context=cirq.DecompositionContext(
+                cirq.GreedyQubitManager(prefix="_a", maximize_reuse=True)
+            ),
+        )
+    )
+    initial_state = [0, 0, 0]
+    final_state = [1, 0, 0]
+    assert_circuit_inp_out_cirqsim(circuit, qubits, initial_state, final_state)
diff --git a/qualtran/conftest.py b/qualtran/conftest.py
index 429c3e10d0..41cfd3c024 100644
--- a/qualtran/conftest.py
+++ b/qualtran/conftest.py
@@ -97,6 +97,9 @@ def assert_bloq_example_serializes_for_pytest(bloq_ex: BloqExample):
         'tensor_product_block_encoding',
         'tensor_product_block_encoding_override',
         'tensor_product_block_encoding_symb',
+        'product_block_encoding',
+        'product_block_encoding_override',
+        'product_block_encoding_symb',
     ]:
         pytest.xfail("Skipping serialization test for bloq examples that cannot yet be serialized.")
 
diff --git a/qualtran/serialization/resolver_dict.py b/qualtran/serialization/resolver_dict.py
index 8a039515fc..c90c526522 100644
--- a/qualtran/serialization/resolver_dict.py
+++ b/qualtran/serialization/resolver_dict.py
@@ -34,6 +34,7 @@
 import qualtran.bloqs.basic_gates.z_basis
 import qualtran.bloqs.block_encoding
 import qualtran.bloqs.block_encoding.lcu_select_and_prepare
+import qualtran.bloqs.block_encoding.product
 import qualtran.bloqs.block_encoding.tensor_product
 import qualtran.bloqs.block_encoding.unitary
 import qualtran.bloqs.bookkeeping
@@ -192,6 +193,7 @@
     "qualtran.bloqs.block_encoding.chebyshev_polynomial.ChebyshevPolynomial": qualtran.bloqs.block_encoding.chebyshev_polynomial.ChebyshevPolynomial,
     "qualtran.bloqs.block_encoding.unitary.Unitary": qualtran.bloqs.block_encoding.unitary.Unitary,
     "qualtran.bloqs.block_encoding.tensor_product.TensorProduct": qualtran.bloqs.block_encoding.tensor_product.TensorProduct,
+    "qualtran.bloqs.block_encoding.product.Product": qualtran.bloqs.block_encoding.product.Product,
     "qualtran.bloqs.bookkeeping.allocate.Allocate": qualtran.bloqs.bookkeeping.allocate.Allocate,
     "qualtran.bloqs.bookkeeping.arbitrary_clifford.ArbitraryClifford": qualtran.bloqs.bookkeeping.arbitrary_clifford.ArbitraryClifford,
     "qualtran.bloqs.bookkeeping.auto_partition.AutoPartition": qualtran.bloqs.bookkeeping.auto_partition.AutoPartition,
diff --git a/qualtran/symbolics/math_funcs.py b/qualtran/symbolics/math_funcs.py
index 2add3071c4..cdd4cbb664 100644
--- a/qualtran/symbolics/math_funcs.py
+++ b/qualtran/symbolics/math_funcs.py
@@ -71,12 +71,34 @@ def bit_length(x: SymbolicFloat) -> SymbolicInt:
 
 
 def smax(*args):
-    if any(isinstance(arg, sympy.Basic) for arg in args):
+    if len(args) == 0:
+        raise ValueError("smax expected at least 1 argument, got 0")
+    if len(args) == 1:
+        (it,) = args
+        if isinstance(it, Iterable):
+            args = tuple(arg for arg in it)
+        if len(args) == 0:
+            raise ValueError("smax() arg is an empty sequence")
+    if len(args) == 1:
+        (arg,) = args
+        return arg
+    if is_symbolic(*args):
         return sympy.Max(*args)
     return max(*args)
 
 
 def smin(*args):
+    if len(args) == 0:
+        raise ValueError("smin expected at least 1 argument, got 0")
+    if len(args) == 1:
+        (it,) = args
+        if isinstance(it, Iterable):
+            args = tuple(arg for arg in it)
+        if len(args) == 0:
+            raise ValueError("smin() arg is an empty sequence")
+    if len(args) == 1:
+        (arg,) = args
+        return arg
     if is_symbolic(*args):
         return sympy.Min(*args)
     return min(*args)