diff --git a/dev_tools/autogenerate-bloqs-notebooks-v2.py b/dev_tools/autogenerate-bloqs-notebooks-v2.py
index 4d7f79f31..f56a3be69 100644
--- a/dev_tools/autogenerate-bloqs-notebooks-v2.py
+++ b/dev_tools/autogenerate-bloqs-notebooks-v2.py
@@ -578,6 +578,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 7b0523790..ae05c6c9c 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 d5400fcbc..7f00f5d4a 100644
--- a/qualtran/bloqs/block_encoding/block_encoding.ipynb
+++ b/qualtran/bloqs/block_encoding/block_encoding.ipynb
@@ -890,6 +890,181 @@
     "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\n",
+    "$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 block encoding of $U_1 * \\cdots * U_n$ with normalization\n",
+    "constant $\\prod_i \\alpha_i$, ancilla bitsize $n - 1 + \\max_i a_i$, and precision\n",
+    "$\\sum_i \\alpha_i \\epsilon_i$.\n",
+    "\n",
+    "Following Fig. 2 in Dalzell et al. (2023), Ch. 10.2, the product is encoded by concatenating\n",
+    "each constituent block encoding, using a shared ancilla register and a set of flag qubits to\n",
+    "verify that the ancilla is left as zero after each use:\n",
+    "```\n",
+    "       ┌────────┐\n",
+    "  |0> ─┤        ├─     |0> ───────────X──────X────\n",
+    "       │        │                     │\n",
+    "       │ U_(AB) │  =        ┌─────┐   │   ┌─────┐\n",
+    "  |0> ─┤        ├─     |0> ─┤     ├──(0)──┤     ├─\n",
+    "       │        │           │ U_B │       │ U_A │\n",
+    "|Psi> ─┤        ├─   |Psi> ─┤     ├───────┤     ├─\n",
+    "       └────────┘           └─────┘       └─────┘\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": "feefd8ac",
+   "metadata": {
+    "cq.autogen": "Product.example_instances.md"
+   },
+   "source": [
+    "### Example Instances"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1067eef2",
+   "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": "a3c9bef3",
+   "metadata": {
+    "cq.autogen": "Product.product_block_encoding_properties"
+   },
+   "outputs": [],
+   "source": [
+    "from qualtran.bloqs.basic_gates import Hadamard, TGate\n",
+    "from qualtran.bloqs.block_encoding.unitary import Unitary\n",
+    "\n",
+    "u1 = Unitary(TGate(), alpha=0.5, ancilla_bitsize=2, resource_bitsize=1, epsilon=0.01)\n",
+    "u2 = Unitary(Hadamard(), alpha=0.5, ancilla_bitsize=1, resource_bitsize=1, epsilon=0.1)\n",
+    "product_block_encoding_properties = Product((u1, u2))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a36f7cd9",
+   "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": "f372c55c",
+   "metadata": {
+    "cq.autogen": "Product.graphical_signature.md"
+   },
+   "source": [
+    "#### Graphical Signature"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "189c579f",
+   "metadata": {
+    "cq.autogen": "Product.graphical_signature.py"
+   },
+   "outputs": [],
+   "source": [
+    "from qualtran.drawing import show_bloqs\n",
+    "show_bloqs([product_block_encoding, product_block_encoding_properties, product_block_encoding_symb],\n",
+    "           ['`product_block_encoding`', '`product_block_encoding_properties`', '`product_block_encoding_symb`'])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "442b1e19",
+   "metadata": {
+    "cq.autogen": "Product.call_graph.md"
+   },
+   "source": [
+    "### Call Graph"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "399ee3dd",
+   "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 000000000..ff067e2d4
--- /dev/null
+++ b/qualtran/bloqs/block_encoding/product.py
@@ -0,0 +1,284 @@
+#  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
+
+import cirq
+from attrs import evolve, field, frozen, validators
+from numpy.typing import NDArray
+
+from qualtran import (
+    bloq_example,
+    BloqBuilder,
+    BloqDocSpec,
+    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.bloqs.mcmt.multi_control_multi_target_pauli import MultiControlPauli
+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 block encoding of $U_1 * \cdots * U_n$ with normalization
+    constant $\prod_i \alpha_i$, ancilla bitsize $n - 1 + \max_i a_i$, and precision
+    $\sum_i \alpha_i \epsilon_i$.
+
+    Following Fig. 2 in Dalzell et al. (2023), Ch. 10.2, the product is encoded by concatenating
+    each constituent block encoding, using a shared ancilla register and a set of flag qubits to
+    verify that the ancilla is left as zero after each use:
+    ```
+           ┌────────┐
+      |0> ─┤        ├─     |0> ───────────X──────X────
+           │        │                     │
+           │ U_(AB) │  =        ┌─────┐   │   ┌─────┐
+      |0> ─┤        ├─     |0> ─┤     ├──(0)──┤     ├─
+           │        │           │ U_B │       │ U_A │
+    |Psi> ─┤        ├─   |Psi> ─┤     ├───────┤     ├─
+           └────────┘           └─────┘       └─────┘
+    ```
+
+    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.
+        # Github issue: https://github.com/quantumlib/Qualtran/issues/1104
+        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(
+                            MultiControlPauli(tuple([0] * anc_bits), cirq.X),
+                            controls=bb.split(anc_used_soq),
+                            target=flag_bits_soq[i],
+                        )
+                        flag_bits_soq[i] = bb.add(XGate(), q=flag_bits_soq[i])
+                        anc_part_soqs["anc_used"] = bb.join(cast(NDArray, 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_properties() -> Product:
+    from qualtran.bloqs.basic_gates import Hadamard, TGate
+    from qualtran.bloqs.block_encoding.unitary import Unitary
+
+    u1 = Unitary(TGate(), alpha=0.5, ancilla_bitsize=2, resource_bitsize=1, epsilon=0.01)
+    u2 = Unitary(Hadamard(), alpha=0.5, ancilla_bitsize=1, resource_bitsize=1, epsilon=0.1)
+    product_block_encoding_properties = Product((u1, u2))
+    return product_block_encoding_properties
+
+
+@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_properties,
+        _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 000000000..0eceddde8
--- /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_properties,
+    _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_properties().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_properties()
+    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_properties_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_properties(), 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 60f36510c..23d8bfb7f 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',
         'apply_lth_bloq',
     ]:
         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 9ca40696a..016f9a35a 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
@@ -193,6 +194,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 fefa8c37a..743dd4e14 100644
--- a/qualtran/symbolics/math_funcs.py
+++ b/qualtran/symbolics/math_funcs.py
@@ -84,12 +84,52 @@ def bit_length(x: SymbolicFloat) -> SymbolicInt:
 
 
 def smax(*args):
-    if any(isinstance(arg, sympy.Basic) for arg in args):
+    """Returns the maximum of the given arguments, which may be symbolic.
+
+    Args:
+        args: Either a pack of arguments or a single Iterable of arguments.
+              At least one argument must be provided in this pack or Iterable.
+
+    Returns:
+        The maximum of the given arguments.
+    """
+    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):
+    """Returns the minimum of the given arguments, which may be symbolic.
+
+    Args:
+        args: Either a pack of arguments or a single Iterable of arguments.
+              At least one argument must be provided in this pack or Iterable.
+
+    Returns:
+        The minimum of the given arguments.
+    """
+    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)
diff --git a/qualtran/symbolics/math_funcs_test.py b/qualtran/symbolics/math_funcs_test.py
index 9aad7a0f8..60c8727a1 100644
--- a/qualtran/symbolics/math_funcs_test.py
+++ b/qualtran/symbolics/math_funcs_test.py
@@ -11,6 +11,8 @@
 #  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 numpy as np
 import pytest
 import sympy
@@ -62,12 +64,46 @@ def test_smax():
     assert smax(1, 2) == 2
     assert smax(1.1, 2.2) == 2.2
     assert smax(1, sympy.Symbol('x')) == sympy.Max(1, sympy.Symbol('x'))
+    with pytest.raises(ValueError):
+        _ = smax()
+    with pytest.raises(ValueError):
+        _ = smax(())
+    with pytest.raises(ValueError):
+        _ = smax(x for x in cast(list[int], []))
+    assert smax(1) == 1
+    assert smax([1]) == 1
+    assert smax((1, 2)) == 2
+    assert smax(None) == None
+    with pytest.raises(TypeError):
+        _ = smax(None, None)
+    with pytest.raises(TypeError):
+        _ = smax((None, None))
+    assert smax(x for x in [1]) == 1
+    assert smax(x for x in [1, 2]) == 2
+    assert smax(x for x in [1, sympy.Symbol('x')]) == sympy.Max(1, sympy.Symbol('x'))
 
 
 def test_smin():
     assert smin(1, 2) == 1
     assert smin(1.1, 2.2) == 1.1
     assert smin(1, sympy.Symbol('x')) == sympy.Min(1, sympy.Symbol('x'))
+    with pytest.raises(ValueError):
+        _ = smin()
+    with pytest.raises(ValueError):
+        _ = smin(())
+    with pytest.raises(ValueError):
+        _ = smin(x for x in cast(list[int], []))
+    assert smin(1) == 1
+    assert smin([1]) == 1
+    assert smin((1, 2)) == 1
+    assert smin(None) == None
+    with pytest.raises(TypeError):
+        _ = smin(None, None)
+    with pytest.raises(TypeError):
+        _ = smin((None, None))
+    assert smin(x for x in [1]) == 1
+    assert smin(x for x in [2, 1]) == 1
+    assert smin(x for x in [1, sympy.Symbol('x')]) == sympy.Min(1, sympy.Symbol('x'))
 
 
 def test_bit_length():