--- a/liboctave/dMatrix.cc
+++ b/liboctave/dMatrix.cc
@@ -1991,203 +1991,42 @@
 Matrix
 Matrix::cumprod (int dim) const
 {
-  int nr = rows ();
-  int nc = cols ();
-  Matrix retval (nr, nc);
-
-  if (nr > 0 && nc >0)
-    {
-      if ((nr == 1 && dim == 0) || dim == 1)
-	{
-	  for (int i = 0; i < nr; i++)
-	    {
-	      double prod = elem (i, 0);
-	      for (int j = 0; j < nc; j++)
-		{
-		  retval.elem (i, j) = prod;
-		  if (j < nc - 1)
-		    prod *= elem (i, j+1);
-		}
-	    }
-	}
-      else
-	{
-	  for (int j = 0; j < nc; j++)
-	    {
-	      double prod = elem (0, j);
-	      for (int i = 0; i < nr; i++)
-		{
-		  retval.elem (i, j) = prod;
-		  if (i < nr - 1)
-		    prod *= elem (i+1, j);
-		}
-	    }
-	}
-    }
-  return retval;
+  MX_CUMMULATIVE_OP (Matrix, double, *=);
 }
 
 Matrix
 Matrix::cumsum (int dim) const
 {
-  int nr = rows ();
-  int nc = cols ();
-  Matrix retval (nr, nc);
-
-  if (nr > 0 && nc > 0)
-    {
-      if ((nr == 1 && dim == 0) || dim == 1)
-	{
-	  for (int i = 0; i < nr; i++)
-	    {
-	      double sum = elem (i, 0);
-	      for (int j = 0; j < nc; j++)
-		{
-		  retval.elem (i, j) = sum;
-		  if (j < nc - 1)
-		    sum += elem (i, j+1);
-		}
-	    }
-	}
-      else
-	{
-	  for (int j = 0; j < nc; j++)
-	    {
-	      double sum = elem (0, j);
-	      for (int i = 0; i < nr; i++)
-		{
-		  retval.elem (i, j) = sum;
-		  if (i < nr - 1)
-		    sum += elem (i+1, j);
-		}
-	    }
-	}
-    }
-  return retval;
+  MX_CUMMULATIVE_OP (Matrix, double, +=);
 }
 
 Matrix
 Matrix::prod (int dim) const
 {
-  Matrix retval;
-
-  int nr = rows ();
-  int nc = cols ();
-
-  if (nr > 0 && nc > 0)
-    {
-      if ((nr == 1 && dim == 0) || dim == 1)
-	{
-	  retval.resize(nr, 1);
-	  for (int i = 0; i < nr; i++)
-	    {
-	      retval.elem (i, 0) = 1.0;
-	      for (int j = 0; j < nc; j++)
-		retval.elem (i, 0) *= elem (i, j);
-	    }
-	}
-      else
-	{
-	  retval.resize (1, nc);
-	  for (int j = 0; j < nc; j++)
-	    {
-	      retval.elem (0, j) = 1.0;
-	      for (int i = 0; i < nr; i++)
-		retval.elem (0, j) *= elem (i, j);
-	    }
-	}
-    }
-  else
-    {
-      retval.resize (1, 1);
-      retval.elem (0, 0) = 1.0;
-    }
-
-  return retval;
+  MX_REDUCTION_OP (Matrix, *=, 1.0, 1.0);
 }
 
 Matrix
 Matrix::sum (int dim) const
 {
-  Matrix retval;
-
-  int nr = rows ();
-  int nc = cols ();
-
-  if (nr > 0 && nc > 0)
-    {
-      if ((nr == 1 && dim == 0) || dim == 1)
-	{
-	  retval.resize (nr, 1);
-	  for (int i = 0; i < nr; i++)
-	    {
-	      retval.elem (i, 0) = 0.0;
-	      for (int j = 0; j < nc; j++)
-		retval.elem (i, 0) += elem (i, j);
-	    }
-	}
-      else
-	{
-	  retval.resize (1, nc);
-	  for (int j = 0; j < nc; j++)
-	    {
-	      retval.elem (0, j) = 0.0;
-	      for (int i = 0; i < nr; i++)
-		retval.elem (0, j) += elem (i, j);
-	    }
-	}
-    }
-  else
-    {
-      retval.resize (1, 1);
-      retval.elem (0, 0) = 0.0;
-    }
-  return retval;
+  MX_REDUCTION_OP (Matrix, +=, 0.0, 0.0);
 }
 
 Matrix
 Matrix::sumsq (int dim) const
 {
-  Matrix retval;
-
-  int nr = rows ();
-  int nc = cols ();
-
-  if (nr > 0 && nc > 0)
-    {
-      if ((nr == 1 && dim == 0) || dim == 1)
-	{
-	  retval.resize (nr, 1);
-	  for (int i = 0; i < nr; i++)
-	    {
-	      retval.elem (i, 0) = 0.0;
-	      for (int j = 0; j < nc; j++)
-		{
-		  double d = elem (i, j);
-		  retval.elem (i, 0) += d * d;
-		}
-	    }
-	}
-      else
-	{
-	  retval.resize (1, nc);
-	  for (int j = 0; j < nc; j++)
-	    {
-	      retval.elem (0, j) = 0.0;
-	      for (int i = 0; i < nr; i++)
-		{
-		  double d = elem (i, j);
-		  retval.elem (0, j) += d * d;
-		}
-	    }
-	}
-    }
-  else
-    {
-      retval.resize (1, 1);
-      retval.elem (0, 0) = 0.0;
-    }
-  return retval;
+#define ROW_EXPR \
+  double d = elem (i, j); \
+  retval.elem (i, 0) += d * d
+
+#define COL_EXPR \
+  double d = elem (i, j); \
+  retval.elem (0, j) += d * d
+
+  MX_BASE_REDUCTION_OP (Matrix, ROW_EXPR, COL_EXPR, 0.0, 0.0);
+
+#undef ROW_EXPR
+#undef COL_EXPR
 }
 
 Matrix