diff --git a/pyproject.toml b/pyproject.toml
index f6c623ee5..9054fc773 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -88,7 +88,6 @@ disable = [
     "invalid-unary-operand-type",
     "missing-function-docstring",
     "missing-module-docstring",
-    "no-else-return",
     "not-callable",  # pylint doesn't understand metaclass properties
     "protected-access",
     "too-many-ancestors",
diff --git a/src/galois/_codes/_linear.py b/src/galois/_codes/_linear.py
index 6c7252565..a96f749ea 100644
--- a/src/galois/_codes/_linear.py
+++ b/src/galois/_codes/_linear.py
@@ -91,9 +91,9 @@ def encode(self, message: ArrayLike, output: Literal["codeword", "parity"] = "co
 
         if output == "codeword":
             return codeword
-        else:
-            parity = self._convert_codeword_to_parity(codeword)
-            return parity
+
+        parity = self._convert_codeword_to_parity(codeword)
+        return parity
 
     def detect(self, codeword: ArrayLike) -> bool | np.ndarray:
         r"""
@@ -187,11 +187,11 @@ def decode(self, codeword, output="message", errors=False):
         if is_codeword_1d:
             decoded, N_errors = decoded[0, :], int(N_errors[0])
 
-        if not errors:
-            return decoded
-        else:
+        if errors:
             return decoded, N_errors
 
+        return decoded
+
     # def dual_code(self) -> _LinearCode:
     #     n = self.n
     #     k = self.n - self.k
diff --git a/src/galois/_domains/_calculate.py b/src/galois/_domains/_calculate.py
index 8d4c11cf7..1c7892145 100644
--- a/src/galois/_domains/_calculate.py
+++ b/src/galois/_domains/_calculate.py
@@ -162,13 +162,13 @@ def __call__(self, ufunc, method, inputs, kwargs, meta):
         if self.field.ufunc_mode == "jit-lookup" or method != "__call__":
             # Use the lookup ufunc on each array entry
             return super().__call__(ufunc, method, inputs, kwargs, meta)
-        else:
-            # Convert entire array to polynomial/vector representation, perform array operation in GF(p), and convert back to GF(p^m)
-            self._verify_operands_in_same_field(ufunc, inputs, meta)
-            inputs, kwargs = self._convert_inputs_to_vector(inputs, kwargs)
-            output = getattr(ufunc, method)(*inputs, **kwargs)
-            output = self._convert_output_from_vector(output, meta["dtype"])
-            return output
+
+        # Convert entire array to polynomial/vector representation, perform array operation in GF(p), and convert back to GF(p^m)
+        self._verify_operands_in_same_field(ufunc, inputs, meta)
+        inputs, kwargs = self._convert_inputs_to_vector(inputs, kwargs)
+        output = getattr(ufunc, method)(*inputs, **kwargs)
+        output = self._convert_output_from_vector(output, meta["dtype"])
+        return output
 
     def set_calculate_globals(self):
         global CHARACTERISTIC, DEGREE
@@ -213,13 +213,13 @@ def __call__(self, ufunc, method, inputs, kwargs, meta):
         if self.field.ufunc_mode == "jit-lookup" or method != "__call__":
             # Use the lookup ufunc on each array entry
             return super().__call__(ufunc, method, inputs, kwargs, meta)
-        else:
-            # Convert entire array to polynomial/vector representation, perform array operation in GF(p), and convert back to GF(p^m)
-            self._verify_operands_in_same_field(ufunc, inputs, meta)
-            inputs, kwargs = self._convert_inputs_to_vector(inputs, kwargs)
-            output = getattr(ufunc, method)(*inputs, **kwargs)
-            output = self._convert_output_from_vector(output, meta["dtype"])
-            return output
+
+        # Convert entire array to polynomial/vector representation, perform array operation in GF(p), and convert back to GF(p^m)
+        self._verify_operands_in_same_field(ufunc, inputs, meta)
+        inputs, kwargs = self._convert_inputs_to_vector(inputs, kwargs)
+        output = getattr(ufunc, method)(*inputs, **kwargs)
+        output = self._convert_output_from_vector(output, meta["dtype"])
+        return output
 
     def set_calculate_globals(self):
         global CHARACTERISTIC, DEGREE
@@ -264,13 +264,13 @@ def __call__(self, ufunc, method, inputs, kwargs, meta):
         if self.field.ufunc_mode == "jit-lookup" or method != "__call__":
             # Use the lookup ufunc on each array entry
             return super().__call__(ufunc, method, inputs, kwargs, meta)
-        else:
-            # Convert entire array to polynomial/vector representation, perform array operation in GF(p), and convert back to GF(p^m)
-            self._verify_operands_in_same_field(ufunc, inputs, meta)
-            inputs, kwargs = self._convert_inputs_to_vector(inputs, kwargs)
-            output = getattr(ufunc, method)(*inputs, **kwargs)
-            output = self._convert_output_from_vector(output, meta["dtype"])
-            return output
+
+        # Convert entire array to polynomial/vector representation, perform array operation in GF(p), and convert back to GF(p^m)
+        self._verify_operands_in_same_field(ufunc, inputs, meta)
+        inputs, kwargs = self._convert_inputs_to_vector(inputs, kwargs)
+        output = getattr(ufunc, method)(*inputs, **kwargs)
+        output = self._convert_output_from_vector(output, meta["dtype"])
+        return output
 
     def set_calculate_globals(self):
         global CHARACTERISTIC, DEGREE
@@ -658,28 +658,25 @@ def compute_x(x):
             # Equation 3.2
             if x % 3 == 1:
                 return MULTIPLY(beta, x)
-            elif x % 3 == 2:
+            if x % 3 == 2:
                 return MULTIPLY(x, x)
-            else:
-                return MULTIPLY(alpha, x)
+            return MULTIPLY(alpha, x)
 
         def compute_a(a, x):
             # Equation 3.3
             if x % 3 == 1:
                 return a
-            elif x % 3 == 2:
+            if x % 3 == 2:
                 return (2 * a) % n
-            else:
-                return (a + 1) % n
+            return (a + 1) % n
 
         def compute_b(b, x):
             # Equation 3.4
             if x % 3 == 1:
                 return (b + 1) % n
-            elif x % 3 == 2:
+            if x % 3 == 2:
                 return (2 * b) % n
-            else:
-                return b
+            return b
 
         while True:
             xi, ai, bi = compute_x(xi), compute_a(ai, xi), compute_b(bi, xi)
@@ -693,14 +690,14 @@ def compute_b(b, x):
                     d, r_inv = EGCD(r, n)[0:2]
                     assert d == 1
                     return (r_inv * (a2i - ai)) % n
-                else:
-                    # Re-try with different x0, a0, and b0
-                    a0 += 1
-                    b0 += 1
-                    x0 = MULTIPLY(x0, beta)
-                    x0 = MULTIPLY(x0, alpha)
-                    xi, ai, bi = x0, a0, b0
-                    x2i, a2i, b2i = xi, ai, bi
+
+                # Re-try with different x0, a0, and b0
+                a0 += 1
+                b0 += 1
+                x0 = MULTIPLY(x0, beta)
+                x0 = MULTIPLY(x0, alpha)
+                xi, ai, bi = x0, a0, b0
+                x2i, a2i, b2i = xi, ai, bi
 
 
 class log_pohlig_hellman(_lookup.log_ufunc):
diff --git a/src/galois/_domains/_function.py b/src/galois/_domains/_function.py
index caf2d2d50..61a93b27f 100644
--- a/src/galois/_domains/_function.py
+++ b/src/galois/_domains/_function.py
@@ -60,10 +60,9 @@ def function(self):
         """
         Returns a JIT-compiled or pure-Python function based on field size.
         """
-        if self.field.ufunc_mode != "python-calculate":
-            return self.jit
-        else:
+        if self.field.ufunc_mode == "python-calculate":
             return self.python
+        return self.jit
 
     @property
     def jit(self) -> numba.types.FunctionType:
diff --git a/src/galois/_domains/_linalg.py b/src/galois/_domains/_linalg.py
index e3abf3f40..5d21702bc 100644
--- a/src/galois/_domains/_linalg.py
+++ b/src/galois/_domains/_linalg.py
@@ -80,13 +80,13 @@ def __call__(self, a: Array, b: Array, out=None) -> Array:
             return _lapack_linalg(self.field, a, b, np.dot, out=out)
 
         if a.ndim == 0 or b.ndim == 0:
-            return a * b
+            dot = a * b
         elif a.ndim == 1 and b.ndim == 1:
-            return np.sum(a * b)
+            dot = np.sum(a * b)
         elif a.ndim == 2 and b.ndim == 2:
-            return np.matmul(a, b, out=out)
+            dot = np.matmul(a, b, out=out)
         elif a.ndim >= 2 and b.ndim == 1:
-            return np.sum(a * b, axis=-1, out=out)
+            dot = np.sum(a * b, axis=-1, out=out)
         # elif a.dnim >= 2 and b.ndim >= 2:
         else:
             raise NotImplementedError(
@@ -94,6 +94,8 @@ def __call__(self, a: Array, b: Array, out=None) -> Array:
                 "Please open a GitHub issue at https://github.com/mhostetter/galois/issues."
             )
 
+        return dot
+
 
 class vdot_jit(Function):
     """
@@ -150,8 +152,8 @@ def __call__(self, a: Array, b: Array, out=None) -> Array:
 
         if self.field._is_prime_field:
             return _lapack_linalg(self.field, a, b, np.outer, out=out, n_sum=1)
-        else:
-            return np.multiply.outer(a.ravel(), b.ravel(), out=out)
+
+        return np.multiply.outer(a.ravel(), b.ravel(), out=out)
 
 
 class matmul_jit(Function):
@@ -395,9 +397,9 @@ def __call__(self, A: Array) -> Array:
         n = A.shape[0]
 
         if n == 2:
-            return A[0, 0] * A[1, 1] - A[0, 1] * A[1, 0]
+            det = A[0, 0] * A[1, 1] - A[0, 1] * A[1, 0]
         elif n == 3:
-            return (
+            det = (
                 A[0, 0] * (A[1, 1] * A[2, 2] - A[1, 2] * A[2, 1])
                 - A[0, 1] * (A[1, 0] * A[2, 2] - A[1, 2] * A[2, 0])
                 + A[0, 2] * (A[1, 0] * A[2, 1] - A[1, 1] * A[2, 0])
@@ -408,7 +410,9 @@ def __call__(self, A: Array) -> Array:
             det_P = (-self.field(1)) ** N_permutations
             det_L = triangular_det_jit(self.field)(L)
             det_U = triangular_det_jit(self.field)(U)
-            return det_P * det_L * det_U
+            det = det_P * det_L * det_U
+
+        return det
 
 
 ###############################################################################
diff --git a/src/galois/_domains/_lookup.py b/src/galois/_domains/_lookup.py
index 6cfc2d06e..20fa9b25b 100644
--- a/src/galois/_domains/_lookup.py
+++ b/src/galois/_domains/_lookup.py
@@ -42,7 +42,7 @@ def lookup(a: int, b: int) -> int:  # pragma: no cover
         """
         if a == 0:
             return b
-        elif b == 0:
+        if b == 0:
             return a
 
         m = LOG[a]
@@ -56,8 +56,8 @@ def lookup(a: int, b: int) -> int:  # pragma: no cover
         if n - m == ZECH_E:
             # zech_log(zech_e) = -Inf and α^(-Inf) = 0
             return 0
-        else:
-            return EXP[m + ZECH_LOG[n - m]]
+
+        return EXP[m + ZECH_LOG[n - m]]
 
 
 class negative_ufunc(_ufunc.negative_ufunc):
@@ -85,9 +85,9 @@ def lookup(a: int) -> int:  # pragma: no cover
         """
         if a == 0:
             return 0
-        else:
-            m = LOG[a]
-            return EXP[ZECH_E + m]
+
+        m = LOG[a]
+        return EXP[ZECH_E + m]
 
 
 class subtract_ufunc(_ufunc.subtract_ufunc):
@@ -123,7 +123,7 @@ def lookup(a: int, b: int) -> int:  # pragma: no cover
 
         if b == 0:
             return a
-        elif a == 0:
+        if a == 0:
             return EXP[n]
 
         if m > n:
@@ -165,10 +165,10 @@ def lookup(a: int, b: int) -> int:  # pragma: no cover
         """
         if a == 0 or b == 0:
             return 0
-        else:
-            m = LOG[a]
-            n = LOG[b]
-            return EXP[m + n]
+
+        m = LOG[a]
+        n = LOG[b]
+        return EXP[m + n]
 
 
 class reciprocal_ufunc(_ufunc.reciprocal_ufunc):
@@ -232,10 +232,10 @@ def lookup(a: int, b: int) -> int:  # pragma: no cover
 
         if a == 0:
             return 0
-        else:
-            m = LOG[a]
-            n = LOG[b]
-            return EXP[(ORDER - 1) + m - n]  # We add `ORDER - 1` to guarantee the index is non-negative
+
+        m = LOG[a]
+        n = LOG[b]
+        return EXP[(ORDER - 1) + m - n]  # We add `ORDER - 1` to guarantee the index is non-negative
 
 
 class power_ufunc(_ufunc.power_ufunc):
@@ -269,11 +269,11 @@ def lookup(a: int, b: int) -> int:  # pragma: no cover
 
         if b == 0:
             return 1
-        elif a == 0:
+        if a == 0:
             return 0
-        else:
-            m = LOG[a]
-            return EXP[(m * b) % (ORDER - 1)]  # TODO: Do b % (ORDER - 1) first? b could be very large and overflow int64
+
+        m = LOG[a]
+        return EXP[(m * b) % (ORDER - 1)]  # TODO: Do b % (ORDER - 1) first? b could be very large and overflow int64
 
 
 class log_ufunc(_ufunc.log_ufunc):
diff --git a/src/galois/_domains/_meta.py b/src/galois/_domains/_meta.py
index 7f0789ec3..fc2fa7f21 100644
--- a/src/galois/_domains/_meta.py
+++ b/src/galois/_domains/_meta.py
@@ -70,13 +70,13 @@ def __init__(cls, name, bases, namespace, **kwargs):
     def __repr__(cls) -> str:
         if cls.order == 0:
             # This is not a runtime-created subclass, so return the base class name.
-            return f"<class '{cls.__module__}.{cls.__name__}'>"
-        else:
-            # When FieldArray instances are created they are added to the `galois._fields._factory` module with a name
-            # like `FieldArray_<p>_<primitive_element>` or `FieldArray_<p>_<m>_<primitive_element>_<irreducible_poly>`.
-            # This is visually unappealing. So here we override the repr() to be more succinct and indicate how the class
-            # was created. So galois._fields._factory.FieldArray_31_3 is converted to galois.GF(31).
-            return f"<class 'galois.{cls.name}'>"
+            return f"<class 'galois.{cls.__name__}'>"
+
+        # When FieldArray instances are created they are added to the `galois._fields._factory` module with a name
+        # like `FieldArray_<p>_<primitive_element>` or `FieldArray_<p>_<m>_<primitive_element>_<irreducible_poly>`.
+        # This is visually unappealing. So here we override the repr() to be more succinct and indicate how the class
+        # was created. So galois._fields._factory.FieldArray_31_3 is converted to galois.GF(31).
+        return f"<class 'galois.{cls.name}'>"
 
     def __dir__(cls) -> list[str]:
         """
diff --git a/src/galois/_domains/_ufunc.py b/src/galois/_domains/_ufunc.py
index bd167e9d8..8dd7af8a5 100644
--- a/src/galois/_domains/_ufunc.py
+++ b/src/galois/_domains/_ufunc.py
@@ -70,12 +70,11 @@ def ufunc(self):
         """
         A ufunc based on the current state of `ufunc_mode`.
         """
+        if self.field.ufunc_mode == "python-calculate":
+            return self.python_calculate
         if self.field.ufunc_mode == "jit-lookup" and not self.always_calculate:
             return self.jit_lookup
-        elif self.field.ufunc_mode == "python-calculate":
-            return self.python_calculate
-        else:
-            return self.jit_calculate
+        return self.jit_calculate
 
     @property
     def ufunc_call_only(self):
@@ -86,14 +85,13 @@ def ufunc_call_only(self):
         (see https://github.com/mhostetter/galois/issues/358). This ufunc should be used in JIT function
         implementations (`Function.implementation`) to ensure bug #358 does not manifest.
         """
-        if self.field.ufunc_mode == "jit-lookup" and not self.always_calculate:
-            return self.jit_lookup
-        elif self.field.ufunc_mode == "python-calculate":
+        if self.field.ufunc_mode == "python-calculate":
             # Specify `dtype=np.object_` for overridden ufuncs so Python int objects are returned, not np.intc (which
             # will eventually overflow and produce incorrect results).
             return self.python_calculate_call_only
-        else:
-            return self.jit_calculate
+        if self.field.ufunc_mode == "jit-lookup" and not self.always_calculate:
+            return self.jit_lookup
+        return self.jit_calculate
 
     @property
     def jit_calculate(self) -> numba.types.FunctionType:
@@ -155,8 +153,8 @@ def python_calculate(self) -> Callable:
 
         if self.type == "unary":
             return np.frompyfunc(self.calculate, 1, 1)
-        else:
-            return np.frompyfunc(self.calculate, 2, 1)
+
+        return np.frompyfunc(self.calculate, 2, 1)
 
     @property
     def python_calculate_call_only(self) -> Callable:
@@ -170,8 +168,8 @@ def python_calculate_call_only(self) -> Callable:
             # Specify `dtype=np.object_` for overridden ufuncs so Python int objects are returned, not np.intc (which
             # will eventually overflow and produce incorrect results).
             return lambda *args, **kwargs: self.override(*args, **kwargs, dtype=np.object_)
-        else:
-            return self.python_calculate
+
+        return self.python_calculate
 
     ###############################################################################
     # Input/output verification
@@ -242,13 +240,14 @@ def _verify_operands_first_field_second_int(self, ufunc, inputs, meta):
         second = inputs[meta["operands"][1]]
         if isinstance(second, (int, np.integer)):
             return
-        # elif type(second) is np.ndarray:
+
+        # if type(second) is np.ndarray:
         #     if not np.issubdtype(second.dtype, np.integer):
         #         raise ValueError(
         #             f"Operation {ufunc.__name__!r} requires the second operand with type np.ndarray to have integer dtype, "
         #             f"not {second.dtype}."
         #         )
-        elif isinstance(second, np.ndarray):
+        if isinstance(second, np.ndarray):
             if self.field.dtypes == [np.object_]:
                 if not (second.dtype == np.object_ or np.issubdtype(second.dtype, np.integer)):
                     raise ValueError(
@@ -291,8 +290,8 @@ def _convert_inputs_to_vector(self, inputs, kwargs):
     def _convert_output_from_vector(self, output, dtype):
         if output is None:
             return None
-        else:
-            return self.field.Vector(output, dtype=dtype)
+
+        return self.field.Vector(output, dtype=dtype)
 
     ###############################################################################
     # Input/output type viewing
@@ -322,12 +321,11 @@ def _view_inputs_as_ndarray(self, inputs, kwargs, dtype=None):
     def _view_output_as_field(self, output, field, dtype):
         if isinstance(type(output), field):
             return output
-        elif isinstance(output, np.ndarray):
+        if isinstance(output, np.ndarray):
             return field._view(output.astype(dtype))
-        elif output is None:
+        if output is None:
             return None
-        else:
-            return field(output, dtype=dtype)
+        return field(output, dtype=dtype)
 
 
 ###############################################################################
@@ -681,24 +679,24 @@ def __array_ufunc__(self, ufunc, method, *inputs, **kwargs):
 
             return getattr(field, field._OVERRIDDEN_UFUNCS[ufunc])(ufunc, method, inputs, kwargs, meta)
 
-        elif ufunc in field._UNSUPPORTED_UFUNCS:
+        if ufunc in field._UNSUPPORTED_UFUNCS:
             raise NotImplementedError(
                 f"The NumPy ufunc {ufunc.__name__!r} is not supported on {field.name} arrays. "
                 "If you believe this ufunc should be supported, "
                 "please submit a GitHub issue at https://github.com/mhostetter/galois/issues."
             )
 
-        else:
-            if ufunc in [np.bitwise_and, np.bitwise_or, np.bitwise_xor] and method not in ["reduce", "accumulate", "at", "reduceat"]:
-                kwargs["casting"] = "unsafe"
+        # Process with our custom ufuncs
+        if ufunc in [np.bitwise_and, np.bitwise_or, np.bitwise_xor] and method not in ["reduce", "accumulate", "at", "reduceat"]:
+            kwargs["casting"] = "unsafe"
 
-            inputs, kwargs = UFunc(field)._view_inputs_as_ndarray(inputs, kwargs)
-            output = super().__array_ufunc__(ufunc, method, *inputs, **kwargs)  # pylint: disable=no-member
+        inputs, kwargs = UFunc(field)._view_inputs_as_ndarray(inputs, kwargs)
+        output = super().__array_ufunc__(ufunc, method, *inputs, **kwargs)  # pylint: disable=no-member
 
-            if ufunc in field._UFUNCS_REQUIRING_VIEW and output is not None:
-                output = field._view(output) if not np.isscalar(output) else field(output, dtype=self.dtype)
+        if ufunc in field._UFUNCS_REQUIRING_VIEW and output is not None:
+            output = field._view(output) if not np.isscalar(output) else field(output, dtype=self.dtype)
 
-            return output
+        return output
 
     def __pow__(self, other):
         # We call power here instead of `super().__pow__(other)` because when doing so `x ** GF(2)` will invoke `np.square(x)`
diff --git a/src/galois/_fields/_array.py b/src/galois/_fields/_array.py
index 241099149..700d4a263 100644
--- a/src/galois/_fields/_array.py
+++ b/src/galois/_fields/_array.py
@@ -159,10 +159,9 @@ def _verify_array_like_types_and_values(cls, x: ElementLike | ArrayLike) -> Elem
     def _verify_element_types_and_convert(cls, array: np.ndarray, object_=False) -> np.ndarray:
         if array.size == 0:
             return array
-        elif object_:
+        if object_:
             return np.vectorize(cls._convert_to_element, otypes=[object])(array)
-        else:
-            return np.vectorize(cls._convert_to_element)(array)
+        return np.vectorize(cls._convert_to_element)(array)
 
     @classmethod
     def _verify_scalar_value(cls, scalar: int):
@@ -1355,14 +1354,16 @@ def field_trace(self) -> FieldArray:
         x = self
 
         if field.is_prime_field:
-            return x.copy()
+            trace = x.copy()
         else:
             subfield = field.prime_subfield
             p = field.characteristic
             m = field.degree
             conjugates = np.power.outer(x, p ** np.arange(0, m, dtype=field.dtypes[-1]))
             trace = np.add.reduce(conjugates, axis=-1)
-            return subfield._view(trace)
+            trace = subfield._view(trace)
+
+        return trace
 
     def field_norm(self) -> FieldArray:
         r"""
@@ -1402,13 +1403,15 @@ def field_norm(self) -> FieldArray:
         x = self
 
         if field.is_prime_field:
-            return x.copy()
+            norm = x.copy()
         else:
             subfield = field.prime_subfield
             p = field.characteristic
             m = field.degree
             norm = x ** ((p**m - 1) // (p - 1))
-            return subfield._view(norm)
+            norm = subfield._view(norm)
+
+        return norm
 
     def characteristic_poly(self) -> Poly:
         r"""
@@ -1470,13 +1473,12 @@ def characteristic_poly(self) -> Poly:
         """
         if self.ndim == 0:
             return _characteristic_poly_element(self)
-        elif self.ndim == 2:
+        if self.ndim == 2:
             return _characteristic_poly_matrix(self)
-        else:
-            raise ValueError(
-                f"The array must be either 0-D to return the characteristic polynomial of a single element "
-                f"or 2-D to return the characteristic polynomial of a square matrix, not have shape {self.shape}."
-            )
+        raise ValueError(
+            f"The array must be either 0-D to return the characteristic polynomial of a single element "
+            f"or 2-D to return the characteristic polynomial of a square matrix, not have shape {self.shape}."
+        )
 
     def minimal_poly(self) -> Poly:
         r"""
@@ -1519,13 +1521,12 @@ def minimal_poly(self) -> Poly:
         """
         if self.ndim == 0:
             return _minimal_poly_element(self)
-        # elif self.ndim == 2:
+        # if self.ndim == 2:
         #     return _minimal_poly_matrix(self)
-        else:
-            raise ValueError(
-                f"The array must be either 0-D to return the minimal polynomial of a single element "
-                f"or 2-D to return the minimal polynomial of a square matrix, not have shape {self.shape}."
-            )
+        raise ValueError(
+            f"The array must be either 0-D to return the minimal polynomial of a single element "
+            f"or 2-D to return the minimal polynomial of a square matrix, not have shape {self.shape}."
+        )
 
     def log(self, base: ElementLike | ArrayLike | None = None) -> int | np.ndarray:
         r"""
@@ -1679,11 +1680,13 @@ def _display(self, mode: Literal["repr", "str"], separator=", ") -> str:
             last_line_width = len(string[idx:] + ", " + order + suffix)
 
             if last_line_width <= np.get_printoptions()["linewidth"]:
-                return prefix + string + ", " + order + suffix
+                output = prefix + string + ", " + order + suffix
             else:
-                return prefix + string + ",\n" + " " * len(prefix) + order + suffix
+                output = prefix + string + ",\n" + " " * len(prefix) + order + suffix
         else:
-            return prefix + string + suffix
+            output = prefix + string + suffix
+
+        return output
 
     @classmethod
     def _formatter(cls, array):
@@ -1784,12 +1787,13 @@ def _characteristic_poly_element(a: FieldArray) -> Poly:
     x = Poly.Identity(field)
 
     if field.is_prime_field:
-        return x - a
+        poly = x - a
     else:
         powers = a ** (field.characteristic ** np.arange(0, field.degree, dtype=field.dtypes[-1]))
         poly = Poly.Roots(powers, field=field)
         poly = Poly(poly.coeffs, field=field.prime_subfield)
-        return poly
+
+    return poly
 
 
 def _characteristic_poly_matrix(A: FieldArray) -> Poly:
@@ -1821,9 +1825,10 @@ def _minimal_poly_element(a: FieldArray) -> Poly:
     x = Poly.Identity(field)
 
     if field.is_prime_field:
-        return x - a
+        poly = x - a
     else:
         conjugates = np.unique(a ** (field.characteristic ** np.arange(0, field.degree, dtype=field.dtypes[-1])))
         poly = Poly.Roots(conjugates, field=field)
         poly = Poly(poly.coeffs, field=field.prime_subfield)
-        return poly
+
+    return poly
diff --git a/src/galois/_fields/_factory.py b/src/galois/_fields/_factory.py
index fd2348e97..b60e7d90e 100644
--- a/src/galois/_fields/_factory.py
+++ b/src/galois/_fields/_factory.py
@@ -274,7 +274,7 @@ def GF(
             raise ValueError(
                 f"Argument 'irreducible_poly' can only be specified for extension fields, not the prime field GF({characteristic})."
             )
-        return _GF_prime(
+        field = _GF_prime(
             characteristic,
             alpha=primitive_element,
             verify=verify,
@@ -282,7 +282,7 @@ def GF(
             repr=repr,
         )
     else:
-        return _GF_extension(
+        field = _GF_extension(
             characteristic,
             degree,
             irreducible_poly_=irreducible_poly,
@@ -292,6 +292,8 @@ def GF(
             repr=repr,
         )
 
+    return field
+
 
 @overload
 def Field(
diff --git a/src/galois/_lfsr.py b/src/galois/_lfsr.py
index ed7dcff12..bacbd63d4 100644
--- a/src/galois/_lfsr.py
+++ b/src/galois/_lfsr.py
@@ -98,11 +98,13 @@ def step(self, steps: int = 1) -> FieldArray:
         verify_isinstance(steps, int)
 
         if steps == 0:
-            return self.field([])
+            y = self.field([])
         elif steps > 0:
-            return self._step_forward(steps)
+            y = self._step_forward(steps)
         else:
-            return self._step_backward(abs(steps))
+            y = self._step_backward(abs(steps))
+
+        return y
 
     def _step_forward(self, steps):
         assert steps > 0
@@ -1601,16 +1603,16 @@ def berlekamp_massey(sequence, output="minimal"):
 
     if output == "minimal":
         return characteristic_poly
-    else:
-        # The first n outputs are the Fibonacci state reversed
-        feedback_poly = characteristic_poly.reverse()
-        state_ = sequence[0 : feedback_poly.degree][::-1]
-        fibonacci_lfsr = FLFSR(feedback_poly, state=state_)
-
-        if output == "fibonacci":
-            return fibonacci_lfsr
-        else:
-            return fibonacci_lfsr.to_galois_lfsr()
+
+    # The first n outputs are the Fibonacci state reversed
+    feedback_poly = characteristic_poly.reverse()
+    state_ = sequence[0 : feedback_poly.degree][::-1]
+    fibonacci_lfsr = FLFSR(feedback_poly, state=state_)
+
+    if output == "fibonacci":
+        return fibonacci_lfsr
+
+    return fibonacci_lfsr.to_galois_lfsr()
 
 
 class berlekamp_massey_jit(Function):
diff --git a/src/galois/_math.py b/src/galois/_math.py
index 7a0b7867e..faef294a0 100644
--- a/src/galois/_math.py
+++ b/src/galois/_math.py
@@ -119,10 +119,10 @@ def isqrt(n: int) -> int:
     # Recursively compute the integer square root
     x = isqrt(n >> 2) << 1
 
-    if (x + 1) ** 2 > n:
-        return x
-    else:
-        return x + 1
+    if (x + 1) ** 2 <= n:
+        x += 1
+
+    return x
 
 
 @export
diff --git a/src/galois/_modular.py b/src/galois/_modular.py
index c396acb4f..d6d0e43bf 100644
--- a/src/galois/_modular.py
+++ b/src/galois/_modular.py
@@ -62,8 +62,8 @@ def totatives(n: int) -> list[int]:
 
     if n == 1:
         return [0]
-    else:
-        return [t for t in range(1, n) if math.gcd(n, t) == 1]
+
+    return [t for t in range(1, n) if math.gcd(n, t) == 1]
 
 
 @export
@@ -391,15 +391,17 @@ def is_cyclic(n: int) -> bool:
 
     if n in [2, 4]:
         return True
-    elif len(p) == 2 and 2 in p and e[p.index(2)] == 1:
+
+    if len(p) == 2 and 2 in p and e[p.index(2)] == 1:
         # n = 2 * p^e
         return True
-    elif len(p) == 1 and p[0] != 2:
+
+    if len(p) == 1 and p[0] != 2:
         # n = p^e
         return True
-    else:
-        # n does not represent a cyclic group
-        return False
+
+    # n does not represent a cyclic group
+    return False
 
 
 @export
@@ -575,11 +577,12 @@ def primitive_root(
 
     try:
         if method == "min":
-            return next(primitive_roots(n, start, stop=stop))
+            root = next(primitive_roots(n, start, stop=stop))
         elif method == "max":
-            return next(primitive_roots(n, start, stop=stop, reverse=True))
+            root = next(primitive_roots(n, start, stop=stop, reverse=True))
         else:
-            return _primitive_root_random_search(n, start, stop)
+            root = _primitive_root_random_search(n, start, stop)
+        return root
     except StopIteration as e:
         raise RuntimeError(f"No primitive roots modulo {n} exist in the range [{start}, {stop}).") from e
 
diff --git a/src/galois/_polymorphic.py b/src/galois/_polymorphic.py
index 85fb0f516..f7b5aab3b 100644
--- a/src/galois/_polymorphic.py
+++ b/src/galois/_polymorphic.py
@@ -98,10 +98,9 @@ def gcd(a, b):
     """
     if isinstance(a, (int, np.integer)) and isinstance(b, (int, np.integer)):
         return int_gcd(a, b)
-    elif isinstance(a, Poly) and isinstance(b, Poly):
+    if isinstance(a, Poly) and isinstance(b, Poly):
         return poly_gcd(a, b)
-    else:
-        raise TypeError(f"Arguments `a` and `b` must both be either int or galois.Poly, not {type(a)} and {type(b)}.")
+    raise TypeError(f"Arguments `a` and `b` must both be either int or galois.Poly, not {type(a)} and {type(b)}.")
 
 
 @overload
@@ -189,10 +188,9 @@ def egcd(a, b):
     """
     if isinstance(a, (int, np.integer)) and isinstance(b, (int, np.integer)):
         return int_egcd(a, b)
-    elif isinstance(a, Poly) and isinstance(b, Poly):
+    if isinstance(a, Poly) and isinstance(b, Poly):
         return poly_egcd(a, b)
-    else:
-        raise TypeError(f"Arguments `a` and `b` must both be either int or galois.Poly, not {type(a)} and {type(b)}.")
+    raise TypeError(f"Arguments `a` and `b` must both be either int or galois.Poly, not {type(a)} and {type(b)}.")
 
 
 @overload
@@ -262,10 +260,9 @@ def lcm(*values):
 
     if all(isinstance(value, (int, np.integer)) for value in values):
         return int_lcm(*values)
-    elif all(isinstance(value, Poly) for value in values):
+    if all(isinstance(value, Poly) for value in values):
         return poly_lcm(*values)
-    else:
-        raise TypeError(f"All arguments must be either int or galois.Poly, not {[type(value) for value in values]}.")
+    raise TypeError(f"All arguments must be either int or galois.Poly, not {[type(value) for value in values]}.")
 
 
 @overload
@@ -334,10 +331,9 @@ def prod(*values):
 
     if all(isinstance(value, (int, np.integer)) for value in values):
         return int_prod(*values)
-    elif all(isinstance(value, Poly) for value in values):
+    if all(isinstance(value, Poly) for value in values):
         return poly_prod(*values)
-    else:
-        raise TypeError(f"All arguments must be either int or galois.Poly, not {[type(value) for value in values]}.")
+    raise TypeError(f"All arguments must be either int or galois.Poly, not {[type(value) for value in values]}.")
 
 
 @overload
@@ -687,10 +683,9 @@ def factors(value):
     """
     if isinstance(value, (int, np.integer)):
         return int_factors(value)
-    elif isinstance(value, Poly):
+    if isinstance(value, Poly):
         return value.factors()
-    else:
-        raise TypeError(f"Argument 'value' must be either int or Poly, not {type(value)}.")
+    raise TypeError(f"Argument 'value' must be either int or Poly, not {type(value)}.")
 
 
 @overload
@@ -778,7 +773,6 @@ def is_square_free(value):
     """
     if isinstance(value, (int, np.integer)):
         return int_is_square_free(value)
-    elif isinstance(value, Poly):
+    if isinstance(value, Poly):
         return value.is_square_free()
-    else:
-        raise TypeError(f"Argument 'value' must be either int or Poly, not {type(value)}.")
+    raise TypeError(f"Argument 'value' must be either int or Poly, not {type(value)}.")
diff --git a/src/galois/_polys/_conversions.py b/src/galois/_polys/_conversions.py
index da5438578..aaf5b1e8d 100644
--- a/src/galois/_polys/_conversions.py
+++ b/src/galois/_polys/_conversions.py
@@ -15,8 +15,8 @@ def integer_to_degree(integer: int, order: int) -> int:
     """
     if order == 2:
         return max(integer.bit_length() - 1, 0)
-    else:
-        return ilog(integer, order)
+
+    return ilog(integer, order)
 
 
 def integer_to_poly(integer: int, order: int, degree: int | None = None) -> list[int]:
diff --git a/src/galois/_polys/_dense.py b/src/galois/_polys/_dense.py
index 60496ce9c..d90f88c72 100644
--- a/src/galois/_polys/_dense.py
+++ b/src/galois/_polys/_dense.py
@@ -122,8 +122,8 @@ def multiply(a: Array, b: Array) -> Array:
     # c(x) = a(x) * b(x)
     if a.ndim == 0 or b.ndim == 0:
         return a * b
-    else:
-        return np.convolve(a, b)
+
+    return np.convolve(a, b)
 
 
 class divmod_jit(Function):
diff --git a/src/galois/_polys/_irreducible.py b/src/galois/_polys/_irreducible.py
index ae50d44ba..6c837f92f 100644
--- a/src/galois/_polys/_irreducible.py
+++ b/src/galois/_polys/_irreducible.py
@@ -82,11 +82,13 @@ def irreducible_poly(order: int, degree: int, method: Literal["min", "max", "ran
         raise ValueError(f"Argument 'method' must be in ['min', 'max', 'random'], not {method!r}.")
 
     if method == "min":
-        return next(irreducible_polys(order, degree))
+        poly = next(irreducible_polys(order, degree))
     elif method == "max":
-        return next(irreducible_polys(order, degree, reverse=True))
+        poly = next(irreducible_polys(order, degree, reverse=True))
     else:
-        return _random_search(order, degree)
+        poly = _random_search(order, degree)
+
+    return poly
 
 
 @export
diff --git a/src/galois/_polys/_poly.py b/src/galois/_polys/_poly.py
index 2f4054abc..827e5ed3c 100644
--- a/src/galois/_polys/_poly.py
+++ b/src/galois/_polys/_poly.py
@@ -721,8 +721,7 @@ def reverse(self) -> Poly:
         """
         if self._type == "sparse":
             return Poly.Degrees(self.degree - self._nonzero_degrees, self._nonzero_coeffs)
-        else:
-            return Poly(self.coeffs[::-1])
+        return Poly(self.coeffs[::-1])
 
     @overload
     def roots(self, multiplicity: Literal[False] = False) -> Array:
@@ -809,9 +808,9 @@ def roots(self, multiplicity=False):
 
         if not multiplicity:
             return roots
-        else:
-            multiplicities = np.array([_root_multiplicity(self, root) for root in roots])
-            return roots, multiplicities
+
+        multiplicities = np.array([_root_multiplicity(self, root) for root in roots])
+        return roots, multiplicities
 
     def square_free_factors(self) -> tuple[list[Poly], list[int]]:
         r"""
@@ -1276,8 +1275,8 @@ def derivative(self, k: int = 1) -> Poly:
         k -= 1
         if k > 0:
             return p_prime.derivative(k)
-        else:
-            return p_prime
+
+        return p_prime
 
     def is_irreducible(self) -> bool:
         r"""
@@ -1553,8 +1552,7 @@ def __str__(self) -> str:
         """
         if self._type == "sparse":
             return sparse_poly_to_str(self._nonzero_degrees, self._nonzero_coeffs)
-        else:
-            return poly_to_str(self.coeffs)
+        return poly_to_str(self.coeffs)
 
     def __index__(self) -> int:
         # Define __index__ to enable use of bin(), oct(), and hex()
@@ -1676,13 +1674,15 @@ def __call__(self, at, field=None, elementwise=True):
         x = field(at)  # An array of finite field elements
 
         if elementwise:
-            return _dense.evaluate_elementwise_jit(field)(coeffs, x)
+            output = _dense.evaluate_elementwise_jit(field)(coeffs, x)
         else:
             if not (x.ndim == 2 and x.shape[0] == x.shape[1]):
                 raise ValueError(
                     f"Argument 'x' must be a square matrix when evaluating the polynomial not element-wise, not have shape {x.shape}."
                 )
-            return _evaluate_matrix(coeffs, x)
+            output = _evaluate_matrix(coeffs, x)
+
+        return output
 
     def __len__(self) -> int:
         """
@@ -1775,17 +1775,19 @@ def __add__(self, other: Poly | Array) -> Poly:
             a = _convert_to_integer(self, self.field)
             b = _convert_to_integer(other, self.field)
             c = _binary.add(a, b)
-            return Poly.Int(c, field=self.field)
+            output = Poly.Int(c, field=self.field)
         elif "sparse" in types:
             a_degrees, a_coeffs = _convert_to_sparse_coeffs(self, self.field)
             b_degrees, b_coeffs = _convert_to_sparse_coeffs(other, self.field)
             c_degrees, c_coeffs = _sparse.add(a_degrees, a_coeffs, b_degrees, b_coeffs)
-            return Poly.Degrees(c_degrees, c_coeffs, field=self.field)
+            output = Poly.Degrees(c_degrees, c_coeffs, field=self.field)
         else:
             a = _convert_to_coeffs(self, self.field)
             b = _convert_to_coeffs(other, self.field)
             c = _dense.add_jit(self.field)(a, b)
-            return Poly(c, field=self.field)
+            output = Poly(c, field=self.field)
+
+        return output
 
     def __radd__(self, other: Poly | Array) -> Poly:
         _check_input_is_poly(other, self.field)
@@ -1795,31 +1797,35 @@ def __radd__(self, other: Poly | Array) -> Poly:
             a = _convert_to_integer(other, self.field)
             b = _convert_to_integer(self, self.field)
             c = _binary.add(a, b)
-            return Poly.Int(c, field=self.field)
+            output = Poly.Int(c, field=self.field)
         elif "sparse" in types:
             a_degrees, a_coeffs = _convert_to_sparse_coeffs(other, self.field)
             b_degrees, b_coeffs = _convert_to_sparse_coeffs(self, self.field)
             c_degrees, c_coeffs = _sparse.add(a_degrees, a_coeffs, b_degrees, b_coeffs)
-            return Poly.Degrees(c_degrees, c_coeffs, field=self.field)
+            output = Poly.Degrees(c_degrees, c_coeffs, field=self.field)
         else:
             a = _convert_to_coeffs(other, self.field)
             b = _convert_to_coeffs(self, self.field)
             c = _dense.add_jit(self.field)(a, b)
-            return Poly(c, field=self.field)
+            output = Poly(c, field=self.field)
+
+        return output
 
     def __neg__(self):
         if self._type == "binary":
             a = _convert_to_integer(self, self.field)
             c = _binary.negative(a)
-            return Poly.Int(c, field=self.field)
+            output = Poly.Int(c, field=self.field)
         elif self._type == "sparse":
             a_degrees, a_coeffs = _convert_to_sparse_coeffs(self, self.field)
             c_degrees, c_coeffs = _sparse.negative(a_degrees, a_coeffs)
-            return Poly.Degrees(c_degrees, c_coeffs, field=self.field)
+            output = Poly.Degrees(c_degrees, c_coeffs, field=self.field)
         else:
             a = _convert_to_coeffs(self, self.field)
             c = _dense.negative(a)
-            return Poly(c, field=self.field)
+            output = Poly(c, field=self.field)
+
+        return output
 
     def __sub__(self, other: Poly | Array) -> Poly:
         _check_input_is_poly(other, self.field)
@@ -1829,17 +1835,19 @@ def __sub__(self, other: Poly | Array) -> Poly:
             a = _convert_to_integer(self, self.field)
             b = _convert_to_integer(other, self.field)
             c = _binary.subtract(a, b)
-            return Poly.Int(c, field=self.field)
+            output = Poly.Int(c, field=self.field)
         elif "sparse" in types:
             a_degrees, a_coeffs = _convert_to_sparse_coeffs(self, self.field)
             b_degrees, b_coeffs = _convert_to_sparse_coeffs(other, self.field)
             c_degrees, c_coeffs = _sparse.subtract(a_degrees, a_coeffs, b_degrees, b_coeffs)
-            return Poly.Degrees(c_degrees, c_coeffs, field=self.field)
+            output = Poly.Degrees(c_degrees, c_coeffs, field=self.field)
         else:
             a = _convert_to_coeffs(self, self.field)
             b = _convert_to_coeffs(other, self.field)
             c = _dense.subtract_jit(self.field)(a, b)
-            return Poly(c, field=self.field)
+            output = Poly(c, field=self.field)
+
+        return output
 
     def __rsub__(self, other: Poly | Array) -> Poly:
         _check_input_is_poly(other, self.field)
@@ -1849,17 +1857,19 @@ def __rsub__(self, other: Poly | Array) -> Poly:
             a = _convert_to_integer(other, self.field)
             b = _convert_to_integer(self, self.field)
             c = _binary.subtract(a, b)
-            return Poly.Int(c, field=self.field)
+            output = Poly.Int(c, field=self.field)
         elif "sparse" in types:
             a_degrees, a_coeffs = _convert_to_sparse_coeffs(other, self.field)
             b_degrees, b_coeffs = _convert_to_sparse_coeffs(self, self.field)
             c_degrees, c_coeffs = _sparse.subtract(a_degrees, a_coeffs, b_degrees, b_coeffs)
-            return Poly.Degrees(c_degrees, c_coeffs, field=self.field)
+            output = Poly.Degrees(c_degrees, c_coeffs, field=self.field)
         else:
             a = _convert_to_coeffs(other, self.field)
             b = _convert_to_coeffs(self, self.field)
             c = _dense.subtract_jit(self.field)(a, b)
-            return Poly(c, field=self.field)
+            output = Poly(c, field=self.field)
+
+        return output
 
     def __mul__(self, other: Poly | Array | int) -> Poly:
         _check_input_is_poly_or_int(other, self.field)
@@ -1869,17 +1879,19 @@ def __mul__(self, other: Poly | Array | int) -> Poly:
             a = _convert_to_integer(self, self.field)
             b = _convert_to_integer(other, self.field)
             c = _binary.multiply(a, b)
-            return Poly.Int(c, field=self.field)
+            output = Poly.Int(c, field=self.field)
         elif "sparse" in types:
             a_degrees, a_coeffs = _convert_to_sparse_coeffs(self, self.field)
             b_degrees, b_coeffs = _convert_to_sparse_coeffs(other, self.field)
             c_degrees, c_coeffs = _sparse.multiply(a_degrees, a_coeffs, b_degrees, b_coeffs)
-            return Poly.Degrees(c_degrees, c_coeffs, field=self.field)
+            output = Poly.Degrees(c_degrees, c_coeffs, field=self.field)
         else:
             a = _convert_to_coeffs(self, self.field)
             b = _convert_to_coeffs(other, self.field)
             c = _dense.multiply(a, b)
-            return Poly(c, field=self.field)
+            output = Poly(c, field=self.field)
+
+        return output
 
     def __rmul__(self, other: Poly | Array | int) -> Poly:
         _check_input_is_poly_or_int(other, self.field)
@@ -1889,17 +1901,19 @@ def __rmul__(self, other: Poly | Array | int) -> Poly:
             a = _convert_to_integer(other, self.field)
             b = _convert_to_integer(self, self.field)
             c = _binary.multiply(a, b)
-            return Poly.Int(c, field=self.field)
+            output = Poly.Int(c, field=self.field)
         elif "sparse" in types:
             a_degrees, a_coeffs = _convert_to_sparse_coeffs(other, self.field)
             b_degrees, b_coeffs = _convert_to_sparse_coeffs(self, self.field)
             c_degrees, c_coeffs = _sparse.multiply(a_degrees, a_coeffs, b_degrees, b_coeffs)
-            return Poly.Degrees(c_degrees, c_coeffs, field=self.field)
+            output = Poly.Degrees(c_degrees, c_coeffs, field=self.field)
         else:
             a = _convert_to_coeffs(other, self.field)
             b = _convert_to_coeffs(self, self.field)
             c = _dense.multiply(a, b)
-            return Poly(c, field=self.field)
+            output = Poly(c, field=self.field)
+
+        return output
 
     def __divmod__(self, other: Poly | Array) -> tuple[Poly, Poly]:
         _check_input_is_poly(other, self.field)
@@ -1909,12 +1923,14 @@ def __divmod__(self, other: Poly | Array) -> tuple[Poly, Poly]:
             a = _convert_to_integer(self, self.field)
             b = _convert_to_integer(other, self.field)
             q, r = _binary.divmod(a, b)
-            return Poly.Int(q, field=self.field), Poly.Int(r, field=self.field)
+            output = Poly.Int(q, field=self.field), Poly.Int(r, field=self.field)
         else:
             a = _convert_to_coeffs(self, self.field)
             b = _convert_to_coeffs(other, self.field)
             q, r = _dense.divmod_jit(self.field)(a, b)
-            return Poly(q, field=self.field), Poly(r, field=self.field)
+            output = Poly(q, field=self.field), Poly(r, field=self.field)
+
+        return output
 
     def __rdivmod__(self, other: Poly | Array) -> tuple[Poly, Poly]:
         _check_input_is_poly(other, self.field)
@@ -1924,12 +1940,14 @@ def __rdivmod__(self, other: Poly | Array) -> tuple[Poly, Poly]:
             a = _convert_to_integer(other, self.field)
             b = _convert_to_integer(self, self.field)
             q, r = _binary.divmod(a, b)
-            return Poly.Int(q, field=self.field), Poly.Int(r, field=self.field)
+            output = Poly.Int(q, field=self.field), Poly.Int(r, field=self.field)
         else:
             a = _convert_to_coeffs(other, self.field)
             b = _convert_to_coeffs(self, self.field)
             q, r = _dense.divmod_jit(self.field)(a, b)
-            return Poly(q, field=self.field), Poly(r, field=self.field)
+            output = Poly(q, field=self.field), Poly(r, field=self.field)
+
+        return output
 
     def __truediv__(self, other):
         raise NotImplementedError(
@@ -1951,12 +1969,14 @@ def __floordiv__(self, other: Poly | Array) -> Poly:
             a = _convert_to_integer(self, self.field)
             b = _convert_to_integer(other, self.field)
             q = _binary.floordiv(a, b)
-            return Poly.Int(q, field=self.field)
+            output = Poly.Int(q, field=self.field)
         else:
             a = _convert_to_coeffs(self, self.field)
             b = _convert_to_coeffs(other, self.field)
             q = _dense.floordiv_jit(self.field)(a, b)
-            return Poly(q, field=self.field)
+            output = Poly(q, field=self.field)
+
+        return output
 
     def __rfloordiv__(self, other: Poly | Array) -> Poly:
         _check_input_is_poly(other, self.field)
@@ -1966,12 +1986,14 @@ def __rfloordiv__(self, other: Poly | Array) -> Poly:
             a = _convert_to_integer(other, self.field)
             b = _convert_to_integer(self, self.field)
             q = _binary.floordiv(a, b)
-            return Poly.Int(q, field=self.field)
+            output = Poly.Int(q, field=self.field)
         else:
             a = _convert_to_coeffs(other, self.field)
             b = _convert_to_coeffs(self, self.field)
             q = _dense.floordiv_jit(self.field)(a, b)
-            return Poly(q, field=self.field)
+            output = Poly(q, field=self.field)
+
+        return output
 
     def __mod__(self, other: Poly | Array) -> Poly:
         _check_input_is_poly(other, self.field)
@@ -1981,12 +2003,14 @@ def __mod__(self, other: Poly | Array) -> Poly:
             a = _convert_to_integer(self, self.field)
             b = _convert_to_integer(other, self.field)
             r = _binary.mod(a, b)
-            return Poly.Int(r, field=self.field)
+            output = Poly.Int(r, field=self.field)
         else:
             a = _convert_to_coeffs(self, self.field)
             b = _convert_to_coeffs(other, self.field)
             r = _dense.mod_jit(self.field)(a, b)
-            return Poly(r, field=self.field)
+            output = Poly(r, field=self.field)
+
+        return output
 
     def __rmod__(self, other: Poly | Array) -> Poly:
         _check_input_is_poly(other, self.field)
@@ -1996,12 +2020,14 @@ def __rmod__(self, other: Poly | Array) -> Poly:
             a = _convert_to_integer(other, self.field)
             b = _convert_to_integer(self, self.field)
             r = _binary.mod(a, b)
-            return Poly.Int(r, field=self.field)
+            output = Poly.Int(r, field=self.field)
         else:
             a = _convert_to_coeffs(other, self.field)
             b = _convert_to_coeffs(self, self.field)
             r = _dense.mod_jit(self.field)(a, b)
-            return Poly(r, field=self.field)
+            output = Poly(r, field=self.field)
+
+        return output
 
     def __pow__(self, exponent: int, modulus: Poly | None = None) -> Poly:
         _check_input_is_poly_or_none(modulus, self.field)
@@ -2015,12 +2041,14 @@ def __pow__(self, exponent: int, modulus: Poly | None = None) -> Poly:
             a = _convert_to_integer(self, self.field)
             b = _convert_to_integer(modulus, self.field) if modulus is not None else None
             q = _binary.pow(a, exponent, b)
-            return Poly.Int(q, field=self.field)
+            output = Poly.Int(q, field=self.field)
         else:
             a = _convert_to_coeffs(self, self.field)
             b = _convert_to_coeffs(modulus, self.field) if modulus is not None else None
             q = _dense.pow_jit(self.field)(a, exponent, b)
-            return Poly(q, field=self.field)
+            output = Poly(q, field=self.field)
+
+        return output
 
     ###############################################################################
     # Instance properties
@@ -2305,12 +2333,14 @@ def _convert_to_coeffs(a: Poly | Array | int, field: Type[Array]) -> Array:
     Convert the polynomial or finite field scalar into a coefficient array.
     """
     if isinstance(a, Poly):
-        return a.coeffs
+        coeffs = a.coeffs
     elif isinstance(a, int):
         # Scalar multiplication
-        return np.atleast_1d(field(a % field.characteristic))
+        coeffs = np.atleast_1d(field(a % field.characteristic))
     else:
-        return np.atleast_1d(a)
+        coeffs = np.atleast_1d(a)
+
+    return coeffs
 
 
 def _convert_to_integer(a: Poly | Array | int, field: Type[Array]) -> int:
@@ -2319,9 +2349,11 @@ def _convert_to_integer(a: Poly | Array | int, field: Type[Array]) -> int:
     """
     if isinstance(a, int):
         # Scalar multiplication
-        return a % field.characteristic
+        integer = a % field.characteristic
     else:
-        return int(a)
+        integer = int(a)
+
+    return integer
 
 
 def _convert_to_sparse_coeffs(a: Poly | Array | int, field: Type[Array]) -> tuple[np.ndarray, Array]:
@@ -2329,9 +2361,14 @@ def _convert_to_sparse_coeffs(a: Poly | Array | int, field: Type[Array]) -> tupl
     Convert the polynomial or finite field scalar into its non-zero degrees and coefficients.
     """
     if isinstance(a, Poly):
-        return a.nonzero_degrees, a.nonzero_coeffs
+        degrees = a.nonzero_degrees
+        coeffs = a.nonzero_coeffs
     elif isinstance(a, int):
         # Scalar multiplication
-        return np.array([0]), np.atleast_1d(field(a % field.characteristic))
+        degrees = np.array([0])
+        coeffs = np.atleast_1d(field(a % field.characteristic))
     else:
-        return np.array([0]), np.atleast_1d(a)
+        degrees = np.array([0])
+        coeffs = np.atleast_1d(a)
+
+    return degrees, coeffs
diff --git a/src/galois/_polys/_primitive.py b/src/galois/_polys/_primitive.py
index 2cd0464dd..fe25af86a 100644
--- a/src/galois/_polys/_primitive.py
+++ b/src/galois/_polys/_primitive.py
@@ -102,11 +102,13 @@ def primitive_poly(order: int, degree: int, method: Literal["min", "max", "rando
         raise ValueError(f"Argument 'method' must be in ['min', 'max', 'random'], not {method!r}.")
 
     if method == "min":
-        return next(primitive_polys(order, degree))
+        poly = next(primitive_polys(order, degree))
     elif method == "max":
-        return next(primitive_polys(order, degree, reverse=True))
+        poly = next(primitive_polys(order, degree, reverse=True))
     else:
-        return _random_search(order, degree)
+        poly = _random_search(order, degree)
+
+    return poly
 
 
 @export
@@ -372,11 +374,13 @@ def matlab_primitive_poly(characteristic: int, degree: int) -> Poly:
     if characteristic == 2 and degree == 7:
         # Not the lexicographically-minimal of `x^7 + x + 1`
         return Poly.Degrees([7, 3, 0])
-    elif characteristic == 2 and degree == 14:
+
+    if characteristic == 2 and degree == 14:
         # Not the lexicographically-minimal of `x^14 + x^5 + x^3 + x + 1`
         return Poly.Degrees([14, 10, 6, 1, 0])
-    elif characteristic == 2 and degree == 16:
+
+    if characteristic == 2 and degree == 16:
         # Not the lexicographically-minimal of `x^16 + x^5 + x^3 + x^2 + 1`
         return Poly.Degrees([16, 12, 3, 1, 0])
-    else:
-        return primitive_poly(characteristic, degree)
+
+    return primitive_poly(characteristic, degree)
diff --git a/src/galois/_prime.py b/src/galois/_prime.py
index ef2eded4a..2f1e2c1d4 100644
--- a/src/galois/_prime.py
+++ b/src/galois/_prime.py
@@ -798,10 +798,10 @@ def jacobi_symbol(a: int, n: int) -> int:
         s = -s
 
     n1 = n % a1
-    if a1 == 1:
-        return s
-    else:
-        return s * jacobi_symbol(n1, a1)
+    if a1 != 1:
+        s *= jacobi_symbol(n1, a1)
+
+    return s
 
 
 @export
@@ -846,10 +846,9 @@ def kronecker_symbol(a: int, n: int) -> int:
     if n == 2:
         if a % 2 == 0:
             return 0
-        elif a % 8 in [1, 7]:
+        if a % 8 in [1, 7]:
             return 1
-        else:
-            return -1
+        return -1
 
     # Factor out the unit +/- 1
     u = -1 if n < 0 else 1
@@ -860,11 +859,12 @@ def kronecker_symbol(a: int, n: int) -> int:
     while n % 2 == 0:
         n, e = n // 2, e + 1
 
+    s = kronecker_symbol(a, u) * kronecker_symbol(a, 2) ** e
     if n >= 3:
         # Handle the remaining odd n using the Jacobi symbol
-        return kronecker_symbol(a, u) * kronecker_symbol(a, 2) ** e * jacobi_symbol(a, n)
-    else:
-        return kronecker_symbol(a, u) * kronecker_symbol(a, 2) ** e
+        s *= jacobi_symbol(a, n)
+
+    return s
 
 
 ###############################################################################
@@ -1454,8 +1454,8 @@ def divisor_sigma(n: int, k: int = 1) -> int:
 
     if n == 0:
         return len(d)
-    else:
-        return sum(di**k for di in d)
+
+    return sum(di**k for di in d)
 
 
 ###############################################################################
@@ -1516,7 +1516,7 @@ def is_prime(n: int) -> bool:
     for p in PRIMES[0:250]:
         if n == p:
             return True
-        elif n % p == 0:
+        if n % p == 0:
             return False
 
     if not fermat_primality_test(n):