diff --git a/numpower.c b/numpower.c
index fcc4989..bede2d7 100644
--- a/numpower.c
+++ b/numpower.c
@@ -4750,6 +4750,60 @@ PHP_METHOD(NDArray, prod) {
     RETURN_NDARRAY(rtn, return_value);
 }
 
+/**
+ * NDArray::cumProd
+ */
+ZEND_BEGIN_ARG_INFO_EX(arginfo_ndarray_cumProd, 0, 0, 1)
+ZEND_ARG_INFO(0, a)
+ZEND_ARG_INFO(0, axis)
+ZEND_END_ARG_INFO()
+PHP_METHOD(NDArray, cumProd) {
+    NDArray *rtn = NULL;
+    zval *a;
+    long axis = -1;
+    int axis_i;
+    ZEND_PARSE_PARAMETERS_START(1, 2)
+    Z_PARAM_ZVAL(a)
+    Z_PARAM_OPTIONAL
+    Z_PARAM_LONG(axis)
+    ZEND_PARSE_PARAMETERS_END();
+    axis_i = (int)axis;
+    NDArray *nda = ZVAL_TO_NDARRAY(a);
+    if (nda == NULL) {
+        return;
+    }
+    rtn = NDArray_Cum_Prod(nda, &axis_i);
+    CHECK_INPUT_AND_FREE(a, nda);
+    RETURN_NDARRAY(rtn, return_value);
+}
+
+/**
+ * NDArray::cumSum
+ */
+ZEND_BEGIN_ARG_INFO_EX(arginfo_ndarray_cumSum, 0, 0, 1)
+ZEND_ARG_INFO(0, a)
+ZEND_ARG_INFO(0, axis)
+ZEND_END_ARG_INFO()
+PHP_METHOD(NDArray, cumSum) {
+    NDArray *rtn = NULL;
+    zval *a;
+    long axis = -1;
+    int axis_i;
+    ZEND_PARSE_PARAMETERS_START(1, 2)
+    Z_PARAM_ZVAL(a)
+    Z_PARAM_OPTIONAL
+    Z_PARAM_LONG(axis)
+    ZEND_PARSE_PARAMETERS_END();
+    axis_i = (int)axis;
+    NDArray *nda = ZVAL_TO_NDARRAY(a);
+    if (nda == NULL) {
+        return;
+    }
+    rtn = NDArray_Cum_Sum(nda, &axis_i);
+    CHECK_INPUT_AND_FREE(a, nda);
+    RETURN_NDARRAY(rtn, return_value);
+}
+
 ZEND_BEGIN_ARG_INFO(arginfo_ndarray_array, 0)
 ZEND_ARG_INFO(0, a)
 ZEND_END_ARG_INFO()
@@ -5190,6 +5244,9 @@ static const zend_function_entry class_NDArray_methods[] = {
     ZEND_ME(NDArray, sum, arginfo_ndarray_sum, ZEND_ACC_PUBLIC | ZEND_ACC_STATIC)
     ZEND_ME(NDArray, prod, arginfo_ndarray_prod, ZEND_ACC_PUBLIC | ZEND_ACC_STATIC)
     ZEND_ME(NDArray, mod, arginfo_ndarray_mod, ZEND_ACC_PUBLIC | ZEND_ACC_STATIC)
+    ZEND_ME(NDArray, cumProd, arginfo_ndarray_cumProd, ZEND_ACC_PUBLIC | ZEND_ACC_STATIC)
+    ZEND_ME(NDArray, cumSum, arginfo_ndarray_cumSum, ZEND_ACC_PUBLIC | ZEND_ACC_STATIC)
+
 
     ZEND_ME(NDArray, size, arginfo_size, ZEND_ACC_PUBLIC)
     ZEND_ME(NDArray, count, arginfo_count, ZEND_ACC_PUBLIC)
diff --git a/src/ndmath/arithmetics.c b/src/ndmath/arithmetics.c
index 158c5f0..333ac75 100644
--- a/src/ndmath/arithmetics.c
+++ b/src/ndmath/arithmetics.c
@@ -945,3 +945,105 @@ NDArray_Abs(NDArray *nda) {
     }
     return rtn;
 }
+
+NDArray*
+NDArray_Cum_Flat(NDArray *a, ElementWiseFloatOperation1F op) {
+    NDArray *rtn = NDArray_Copy(a, NDArray_DEVICE(a));
+    int num_elements = NDArray_NUMELEMENTS(rtn);
+    int *flat_shape = emalloc(sizeof(int) * 2);
+    flat_shape[0] = 1;
+    flat_shape[1] = num_elements;
+    rtn = NDArray_Reshape(rtn, flat_shape, 2);
+    for (int i = 1; i < num_elements; i++) {
+        NDArray_FDATA(rtn)[i] = op(NDArray_FDATA(rtn)[i], NDArray_FDATA(rtn)[i - 1]);
+    }
+    return rtn;
+}
+
+NDArray*
+NDArray_Cum_Axis(NDArray *a, int *axis, NDArray *(*operation)(NDArray *, NDArray *)) {
+    if (*axis == 1) {
+        a = NDArray_Transpose(a, NULL);
+    }
+    
+    int rows = NDArray_SHAPE(a)[0];
+    int cols = NDArray_SHAPE(a)[1];
+    // setup indices for slicing
+    int *indices_shape = emalloc(sizeof(int) * 2);
+    indices_shape[0] = 2;
+    indices_shape[1] = 1;
+    NDArray** indices_axis = emalloc(sizeof(NDArray*) * 2);
+    indices_axis[0] =  NDArray_Zeros(indices_shape, 1, NDARRAY_TYPE_FLOAT32, NDARRAY_DEVICE_CPU);
+    indices_axis[1] =  NDArray_Zeros(indices_shape, 1, NDARRAY_TYPE_FLOAT32, NDARRAY_DEVICE_CPU);
+    NDArray **row_array = emalloc(sizeof(NDArray*) * rows);
+    // Set Column
+    NDArray_FDATA(indices_axis[1])[0] = 0;
+    NDArray_FDATA(indices_axis[1])[1] = cols;
+    for (int i = 0; i < rows; i++) {
+        // Set Row
+        NDArray_FDATA(indices_axis[0])[0] = i;
+        NDArray_FDATA(indices_axis[0])[1] = i + 1;
+        if (i == 0) {
+            row_array[i] = NDArray_Slice(a, indices_axis, 2);
+        } else {
+            NDArray *curr_row = NDArray_Slice(a, indices_axis, 2);
+            row_array[i] = operation(row_array[i-1], curr_row);
+            NDArray_FREE(curr_row);
+        }
+    }
+    efree(indices_axis[0]);
+    efree(indices_axis[1]);
+    efree(indices_axis);
+    efree(indices_shape);
+    NDArray *combined = NDArray_Concatenate(row_array, rows, 0);
+    for (int i = 0; i < rows; i++) {
+        NDArray_FREE(row_array[i]);
+    }
+    efree(row_array);
+    if (*axis == 1) {
+        combined = NDArray_Transpose(combined, NULL);
+    }
+    return combined;
+}
+
+/**
+ * NDArray::cumprod
+ *
+ * @param a
+ * @param axis
+ * @return
+ */
+NDArray*
+NDArray_Cum_Prod(NDArray *a, int *axis) {
+    NDArray *rtn = NULL;
+    
+    if (*axis == -1) {
+        rtn = NDArray_Cum_Flat(a, float_product);
+    } else {
+        rtn = NDArray_Cum_Axis(a, axis, NDArray_Multiply_Float);
+    }
+
+    return rtn;
+}
+
+/**
+ * NDArray::cumsum
+ *
+ * @param a
+ * @param axis
+ * @return
+ */
+NDArray*
+NDArray_Cum_Sum(NDArray *a, int *axis) {
+    NDArray *rtn = NULL;
+
+    if (*axis == -1) {
+        rtn = NDArray_Cum_Flat(a, float_sum);
+    } else {
+        rtn = NDArray_Cum_Axis(a, axis, NDArray_Add_Float);
+    }
+
+    return rtn;
+}
+
+
diff --git a/src/ndmath/arithmetics.h b/src/ndmath/arithmetics.h
index 5aeb4da..332a4f9 100644
--- a/src/ndmath/arithmetics.h
+++ b/src/ndmath/arithmetics.h
@@ -15,4 +15,6 @@ float NDArray_Mean_Float(NDArray* a);
 float NDArray_Mean_Float_Axis(NDArray* a, NDArray *b);
 NDArray* NDArray_Abs(NDArray *nda);
 float NDArray_Median_Float(NDArray* a);
+NDArray* NDArray_Cum_Prod(NDArray* a, int *axis);
+NDArray* NDArray_Cum_Sum(NDArray* a, int *axis);
 #endif //PHPSCI_NDARRAY_ARITHMETICS_H
diff --git a/src/ndmath/double_math.c b/src/ndmath/double_math.c
index 80a42b8..51d172f 100644
--- a/src/ndmath/double_math.c
+++ b/src/ndmath/double_math.c
@@ -262,4 +262,12 @@ float float_arctan2(float x, float y) {
 
 float float_reciprocal(float val) {
     return 1 / val;
+}
+
+float float_sum(float val1, float val2) {
+    return val1 + val2;
+}
+
+float float_product(float val1, float val2) {
+    return val1 * val2;
 }
\ No newline at end of file
diff --git a/src/ndmath/double_math.h b/src/ndmath/double_math.h
index e577214..5e85639 100644
--- a/src/ndmath/double_math.h
+++ b/src/ndmath/double_math.h
@@ -42,4 +42,6 @@ float float_rsqrt(float val);
 float float_arctan2(float x, float y);
 float float_positive(float val);
 float float_reciprocal(float val);
+float float_sum(float val1, float val2);
+float float_product(float val1, float val2);
 #endif //PHPSCI_NDARRAY_DOUBLE_MATH_H
diff --git a/tests/math/045-ndarray-cumsum.phpt b/tests/math/045-ndarray-cumsum.phpt
new file mode 100644
index 0000000..f02d4c5
--- /dev/null
+++ b/tests/math/045-ndarray-cumsum.phpt
@@ -0,0 +1,76 @@
+--TEST--
+NDArray::cumsum
+--FILE--
+<?php
+$empty = \NDArray::array([[]]);
+$single = \NDArray::array([[3]]);
+$a = \NDArray::array([[1, 2, 3], [4, 5, 6]]);
+print_r(\NDArray::cumsum($empty)->toArray());
+print_r(\NDArray::cumsum($single)->toArray());
+print_r(\NDArray::cumsum($a)->toArray());
+print_r(\NDArray::cumsum($a, 0)->toArray());
+print_r(\NDArray::cumsum($a, 1)->toArray());
+?>
+--EXPECT--
+Array
+(
+    [0] => Array
+        (
+        )
+
+)
+Array
+(
+    [0] => Array
+        (
+            [0] => 3
+        )
+
+)
+Array
+(
+    [0] => Array
+        (
+            [0] => 1
+            [1] => 3
+            [2] => 6
+            [3] => 10
+            [4] => 15
+            [5] => 21
+        )
+
+)
+Array
+(
+    [0] => Array
+        (
+            [0] => 1
+            [1] => 2
+            [2] => 3
+        )
+
+    [1] => Array
+        (
+            [0] => 5
+            [1] => 7
+            [2] => 9
+        )
+
+)
+Array
+(
+    [0] => Array
+        (
+            [0] => 1
+            [1] => 3
+            [2] => 6
+        )
+
+    [1] => Array
+        (
+            [0] => 4
+            [1] => 9
+            [2] => 15
+        )
+
+)
\ No newline at end of file
diff --git a/tests/math/046-ndarray-cumprod.phpt b/tests/math/046-ndarray-cumprod.phpt
new file mode 100644
index 0000000..4d0587b
--- /dev/null
+++ b/tests/math/046-ndarray-cumprod.phpt
@@ -0,0 +1,76 @@
+--TEST--
+NDArray::cumprod
+--FILE--
+<?php
+$empty = \NDArray::array([[]]);
+$single = \NDArray::array([[3]]);
+$a = \NDArray::array([[1, 2, 3], [4, 5, 6]]);
+print_r(\NDArray::cumprod($empty)->toArray());
+print_r(\NDArray::cumprod($single)->toArray());
+print_r(\NDArray::cumprod($a)->toArray());
+print_r(\NDArray::cumprod($a, 0)->toArray());
+print_r(\NDArray::cumprod($a, 1)->toArray());
+?>
+--EXPECT--
+Array
+(
+    [0] => Array
+        (
+        )
+
+)
+Array
+(
+    [0] => Array
+        (
+            [0] => 3
+        )
+
+)
+Array
+(
+    [0] => Array
+        (
+            [0] => 1
+            [1] => 2
+            [2] => 6
+            [3] => 24
+            [4] => 120
+            [5] => 720
+        )
+
+)
+Array
+(
+    [0] => Array
+        (
+            [0] => 1
+            [1] => 2
+            [2] => 3
+        )
+
+    [1] => Array
+        (
+            [0] => 4
+            [1] => 10
+            [2] => 18
+        )
+
+)
+Array
+(
+    [0] => Array
+        (
+            [0] => 1
+            [1] => 2
+            [2] => 6
+        )
+
+    [1] => Array
+        (
+            [0] => 4
+            [1] => 20
+            [2] => 120
+        )
+
+)
\ No newline at end of file