diff --git a/numpy_groupies/aggregate_numpy.py b/numpy_groupies/aggregate_numpy.py
index f971176..1eff227 100644
--- a/numpy_groupies/aggregate_numpy.py
+++ b/numpy_groupies/aggregate_numpy.py
@@ -250,7 +250,10 @@ def _cumsum(group_idx, a, size, fill_value=None, dtype=None):
 
     increasing = np.arange(len(a), dtype=int)
     group_starts = _min(group_idx_srt, increasing, size, fill_value=0)[group_idx_srt]
-    a_srt_cumsum += -a_srt_cumsum[group_starts] + a_srt[group_starts]
+    # First subtract large numbers
+    a_srt_cumsum -= a_srt_cumsum[group_starts]
+    # Then add potentially small numbers
+    a_srt_cumsum += a_srt[group_starts]
     return a_srt_cumsum[invsortidx]
 
 
diff --git a/numpy_groupies/tests/test_generic.py b/numpy_groupies/tests/test_generic.py
index 2a2811d..7a8b772 100644
--- a/numpy_groupies/tests/test_generic.py
+++ b/numpy_groupies/tests/test_generic.py
@@ -570,3 +570,14 @@ def test_var_with_nan_fill_value(aggregate_all, ddof, nan_inds, func):
         group_idx, a, axis=-1, fill_value=np.nan, func=func, ddof=ddof
     )
     np.testing.assert_equal(actual, expected)
+
+
+def test_cumsum_accuracy(aggregate_all):
+    array = np.array(
+        [0.00000000e00, 0.00000000e00, 0.00000000e00, 3.27680000e04, 9.99999975e-06]
+    )
+    group_idx = np.array([0, 0, 0, 0, 1])
+
+    actual = aggregate_all(group_idx, array, axis=-1, func="cumsum")
+    expected = array
+    np.testing.assert_allclose(actual, expected)