diff --git a/include/pal_math.h b/include/pal_math.h
index 46713b7..bb12005 100644
--- a/include/pal_math.h
+++ b/include/pal_math.h
@@ -240,6 +240,11 @@ void p_min_f32(const float *a, float *c, int *index, int n);
  *
  ****************************************************************
  */
+/*means of columns of matrix: c = mean ( a[, n-1:0] ) */
+void p_colmeans_f32(const float *a, float *c, size_t m, size_t n);
+
+/*sum of columns of matrix: c = sum ( a[, n-1:0] ) */
+void p_colsums_f32(const float *a, float *c, size_t m, size_t n);
 
 /*sort an array*/
 void p_sort_f32(const float *a, float *c, int n);
diff --git a/src/math/Makefile.am b/src/math/Makefile.am
index 9232fbf..eb3edef 100644
--- a/src/math/Makefile.am
+++ b/src/math/Makefile.am
@@ -16,6 +16,8 @@ libpal_math_la_SOURCES = \
     p_atan.c \
     p_atanh.c \
     p_cbrt.c \
+    p_colmeans.c \
+    p_colsums.c \
     p_cos.c \
     p_cosh.c \
     p_div.c \
diff --git a/src/math/p_colmeans.c b/src/math/p_colmeans.c
new file mode 100644
index 0000000..0450fc8
--- /dev/null
+++ b/src/math/p_colmeans.c
@@ -0,0 +1,33 @@
+#include <pal.h>
+
+/**
+ *
+ * Calculates the means vector 'c' of matrix 'a' columns.
+ *
+ * @param a     Pointer to input matrix
+ *
+ * @param c     Pointer to output vector
+ *
+ * @param m     Number of matrix rows
+ *
+ * @param n     Number of matrix columns and size of output vector
+ *
+ * @return      None
+ *
+ */
+
+void p_colmeans_f32(const float *a, float *c, size_t m, size_t n)
+{
+    const float *pa;
+    float *pc;
+    size_t i;
+
+    pa = a;
+    pc = c;
+
+    p_colsums_f32(pa, pc, m, n);
+
+    for (i = 0; i < n; ++i) {
+        *(pc + i) /= n;
+    }
+}
diff --git a/src/math/p_colsums.c b/src/math/p_colsums.c
new file mode 100644
index 0000000..5807143
--- /dev/null
+++ b/src/math/p_colsums.c
@@ -0,0 +1,31 @@
+#include <pal.h>
+
+/**
+ *
+ * Calculates the sums vector 'c' of matrix 'a' columns.
+ *
+ * @param a     Pointer to input matrix
+ *
+ * @param c     Pointer to output vector
+ *
+ * @param m     Number of matrix rows
+ *
+ * @param n     Number of matrix columns and size of output vector
+ *
+ * @return      None
+ *
+ */
+
+void p_colsums_f32(const float *a, float *c, size_t m, size_t n)
+{
+    const float *pa;
+    float *pc;
+    size_t i;
+
+    pa = a;
+    pc = c;
+
+    for (i = 0; i < m * n; ++i) {
+        *(pc + i % n) += *(pa + i);
+    }
+}
diff --git a/tests/math/Makefile.am b/tests/math/Makefile.am
index 652f282..994ee64 100644
--- a/tests/math/Makefile.am
+++ b/tests/math/Makefile.am
@@ -81,6 +81,8 @@ check_PROGRAMS = \
 	check_p_atan_f32 \
 	check_p_atanh_f32 \
 	check_p_cbrt_f32 \
+	check_p_colmeans_f32 \
+	check_p_colsums_f32 \
 	check_p_cos_f32 \
 	check_p_cosh_f32 \
 	check_p_div_f32 \
@@ -126,6 +128,8 @@ check_p_atan2_f32_SOURCES       = $(SIMPLE) p_atan2.c
 check_p_atan_f32_SOURCES        = $(SIMPLE) p_atan.c
 check_p_atanh_f32_SOURCES       = $(SIMPLE) p_atanh.c
 check_p_cbrt_f32_SOURCES        = $(SIMPLE) p_cbrt.c
+check_p_colmeans_f32_SOURCES    = notest.c
+check_p_colsums_f32_SOURCES     = notest.c
 check_p_cos_f32_SOURCES         = $(SIMPLE) p_cos.c
 check_p_cosh_f32_SOURCES        = $(SIMPLE) p_cosh.c
 check_p_div_f32_SOURCES         = $(SIMPLE)
@@ -172,6 +176,8 @@ check_p_atan2_f32_CFLAGS        = -DFUNCTION=p_atan2_f32   -DIS_BINARY
 check_p_atan_f32_CFLAGS         = -DFUNCTION=p_atan_f32    -DIS_UNARY
 check_p_atanh_f32_CFLAGS        = -DFUNCTION=p_atanh_f32   -DIS_UNARY
 check_p_cbrt_f32_CFLAGS         = -DFUNCTION=p_cbrt_f32    -DIS_UNARY
+check_p_colmeans_f32_CFLAGS     = -DFUNCTION=p_colmeans_f32 -DIS_UNARY
+check_p_colsums_f32_CFLAGS      = -DFUNCTION=p_colsums_f32 -DIS_UNARY
 check_p_cos_f32_CFLAGS          = -DFUNCTION=p_cos_f32     -DIS_UNARY
 check_p_cosh_f32_CFLAGS         = -DFUNCTION=p_cosh_f32    -DIS_UNARY
 check_p_div_f32_CFLAGS          = -DFUNCTION=p_div_f32     -DIS_BINARY