diff --git a/pyzx/circuit/__init__.py b/pyzx/circuit/__init__.py
index 8272be186..a9ba0a24f 100644
--- a/pyzx/circuit/__init__.py
+++ b/pyzx/circuit/__init__.py
@@ -403,7 +403,8 @@ def from_quipper_file(fname: str) -> 'Circuit':
     @staticmethod
     def from_qasm(s: str) -> 'Circuit':
         """Produces a :class:`Circuit` based on a QASM input string.
-        It ignores all the non-unitary instructions like measurements in the file.
+        Supports OpenQASM 2 and 3, including ``reset`` and ``measure``
+        statements.
         It currently doesn't support custom gates that have parameters."""
         from .qasmparser import QASMParser
         p = QASMParser()
@@ -412,7 +413,8 @@ def from_qasm(s: str) -> 'Circuit':
     @staticmethod
     def from_qasm_file(fname: str) -> 'Circuit':
         """Produces a :class:`Circuit` based on a QASM description of a circuit.
-        It ignores all the non-unitary instructions like measurements in the file.
+        Supports OpenQASM 2 and 3, including ``reset`` and ``measure``
+        statements.
         It currently doesn't support custom gates that have parameters."""
         from .qasmparser import QASMParser
         p = QASMParser()
@@ -438,6 +440,21 @@ def to_qasm(self, version: int = 2) -> str:
         else:
             s = """OPENQASM 2.0;\ninclude "qelib1.inc";\n"""
             s += "qreg q[{!s}];\n".format(self.qubits)
+        # Collect classical register declarations from measurement gates.
+        cregs: Dict[str, int] = {}
+        for g in self.gates:
+            if isinstance(g, Measurement) and g.result_symbol is not None:
+                # result_symbol is e.g. "c[2]"; extract register name and index.
+                if '[' in g.result_symbol:
+                    regname, idx_str = g.result_symbol.split('[', 1)
+                    idx = int(idx_str.rstrip(']')) + 1
+                    cregs[regname] = max(cregs.get(regname, 0), idx)
+        if version == 3:
+            for regname, size in cregs.items():
+                s += "bit[{!s}] {};\n".format(size, regname)
+        else:
+            for regname, size in cregs.items():
+                s += "creg {}[{!s}];\n".format(regname, size)
         for g in self.gates:
             s += g.to_qasm() + "\n"
         return s
diff --git a/pyzx/circuit/gates.py b/pyzx/circuit/gates.py
index 609c50e14..8cf2bf1d4 100644
--- a/pyzx/circuit/gates.py
+++ b/pyzx/circuit/gates.py
@@ -1295,6 +1295,42 @@ def reposition(self, mask, bit_mask=None):
         return g
 
 
+class Reset(Gate):
+    """Reset an existing qubit to |0⟩.
+
+    Corresponds to the OpenQASM ``reset`` instruction, which discards
+    the current qubit state and unconditionally prepares ``|0⟩``.
+
+    In the ZX-diagram this is represented as a Z spider connected to
+    ground (tracing out / discarding the qubit) followed by a
+    disconnected X spider with phase 0 (state preparation ``|0⟩``).
+    This mirrors the ``DiscardBit`` pattern and models reset as a CPTP
+    map.
+    """
+    name = 'Reset'
+
+    def __init__(self, target: int):
+        self.target = target
+        self.label = target
+
+    def __str__(self) -> str:
+        return "Reset({})".format(self.target)
+
+    def __eq__(self, other: object) -> bool:
+        if not isinstance(other, Reset):
+            return False
+        return self.target == other.target
+
+    def to_qasm(self) -> str:
+        return "reset q[{:d}];".format(self.target)
+
+    def reposition(self, mask, bit_mask=None):
+        g = self.copy()
+        g.target = mask[self.target]
+        g.label = g.target
+        return g
+
+
 class PostSelect(Gate):
     """Post-select a qubit in a specified state.
 
@@ -1384,8 +1420,17 @@ def __eq__(self, other: object) -> bool:
         if not isinstance(other, Measurement): return False
         if self.target != other.target: return False
         if self.result_bit != other.result_bit: return False
+        if self.result_symbol != other.result_symbol: return False
         return True
 
+    def to_qasm(self) -> str:
+        if self.result_symbol is not None:
+            return "measure q[{:d}] -> {};".format(self.target, self.result_symbol)
+        if self.result_bit is not None:
+            return "measure q[{:d}] -> c[{:d}];".format(self.target, self.result_bit)
+        raise TypeError("Measurement on qubit {} has no result destination".format(
+            self.target))
+
     def reposition(self, mask, bit_mask = None):
         g = self.copy()
         g.target = mask[self.target]
@@ -1476,6 +1521,7 @@ def to_graph(self, g, q_mapper, c_mapper, ground=False):
     "RXX": RXX,
     "FSim": FSim,
     "InitAncilla": InitAncilla,
+    "Reset": Reset,
     "PostSelect": PostSelect,
     "DiscardBit": DiscardBit,
     "Measurement": Measurement,
diff --git a/pyzx/circuit/graphparser.py b/pyzx/circuit/graphparser.py
index 3dc507623..68e54a7af 100644
--- a/pyzx/circuit/graphparser.py
+++ b/pyzx/circuit/graphparser.py
@@ -17,7 +17,7 @@
 from typing import Dict, List, Optional
 
 from . import Circuit
-from .gates import Measurement, TargetMapper
+from .gates import InitAncilla, Measurement, Reset, TargetMapper
 from ..utils import EdgeType, VertexType, FloatInt, FractionLike
 from ..graph import Graph
 from ..graph.base import BaseGraph, VT, ET
@@ -30,9 +30,17 @@ def graph_to_circuit(g:BaseGraph[VT,ET], split_phases:bool=True) -> Circuit:
     ty = g.types()
     phases = g.phases()
     rows: Dict[FloatInt,List[VT]] = {}
+    # Map (vertex_type, phase) → InitAncilla state string.
+    _reverse_state_map: Dict[tuple, str] = {
+        (VertexType.Z, 0): '+',
+        (VertexType.Z, 1): '-',
+        (VertexType.X, 0): '0',
+        (VertexType.X, 1): '1',
+    }
+    input_qubits = set(qs[v] for v in inputs)
 
     c = Circuit(len(inputs))
-    
+
     for v in g.vertices():
         if v in inputs: continue
         r = g.row(v)
@@ -44,6 +52,23 @@ def graph_to_circuit(g:BaseGraph[VT,ET], split_phases:bool=True) -> Circuit:
             phase = phases[v]
             t = ty[v]
             neigh = [w for w in g.neighbors(v) if rs[w]<r]
+            if len(neigh) == 0:
+                # No backward neighbour: state preparation vertex.
+                qi = int(q)
+                if t == VertexType.X and phase == 0 and q in input_qubits:
+                    # Qubit already has an input boundary → mid-circuit reset.
+                    # NOTE: this heuristic assumes any disconnected X(0) spider
+                    # on an input qubit is a reset.  It could misidentify a
+                    # manually constructed graph that uses X(0) preparation on
+                    # an input qubit for a different purpose.
+                    c.add_gate(Reset(qi))
+                    continue
+                state = _reverse_state_map.get((t, phase))
+                if state is not None:
+                    c.add_gate(InitAncilla(qi, state))
+                    continue
+                raise TypeError("Graph doesn't seem circuit like: "
+                                "vertex {} has no parents".format(v))
             if len(neigh) != 1:
                 raise TypeError("Graph doesn't seem circuit like: multiple parents")
             n = neigh[0]
@@ -53,6 +78,11 @@ def graph_to_circuit(g:BaseGraph[VT,ET], split_phases:bool=True) -> Circuit:
                 c.add_gate("HAD", q)
             if t == VertexType.BOUNDARY: #vertex is an output
                 continue
+            if g.is_ground(v):
+                # Ground vertex: discard/trace-out, part of a Reset pair.
+                # (The corresponding Reset gate is emitted when the state
+                # preparation vertex on this qubit is processed.)
+                continue
             if phase!=0 and not split_phases:
                 if t == VertexType.Z: c.add_gate("ZPhase", q, phase=phase)
                 else: c.add_gate("XPhase", q, phase=phase)
@@ -118,8 +148,9 @@ def circuit_to_graph(
         qubit = i+c.qubits
         v = g.add_vertex(VertexType.BOUNDARY, qubit, 0)
         inputs.append(v)
-        q_mapper.add_label(qubit, 1)
-        c_mapper.set_prev_vertex(qubit, v)
+        c_mapper.add_label(i, 1)
+        c_mapper.set_qubit(i, qubit)
+        c_mapper.set_prev_vertex(i, v)
 
 
     for gate in c.gates:
@@ -128,24 +159,41 @@ def circuit_to_graph(
             measure_targets.add(gate.target)
         if gate.name == 'InitAncilla':
             l = gate.label # type: ignore
-            try:
-                q_mapper.add_label(l, q_mapper.next_row_or_default(l, q_mapper.max_row() - 1))
-            except ValueError:
+            if l in q_mapper.labels():
                 raise ValueError("Ancilla label {} already in use".format(str(l)))
             vtype, phase = gate.get_vertex_info() # type: ignore
+            q_mapper.add_label(l, q_mapper.next_row_or_default(l, q_mapper.max_row() - 1))
             v = g.add_vertex(vtype, q_mapper.to_qubit(l), q_mapper.next_row(l), phase)
             q_mapper.set_prev_vertex(l, v)
             q_mapper.advance_next_row(l)
+        elif gate.name == 'Reset':
+            # Model reset as discard (ground) then preparation |0⟩,
+            # creating a wire break.  The ground connection traces out
+            # the qubit, modelling reset as a CPTP map.
+            l = gate.label # type: ignore
+            q = q_mapper.to_qubit(l)
+            # The qubit is being reused, so clear any pending measurement
+            # marker.  If the qubit is measured again later without a
+            # subsequent reset, it will be re-added.
+            measure_targets.discard(q)
+            r = q_mapper.next_row(l)
+            u = q_mapper.prev_vertex(l)
+            # Effect vertex: discard the old state by tracing out.
+            effect_v = g.add_vertex(VertexType.Z, q, r, 0, ground=True)
+            g.add_edge((u, effect_v), EdgeType.SIMPLE)
+            # State vertex: prepare the new |0⟩ state.
+            state_v = g.add_vertex(VertexType.X, q, r + 1, 0)
+            q_mapper.set_prev_vertex(l, state_v)
+            q_mapper.set_next_row(l, r + 2)
         elif gate.name == 'PostSelect':
             l = gate.label # type: ignore
-            try:
-                q = q_mapper.to_qubit(l)
-                r = q_mapper.next_row(l)
-                u = q_mapper.prev_vertex(l)
-                q_mapper.set_next_row(l, r + 1)
-                q_mapper.remove_label(l)
-            except ValueError:
+            if l not in q_mapper.labels():
                 raise ValueError("PostSelect label {} is not in use".format(str(l)))
+            q = q_mapper.to_qubit(l)
+            r = q_mapper.next_row(l)
+            u = q_mapper.prev_vertex(l)
+            q_mapper.set_next_row(l, r + 1)
+            q_mapper.remove_label(l)
             vtype, phase = gate.get_vertex_info() # type: ignore
             v = g.add_vertex(vtype, q, r, phase)
             g.add_edge((u,v),EdgeType.SIMPLE)
diff --git a/pyzx/circuit/qasmparser.py b/pyzx/circuit/qasmparser.py
index 6cf36ea31..3c33e6324 100644
--- a/pyzx/circuit/qasmparser.py
+++ b/pyzx/circuit/qasmparser.py
@@ -21,7 +21,7 @@
 from typing import List, Dict, Tuple, Optional
 
 from . import Circuit
-from .gates import Gate, qasm_gate_table, Measurement
+from .gates import Gate, qasm_gate_table, Measurement, Reset
 from ..utils import settings
 
 
@@ -33,6 +33,7 @@ def __init__(self) -> None:
         self.gates: List[Gate] = []
         self.custom_gates: Dict[str,Circuit] = {}
         self.registers: Dict[str,Tuple[int,int]] = {}
+        self.cregisters: Dict[str,int] = {}
         self.qubit_count: int = 0
         self.bit_count: int = 0
         self.circuit: Optional[Circuit] = None
@@ -41,7 +42,9 @@ def parse(self, s: str, strict:bool=True) -> Circuit:
         self.gates = []
         self.custom_gates = {}
         self.registers = {}
+        self.cregisters = {}
         self.qubit_count = 0
+        self.bit_count = 0
         self.circuit = None
         lines = s.splitlines()
         r = []
@@ -78,7 +81,8 @@ def parse(self, s: str, strict:bool=True) -> Circuit:
         for c in commands:
             self.gates.extend(self.parse_command(c, self.registers))
 
-        circ = Circuit(self.qubit_count)
+        self.bit_count = sum(self.cregisters.values())
+        circ = Circuit(self.qubit_count, bit_amount=self.bit_count)
         circ.gates = self.gates
         self.circuit = circ
         return self.circuit
@@ -141,16 +145,70 @@ def extract_command_parts(self, c: str) -> Tuple[str,List[Fraction],List[str]]:
     def parse_command(self, c: str, registers: Dict[str,Tuple[int,int]]) -> List[Gate]:
         gates: List[Gate] = []
         name, phases, args = self.extract_command_parts(c)
-        if name in ("barrier","creg", "id"): return gates
+        if name in ("barrier", "id"): return gates
+        if name == "creg":
+            regname, sizep = args[0].split("[", 1)
+            size = int(sizep[:-1])
+            self.cregisters[regname] = size
+            return gates
         if name == "measure":
-            target, result_bit = args[0].split(' -> ')
-            # Extract the register name and index separately for both target and result
-            _, target_idx = target.split('[')
-            result_reg_name, result_idx = result_bit.split('[')
-            # Remove the trailing ']' and convert to int
-            target_qbit = int(target_idx[:-1])
-            gate = Measurement(target_qbit, result_symbol=f"{result_reg_name}[{result_idx[:-1]}]")
-            gates.append(gate)
+            # Strip all spaces so we can split on '->' regardless of spacing.
+            raw = args[0].replace(' ', '')
+            target, result_bit = raw.split('->')
+            if '[' in target:
+                # Single qubit: measure q[0] -> c[0];
+                regname, target_idx = target.split('[', 1)
+                if regname not in registers:
+                    raise TypeError("Invalid register {}".format(regname))
+                target_qbit = registers[regname][0] + int(target_idx[:-1])
+                result_reg_name, result_idx = result_bit.split('[', 1)
+                result_index = int(result_idx[:-1])
+                if result_reg_name not in self.cregisters:
+                    raise TypeError("Undeclared classical register {}".format(
+                        result_reg_name))
+                if result_index >= self.cregisters[result_reg_name]:
+                    raise TypeError(
+                        "Index {} out of range for classical register {} "
+                        "of size {}".format(
+                            result_index, result_reg_name,
+                            self.cregisters[result_reg_name]))
+                gate = Measurement(target_qbit, result_symbol=f"{result_reg_name}[{result_index}]")
+                gates.append(gate)
+            else:
+                # Register broadcast: measure q -> c.
+                if target not in registers:
+                    raise TypeError("Invalid register {}".format(target))
+                if result_bit not in self.cregisters:
+                    raise TypeError("Undeclared classical register {}".format(
+                        result_bit))
+                start, size = registers[target]
+                if size > self.cregisters[result_bit]:
+                    raise TypeError(
+                        "Quantum register {} (size {}) is larger than "
+                        "classical register {} (size {})".format(
+                            target, size, result_bit,
+                            self.cregisters[result_bit]))
+                for i in range(size):
+                    gate = Measurement(start + i, result_symbol=f"{result_bit}[{i}]")
+                    gates.append(gate)
+            return gates
+        if name == "reset":
+            # Reset initializes a qubit to |0⟩ state.
+            # In OpenQASM: `reset q[0];` or `reset q;` (for entire register)
+            for a in args:
+                if "[" in a:
+                    regname, valp = a.split("[", 1)
+                    val = int(valp[:-1])
+                    if regname not in registers:
+                        raise TypeError("Invalid register {}".format(regname))
+                    qubit = registers[regname][0] + val
+                    gates.append(Reset(qubit))
+                else:
+                    if a not in registers:
+                        raise TypeError("Invalid register {}".format(a))
+                    start, size = registers[a]
+                    for i in range(size):
+                        gates.append(Reset(start + i))
             return gates
         if name in ("opaque", "if"):
             raise TypeError("Unsupported operation {}".format(c))
diff --git a/pyzx/symbolic.py b/pyzx/symbolic.py
index cc82b17a7..6b9792dfd 100644
--- a/pyzx/symbolic.py
+++ b/pyzx/symbolic.py
@@ -405,7 +405,7 @@ def new_const(coeff: Union[int, Fraction]) -> Poly:
     ?factor    : base ("^" exponent)?
     base       : intf | frac | decimal | pi | pifrac | var | "(" expr ")"
     exponent   : intf
-    var        : CNAME
+    var        : CNAME ("[" INT "]")?
     intf       : INT
     decimal    : DECIMAL
     pi         : "\\pi" | "pi" | "π"
@@ -464,6 +464,8 @@ def exponent(self, items: List[Any]) -> int:
 
     def var(self, items: List[Any]) -> Poly:
         v = str(items[0])
+        if len(items) > 1 and items[1] is not None:
+            v += "[{}]".format(int(items[1]))
         return self._new_var(v)
 
     def pi(self, _: List[Any]) -> Poly:
diff --git a/tests/test_init_postselect.py b/tests/test_init_postselect.py
index 9733623e1..15be25995 100644
--- a/tests/test_init_postselect.py
+++ b/tests/test_init_postselect.py
@@ -23,7 +23,7 @@
     sys.path.append('.')
 
 from pyzx.circuit import Circuit
-from pyzx.circuit.gates import InitAncilla, PostSelect
+from pyzx.circuit.gates import InitAncilla, Measurement, PostSelect, Reset
 from pyzx.utils import VertexType
 
 
@@ -102,6 +102,66 @@ def test_equality(self):
         self.assertNotEqual(PostSelect(0), InitAncilla(0))
 
 
+class TestReset(unittest.TestCase):
+    """Tests for Reset gate."""
+
+    def test_str(self):
+        self.assertEqual(str(Reset(0)), "Reset(0)")
+        self.assertEqual(str(Reset(3)), "Reset(3)")
+
+    def test_equality(self):
+        self.assertEqual(Reset(0), Reset(0))
+        self.assertNotEqual(Reset(0), Reset(1))
+        self.assertNotEqual(Reset(0), InitAncilla(0, '0'))
+
+    def test_to_qasm(self):
+        self.assertEqual(Reset(2).to_qasm(), "reset q[2];")
+
+    def test_reposition(self):
+        mask = [0, 5, 2, 3]
+        g1 = Reset(1)
+        g2 = g1.reposition(mask)
+        self.assertEqual(g2.target, 5)
+        self.assertEqual(g2.label, 5)
+
+    def test_to_graph(self):
+        """Test that Reset on an existing qubit produces ground + state vertices."""
+        c = Circuit(1)
+        c.add_gate(Reset(0))
+        g = c.to_graph()
+
+        # Effect vertex: Z spider connected to ground (discard).
+        ground_verts = list(g.grounds())
+        self.assertEqual(len(ground_verts), 1)
+        self.assertEqual(g.type(ground_verts[0]), VertexType.Z)
+        self.assertEqual(g.phase(ground_verts[0]), 0)
+
+        # State vertex: X spider phase 0 (|0⟩ preparation).
+        x_verts = [v for v in g.vertices() if g.type(v) == VertexType.X]
+        self.assertEqual(len(x_verts), 1)
+        self.assertEqual(g.phase(x_verts[0]), 0)
+
+
+class TestMeasurementEquality(unittest.TestCase):
+    """Tests for Measurement equality."""
+
+    def test_same(self):
+        self.assertEqual(Measurement(0, result_symbol="c[0]"),
+                         Measurement(0, result_symbol="c[0]"))
+
+    def test_different_target(self):
+        self.assertNotEqual(Measurement(0, result_symbol="c[0]"),
+                            Measurement(1, result_symbol="c[0]"))
+
+    def test_different_result_symbol(self):
+        self.assertNotEqual(Measurement(0, result_symbol="c[0]"),
+                            Measurement(0, result_symbol="d[0]"))
+
+    def test_different_result_bit(self):
+        self.assertNotEqual(Measurement(0, result_bit=0),
+                            Measurement(0, result_bit=1))
+
+
 class TestReposition(unittest.TestCase):
     """Tests for repositioning InitAncilla/PostSelect gates."""
 
diff --git a/tests/test_qasm.py b/tests/test_qasm.py
index eb0639f44..1152680e1 100644
--- a/tests/test_qasm.py
+++ b/tests/test_qasm.py
@@ -45,6 +45,12 @@
 except ImportError:
     QuantumCircuit = None
 
+stim: Optional[ModuleType]
+try:
+    import stim  # type: ignore[no-redef,import]
+except ImportError:
+    stim = None
+
 
 @unittest.skipUnless(np, "numpy needs to be installed for this to run")
 class TestQASM(unittest.TestCase):
@@ -359,6 +365,703 @@ def test_qiskit_transpile_pyzx_optimization_round_trip(self):
 
         self.assertTrue(compare_tensors(t1, t3))
 
+    def test_reset_single_qubit(self):
+        """Test that 'reset' command is parsed for a single qubit."""
+        from pyzx.circuit.gates import Reset
+        c = Circuit.from_qasm("""
+        OPENQASM 2.0;
+        include "qelib1.inc";
+        qreg q[2];
+        reset q[0];
+        h q[0];
+        """)
+        self.assertEqual(c.qubits, 2)
+        self.assertEqual(len(c.gates), 2)
+        self.assertIsInstance(c.gates[0], Reset)
+        self.assertEqual(c.gates[0].target, 0)
+
+    def test_reset_entire_register(self):
+        """Test that 'reset' command broadcasts over entire register."""
+        from pyzx.circuit.gates import Reset
+        c = Circuit.from_qasm("""
+        OPENQASM 2.0;
+        include "qelib1.inc";
+        qreg q[3];
+        reset q;
+        """)
+        self.assertEqual(c.qubits, 3)
+        self.assertEqual(len(c.gates), 3)
+        for i, gate in enumerate(c.gates):
+            self.assertIsInstance(gate, Reset)
+            self.assertEqual(gate.target, i)
+
+    def test_reset_openqasm3(self):
+        """Test that 'reset' command works with OpenQASM 3."""
+        from pyzx.circuit.gates import Reset
+        c = Circuit.from_qasm("""
+        OPENQASM 3;
+        include "stdgates.inc";
+        qubit[2] q;
+        reset q[1];
+        x q[1];
+        """)
+        self.assertEqual(c.qubits, 2)
+        self.assertEqual(len(c.gates), 2)
+        self.assertIsInstance(c.gates[0], Reset)
+        self.assertEqual(c.gates[0].target, 1)
+
+    def test_reset_to_graph(self):
+        """Test that reset creates appropriate vertices when converted to graph."""
+        from pyzx.utils import VertexType
+        c = Circuit.from_qasm("""
+        OPENQASM 2.0;
+        include "qelib1.inc";
+        qreg q[1];
+        reset q[0];
+        """)
+        g = c.to_graph()
+
+        # Effect vertex: Z spider connected to ground (discard).
+        ground_verts = list(g.grounds())
+        self.assertEqual(len(ground_verts), 1)
+        self.assertEqual(g.type(ground_verts[0]), VertexType.Z)
+
+        # State vertex: X spider phase 0 (|0⟩ preparation).
+        x_vertices = [v for v in g.vertices() if g.type(v) == VertexType.X]
+        self.assertEqual(len(x_vertices), 1)
+        self.assertEqual(g.phase(x_vertices[0]), 0)
+
+    def test_reset_to_graph_has_output(self):
+        """Test that a qubit has an output boundary after reset."""
+        c = Circuit.from_qasm("""
+        OPENQASM 2.0;
+        include "qelib1.inc";
+        qreg q[1];
+        reset q[0];
+        h q[0];
+        """)
+        g = c.to_graph()
+        self.assertEqual(len(g.inputs()), 1)
+        self.assertEqual(len(g.outputs()), 1)
+
+    def test_measure_reset_has_output(self):
+        """Test that a qubit has an output after measure followed by reset."""
+        c = Circuit.from_qasm("""
+        OPENQASM 2.0;
+        include "qelib1.inc";
+        qreg q[1];
+        creg c[1];
+        h q[0];
+        measure q[0] -> c[0];
+        reset q[0];
+        x q[0];
+        """)
+        g = c.to_graph()
+        # 1 qubit + 1 classical bit.
+        self.assertEqual(len(g.inputs()), 2)
+        # q[0] after reset+x gets output, c[0] gets pass-through output.
+        self.assertEqual(len(g.outputs()), 2)
+
+    def test_measure_reset_measure_no_output(self):
+        """Test that a qubit ending with measurement has no output boundary."""
+        c = Circuit.from_qasm("""
+        OPENQASM 2.0;
+        include "qelib1.inc";
+        qreg q[2];
+        creg c[1];
+        h q[0];
+        measure q[0] -> c[0];
+        reset q[0];
+        h q[0];
+        measure q[0] -> c[0];
+        """)
+        g = c.to_graph()
+        # 2 qubits + 1 classical bit.
+        self.assertEqual(len(g.inputs()), 3)
+        # q[0] ends with a measurement (no output), q[1] and c[0] get outputs.
+        self.assertEqual(len(g.outputs()), 2)
+
+    def test_reset_qasm_round_trip(self):
+        """Test that reset survives a QASM round-trip."""
+        from pyzx.circuit.gates import Reset
+        qasm_in = """
+        OPENQASM 2.0;
+        include "qelib1.inc";
+        qreg q[2];
+        h q[0];
+        reset q[0];
+        cx q[0],q[1];
+        """
+        c1 = Circuit.from_qasm(qasm_in)
+        qasm_out = c1.to_qasm()
+        c2 = Circuit.from_qasm(qasm_out)
+        self.assertEqual(len(c1.gates), len(c2.gates))
+        for g1, g2 in zip(c1.gates, c2.gates):
+            self.assertEqual(type(g1), type(g2))
+            if isinstance(g1, Reset):
+                self.assertEqual(g1.target, g2.target)
+
+
+    def test_measure_register_broadcast(self):
+        """Test that 'measure q -> c' broadcasts over entire register."""
+        from pyzx.circuit.gates import Measurement
+        c = Circuit.from_qasm("""
+        OPENQASM 2.0;
+        include "qelib1.inc";
+        qreg q[3];
+        creg c[3];
+        h q[0];
+        h q[1];
+        h q[2];
+        measure q -> c;
+        """)
+        self.assertEqual(c.qubits, 3)
+        self.assertEqual(len(c.gates), 6)
+        for i in range(3):
+            gate = c.gates[3 + i]
+            self.assertIsInstance(gate, Measurement)
+            self.assertEqual(gate.target, i)
+
+
+    def test_graph_to_circuit_reset(self):
+        """Test that graph_to_circuit recovers Reset gates."""
+        from pyzx.circuit.gates import Reset
+        from pyzx.circuit.graphparser import graph_to_circuit
+        c1 = Circuit.from_qasm("""
+        OPENQASM 2.0;
+        include "qelib1.inc";
+        qreg q[1];
+        h q[0];
+        reset q[0];
+        x q[0];
+        """)
+        g = c1.to_graph()
+        c2 = graph_to_circuit(g)
+
+        # The extracted circuit should contain a Reset gate.
+        reset_gates = [gt for gt in c2.gates if isinstance(gt, Reset)]
+        self.assertEqual(len(reset_gates), 1)
+        self.assertEqual(reset_gates[0].target, 0)
+
+
+    def test_issue_345_circuit1_measure_reset(self):
+        """End-to-end test for issue #345 circuit 1 (Steane code with reset).
+
+        Verifies the full workflow: QASM with measure+reset parses, converts
+        to a valid ZX-graph, and round-trips through QASM.
+        """
+        qasm = """
+        OPENQASM 2.0;
+        include "qelib1.inc";
+        qreg q[8];
+        creg c[1];
+        h q[0];
+        cx q[0],q[1];
+        cx q[0],q[2];
+        cx q[0],q[3];
+        cx q[0],q[4];
+        h q[0];
+        measure q[0] -> c[0];
+        reset q[0];
+        h q[0];
+        cx q[0],q[1];
+        cx q[0],q[2];
+        cx q[0],q[5];
+        cx q[0],q[6];
+        h q[0];
+        measure q[0] -> c[0];
+        reset q[0];
+        h q[0];
+        cx q[0],q[1];
+        cx q[0],q[3];
+        cx q[0],q[5];
+        cx q[0],q[7];
+        h q[0];
+        measure q[0] -> c[0];
+        """
+        # Parse.
+        c = Circuit.from_qasm(qasm)
+        self.assertEqual(c.qubits, 8)
+
+        # Convert to graph.
+        g = c.to_graph()
+        # 8 qubits + 1 classical bit.
+        self.assertEqual(len(g.inputs()), 9)
+        # q[0] ends with measurement (no reset after); 7 qubit + 1 classical output.
+        self.assertEqual(len(g.outputs()), 8)
+
+        # QASM round-trip.
+        c2 = Circuit.from_qasm(c.to_qasm())
+        self.assertEqual(len(c.gates), len(c2.gates))
+        for g1, g2 in zip(c.gates, c2.gates):
+            self.assertEqual(type(g1), type(g2))
+
+    def test_issue_345_circuit2_ancilla_measure(self):
+        """End-to-end test for issue #345 circuit 2 (Steane code, no reset).
+
+        Uses three separate ancilla qubits instead of resetting one qubit.
+        Verifies that QASM with multiple registers and measurements parses
+        and converts to a valid ZX-graph.
+        """
+        qasm = """
+        OPENQASM 2.0;
+        include "qelib1.inc";
+        qreg ancilla[3];
+        qreg data[7];
+        creg bits[3];
+        h ancilla[0];
+        cx ancilla[0],data[0];
+        cx ancilla[0],data[1];
+        cx ancilla[0],data[2];
+        cx ancilla[0],data[3];
+        h ancilla[0];
+        measure ancilla[0] -> bits[0];
+        h ancilla[1];
+        cx ancilla[1],data[0];
+        cx ancilla[1],data[1];
+        cx ancilla[1],data[4];
+        cx ancilla[1],data[5];
+        h ancilla[1];
+        measure ancilla[1] -> bits[1];
+        h ancilla[2];
+        cx ancilla[2],data[0];
+        cx ancilla[2],data[2];
+        cx ancilla[2],data[4];
+        cx ancilla[2],data[6];
+        h ancilla[2];
+        measure ancilla[2] -> bits[2];
+        """
+        # Parse.
+        c = Circuit.from_qasm(qasm)
+        self.assertEqual(c.qubits, 10)
+
+        # Convert to graph.
+        g = c.to_graph()
+        # 10 qubits + 3 classical bits.
+        self.assertEqual(len(g.inputs()), 13)
+        # Ancilla qubits 0-2 are measured (no output), data qubits 3-9 and
+        # 3 classical bits get outputs.
+        self.assertEqual(len(g.outputs()), 10)
+
+        # QASM round-trip.
+        c2 = Circuit.from_qasm(c.to_qasm())
+        self.assertEqual(len(c.gates), len(c2.gates))
+        for g1, g2 in zip(c.gates, c2.gates):
+            self.assertEqual(type(g1), type(g2))
+
+    def test_measure_multi_register_offset(self):
+        """Measure on a register that is not the first declared uses
+        the correct global qubit index."""
+        from pyzx.circuit.gates import Measurement
+        qasm = """
+        OPENQASM 2.0;
+        include "qelib1.inc";
+        qreg data[4];
+        qreg anc[2];
+        creg c[2];
+        measure anc[0] -> c[0];
+        measure anc[1] -> c[1];
+        """
+        c = Circuit.from_qasm(qasm)
+        measurements = [g for g in c.gates if isinstance(g, Measurement)]
+        self.assertEqual(len(measurements), 2)
+        # anc[0] is global qubit 4, anc[1] is global qubit 5.
+        self.assertEqual(measurements[0].target, 4)
+        self.assertEqual(measurements[1].target, 5)
+
+    def test_measure_qasm_round_trip(self):
+        """Test that measure survives a QASM round-trip."""
+        from pyzx.circuit.gates import Measurement
+        qasm_in = """
+        OPENQASM 2.0;
+        include "qelib1.inc";
+        qreg q[2];
+        creg c[2];
+        h q[0];
+        h q[1];
+        measure q[0] -> c[0];
+        measure q[1] -> c[1];
+        """
+        c1 = Circuit.from_qasm(qasm_in)
+        self.assertEqual(c1.bits, 2)
+        qasm_out = c1.to_qasm()
+        self.assertIn("creg c[2]", qasm_out)
+        c2 = Circuit.from_qasm(qasm_out)
+        self.assertEqual(c2.bits, 2)
+        self.assertEqual(len(c1.gates), len(c2.gates))
+        for g1, g2 in zip(c1.gates, c2.gates):
+            self.assertEqual(type(g1), type(g2))
+            if isinstance(g1, Measurement):
+                self.assertEqual(g1.target, g2.target)
+
+    def test_measure_graph_json_round_trip(self):
+        """Test that measurement phases survive a graph JSON round-trip.
+
+        Measurement gates produce symbolic Poly phases like c[0] which
+        must be correctly serialised and deserialised by the JSON layer.
+        """
+        from pyzx.graph.jsonparser import graph_to_json, json_to_graph
+        from pyzx.symbolic import Poly
+
+        qasm = """
+        OPENQASM 2.0;
+        include "qelib1.inc";
+        qreg q[2];
+        creg c[2];
+        h q[0];
+        h q[1];
+        measure q[0] -> c[0];
+        measure q[1] -> c[1];
+        """
+        g1 = Circuit.from_qasm(qasm).to_graph()
+
+        # Collect the symbolic phases before the round-trip.
+        poly_phases_before = {
+            v: g1.phase(v) for v in g1.vertices()
+            if isinstance(g1.phase(v), Poly) and g1.phase(v) != 0
+        }
+        self.assertTrue(len(poly_phases_before) > 0,
+                        "Graph should contain at least one Poly phase")
+
+        # Round-trip through JSON.
+        g2 = json_to_graph(graph_to_json(g1))
+
+        # The round-tripped graph must have the same number of vertices,
+        # edges, and matching Poly phase strings.
+        self.assertEqual(g1.num_vertices(), g2.num_vertices())
+        self.assertEqual(g1.num_edges(), g2.num_edges())
+
+        poly_phases_after = {
+            v: g2.phase(v) for v in g2.vertices()
+            if isinstance(g2.phase(v), Poly) and g2.phase(v) != 0
+        }
+        self.assertEqual(
+            {v: str(p) for v, p in poly_phases_before.items()},
+            {v: str(p) for v, p in poly_phases_after.items()})
+
+    def test_circuit_bits_from_creg(self):
+        """Test that Circuit.bits reflects declared classical registers."""
+        c = Circuit.from_qasm("""
+        OPENQASM 2.0;
+        include "qelib1.inc";
+        qreg q[2];
+        creg a[3];
+        creg b[2];
+        measure q[0] -> a[0];
+        measure q[1] -> b[1];
+        """)
+        self.assertEqual(c.bits, 5)
+
+    def test_measure_undeclared_creg(self):
+        """Test that measure into an undeclared classical register raises."""
+        with self.assertRaises(TypeError) as ctx:
+            Circuit.from_qasm("""
+            OPENQASM 2.0;
+            include "qelib1.inc";
+            qreg q[1];
+            measure q[0] -> c[0];
+            """)
+        self.assertIn("Undeclared classical register", str(ctx.exception))
+
+    def test_measure_creg_index_out_of_range(self):
+        """Test that an out-of-range classical index raises."""
+        with self.assertRaises(TypeError) as ctx:
+            Circuit.from_qasm("""
+            OPENQASM 2.0;
+            include "qelib1.inc";
+            qreg q[1];
+            creg c[1];
+            measure q[0] -> c[1];
+            """)
+        self.assertIn("out of range", str(ctx.exception))
+
+    def test_measure_broadcast_undeclared_creg(self):
+        """Test that broadcast measure into an undeclared creg raises."""
+        with self.assertRaises(TypeError) as ctx:
+            Circuit.from_qasm("""
+            OPENQASM 2.0;
+            include "qelib1.inc";
+            qreg q[2];
+            measure q -> c;
+            """)
+        self.assertIn("Undeclared classical register", str(ctx.exception))
+
+    def test_measure_broadcast_creg_too_small(self):
+        """Test that broadcast measure into a too-small creg raises."""
+        with self.assertRaises(TypeError) as ctx:
+            Circuit.from_qasm("""
+            OPENQASM 2.0;
+            include "qelib1.inc";
+            qreg q[3];
+            creg c[2];
+            measure q -> c;
+            """)
+        self.assertIn("is larger than", str(ctx.exception))
+
+    def test_classical_bits_in_c_mapper(self):
+        """Classical bit labels should be in c_mapper, not q_mapper."""
+        from pyzx.circuit.graphparser import circuit_to_graph
+        c = Circuit(2, bit_amount=1)
+        g = circuit_to_graph(c)
+        # 2 quantum inputs + 1 classical input = 3 inputs.
+        self.assertEqual(len(g.inputs()), 3)
+        # All wires are straight-through, so 3 outputs.
+        self.assertEqual(len(g.outputs()), 3)
+
+    def test_classical_bits_output_positions(self):
+        """Classical output boundaries should have the right qubit positions."""
+        from pyzx.circuit.graphparser import circuit_to_graph
+        c = Circuit(2, bit_amount=2)
+        g = circuit_to_graph(c)
+        self.assertEqual(len(g.inputs()), 4)
+        self.assertEqual(len(g.outputs()), 4)
+        # Output boundary qubits should be 0, 1, 2, 3.
+        output_qubits = sorted(g.qubit(v) for v in g.outputs())
+        self.assertEqual(output_qubits, [0, 1, 2, 3])
+
+    def test_measurement_ground_mode_with_result_bit(self):
+        """Measurement with result_bit in ground mode should use c_mapper."""
+        from pyzx.circuit.gates import Measurement
+        from pyzx.circuit.graphparser import circuit_to_graph
+        c = Circuit(1, bit_amount=1)
+        c.gates = [Measurement(0, result_bit=0)]
+        g = circuit_to_graph(c)
+        # 1 quantum input + 1 classical input = 2 inputs.
+        self.assertEqual(len(g.inputs()), 2)
+        # Measurement consumes the qubit (no quantum output), but
+        # the classical bit still gets an output boundary.
+        self.assertEqual(len(g.outputs()), 1)
+
+    def test_measurement_ground_mode_graph_structure(self):
+        """Ground-mode measurement should create correct graph structure."""
+        from pyzx.circuit.gates import Measurement
+        from pyzx.circuit.graphparser import circuit_to_graph
+        from pyzx.utils import VertexType
+        c = Circuit(1, bit_amount=1)
+        m = Measurement(0, result_bit=0)
+        c.gates = [m]
+        g = circuit_to_graph(c)
+        # Count vertex types.
+        types = [g.type(v) for v in g.vertices()]
+        n_boundary = types.count(VertexType.BOUNDARY)
+        # 1 quantum input + 1 classical input + 1 classical output = 3.
+        self.assertEqual(n_boundary, 3)
+
+    def test_measurement_ground_mode_direct(self):
+        """Directly test to_graph_ground with c_mapper labels."""
+        from pyzx.circuit.gates import Measurement, TargetMapper
+        from pyzx.utils import VertexType
+        from pyzx.graph import Graph
+        g = Graph()
+        q_mapper = TargetMapper()
+        c_mapper = TargetMapper()
+        # Set up one quantum wire and one classical wire.
+        q_in = g.add_vertex(VertexType.BOUNDARY, 0, 0)
+        c_in = g.add_vertex(VertexType.BOUNDARY, 1, 0)
+        q_mapper.add_label(0, 1)
+        q_mapper.set_prev_vertex(0, q_in)
+        c_mapper.add_label(0, 1)
+        c_mapper.set_qubit(0, 1)
+        c_mapper.set_prev_vertex(0, c_in)
+        # Call to_graph_ground directly.
+        m = Measurement(0, result_bit=0)
+        m.to_graph_ground(g, q_mapper, c_mapper)
+        # Should not crash, and the graph should have vertices.
+        self.assertGreater(len(list(g.vertices())), 2)
+
+    def test_hththt_t_injection_feedforward(self):
+        """HTHTHT... circuit with T gates implemented via injection.
+
+        This is the T-injection protocol example from tqec/tqec#708: prepare an
+        ancilla in |T>, CNOT with the data qubit, measure the ancilla, and
+        conditionally apply S-dagger. Three rounds are performed.
+
+        The ``if(c==1) sdg q[0];`` conditional corrections are not yet
+        supported by the parser (they require classical control / feedforward).
+        For now they are omitted, and the test verifies the structural
+        skeleton: reset, gates, measure, and QASM round-trip.
+        """
+        from pyzx.circuit.gates import Measurement, Reset
+
+        # Full circuit would include ``if(c==1) sdg q[0];`` after each
+        # measure, but that requires feedforward support (ConditionalGate).
+        qasm = """
+        OPENQASM 2.0;
+        include "qelib1.inc";
+        qreg q[2];
+        creg c[1];
+        h q[0];
+        reset q[1];
+        h q[1];
+        t q[1];
+        cx q[0],q[1];
+        h q[1];
+        measure q[1] -> c[0];
+        h q[0];
+        reset q[1];
+        h q[1];
+        t q[1];
+        cx q[0],q[1];
+        h q[1];
+        measure q[1] -> c[0];
+        h q[0];
+        reset q[1];
+        h q[1];
+        t q[1];
+        cx q[0],q[1];
+        h q[1];
+        measure q[1] -> c[0];
+        """
+
+        # Parse.
+        c = Circuit.from_qasm(qasm)
+        self.assertEqual(c.qubits, 2)
+
+        # Count gate types.
+        measurements = [g for g in c.gates if isinstance(g, Measurement)]
+        resets = [g for g in c.gates if isinstance(g, Reset)]
+        self.assertEqual(len(measurements), 3)
+        self.assertEqual(len(resets), 3)
+        # TODO: once feedforward is implemented, assert 3 ConditionalGate
+        # instances wrapping sdg on q[0], each conditioned on c==1.
+
+        # Convert to graph.
+        g = c.to_graph()
+        self.assertEqual(len(g.inputs()), 3)   # 2 qubits + 1 classical bit
+        self.assertEqual(len(g.outputs()), 2)  # q[0] + c[0]; q[1] consumed by measurement
+
+        # QASM round-trip.
+        qasm_out = c.to_qasm()
+        c2 = Circuit.from_qasm(qasm_out)
+        self.assertEqual(len(c.gates), len(c2.gates))
+        for g1, g2 in zip(c.gates, c2.gates):
+            self.assertEqual(type(g1), type(g2))
+
+
+@unittest.skipUnless(stim, "stim needs to be installed for this to run")
+class TestStimQASMInterop(unittest.TestCase):
+    """Test loading Stim-generated QASMs into PyZX.
+
+    Stim can export circuits to QASM via circuit.to_qasm().  These tests
+    verify that PyZX can parse the resulting QASM for a range of QEC
+    circuits.  See tqec/tqec#708 for context.
+
+    QASM 2 with skip_dets_and_obs=True is the primary interop path.
+    QASM 3 works when Stim does not use its MR subroutine definition.
+    """
+
+    def _load_stim_qasm2(self, stim_circuit):
+        """Export a Stim circuit to QASM 2, parse it in PyZX, and return
+        the PyZX Circuit."""
+        qasm = stim_circuit.to_qasm(
+            open_qasm_version=2, skip_dets_and_obs=True)
+        return Circuit.from_qasm(qasm)
+
+    def test_bell_pair_measurement(self):
+        c_stim = stim.Circuit("H 0\nCNOT 0 1\nM 0 1")
+        c = self._load_stim_qasm2(c_stim)
+        self.assertEqual(c.qubits, 2)
+        self.assertEqual(len(c.gates), 4)
+
+    def test_measure_reset_reuse(self):
+        c_stim = stim.Circuit("H 0\nM 0\nR 0\nH 0\nM 0")
+        c = self._load_stim_qasm2(c_stim)
+        self.assertEqual(c.qubits, 1)
+        self.assertEqual(len(c.gates), 5)
+
+    def test_entangle_measure_reset_cycle(self):
+        c_stim = stim.Circuit(
+            "H 0\nCNOT 0 1\nM 0 1\nR 0 1\n"
+            "H 0\nCNOT 0 1\nM 0 1")
+        c = self._load_stim_qasm2(c_stim)
+        self.assertEqual(c.qubits, 2)
+        self.assertEqual(len(c.gates), 10)
+
+    def test_mr_combined_gate(self):
+        """Stim's MR (combined measure-reset) decomposes to separate
+        measure + reset on a single line in QASM 2."""
+        c_stim = stim.Circuit("H 0\nCNOT 0 1\nMR 0\nH 0\nM 0 1")
+        c = self._load_stim_qasm2(c_stim)
+        self.assertEqual(c.qubits, 2)
+        self.assertEqual(len(c.gates), 7)
+
+    def test_surface_code_d3(self):
+        c_stim = stim.Circuit.generated(
+            "surface_code:rotated_memory_z",
+            rounds=2, distance=3,
+            after_clifford_depolarization=0,
+            after_reset_flip_probability=0,
+            before_measure_flip_probability=0,
+            before_round_data_depolarization=0)
+        c = self._load_stim_qasm2(c_stim)
+        # d=3 rotated surface code uses 9 data + 8 ancilla = 17 qubits,
+        # but Stim may allocate more. Just check it parses and round-trips.
+        self.assertGreater(c.qubits, 0)
+        g = c.to_graph()
+        self.assertGreater(len(list(g.vertices())), 0)
+        # QASM round-trip.
+        c2 = Circuit.from_qasm(c.to_qasm())
+        self.assertEqual(len(c.gates), len(c2.gates))
+
+    def test_surface_code_d5(self):
+        c_stim = stim.Circuit.generated(
+            "surface_code:rotated_memory_z",
+            rounds=3, distance=5,
+            after_clifford_depolarization=0,
+            after_reset_flip_probability=0,
+            before_measure_flip_probability=0,
+            before_round_data_depolarization=0)
+        c = self._load_stim_qasm2(c_stim)
+        self.assertGreater(c.qubits, 0)
+        g = c.to_graph()
+        self.assertGreater(len(list(g.vertices())), 0)
+
+    def test_repetition_code(self):
+        c_stim = stim.Circuit.generated(
+            "repetition_code:memory",
+            rounds=3, distance=3,
+            after_clifford_depolarization=0,
+            after_reset_flip_probability=0,
+            before_measure_flip_probability=0,
+            before_round_data_depolarization=0)
+        c = self._load_stim_qasm2(c_stim)
+        self.assertGreater(c.qubits, 0)
+        g = c.to_graph()
+        self.assertGreater(len(list(g.vertices())), 0)
+
+    def test_colour_code(self):
+        c_stim = stim.Circuit.generated(
+            "color_code:memory_xyz",
+            rounds=2, distance=3,
+            after_clifford_depolarization=0,
+            after_reset_flip_probability=0,
+            before_measure_flip_probability=0,
+            before_round_data_depolarization=0)
+        c = self._load_stim_qasm2(c_stim)
+        self.assertGreater(c.qubits, 0)
+        g = c.to_graph()
+        self.assertGreater(len(list(g.vertices())), 0)
+
+    def test_qasm3_without_mr(self):
+        """QASM 3 works when Stim uses separate M + R (not MR)."""
+        c_stim = stim.Circuit(
+            "H 0\nCNOT 0 1\nM 0\nR 0\nH 0\nM 0 1")
+        qasm3 = c_stim.to_qasm(
+            open_qasm_version=3, skip_dets_and_obs=True)
+        c = Circuit.from_qasm(qasm3)
+        self.assertEqual(c.qubits, 2)
+        self.assertEqual(len(c.gates), 7)
+
+    def test_qasm3_mr_subroutine_fails(self):
+        """QASM 3 with Stim's MR subroutine definition is not supported."""
+        c_stim = stim.Circuit("H 0\nMR 0\nH 0\nM 0")
+        qasm3 = c_stim.to_qasm(
+            open_qasm_version=3, skip_dets_and_obs=True)
+        # Stim emits `def mr(qubit q0) -> bit { ... }` which PyZX
+        # cannot parse.
+        with self.assertRaises((TypeError, ValueError)):
+            Circuit.from_qasm(qasm3)
+
 
 if __name__ == '__main__':
     unittest.main()
diff --git a/tests/test_symbolic_parsing.py b/tests/test_symbolic_parsing.py
index ef8488feb..44b5fd57b 100644
--- a/tests/test_symbolic_parsing.py
+++ b/tests/test_symbolic_parsing.py
@@ -210,6 +210,27 @@ def test_dot_multiplication(self):
         expected = self.new_var("x") * (self.new_var("y") + self.new_var("z"))
         self.assertEqual(result, expected)
 
+    def test_subscripted_variables(self):
+        """Test that subscripted variable names like c[0] can be parsed."""
+        result = parse("c[0]", self.new_var)
+        expected = self.new_var("c[0]")
+        self.assertEqual(result, expected)
+
+        # Multiple subscripted variables.
+        result = parse("c[0] + c[1]", self.new_var)
+        expected = self.new_var("c[0]") + self.new_var("c[1]")
+        self.assertEqual(result, expected)
+
+        # Mixed plain and subscripted variables.
+        result = parse("x + c[2]", self.new_var)
+        expected = self.new_var("x") + self.new_var("c[2]")
+        self.assertEqual(result, expected)
+
+        # Coefficient with subscripted variable.
+        result = parse("3*c[0]", self.new_var)
+        expected = new_const(3) * self.new_var("c[0]")
+        self.assertEqual(result, expected)
+
     def test_complex_expressions(self):
         """Test complex mathematical expressions combining all features."""
         # Test polynomial expression