diff --git a/stumpy/core.py b/stumpy/core.py
index af3bf7a3b..085397296 100644
--- a/stumpy/core.py
+++ b/stumpy/core.py
@@ -2254,6 +2254,7 @@ def preprocess_diagonal(
 
     M_T, Σ_T = compute_mean_std(T, m)
     Σ_T[T_subseq_isconstant] = 1.0  # Avoid divide by zero in next inversion step
+    Σ_T[~T_subseq_isfinite] = np.nan
     Σ_T_inverse = 1.0 / Σ_T
     M_T_m_1, _ = compute_mean_std(T, m - 1)
 
diff --git a/stumpy/gpu_ostinato.py b/stumpy/gpu_ostinato.py
index cdcda2e5a..5effc933b 100644
--- a/stumpy/gpu_ostinato.py
+++ b/stumpy/gpu_ostinato.py
@@ -113,9 +113,10 @@ def gpu_ostinato(Ts, m, device_id=0, normalize=True, p=2.0, Ts_subseq_isconstant
     M_Ts = [None] * len(Ts)
     Σ_Ts = [None] * len(Ts)
     for i, T in enumerate(Ts):
-        Ts_copy[i], M_Ts[i], Σ_Ts[i], Ts_subseq_isconstant[i] = core.preprocess(
+        _, M_Ts[i], Σ_Ts[i], Ts_subseq_isconstant[i] = core.preprocess(
             T, m, copy=True, T_subseq_isconstant=Ts_subseq_isconstant[i]
         )
+        Ts_copy[i] = T.copy()
 
     bsf_radius, bsf_Ts_idx, bsf_subseq_idx = _ostinato(
         Ts_copy,
diff --git a/stumpy/ostinato.py b/stumpy/ostinato.py
index 424da525f..746c8ea0e 100644
--- a/stumpy/ostinato.py
+++ b/stumpy/ostinato.py
@@ -58,10 +58,9 @@ def _across_series_nearest_neighbors(
     nns_subseq_idx = np.zeros(k, dtype=np.int64)
 
     for i in range(k):
-        QT = core.sliding_dot_product(Ts[Ts_idx][subseq_idx : subseq_idx + m], Ts[i])
-        distance_profile = core._mass(
-            Q,
-            Ts[i],
+        QT = core.sliding_dot_product(Q, Ts[i])
+        distance_profile = core.calculate_distance_profile(
+            m,
             QT,
             M_Ts[Ts_idx][subseq_idx],
             Σ_Ts[Ts_idx][subseq_idx],
@@ -256,9 +255,8 @@ def _ostinato(
                         (
                             radius,
                             np.min(
-                                core._mass(
-                                    Ts[j][q : q + m],
-                                    Ts[i],
+                                core.calculate_distance_profile(
+                                    m,
                                     QT,
                                     M_Ts[j][q],
                                     Σ_Ts[j][q],
@@ -373,9 +371,10 @@ def ostinato(Ts, m, normalize=True, p=2.0, Ts_subseq_isconstant=None):
         Ts_subseq_isconstant = [None] * len(Ts)
 
     for i, T in enumerate(Ts):
-        Ts_copy[i], M_Ts[i], Σ_Ts[i], Ts_subseq_isconstant[i] = core.preprocess(
+        _, M_Ts[i], Σ_Ts[i], Ts_subseq_isconstant[i] = core.preprocess(
             T, m, copy=True, T_subseq_isconstant=Ts_subseq_isconstant[i]
         )
+        Ts_copy[i] = T.copy()
 
     bsf_radius, bsf_Ts_idx, bsf_subseq_idx = _ostinato(
         Ts_copy, m, M_Ts, Σ_Ts, Ts_subseq_isconstant
@@ -516,9 +515,10 @@ def ostinatoed(client, Ts, m, normalize=True, p=2.0, Ts_subseq_isconstant=None):
     if Ts_subseq_isconstant is None:
         Ts_subseq_isconstant = [None] * len(Ts)
     for i, T in enumerate(Ts):
-        Ts_copy[i], M_Ts[i], Σ_Ts[i], Ts_subseq_isconstant[i] = core.preprocess(
+        _, M_Ts[i], Σ_Ts[i], Ts_subseq_isconstant[i] = core.preprocess(
             T, m, copy=True, T_subseq_isconstant=Ts_subseq_isconstant[i]
         )
+        Ts_copy[i] = T.copy()
 
     bsf_radius, bsf_Ts_idx, bsf_subseq_idx = _ostinato(
         Ts_copy,
diff --git a/tests/test_core.py b/tests/test_core.py
index 0e84cf8c3..01267251e 100644
--- a/tests/test_core.py
+++ b/tests/test_core.py
@@ -880,9 +880,11 @@ def test_preprocess_non_normalized():
 def test_preprocess_diagonal():
     T = np.array([0, np.nan, 2, 3, 4, 5, 6, 7, np.inf, 9])
     m = 3
+    T_subseq_isfinite = core.rolling_isfinite(T, m)
 
     ref_T = np.array([0, 0, 2, 3, 4, 5, 6, 7, 0, 9], dtype=float)
     ref_M, ref_Σ = naive.compute_mean_std(ref_T, m)
+    ref_Σ[~T_subseq_isfinite] = np.nan
     ref_Σ_inverse = 1.0 / ref_Σ
     ref_M_m_1, _ = naive.compute_mean_std(ref_T, m - 1)
 
@@ -1753,3 +1755,13 @@ def test_process_isconstant_2d():
     T_subseq_isconstant_comp = core.process_isconstant(T, m, T_subseq_isconstant)
 
     npt.assert_array_equal(T_subseq_isconstant_ref, T_subseq_isconstant_comp)
+
+
+@pytest.mark.filterwarnings(
+    "error:divide by zero encountered in divide", category=RuntimeWarning
+)
+def test_preprocess_diagonal_divide_by_zero():
+    T = np.random.rand(64)
+    m = 3
+    T[:m] = np.nan
+    core.preprocess_diagonal(T, m)
diff --git a/tests/test_ostinato.py b/tests/test_ostinato.py
index 7916a34d8..7e0181de5 100644
--- a/tests/test_ostinato.py
+++ b/tests/test_ostinato.py
@@ -249,3 +249,24 @@ def test_extract_several_consensus_ostinatoed(dask_cluster):
                     Ts_ref[i][np.isfinite(Ts_ref[i])],
                     Ts_comp[i][np.isfinite(Ts_comp[i])],
                 )
+
+
+@pytest.mark.filterwarnings(
+    "error:divide by zero encountered in divide", category=RuntimeWarning
+)
+def test_divide_by_zero_ostinato():
+    Ts = [np.random.rand(n) for n in [64, 128, 256]]
+    m = 5
+    Ts[0][:m] = np.nan
+    stumpy.ostinato(Ts, m)
+
+
+@pytest.mark.filterwarnings(
+    "error:divide by zero encountered in divide", category=RuntimeWarning
+)
+def test_divide_by_zero_ostinatoed(dask_cluster):
+    Ts = [np.random.rand(n) for n in [64, 128, 256]]
+    m = 5
+    Ts[0][:m] = np.nan
+    with Client(dask_cluster) as dask_client:
+        stumpy.ostinatoed(dask_client, Ts, m)