diff --git a/cirq-core/cirq/__init__.py b/cirq-core/cirq/__init__.py
index 46d8f558205..22dcdc8768f 100644
--- a/cirq-core/cirq/__init__.py
+++ b/cirq-core/cirq/__init__.py
@@ -441,6 +441,7 @@
 
 from cirq.sim import (
     CIRCUIT_LIKE,
+    ClassicalStateSimulator,
     CliffordSimulator,
     CliffordState,
     CliffordSimulatorStepResult,
diff --git a/cirq-core/cirq/protocols/json_test_data/spec.py b/cirq-core/cirq/protocols/json_test_data/spec.py
index 05aeaa6be86..5b32c638086 100644
--- a/cirq-core/cirq/protocols/json_test_data/spec.py
+++ b/cirq-core/cirq/protocols/json_test_data/spec.py
@@ -58,6 +58,7 @@
         'ZerosSampler',
     ],
     should_not_be_serialized=[
+        'ClassicalStateSimulator',
         # Heatmaps
         'Heatmap',
         'TwoQubitInteractionHeatmap',
diff --git a/cirq-core/cirq/sim/__init__.py b/cirq-core/cirq/sim/__init__.py
index 57fc6f1d8ae..980b649db7a 100644
--- a/cirq-core/cirq/sim/__init__.py
+++ b/cirq-core/cirq/sim/__init__.py
@@ -68,6 +68,8 @@
 
 from cirq.sim.state_vector_simulation_state import StateVectorSimulationState
 
+from cirq.sim.classical_simulator import ClassicalStateSimulator
+
 from cirq.sim.state_vector_simulator import (
     SimulatesIntermediateStateVector,
     StateVectorStepResult,
diff --git a/cirq-core/cirq/sim/classical_simulator.py b/cirq-core/cirq/sim/classical_simulator.py
new file mode 100644
index 00000000000..80a478f3eb2
--- /dev/null
+++ b/cirq-core/cirq/sim/classical_simulator.py
@@ -0,0 +1,127 @@
+# Copyright 2023 The Cirq Developers
+#
+# 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 Dict
+from collections import defaultdict
+from cirq.sim.simulator import SimulatesSamples
+from cirq import ops, protocols
+from cirq.study.resolver import ParamResolver
+from cirq.circuits.circuit import AbstractCircuit
+from cirq.ops.raw_types import Qid
+import numpy as np
+
+
+class ClassicalStateSimulator(SimulatesSamples):
+    """A simulator that only accepts only gates with classical counterparts.
+
+    This simulator evolves a single state, using only gates that output a single state for each
+    input state. The simulator runs in linear time, at the cost of not supporting superposition.
+    It can be used to estimate costs and simulate circuits for simple non-quantum algorithms using
+    many more qubits than fully capable quantum simulators.
+
+    The supported gates are:
+        - cirq.X
+        - cirq.CNOT
+        - cirq.SWAP
+        - cirq.TOFFOLI
+        - cirq.measure
+
+    Args:
+        circuit: The circuit to simulate.
+        param_resolver: Parameters to run with the program.
+        repetitions: Number of times to repeat the run. It is expected that
+            this is validated greater than zero before calling this method.
+
+    Returns:
+        A dictionary mapping measurement keys to measurement results.
+
+    Raises:
+        ValueError: If one of the gates is not an X, CNOT, SWAP, TOFFOLI or a measurement.
+    """
+
+    def _run(
+        self, circuit: AbstractCircuit, param_resolver: ParamResolver, repetitions: int
+    ) -> Dict[str, np.ndarray]:
+        results_dict: Dict[str, np.ndarray] = {}
+        values_dict: Dict[Qid, int] = defaultdict(int)
+        param_resolver = param_resolver or ParamResolver({})
+        resolved_circuit = protocols.resolve_parameters(circuit, param_resolver)
+
+        for moment in resolved_circuit:
+            for op in moment:
+                gate = op.gate
+                if gate == ops.X:
+                    values_dict[op.qubits[0]] = 1 - values_dict[op.qubits[0]]
+
+                elif (
+                    isinstance(gate, ops.CNotPowGate)
+                    and gate.exponent == 1
+                    and gate.global_shift == 0
+                ):
+                    if values_dict[op.qubits[0]] == 1:
+                        values_dict[op.qubits[1]] = 1 - values_dict[op.qubits[1]]
+
+                elif (
+                    isinstance(gate, ops.SwapPowGate)
+                    and gate.exponent == 1
+                    and gate.global_shift == 0
+                ):
+                    hold_qubit = values_dict[op.qubits[1]]
+                    values_dict[op.qubits[1]] = values_dict[op.qubits[0]]
+                    values_dict[op.qubits[0]] = hold_qubit
+
+                elif (
+                    isinstance(gate, ops.CCXPowGate)
+                    and gate.exponent == 1
+                    and gate.global_shift == 0
+                ):
+                    if (values_dict[op.qubits[0]] == 1) and (values_dict[op.qubits[1]] == 1):
+                        values_dict[op.qubits[2]] = 1 - values_dict[op.qubits[2]]
+
+                elif isinstance(gate, ops.MeasurementGate):
+                    qubits_in_order = op.qubits
+                    # add the new instance of a key to the numpy array in results dictionary
+                    if gate.key in results_dict:
+                        shape = len(qubits_in_order)
+                        current_array = results_dict[gate.key]
+                        new_instance = np.zeros(shape, dtype=np.uint8)
+                        for bits in range(0, len(qubits_in_order)):
+                            new_instance[bits] = values_dict[qubits_in_order[bits]]
+                            results_dict[gate.key] = np.insert(
+                                current_array, len(current_array[0]), new_instance, axis=1
+                            )
+                    else:
+                        # create the array for the results dictionary
+                        new_array_shape = (repetitions, 1, len(qubits_in_order))
+                        new_array = np.zeros(new_array_shape, dtype=np.uint8)
+                        for reps in range(0, repetitions):
+                            for instances in range(1):
+                                for bits in range(0, len(qubits_in_order)):
+                                    new_array[reps][instances][bits] = values_dict[
+                                        qubits_in_order[bits]
+                                    ]
+                        results_dict[gate.key] = new_array
+
+                elif not (
+                    (isinstance(gate, ops.XPowGate) and gate.exponent == 0)
+                    or (isinstance(gate, ops.CCXPowGate) and gate.exponent == 0)
+                    or (isinstance(gate, ops.SwapPowGate) and gate.exponent == 0)
+                    or (isinstance(gate, ops.CNotPowGate) and gate.exponent == 0)
+                ):
+                    raise ValueError(
+                        "Can not simulate gates other than cirq.XGate, "
+                        + "cirq.CNOT, cirq.SWAP, and cirq.CCNOT"
+                    )
+
+        return results_dict
diff --git a/cirq-core/cirq/sim/classical_simulator_test.py b/cirq-core/cirq/sim/classical_simulator_test.py
new file mode 100644
index 00000000000..55f8d72b5c8
--- /dev/null
+++ b/cirq-core/cirq/sim/classical_simulator_test.py
@@ -0,0 +1,196 @@
+# Copyright 2023 The Cirq Developers
+#
+# 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.
+import numpy as np
+import pytest
+import cirq
+import sympy
+
+
+class TestSimulator:
+    def test_x_gate(self):
+        q0, q1 = cirq.LineQubit.range(2)
+        circuit = cirq.Circuit()
+        circuit.append(cirq.X(q0))
+        circuit.append(cirq.X(q1))
+        circuit.append(cirq.X(q1))
+        circuit.append(cirq.measure((q0, q1), key='key'))
+        expected_results = {'key': np.array([[[1, 0]]], dtype=np.uint8)}
+        sim = cirq.ClassicalStateSimulator()
+        results = sim.run(circuit, param_resolver=None, repetitions=1).records
+        np.testing.assert_equal(results, expected_results)
+
+    def test_CNOT(self):
+        q0, q1 = cirq.LineQubit.range(2)
+        circuit = cirq.Circuit()
+        circuit.append(cirq.X(q0))
+        circuit.append(cirq.CNOT(q0, q1))
+        circuit.append(cirq.measure(q1, key='key'))
+        expected_results = {'key': np.array([[[1]]], dtype=np.uint8)}
+        sim = cirq.ClassicalStateSimulator()
+        results = sim.run(circuit, param_resolver=None, repetitions=1).records
+        np.testing.assert_equal(results, expected_results)
+
+    def test_Swap(self):
+        q0, q1 = cirq.LineQubit.range(2)
+        circuit = cirq.Circuit()
+        circuit.append(cirq.X(q0))
+        circuit.append(cirq.SWAP(q0, q1))
+        circuit.append(cirq.measure((q0, q1), key='key'))
+        expected_results = {'key': np.array([[[0, 1]]], dtype=np.uint8)}
+        sim = cirq.ClassicalStateSimulator()
+        results = sim.run(circuit, param_resolver=None, repetitions=1).records
+        np.testing.assert_equal(results, expected_results)
+
+    def test_CCNOT(self):
+        q0, q1, q2 = cirq.LineQubit.range(3)
+        circuit = cirq.Circuit()
+        circuit.append(cirq.CCNOT(q0, q1, q2))
+        circuit.append(cirq.measure((q0, q1, q2), key='key'))
+        circuit.append(cirq.X(q0))
+        circuit.append(cirq.CCNOT(q0, q1, q2))
+        circuit.append(cirq.measure((q0, q1, q2), key='key'))
+        circuit.append(cirq.X(q1))
+        circuit.append(cirq.X(q0))
+        circuit.append(cirq.CCNOT(q0, q1, q2))
+        circuit.append(cirq.measure((q0, q1, q2), key='key'))
+        circuit.append(cirq.X(q0))
+        circuit.append(cirq.CCNOT(q0, q1, q2))
+        circuit.append(cirq.measure((q0, q1, q2), key='key'))
+        expected_results = {
+            'key': np.array([[[0, 0, 0], [1, 0, 0], [0, 1, 0], [1, 1, 1]]], dtype=np.uint8)
+        }
+        sim = cirq.ClassicalStateSimulator()
+        results = sim.run(circuit, param_resolver=None, repetitions=1).records
+        np.testing.assert_equal(results, expected_results)
+
+    def test_measurement_gate(self):
+        q0, q1 = cirq.LineQubit.range(2)
+        circuit = cirq.Circuit()
+        circuit.append(cirq.measure((q0, q1), key='key'))
+        expected_results = {'key': np.array([[[0, 0]]], dtype=np.uint8)}
+        sim = cirq.ClassicalStateSimulator()
+        results = sim.run(circuit, param_resolver=None, repetitions=1).records
+        np.testing.assert_equal(results, expected_results)
+
+    def test_qubit_order(self):
+        q0, q1 = cirq.LineQubit.range(2)
+        circuit = cirq.Circuit()
+        circuit.append(cirq.CNOT(q0, q1))
+        circuit.append(cirq.X(q0))
+        circuit.append(cirq.measure((q0, q1), key='key'))
+        expected_results = {'key': np.array([[[1, 0]]], dtype=np.uint8)}
+        sim = cirq.ClassicalStateSimulator()
+        results = sim.run(circuit, param_resolver=None, repetitions=1).records
+        np.testing.assert_equal(results, expected_results)
+
+    def test_same_key_instances(self):
+        q0, q1 = cirq.LineQubit.range(2)
+        circuit = cirq.Circuit()
+        circuit.append(cirq.measure((q0, q1), key='key'))
+        circuit.append(cirq.X(q0))
+        circuit.append(cirq.measure((q0, q1), key='key'))
+        expected_results = {'key': np.array([[[0, 0], [1, 0]]], dtype=np.uint8)}
+        sim = cirq.ClassicalStateSimulator()
+        results = sim.run(circuit, param_resolver=None, repetitions=1).records
+        np.testing.assert_equal(results, expected_results)
+
+    def test_same_key_instances_order(self):
+        q0, q1 = cirq.LineQubit.range(2)
+        circuit = cirq.Circuit()
+        circuit.append(cirq.X(q0))
+        circuit.append(cirq.measure((q0, q1), key='key'))
+        circuit.append(cirq.X(q0))
+        circuit.append(cirq.measure((q1, q0), key='key'))
+        expected_results = {'key': np.array([[[1, 0], [0, 0]]], dtype=np.uint8)}
+        sim = cirq.ClassicalStateSimulator()
+        results = sim.run(circuit, param_resolver=None, repetitions=1).records
+        np.testing.assert_equal(results, expected_results)
+
+    def test_repetitions(self):
+        q0 = cirq.LineQubit.range(1)
+        circuit = cirq.Circuit()
+        circuit.append(cirq.measure(q0, key='key'))
+        expected_results = {
+            'key': np.array(
+                [[[0]], [[0]], [[0]], [[0]], [[0]], [[0]], [[0]], [[0]], [[0]], [[0]]],
+                dtype=np.uint8,
+            )
+        }
+        sim = cirq.ClassicalStateSimulator()
+        results = sim.run(circuit, param_resolver=None, repetitions=10).records
+        np.testing.assert_equal(results, expected_results)
+
+    def test_multiple_gates(self):
+        q0, q1 = cirq.LineQubit.range(2)
+        circuit = cirq.Circuit()
+        circuit.append(cirq.X(q0))
+        circuit.append(cirq.CNOT(q0, q1))
+        circuit.append(cirq.CNOT(q0, q1))
+        circuit.append(cirq.CNOT(q0, q1))
+        circuit.append(cirq.X(q1))
+        circuit.append(cirq.measure((q0, q1), key='key'))
+        expected_results = {'key': np.array([[[1, 0]]], dtype=np.uint8)}
+        sim = cirq.ClassicalStateSimulator()
+        results = sim.run(circuit, param_resolver=None, repetitions=1).records
+        np.testing.assert_equal(results, expected_results)
+
+    def test_multiple_gates_order(self):
+        q0, q1 = cirq.LineQubit.range(2)
+        circuit = cirq.Circuit()
+        circuit.append(cirq.X(q0))
+        circuit.append(cirq.CNOT(q0, q1))
+        circuit.append(cirq.CNOT(q1, q0))
+        circuit.append(cirq.measure((q0, q1), key='key'))
+        expected_results = {'key': np.array([[[0, 1]]], dtype=np.uint8)}
+        sim = cirq.ClassicalStateSimulator()
+        results = sim.run(circuit, param_resolver=None, repetitions=1).records
+        np.testing.assert_equal(results, expected_results)
+
+    def test_param_resolver(self):
+        gate = cirq.CNOT ** sympy.Symbol('t')
+        q0, q1 = cirq.LineQubit.range(2)
+        circuit = cirq.Circuit()
+        circuit.append(cirq.X(q0))
+        circuit.append(gate(q0, q1))
+        circuit.append(cirq.measure((q1), key='key'))
+        resolver = cirq.ParamResolver({'t': 0})
+        sim = cirq.ClassicalStateSimulator()
+        results_with_parameter_zero = sim.run(
+            circuit, param_resolver=resolver, repetitions=1
+        ).records
+        resolver = cirq.ParamResolver({'t': 1})
+        results_with_parameter_one = sim.run(
+            circuit, param_resolver=resolver, repetitions=1
+        ).records
+        np.testing.assert_equal(
+            results_with_parameter_zero, {'key': np.array([[[0]]], dtype=np.uint8)}
+        )
+        np.testing.assert_equal(
+            results_with_parameter_one, {'key': np.array([[[1]]], dtype=np.uint8)}
+        )
+
+    def test_unknown_gates(self):
+        gate = cirq.CNOT ** sympy.Symbol('t')
+        q0, q1 = cirq.LineQubit.range(2)
+        circuit = cirq.Circuit()
+        circuit.append(gate(q0, q1))
+        circuit.append(cirq.measure((q0), key='key'))
+        resolver = cirq.ParamResolver({'t': 0.5})
+        sim = cirq.ClassicalStateSimulator()
+        with pytest.raises(
+            ValueError,
+            match="Can not simulate gates other than "
+            + "cirq.XGate, cirq.CNOT, cirq.SWAP, and cirq.CCNOT",
+        ):
+            _ = sim.run(circuit, param_resolver=resolver, repetitions=1).records