diff --git a/pyproject.toml b/pyproject.toml
index 3640b344d..0186afb48 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -50,8 +50,8 @@ classifiers = [
 ]
 requires-python = ">=3.7"
 dependencies = [
-    "numpy >= 1.21.0, < 1.27",    # v1.21.0 is needed for dtype support of ufuncs, see https://numpy.org/devdocs/release/1.21.0-notes.html#ufunc-signature-and-dtype-generalization-and-casting
-    "numba >= 0.55, < 0.60",      # v0.55 is needed for support of NumPy 1.21
+    "numpy >= 1.21.0, < 2.1",     # v1.21.0 is needed for dtype support of ufuncs, see https://numpy.org/devdocs/release/1.21.0-notes.html#ufunc-signature-and-dtype-generalization-and-casting
+    "numba >= 0.55, < 0.61",      # v0.55 is needed for support of NumPy 1.21
     "typing_extensions >= 4.0.0", # v4.0.0 is needed for use of Self (Python 3.11+) and Literal (Python 3.8+)
 ]
 dynamic = ["version"]
diff --git a/src/galois/_domains/_function.py b/src/galois/_domains/_function.py
index 0460043cd..eaf7d4a9f 100644
--- a/src/galois/_domains/_function.py
+++ b/src/galois/_domains/_function.py
@@ -308,7 +308,7 @@ class FunctionMixin(np.ndarray, metaclass=ArrayMeta):
         np.unpackbits,
         np.unwrap,
         np.around,
-        np.round_,
+        np.round,
         np.fix,
         np.gradient,
         np.trapz,
diff --git a/src/galois/_domains/_ufunc.py b/src/galois/_domains/_ufunc.py
index 72168b6da..e89d2456e 100644
--- a/src/galois/_domains/_ufunc.py
+++ b/src/galois/_domains/_ufunc.py
@@ -394,7 +394,11 @@ def __call__(self, ufunc, method, inputs, kwargs, meta):
             self._verify_operands_in_field_or_int(ufunc, inputs, meta)
             inputs, kwargs = self._view_inputs_as_ndarray(inputs, kwargs)
             i = meta["non_field_operands"][0]  # Scalar multiplicand
-            inputs[i] = np.mod(inputs[i], self.field.characteristic)
+            if meta["dtype"] == np.object_:
+                # Need to explicitly cast to np.object_ in NumPy v2.0 or the integer will overflow
+                inputs[i] = np.mod(inputs[i], self.field.characteristic, dtype=np.object_)
+            else:
+                inputs[i] = np.mod(inputs[i], self.field.characteristic)
         inputs, kwargs = self._view_inputs_as_ndarray(inputs, kwargs)
         output = getattr(self.ufunc, method)(*inputs, **kwargs)
         output = self._view_output_as_field(output, self.field, meta["dtype"])
diff --git a/src/galois/_fields/_array.py b/src/galois/_fields/_array.py
index 9f1173f32..2895a4563 100644
--- a/src/galois/_fields/_array.py
+++ b/src/galois/_fields/_array.py
@@ -989,7 +989,7 @@ def additive_order(self) -> int | np.ndarray:
         if x.ndim == 0:
             order = 1 if x == 0 else field.characteristic
         else:
-            order = field.characteristic * np.ones(x.shape, dtype=np.int64)
+            order = field.characteristic * np.ones(x.shape, dtype=field.dtypes[-1])
             order[np.where(x == 0)] = 1
 
         return order
diff --git a/src/galois/_polys/_dense.py b/src/galois/_polys/_dense.py
index b1b3e289c..6edd1e80e 100644
--- a/src/galois/_polys/_dense.py
+++ b/src/galois/_polys/_dense.py
@@ -487,7 +487,7 @@ def implementation(nonzero_degrees, nonzero_coeffs, primitive_element):  # pragm
         # Test if 0 is a root
         if nonzero_degrees[-1] != 0:
             roots.append(0)
-            powers.append(-1)
+            powers.append(nonzero_degrees[-1])  # 0 has multiplicity equal to the lowest degree of x that is non-zero
 
         # Test if 1 is a root
         _sum = 0
@@ -510,4 +510,4 @@ def implementation(nonzero_degrees, nonzero_coeffs, primitive_element):  # pragm
             if len(roots) == degree:
                 break
 
-        return np.array([roots, powers], dtype=nonzero_coeffs.dtype)
+        return np.array([roots, powers])