From c111f57e91a58bc7c88bfa0e4567e01ed5ec59bc Mon Sep 17 00:00:00 2001
From: Aaryan Mehta <73230976+blazingphoenix7@users.noreply.github.com>
Date: Fri, 12 Jun 2026 23:53:30 -0400
Subject: [PATCH 1/4] Accept NumPy scalar coefficients across operator and
 representation types

NumPy scalar types (numpy.int64, numpy.float32, numpy.complex64, and so on)
are not subclasses of Python's int/float/complex, so coefficient checks of
the form isinstance(x, (int, float, complex)) reject them. That is the root
of issue #1097. The SymbolicOperator family (FermionOperator, QubitOperator,
...) was already fixed by adding numbers.Number to its accepted coefficient
types; this completes the same fix for the remaining types that validate
coefficients the same way:

  - PolynomialTensor (and InteractionOperator and friends)
  - DOCIHamiltonian
  - MajoranaOperator
  - the majorana_operator builder in special_operators

numbers.Number matches Python and NumPy numeric scalars alike, so an accepted
coefficient now behaves identically whether it comes from Python or NumPy.

Adds a test in each affected module checking that NumPy scalar coefficients
give the same result as their Python equivalents.

Fixes #1097.
---
 .../hamiltonians/special_operators.py         |  3 ++-
 .../hamiltonians/special_operators_test.py    | 10 ++++++++
 .../ops/operators/majorana_operator.py        | 24 +++++++++++--------
 .../ops/operators/majorana_operator_test.py   | 24 +++++++++++++++++++
 .../ops/representations/doci_hamiltonian.py   |  4 +++-
 .../representations/doci_hamiltonian_test.py  | 16 +++++++++++++
 .../ops/representations/polynomial_tensor.py  |  3 ++-
 .../representations/polynomial_tensor_test.py | 17 +++++++++++++
 8 files changed, 88 insertions(+), 13 deletions(-)

diff --git a/src/openfermion/hamiltonians/special_operators.py b/src/openfermion/hamiltonians/special_operators.py
index 35732c658..5893556ef 100644
--- a/src/openfermion/hamiltonians/special_operators.py
+++ b/src/openfermion/hamiltonians/special_operators.py
@@ -11,6 +11,7 @@
 #   limitations under the License.
 """Commonly used operators (mainly instances of SymbolicOperator)."""
 
+import numbers
 from typing import Optional, Union, Tuple
 
 from openfermion.ops.operators import BosonOperator, FermionOperator
@@ -236,7 +237,7 @@ def majorana_operator(
     Returns:
         FermionOperator
     """
-    if not isinstance(coefficient, (int, float, complex)):
+    if not isinstance(coefficient, (int, float, complex, numbers.Number)):
         raise ValueError('Coefficient must be scalar.')
 
     # If term is a string, convert it to a tuple
diff --git a/src/openfermion/hamiltonians/special_operators_test.py b/src/openfermion/hamiltonians/special_operators_test.py
index ba528fcd8..c4d126093 100644
--- a/src/openfermion/hamiltonians/special_operators_test.py
+++ b/src/openfermion/hamiltonians/special_operators_test.py
@@ -12,6 +12,9 @@
 """testing angular momentum generators. _fermion_spin_operators.py"""
 
 import unittest
+
+import numpy
+
 from openfermion.ops.operators import FermionOperator, BosonOperator
 from openfermion.utils import commutator
 from openfermion.transforms.opconversions import normal_ordered
@@ -204,3 +207,10 @@ def test_bad_term(self):
             majorana_operator('a')
         with self.assertRaises(ValueError):
             majorana_operator(2)
+
+
+def test_majorana_operator_builder_numpy_scalar_coefficient():
+    """The majorana_operator builder accepts NumPy scalar coefficients (issue #1097)."""
+    cases = [(numpy.int64(2), 2), (numpy.float32(0.5), 0.5), (numpy.complex64(1 + 2j), 1 + 2j)]
+    for numpy_scalar, python_scalar in cases:
+        assert majorana_operator((1, 0), numpy_scalar) == majorana_operator((1, 0), python_scalar)
diff --git a/src/openfermion/ops/operators/majorana_operator.py b/src/openfermion/ops/operators/majorana_operator.py
index 6619f1ba7..2af867851 100644
--- a/src/openfermion/ops/operators/majorana_operator.py
+++ b/src/openfermion/ops/operators/majorana_operator.py
@@ -13,8 +13,12 @@
 
 import copy
 import itertools
+import numbers
+
 import numpy
 
+COEFFICIENT_TYPES = (int, float, complex, numbers.Number)
+
 
 class MajoranaOperator:
     r"""A linear combination of products of Majorana operators.
@@ -78,7 +82,7 @@ def from_dict(terms):
 
     def commutes_with(self, other):
         """Test commutation with another MajoranaOperator"""
-        if isinstance(other, (int, float, complex)):
+        if isinstance(other, COEFFICIENT_TYPES):
             return True
 
         if not isinstance(other, type(self)):
@@ -145,7 +149,7 @@ def __iadd__(self, other):
                     self.terms[term] += coefficient
                 else:
                     self.terms[term] = coefficient
-        elif isinstance(other, (int, float, complex)):
+        elif isinstance(other, COEFFICIENT_TYPES):
             self.constant += other
         else:
             raise TypeError("Cannot add invalid type to {}".format(type(self)))
@@ -163,7 +167,7 @@ def __isub__(self, other):
                     self.terms[term] -= coefficient
                 else:
                     self.terms[term] = -coefficient
-        elif isinstance(other, (int, float, complex)):
+        elif isinstance(other, COEFFICIENT_TYPES):
             self.constant -= other
         else:
             raise TypeError("Cannot subtract invalid type from {}".format(type(self)))
@@ -175,10 +179,10 @@ def __sub__(self, other):
         return minuend
 
     def __mul__(self, other):
-        if not isinstance(other, (type(self), int, float, complex)):
+        if not isinstance(other, (type(self), *COEFFICIENT_TYPES)):
             return NotImplemented
 
-        if isinstance(other, (int, float, complex)):
+        if isinstance(other, COEFFICIENT_TYPES):
             terms = {term: coefficient * other for term, coefficient in self.terms.items()}
             return MajoranaOperator.from_dict(terms)
 
@@ -194,10 +198,10 @@ def __mul__(self, other):
         return MajoranaOperator.from_dict(terms)
 
     def __imul__(self, other):
-        if not isinstance(other, (type(self), int, float, complex)):
+        if not isinstance(other, (type(self), *COEFFICIENT_TYPES)):
             return NotImplemented
 
-        if isinstance(other, (int, float, complex)):
+        if isinstance(other, COEFFICIENT_TYPES):
             for term in self.terms:
                 self.terms[term] *= other
             return self
@@ -205,19 +209,19 @@ def __imul__(self, other):
         return self * other
 
     def __rmul__(self, other):
-        if not isinstance(other, (int, float, complex)):
+        if not isinstance(other, COEFFICIENT_TYPES):
             return NotImplemented
         return self * other
 
     def __truediv__(self, other):
-        if not isinstance(other, (int, float, complex)):
+        if not isinstance(other, COEFFICIENT_TYPES):
             return NotImplemented
 
         terms = {term: coefficient / other for term, coefficient in self.terms.items()}
         return MajoranaOperator.from_dict(terms)
 
     def __itruediv__(self, other):
-        if not isinstance(other, (int, float, complex)):
+        if not isinstance(other, COEFFICIENT_TYPES):
             return NotImplemented
 
         for term in self.terms:
diff --git a/src/openfermion/ops/operators/majorana_operator_test.py b/src/openfermion/ops/operators/majorana_operator_test.py
index ad572c389..b0e0f737a 100644
--- a/src/openfermion/ops/operators/majorana_operator_test.py
+++ b/src/openfermion/ops/operators/majorana_operator_test.py
@@ -238,3 +238,27 @@ def test_majorana_operator_str():
 def test_majorana_operator_repr():
     a = MajoranaOperator((0, 1, 5), 1.5)
     assert repr(a) == 'MajoranaOperator.from_dict(terms={(0, 1, 5): 1.5})'
+
+
+def test_majorana_operator_numpy_scalar_coefficients():
+    """NumPy scalar coefficients behave like Python scalars (issue #1097)."""
+    op = MajoranaOperator((0, 1), 1.0) + MajoranaOperator((2, 3), 2.0)
+    cases = [(numpy.int64(2), 2), (numpy.float32(0.5), 0.5), (numpy.complex64(1 + 2j), 1 + 2j)]
+    for numpy_scalar, python_scalar in cases:
+        assert op * numpy_scalar == op * python_scalar
+        assert numpy_scalar * op == python_scalar * op
+        assert op + numpy_scalar == op + python_scalar
+        assert op - numpy_scalar == op - python_scalar
+        assert op / numpy_scalar == op / python_scalar
+    assert op.commutes_with(numpy.int64(5))
+
+    # In-place operators, using exact values so the comparison is not subject
+    # to float32 round-off across the sequence of operations.
+    for numpy_scalar, python_scalar in [(numpy.int64(2), 2), (numpy.float32(0.5), 0.5)]:
+        numpy_op = MajoranaOperator((0, 1), 1.0) + MajoranaOperator((2, 3), 2.0)
+        python_op = MajoranaOperator((0, 1), 1.0) + MajoranaOperator((2, 3), 2.0)
+        numpy_op *= numpy_scalar
+        python_op *= python_scalar
+        numpy_op /= numpy_scalar
+        python_op /= python_scalar
+        assert numpy_op == python_op
diff --git a/src/openfermion/ops/representations/doci_hamiltonian.py b/src/openfermion/ops/representations/doci_hamiltonian.py
index 1587aeb6c..6fabf1588 100644
--- a/src/openfermion/ops/representations/doci_hamiltonian.py
+++ b/src/openfermion/ops/representations/doci_hamiltonian.py
@@ -11,12 +11,14 @@
 #   limitations under the License.
 """Class and functions to store interaction operators."""
 
+import numbers
+
 import numpy
 
 from openfermion.ops import QubitOperator
 from openfermion.ops.representations import PolynomialTensor, get_tensors_from_integrals
 
-COEFFICIENT_TYPES = (int, float, complex)
+COEFFICIENT_TYPES = (int, float, complex, numbers.Number)
 
 
 class DOCIHamiltonian(PolynomialTensor):
diff --git a/src/openfermion/ops/representations/doci_hamiltonian_test.py b/src/openfermion/ops/representations/doci_hamiltonian_test.py
index 70fa603b3..44f4ed1fb 100644
--- a/src/openfermion/ops/representations/doci_hamiltonian_test.py
+++ b/src/openfermion/ops/representations/doci_hamiltonian_test.py
@@ -277,3 +277,19 @@ def test_from_integrals_to_qubit(self):
             + "Hamiltonian\n"
             + str(sub_matrix)
         )
+
+
+def test_doci_numpy_scalar_coefficients():
+    """NumPy scalar coefficients behave like Python scalars (issue #1097)."""
+    doci = DOCIHamiltonian(
+        1.0,
+        numpy.array([1.0, 2.0]),
+        numpy.array([[0.0, 0.5], [0.5, 0.0]]),
+        numpy.array([[0.0, 0.3], [0.3, 0.0]]),
+    )
+    cases = [(numpy.int64(2), 2), (numpy.float32(0.5), 0.5)]
+    for numpy_scalar, python_scalar in cases:
+        assert doci * numpy_scalar == doci * python_scalar
+        assert doci + numpy_scalar == doci + python_scalar
+        assert doci - numpy_scalar == doci - python_scalar
+        assert doci / numpy_scalar == doci / python_scalar
diff --git a/src/openfermion/ops/representations/polynomial_tensor.py b/src/openfermion/ops/representations/polynomial_tensor.py
index 523e55c52..ded9f26b6 100644
--- a/src/openfermion/ops/representations/polynomial_tensor.py
+++ b/src/openfermion/ops/representations/polynomial_tensor.py
@@ -14,13 +14,14 @@
 
 import copy
 import itertools
+import numbers
 import operator
 
 import numpy
 
 from openfermion.config import EQ_TOLERANCE
 
-COEFFICIENT_TYPES = (int, float, complex)
+COEFFICIENT_TYPES = (int, float, complex, numbers.Number)
 
 
 class PolynomialTensorError(Exception):
diff --git a/src/openfermion/ops/representations/polynomial_tensor_test.py b/src/openfermion/ops/representations/polynomial_tensor_test.py
index 8b166111b..7a98edbc6 100644
--- a/src/openfermion/ops/representations/polynomial_tensor_test.py
+++ b/src/openfermion/ops/representations/polynomial_tensor_test.py
@@ -536,3 +536,20 @@ def do_rotate_basis_high_order(self, order):
         polynomial_tensor.rotate_basis(numpy.array([[rotation]]))
 
         return polynomial_tensor, want_polynomial_tensor
+
+
+def test_numpy_scalar_coefficients():
+    """NumPy scalar coefficients behave like Python scalars (issue #1097)."""
+    tensor = PolynomialTensor(
+        {
+            (): 1.0,
+            (1, 0): numpy.array([[1.0, 2.0], [3.0, 4.0]]),
+            (1, 1, 0, 0): numpy.arange(16, dtype=float).reshape((2, 2, 2, 2)),
+        }
+    )
+    cases = [(numpy.int64(2), 2), (numpy.int32(3), 3), (numpy.float32(0.5), 0.5)]
+    for numpy_scalar, python_scalar in cases:
+        assert tensor * numpy_scalar == tensor * python_scalar
+        assert tensor + numpy_scalar == tensor + python_scalar
+        assert tensor - numpy_scalar == tensor - python_scalar
+        assert tensor / numpy_scalar == tensor / python_scalar

From 9355ad2459267bde3f8cb62da21cec16dd0bf0ee Mon Sep 17 00:00:00 2001
From: Aaryan Mehta <73230976+blazingphoenix7@users.noreply.github.com>
Date: Sat, 13 Jun 2026 14:54:12 -0400
Subject: [PATCH 2/4] Simplify MajoranaOperator type checks and move new tests
 into TestCase classes

Two changes in response to the automated review on this PR:

- MajoranaOperator.__mul__ and __imul__ unpacked COEFFICIENT_TYPES into a fresh
  tuple on every call via (type(self), *COEFFICIENT_TYPES). Switched to a single
  isinstance with a nested tuple, (type(self), COEFFICIENT_TYPES), which avoids
  the per-call unpacking and stays a single isinstance call (so pylint's
  consider-merging-isinstance check remains satisfied).

- The three new NumPy-scalar tests were module-level functions in files whose
  other tests live in unittest.TestCase classes, so they would be skipped under
  plain unittest discovery. Moved them into PolynomialTensorTest,
  DOCIHamiltonianTest, and MajoranaOperatorTest and switched to self.assertEqual.
---
 .../hamiltonians/special_operators_test.py    | 13 ++++----
 .../ops/operators/majorana_operator.py        |  4 +--
 .../representations/doci_hamiltonian_test.py  | 29 +++++++++--------
 .../representations/polynomial_tensor_test.py | 31 +++++++++----------
 4 files changed, 38 insertions(+), 39 deletions(-)

diff --git a/src/openfermion/hamiltonians/special_operators_test.py b/src/openfermion/hamiltonians/special_operators_test.py
index c4d126093..94aa3e42f 100644
--- a/src/openfermion/hamiltonians/special_operators_test.py
+++ b/src/openfermion/hamiltonians/special_operators_test.py
@@ -208,9 +208,10 @@ def test_bad_term(self):
         with self.assertRaises(ValueError):
             majorana_operator(2)
 
-
-def test_majorana_operator_builder_numpy_scalar_coefficient():
-    """The majorana_operator builder accepts NumPy scalar coefficients (issue #1097)."""
-    cases = [(numpy.int64(2), 2), (numpy.float32(0.5), 0.5), (numpy.complex64(1 + 2j), 1 + 2j)]
-    for numpy_scalar, python_scalar in cases:
-        assert majorana_operator((1, 0), numpy_scalar) == majorana_operator((1, 0), python_scalar)
+    def test_builder_numpy_scalar_coefficient(self):
+        """The majorana_operator builder accepts NumPy scalar coefficients (issue #1097)."""
+        cases = [(numpy.int64(2), 2), (numpy.float32(0.5), 0.5), (numpy.complex64(1 + 2j), 1 + 2j)]
+        for numpy_scalar, python_scalar in cases:
+            self.assertEqual(
+                majorana_operator((1, 0), numpy_scalar), majorana_operator((1, 0), python_scalar)
+            )
diff --git a/src/openfermion/ops/operators/majorana_operator.py b/src/openfermion/ops/operators/majorana_operator.py
index 2af867851..7e6724013 100644
--- a/src/openfermion/ops/operators/majorana_operator.py
+++ b/src/openfermion/ops/operators/majorana_operator.py
@@ -179,7 +179,7 @@ def __sub__(self, other):
         return minuend
 
     def __mul__(self, other):
-        if not isinstance(other, (type(self), *COEFFICIENT_TYPES)):
+        if not isinstance(other, (type(self), COEFFICIENT_TYPES)):
             return NotImplemented
 
         if isinstance(other, COEFFICIENT_TYPES):
@@ -198,7 +198,7 @@ def __mul__(self, other):
         return MajoranaOperator.from_dict(terms)
 
     def __imul__(self, other):
-        if not isinstance(other, (type(self), *COEFFICIENT_TYPES)):
+        if not isinstance(other, (type(self), COEFFICIENT_TYPES)):
             return NotImplemented
 
         if isinstance(other, COEFFICIENT_TYPES):
diff --git a/src/openfermion/ops/representations/doci_hamiltonian_test.py b/src/openfermion/ops/representations/doci_hamiltonian_test.py
index 44f4ed1fb..4f2aa2ded 100644
--- a/src/openfermion/ops/representations/doci_hamiltonian_test.py
+++ b/src/openfermion/ops/representations/doci_hamiltonian_test.py
@@ -278,18 +278,17 @@ def test_from_integrals_to_qubit(self):
             + str(sub_matrix)
         )
 
-
-def test_doci_numpy_scalar_coefficients():
-    """NumPy scalar coefficients behave like Python scalars (issue #1097)."""
-    doci = DOCIHamiltonian(
-        1.0,
-        numpy.array([1.0, 2.0]),
-        numpy.array([[0.0, 0.5], [0.5, 0.0]]),
-        numpy.array([[0.0, 0.3], [0.3, 0.0]]),
-    )
-    cases = [(numpy.int64(2), 2), (numpy.float32(0.5), 0.5)]
-    for numpy_scalar, python_scalar in cases:
-        assert doci * numpy_scalar == doci * python_scalar
-        assert doci + numpy_scalar == doci + python_scalar
-        assert doci - numpy_scalar == doci - python_scalar
-        assert doci / numpy_scalar == doci / python_scalar
+    def test_numpy_scalar_coefficients(self):
+        """NumPy scalar coefficients behave like Python scalars (issue #1097)."""
+        doci = DOCIHamiltonian(
+            1.0,
+            numpy.array([1.0, 2.0]),
+            numpy.array([[0.0, 0.5], [0.5, 0.0]]),
+            numpy.array([[0.0, 0.3], [0.3, 0.0]]),
+        )
+        cases = [(numpy.int64(2), 2), (numpy.float32(0.5), 0.5)]
+        for numpy_scalar, python_scalar in cases:
+            self.assertEqual(doci * numpy_scalar, doci * python_scalar)
+            self.assertEqual(doci + numpy_scalar, doci + python_scalar)
+            self.assertEqual(doci - numpy_scalar, doci - python_scalar)
+            self.assertEqual(doci / numpy_scalar, doci / python_scalar)
diff --git a/src/openfermion/ops/representations/polynomial_tensor_test.py b/src/openfermion/ops/representations/polynomial_tensor_test.py
index 7a98edbc6..7ee4a31b3 100644
--- a/src/openfermion/ops/representations/polynomial_tensor_test.py
+++ b/src/openfermion/ops/representations/polynomial_tensor_test.py
@@ -537,19 +537,18 @@ def do_rotate_basis_high_order(self, order):
 
         return polynomial_tensor, want_polynomial_tensor
 
-
-def test_numpy_scalar_coefficients():
-    """NumPy scalar coefficients behave like Python scalars (issue #1097)."""
-    tensor = PolynomialTensor(
-        {
-            (): 1.0,
-            (1, 0): numpy.array([[1.0, 2.0], [3.0, 4.0]]),
-            (1, 1, 0, 0): numpy.arange(16, dtype=float).reshape((2, 2, 2, 2)),
-        }
-    )
-    cases = [(numpy.int64(2), 2), (numpy.int32(3), 3), (numpy.float32(0.5), 0.5)]
-    for numpy_scalar, python_scalar in cases:
-        assert tensor * numpy_scalar == tensor * python_scalar
-        assert tensor + numpy_scalar == tensor + python_scalar
-        assert tensor - numpy_scalar == tensor - python_scalar
-        assert tensor / numpy_scalar == tensor / python_scalar
+    def test_numpy_scalar_coefficients(self):
+        """NumPy scalar coefficients behave like Python scalars (issue #1097)."""
+        tensor = PolynomialTensor(
+            {
+                (): 1.0,
+                (1, 0): numpy.array([[1.0, 2.0], [3.0, 4.0]]),
+                (1, 1, 0, 0): numpy.arange(16, dtype=float).reshape((2, 2, 2, 2)),
+            }
+        )
+        cases = [(numpy.int64(2), 2), (numpy.int32(3), 3), (numpy.float32(0.5), 0.5)]
+        for numpy_scalar, python_scalar in cases:
+            self.assertEqual(tensor * numpy_scalar, tensor * python_scalar)
+            self.assertEqual(tensor + numpy_scalar, tensor + python_scalar)
+            self.assertEqual(tensor - numpy_scalar, tensor - python_scalar)
+            self.assertEqual(tensor / numpy_scalar, tensor / python_scalar)

From cad8b5eee4cc93760a18cadafcf454031754d34b Mon Sep 17 00:00:00 2001
From: Aaryan Mehta <73230976+blazingphoenix7@users.noreply.github.com>
Date: Tue, 16 Jun 2026 13:18:46 -0400
Subject: [PATCH 3/4] Centralize COEFFICIENT_TYPES in config and precompute
 Majorana mul types

Addresses the review comments on this PR:

- COEFFICIENT_TYPES was defined identically in several files. Moved the
  definition to openfermion/config.py and have the other modules pick it up
  via `import openfermion.config as config` plus a local
  COEFFICIENT_TYPES = config.COEFFICIENT_TYPES binding, so there is a single
  source of truth and no repeated attribute lookups. While editing
  polynomial_tensor.py I routed its EQ_TOLERANCE through the same config import
  so the file does not mix import styles; special_operators.py uses
  config.COEFFICIENT_TYPES directly for its one check.

- MajoranaOperator.__mul__ and __imul__ built (type(self), COEFFICIENT_TYPES)
  on every call. Added a module-level _MAJORANA_MUL_TYPES tuple, built once
  after the class, and use it in both.
---
 src/openfermion/config.py                         |  6 ++++++
 src/openfermion/hamiltonians/special_operators.py |  4 ++--
 .../ops/operators/majorana_operator.py            | 15 +++++++++++----
 .../ops/representations/doci_hamiltonian.py       |  5 ++---
 .../ops/representations/polynomial_tensor.py      |  6 +++---
 5 files changed, 24 insertions(+), 12 deletions(-)

diff --git a/src/openfermion/config.py b/src/openfermion/config.py
index 9b6cf44f2..ec9600cf9 100644
--- a/src/openfermion/config.py
+++ b/src/openfermion/config.py
@@ -10,11 +10,17 @@
 #   See the License for the specific language governing permissions and
 #   limitations under the License.
 
+import numbers
 import os
 
 # Tolerance to consider number zero.
 EQ_TOLERANCE = 1e-8
 
+# Numeric types accepted as operator and tensor coefficients. numbers.Number
+# covers Python and NumPy scalar types alike (NumPy scalars are not subclasses
+# of the built-in int/float/complex types).
+COEFFICIENT_TYPES = (int, float, complex, numbers.Number)
+
 # Molecular data directory.
 THIS_DIRECTORY = os.path.dirname(os.path.realpath(__file__))
 DATA_DIRECTORY = os.path.realpath(os.path.join(THIS_DIRECTORY, 'testing/data'))
diff --git a/src/openfermion/hamiltonians/special_operators.py b/src/openfermion/hamiltonians/special_operators.py
index 5893556ef..f4b4c6d99 100644
--- a/src/openfermion/hamiltonians/special_operators.py
+++ b/src/openfermion/hamiltonians/special_operators.py
@@ -11,9 +11,9 @@
 #   limitations under the License.
 """Commonly used operators (mainly instances of SymbolicOperator)."""
 
-import numbers
 from typing import Optional, Union, Tuple
 
+import openfermion.config as config
 from openfermion.ops.operators import BosonOperator, FermionOperator
 from openfermion.utils.indexing import down_index, up_index
 
@@ -237,7 +237,7 @@ def majorana_operator(
     Returns:
         FermionOperator
     """
-    if not isinstance(coefficient, (int, float, complex, numbers.Number)):
+    if not isinstance(coefficient, config.COEFFICIENT_TYPES):
         raise ValueError('Coefficient must be scalar.')
 
     # If term is a string, convert it to a tuple
diff --git a/src/openfermion/ops/operators/majorana_operator.py b/src/openfermion/ops/operators/majorana_operator.py
index 7e6724013..a8613d7d0 100644
--- a/src/openfermion/ops/operators/majorana_operator.py
+++ b/src/openfermion/ops/operators/majorana_operator.py
@@ -13,11 +13,12 @@
 
 import copy
 import itertools
-import numbers
 
 import numpy
 
-COEFFICIENT_TYPES = (int, float, complex, numbers.Number)
+import openfermion.config as config
+
+COEFFICIENT_TYPES = config.COEFFICIENT_TYPES
 
 
 class MajoranaOperator:
@@ -179,7 +180,7 @@ def __sub__(self, other):
         return minuend
 
     def __mul__(self, other):
-        if not isinstance(other, (type(self), COEFFICIENT_TYPES)):
+        if not isinstance(other, _MAJORANA_MUL_TYPES):
             return NotImplemented
 
         if isinstance(other, COEFFICIENT_TYPES):
@@ -198,7 +199,7 @@ def __mul__(self, other):
         return MajoranaOperator.from_dict(terms)
 
     def __imul__(self, other):
-        if not isinstance(other, (type(self), COEFFICIENT_TYPES)):
+        if not isinstance(other, _MAJORANA_MUL_TYPES):
             return NotImplemented
 
         if isinstance(other, COEFFICIENT_TYPES):
@@ -279,6 +280,12 @@ def __repr__(self):
         return 'MajoranaOperator.from_dict(terms={!r})'.format(self.terms)
 
 
+# Types accepted by MajoranaOperator multiplication: another MajoranaOperator or
+# a scalar coefficient. Defined here, after the class, so the tuple is built once
+# at import rather than on every __mul__/__imul__ call.
+_MAJORANA_MUL_TYPES = (MajoranaOperator,) + COEFFICIENT_TYPES
+
+
 def _sort_majorana_term(term):
     """Sort a Majorana term.
 
diff --git a/src/openfermion/ops/representations/doci_hamiltonian.py b/src/openfermion/ops/representations/doci_hamiltonian.py
index 6fabf1588..151c8ad4d 100644
--- a/src/openfermion/ops/representations/doci_hamiltonian.py
+++ b/src/openfermion/ops/representations/doci_hamiltonian.py
@@ -11,14 +11,13 @@
 #   limitations under the License.
 """Class and functions to store interaction operators."""
 
-import numbers
-
 import numpy
 
+import openfermion.config as config
 from openfermion.ops import QubitOperator
 from openfermion.ops.representations import PolynomialTensor, get_tensors_from_integrals
 
-COEFFICIENT_TYPES = (int, float, complex, numbers.Number)
+COEFFICIENT_TYPES = config.COEFFICIENT_TYPES
 
 
 class DOCIHamiltonian(PolynomialTensor):
diff --git a/src/openfermion/ops/representations/polynomial_tensor.py b/src/openfermion/ops/representations/polynomial_tensor.py
index ded9f26b6..df703abfe 100644
--- a/src/openfermion/ops/representations/polynomial_tensor.py
+++ b/src/openfermion/ops/representations/polynomial_tensor.py
@@ -14,14 +14,14 @@
 
 import copy
 import itertools
-import numbers
 import operator
 
 import numpy
 
-from openfermion.config import EQ_TOLERANCE
+import openfermion.config as config
 
-COEFFICIENT_TYPES = (int, float, complex, numbers.Number)
+EQ_TOLERANCE = config.EQ_TOLERANCE
+COEFFICIENT_TYPES = config.COEFFICIENT_TYPES
 
 
 class PolynomialTensorError(Exception):

From a33759a54000702b1d0e9529cc34239fd14436e3 Mon Sep 17 00:00:00 2001
From: Aaryan Mehta <73230976+blazingphoenix7@users.noreply.github.com>
Date: Tue, 16 Jun 2026 17:41:15 -0400
Subject: [PATCH 4/4] Bind COEFFICIENT_TYPES at module level in
 special_operators too

special_operators.py still referenced config.COEFFICIENT_TYPES inline. Pulled
it up to a module-level COEFFICIENT_TYPES = config.COEFFICIENT_TYPES binding to
match majorana_operator.py and the representation modules.
---
 src/openfermion/hamiltonians/special_operators.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/openfermion/hamiltonians/special_operators.py b/src/openfermion/hamiltonians/special_operators.py
index f4b4c6d99..92f2f0cb2 100644
--- a/src/openfermion/hamiltonians/special_operators.py
+++ b/src/openfermion/hamiltonians/special_operators.py
@@ -17,6 +17,8 @@
 from openfermion.ops.operators import BosonOperator, FermionOperator
 from openfermion.utils.indexing import down_index, up_index
 
+COEFFICIENT_TYPES = config.COEFFICIENT_TYPES
+
 
 def s_plus_operator(n_spatial_orbitals: int) -> FermionOperator:
     r"""Return the s+ operator.
@@ -237,7 +239,7 @@ def majorana_operator(
     Returns:
         FermionOperator
     """
-    if not isinstance(coefficient, config.COEFFICIENT_TYPES):
+    if not isinstance(coefficient, COEFFICIENT_TYPES):
         raise ValueError('Coefficient must be scalar.')
 
     # If term is a string, convert it to a tuple