Mercurial > octave-nkf

changeset 9513:9f870f73ab7d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
implement built-in diff
author Jaroslav Hajek <highegg@gmail.com>
date Tue, 11 Aug 2009 09:08:12 +0200
parents 56e850e3b06f
children af86991d8d47
files liboctave/CNDArray.cc liboctave/CNDArray.h liboctave/ChangeLog liboctave/dNDArray.cc liboctave/dNDArray.h liboctave/fCNDArray.cc liboctave/fCNDArray.h liboctave/fNDArray.cc liboctave/fNDArray.h liboctave/intNDArray.cc liboctave/intNDArray.h liboctave/mx-inlines.cc scripts/ChangeLog scripts/general/Makefile.in scripts/general/diff.m src/ChangeLog src/data.cc
diffstat 17 files changed, 420 insertions(+), 151 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/CNDArray.cc	Mon Aug 10 15:18:41 2009 -0400
+++ b/liboctave/CNDArray.cc	Tue Aug 11 09:08:12 2009 +0200
@@ -672,6 +672,12 @@
 }
 
 ComplexNDArray
+ComplexNDArray::diff (octave_idx_type order, int dim) const
+{
+  return do_mx_diff_op<ComplexNDArray> (*this, dim, order, mx_inline_diff);
+}
+
+ComplexNDArray
 ComplexNDArray::concat (const ComplexNDArray& rb, const Array<octave_idx_type>& ra_idx)
 {
   if (rb.numel () > 0)
--- a/liboctave/CNDArray.h	Mon Aug 10 15:18:41 2009 -0400
+++ b/liboctave/CNDArray.h	Tue Aug 11 09:08:12 2009 +0200
@@ -93,6 +93,8 @@
   ComplexNDArray cummin (int dim = 0) const;
   ComplexNDArray cummin (ArrayN<octave_idx_type>& index, int dim = 0) const;
 
+  ComplexNDArray diff (octave_idx_type order = 1, int dim = 0) const;
+
   ComplexNDArray& insert (const NDArray& a, octave_idx_type r, octave_idx_type c);
   ComplexNDArray& insert (const ComplexNDArray& a, octave_idx_type r, octave_idx_type c);
   ComplexNDArray& insert (const ComplexNDArray& a, const Array<octave_idx_type>& ra_idx);
--- a/liboctave/ChangeLog	Mon Aug 10 15:18:41 2009 -0400
+++ b/liboctave/ChangeLog	Tue Aug 11 09:08:12 2009 +0200
@@ -1,3 +1,20 @@
+2009-08-11  Jaroslav Hajek  <highegg@gmail.com>
+
+	* mx-inlines.cc (mx_inline_diff<T>): New overloaded template
+	function.
+	(get_extent_triplet): Use dim_vector::first_non_singleton.
+	(do_mx_diff_op): New template function.
+	* dNDArray.cc (NDArray::diff): New method.
+	* dNDArray.h: Declare it.
+	* fNDArray.cc (FloatNDArray::diff): New method.
+	* fNDArray.h: Declare it.
+	* CNDArray.cc (ComplexNDArray::diff): New method.
+	* CNDArray.h: Declare it.
+	* fCNDArray.cc (FloatComplexNDArray::diff): New method.
+	* fCNDArray.h: Declare it.
+	* intNDArray.cc (intNDArray<T>::diff): New method.
+	* intNDArray.h: Declare it.
+
 2009-08-10  Jaroslav Hajek  <highegg@gmail.com>
 
 	* dim-vector.h (dim_vector::first_non_singleton): New method.
--- a/liboctave/dNDArray.cc	Mon Aug 10 15:18:41 2009 -0400
+++ b/liboctave/dNDArray.cc	Tue Aug 11 09:08:12 2009 +0200
@@ -787,6 +787,12 @@
 }
 
 NDArray
+NDArray::diff (octave_idx_type order, int dim) const
+{
+  return do_mx_diff_op<NDArray> (*this, dim, order, mx_inline_diff);
+}
+
+NDArray
 NDArray::concat (const NDArray& rb, const Array<octave_idx_type>& ra_idx)
 {
   if (rb.numel () > 0)
--- a/liboctave/dNDArray.h	Mon Aug 10 15:18:41 2009 -0400
+++ b/liboctave/dNDArray.h	Tue Aug 11 09:08:12 2009 +0200
@@ -105,6 +105,8 @@
   NDArray cummin (int dim = 0) const;
   NDArray cummin (ArrayN<octave_idx_type>& index, int dim = 0) const;
 
+  NDArray diff (octave_idx_type order = 1, int dim = 0) const;
+
   NDArray& insert (const NDArray& a, octave_idx_type r, octave_idx_type c);
   NDArray& insert (const NDArray& a, const Array<octave_idx_type>& ra_idx);
 
--- a/liboctave/fCNDArray.cc	Mon Aug 10 15:18:41 2009 -0400
+++ b/liboctave/fCNDArray.cc	Tue Aug 11 09:08:12 2009 +0200
@@ -667,6 +667,12 @@
 }
 
 FloatComplexNDArray
+FloatComplexNDArray::diff (octave_idx_type order, int dim) const
+{
+  return do_mx_diff_op<FloatComplexNDArray> (*this, dim, order, mx_inline_diff);
+}
+
+FloatComplexNDArray
 FloatComplexNDArray::concat (const FloatComplexNDArray& rb, const Array<octave_idx_type>& ra_idx)
 {
   if (rb.numel () > 0)
--- a/liboctave/fCNDArray.h	Mon Aug 10 15:18:41 2009 -0400
+++ b/liboctave/fCNDArray.h	Tue Aug 11 09:08:12 2009 +0200
@@ -93,6 +93,8 @@
   FloatComplexNDArray cummin (int dim = 0) const;
   FloatComplexNDArray cummin (ArrayN<octave_idx_type>& index, int dim = 0) const;
 
+  FloatComplexNDArray diff (octave_idx_type order = 1, int dim = 0) const;
+
   FloatComplexNDArray& insert (const NDArray& a, octave_idx_type r, octave_idx_type c);
   FloatComplexNDArray& insert (const FloatComplexNDArray& a, octave_idx_type r, octave_idx_type c);
   FloatComplexNDArray& insert (const FloatComplexNDArray& a, const Array<octave_idx_type>& ra_idx);
--- a/liboctave/fNDArray.cc	Mon Aug 10 15:18:41 2009 -0400
+++ b/liboctave/fNDArray.cc	Tue Aug 11 09:08:12 2009 +0200
@@ -742,6 +742,12 @@
 }
 
 FloatNDArray
+FloatNDArray::diff (octave_idx_type order, int dim) const
+{
+  return do_mx_diff_op<FloatNDArray> (*this, dim, order, mx_inline_diff);
+}
+
+FloatNDArray
 FloatNDArray::concat (const FloatNDArray& rb, const Array<octave_idx_type>& ra_idx)
 {
   if (rb.numel () > 0)
--- a/liboctave/fNDArray.h	Mon Aug 10 15:18:41 2009 -0400
+++ b/liboctave/fNDArray.h	Tue Aug 11 09:08:12 2009 +0200
@@ -102,6 +102,8 @@
   FloatNDArray cummin (int dim = 0) const;
   FloatNDArray cummin (ArrayN<octave_idx_type>& index, int dim = 0) const;
 
+  FloatNDArray diff (octave_idx_type order = 1, int dim = 0) const;
+
   FloatNDArray& insert (const FloatNDArray& a, octave_idx_type r, octave_idx_type c);
   FloatNDArray& insert (const FloatNDArray& a, const Array<octave_idx_type>& ra_idx);
 
--- a/liboctave/intNDArray.cc	Mon Aug 10 15:18:41 2009 -0400
+++ b/liboctave/intNDArray.cc	Tue Aug 11 09:08:12 2009 +0200
@@ -270,6 +270,13 @@
   return do_mx_cumminmax_op<intNDArray<T> > (*this, idx_arg, dim, mx_inline_cummin);
 }
 
+template <class T>
+intNDArray<T>
+intNDArray<T>::diff (octave_idx_type order, int dim) const
+{
+  return do_mx_diff_op<intNDArray<T> > (*this, dim, order, mx_inline_diff);
+}
+
 /*
 ;;; Local Variables: ***
 ;;; mode: C++ ***
--- a/liboctave/intNDArray.h	Mon Aug 10 15:18:41 2009 -0400
+++ b/liboctave/intNDArray.h	Tue Aug 11 09:08:12 2009 +0200
@@ -86,6 +86,8 @@
   intNDArray sum (int dim) const;
   intNDArray cumsum (int dim) const;
 
+  intNDArray diff (octave_idx_type order = 1, int dim = 0) const;
+
   intNDArray abs (void) const;
   intNDArray signum (void) const;
 
--- a/liboctave/mx-inlines.cc	Mon Aug 10 15:18:41 2009 -0400
+++ b/liboctave/mx-inlines.cc	Tue Aug 11 09:08:12 2009 +0200
@@ -844,6 +844,114 @@
 OP_CUMMINMAX_FCNN (mx_inline_cummin)
 OP_CUMMINMAX_FCNN (mx_inline_cummax)
 
+template <class T>
+void mx_inline_diff (const T *v, T *r, octave_idx_type n,
+                     octave_idx_type order)
+{
+  switch (order)
+    {
+    case 1:
+      for (octave_idx_type i = 0; i < n-1; i++)
+        r[i] = v[i+1] - v[i];
+      break;
+    case 2:
+        {
+          T lst;
+          if (n > 1)
+            lst = v[1] - v[0];
+          for (octave_idx_type i = 0; i < n-2; i++)
+            {
+              T dif = v[i+2] - v[i+1];
+              r[i] = dif - lst;
+              lst = dif;
+            }
+        }
+      break;
+    default:
+        {
+          OCTAVE_LOCAL_BUFFER (T, buf, n-1);
+
+          for (octave_idx_type i = 0; i < n-1; i++)
+            buf[i] = v[i+1] - v[i];
+
+          for (octave_idx_type o = 2; o <= order; o++)
+            {
+              for (octave_idx_type i = 0; i < n-o; i++)
+                buf[i] = buf[i+1] - buf[i];
+            }
+
+          for (octave_idx_type i = 0; i < n-order; i++)
+            r[i] = buf[i];
+        }
+    }
+}
+
+template <class T>
+void mx_inline_diff (const T *v, T *r, 
+                     octave_idx_type m, octave_idx_type n,
+                     octave_idx_type order)
+{
+  switch (order)
+    {
+    case 1:
+      for (octave_idx_type i = 0; i < m*(n-1); i++)
+        r[i] = v[i+m] - v[i];
+      break;
+    case 2:
+      for (octave_idx_type i = 0; i < n-2; i++)
+        {
+          for (octave_idx_type j = i*m; j < i*m+m; j++)
+            r[j] = (v[j+m+m] - v[j+m]) + (v[j+m] - v[j]);
+        }
+      break;
+    default:
+        {
+          OCTAVE_LOCAL_BUFFER (T, buf, n-1);
+
+          for (octave_idx_type j = 0; j < m; j++)
+            {
+              for (octave_idx_type i = 0; i < n-1; i++)
+                buf[i] = v[i*m+j+m] - v[i*m+j];
+
+              for (octave_idx_type o = 2; o <= order; o++)
+                {
+                  for (octave_idx_type i = 0; i < n-o; i++)
+                    buf[i] = buf[i+1] - buf[i];
+                }
+
+              for (octave_idx_type i = 0; i < n-order; i++)
+                r[i*m+j] = buf[i];
+            }
+        }
+    }
+}
+
+template <class T>
+inline void
+mx_inline_diff (const T *v, T *r,
+                octave_idx_type l, octave_idx_type n, octave_idx_type u,
+                octave_idx_type order)
+{
+  if (! n) return;
+  if (l == 1)
+    {
+      for (octave_idx_type i = 0; i < u; i++)
+        {
+          mx_inline_diff (v, r, n, order);
+          v += n; r += n-order;
+        }
+    }
+  else
+    {
+      for (octave_idx_type i = 0; i < u; i++)
+        {
+          mx_inline_diff (v, r, l, n, order);
+          v += l*n;
+          r += l*(n-order);
+        }
+    }
+}
+
 // Assistant function
 
 inline void
@@ -861,10 +969,8 @@
   else
     {
       if (dim < 0)
-        {
-          // find first non-singleton dim
-          for (dim = 0; dims(dim) == 1 && dim < ndims - 1; dim++) ;
-        }
+        dim = dims.first_non_singleton ();
+
       // calculate extent triplet.
       l = 1, n = dims(dim), u = 1;
       for (octave_idx_type i = 0; i < dim; i++) 
@@ -1003,6 +1109,40 @@
   return ret;
 }
 
+template <class ArrayType>
+inline ArrayType
+do_mx_diff_op (const ArrayType& src, int dim, octave_idx_type order,
+               void (*mx_diff_op) (const typename ArrayType::element_type *, 
+                                   typename ArrayType::element_type *,
+                                   octave_idx_type, octave_idx_type, octave_idx_type,
+                                   octave_idx_type))
+{
+  octave_idx_type l, n, u;
+  if (order <= 0)
+    return src;
+
+  dim_vector dims = src.dims ();
+
+  get_extent_triplet (dims, dim, l, n, u);
+  if (dim >= dims.length ())
+    dims.resize (dim+1, 1);
+
+  if (dims(dim) <= order)
+    {
+      dims (dim) = 0;
+      return ArrayType (dims);
+    }
+  else
+    {
+      dims(dim) -= order;
+    }
+
+  ArrayType ret (dims);
+  mx_diff_op (src.data (), ret.fortran_vec (), l, n, u, order);
+
+  return ret;
+}
+
 #endif
 
 /*
--- a/scripts/ChangeLog	Mon Aug 10 15:18:41 2009 -0400
+++ b/scripts/ChangeLog	Tue Aug 11 09:08:12 2009 +0200
@@ -1,3 +1,8 @@
+2009-08-11  Jaroslav Hajek  <highegg@gmail.com>
+
+	* general/diff.m: Remove.
+	* general/Makefile.in: Update.
+
 2009-08-07  Jaroslav Hajek  <highegg@gmail.com>
 
 	* general/flipdim.m: Fix omitted check.
--- a/scripts/general/Makefile.in	Mon Aug 10 15:18:41 2009 -0400
+++ b/scripts/general/Makefile.in	Tue Aug 11 09:08:12 2009 +0200
@@ -36,7 +36,7 @@
 SOURCES = __isequal__.m __splinen__.m accumarray.m arrayfun.m \
   bicubic.m bitcmp.m bitget.m bitset.m blkdiag.m cart2pol.m \
   cart2sph.m cellidx.m cell2mat.m celldisp.m circshift.m colon.m common_size.m \
-  cplxpair.m cumtrapz.m dblquad.m deal.m del2.m diff.m display.m flipdim.m \
+  cplxpair.m cumtrapz.m dblquad.m deal.m del2.m display.m flipdim.m \
   fliplr.m flipud.m genvarname.m gradient.m idivide.m int2str.m \
   interp1.m interp1q.m interp2.m interp3.m interpn.m interpft.m \
   is_duplicate_entry.m isa.m isdefinite.m isdir.m isequal.m \
--- a/scripts/general/diff.m	Mon Aug 10 15:18:41 2009 -0400
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,146 +0,0 @@
-## Copyright (C) 1995, 1996, 1999, 2000, 2002, 2004, 2005, 2006, 2007,
-##               2008, 2009 Kurt Hornik
-##
-## This file is part of Octave.
-##
-## Octave is free software; you can redistribute it and/or modify it
-## under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 3 of the License, or (at
-## your option) any later version.
-##
-## Octave is distributed in the hope that it will be useful, but
-## WITHOUT ANY WARRANTY; without even the implied warranty of
-## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-## General Public License for more details.
-##
-## You should have received a copy of the GNU General Public License
-## along with Octave; see the file COPYING.  If not, see
-## <http://www.gnu.org/licenses/>.
-
-## -*- texinfo -*-
-## @deftypefn {Function File} {} diff (@var{x}, @var{k}, @var{dim})
-## If @var{x} is a vector of length @var{n}, @code{diff (@var{x})} is the
-## vector of first differences
-## @tex
-##  $x_2 - x_1, \ldots{}, x_n - x_{n-1}$.
-## @end tex
-## @ifnottex
-##  @var{x}(2) - @var{x}(1), @dots{}, @var{x}(n) - @var{x}(n-1).
-## @end ifnottex
-##
-## If @var{x} is a matrix, @code{diff (@var{x})} is the matrix of column
-## differences along the first non-singleton dimension.
-##
-## The second argument is optional.  If supplied, @code{diff (@var{x},
-## @var{k})}, where @var{k} is a non-negative integer, returns the
-## @var{k}-th differences.  It is possible that @var{k} is larger than
-## then first non-singleton dimension of the matrix.  In this case,
-## @code{diff} continues to take the differences along the next
-## non-singleton dimension.
-##
-## The dimension along which to take the difference can be explicitly
-## stated with the optional variable @var{dim}.  In this case the 
-## @var{k}-th order differences are calculated along this dimension.
-## In the case where @var{k} exceeds @code{size (@var{x}, @var{dim})}
-## then an empty matrix is returned.
-## @end deftypefn
-
-## Author: KH <Kurt.Hornik@wu-wien.ac.at>
-## Created: 2 February 1995
-## Adapted-By: jwe
-
-function x = diff (x, k, dim)
-
-  if (nargin < 1 || nargin > 3)
-    print_usage ();
-  endif
-
-  if (nargin < 2 || isempty(k))
-    k = 1;
-  else
-    if (! (isscalar (k) && k == round (k) && k >= 0))
-      error ("diff: k must be a nonnegative integer");
-    elseif (k == 0)
-      return;
-    endif
-  endif
-
-  nd = ndims (x);
-  sz = size (x);
-  if (nargin != 3)
-    %% Find the first non-singleton dimension
-    dim  = 1;
-    while (dim < nd + 1 && sz (dim) == 1)
-      dim = dim + 1;
-    endwhile
-    if (dim > nd)
-      dim = 1;
-    endif
-  else
-    if (! (isscalar (dim) && dim == round (dim)) && dim > 0 && 
-	dim < (nd + 1))
-      error ("diff: dim must be an integer and valid dimension");
-    endif
-  endif
-
-  if (ischar (x))
-    error ("diff: symbolic differentiation not (yet) supported");
-  endif
-
-
-  if (nargin == 3)
-    if (sz (dim) <= k)
-      sz(dim) = 0;
-      x = zeros (sz);
-    else
-      n = sz (dim);
-      idx1 = cell ();
-      for i = 1:nd
-	idx1{i} = 1:sz(i);
-      endfor
-      idx2 = idx1;
-      for i = 1 : k;
-	idx1{dim} = 2 : (n - i + 1);	
-	idx2{dim} = 1 : (n - i);	
-	x = x(idx1{:}) - x(idx2{:});
-      endfor
-    endif
-  else
-    if (sum (sz - 1) < k)
-      x = [];
-    else
-      idx1 = cell ();
-      for i = 1:nd
-	idx1{i} = 1:sz(i);
-      endfor
-      idx2 = idx1;
-      while (k)
-	n = sz (dim);
-	for i = 1 : min (k, n - 1)
-	  idx1{dim} = 2 : (n - i + 1);	
-	  idx2{dim} = 1 : (n - i);	
-	  x = x(idx1{:}) - x(idx2{:});
-	endfor
-	idx1{dim} = idx2{dim} = 1;
-	k = k - min (k, n - 1);
-	dim = dim + 1;
-      endwhile
-    endif
-  endif
-
-endfunction
-
-%!assert((diff ([1, 2, 3, 4]) == [1, 1, 1]
-%! && diff ([1, 3, 7, 19], 2) == [2, 8]
-%! && diff ([1, 2; 5, 4; 8, 7; 9, 6; 3, 1]) == [4, 2; 3, 3; 1, -1; -6, -5]
-%! && diff ([1, 2; 5, 4; 8, 7; 9, 6; 3, 1], 3) == [-1, -5; -5, 0]
-%! && isempty (diff (1))));
-
-%!error diff ([1, 2; 3, 4], -1);
-
-%!error diff ("foo");
-
-%!error diff ();
-
-%!error diff (1, 2, 3, 4);
-
--- a/src/ChangeLog	Mon Aug 10 15:18:41 2009 -0400
+++ b/src/ChangeLog	Tue Aug 11 09:08:12 2009 +0200
@@ -1,3 +1,8 @@
+2009-08-11  Jaroslav Hajek  <highegg@gmail.com>
+
+	* data.cc (Fdiff): New built-in function.
+	(do_diff): New assistant function.
+
 2009-08-10  John W. Eaton  <jwe@octave.org>
 
 	* DLD-FUNCTIONS/dlmread.cc (Fdlmread): Perform tilde expansion on
--- a/src/data.cc	Mon Aug 10 15:18:41 2009 -0400
+++ b/src/data.cc	Tue Aug 11 09:08:12 2009 +0200
@@ -3,6 +3,7 @@
 Copyright (C) 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002,
               2003, 2004, 2005, 2006, 2007, 2008, 2009 John W. Eaton
 Copyright (C) 2009 Jaroslav Hajek
+Copyright (C) 2009 VZLU Prague
 
 This file is part of Octave.
 
@@ -6159,6 +6160,212 @@
 
 #undef MAKE_INT_BRANCH
 
+template <class SparseT>
+static SparseT
+do_sparse_diff (const SparseT& array, octave_idx_type order,
+                int dim)
+{
+  SparseT retval = array;
+  if (dim == 1)
+    {
+      octave_idx_type k = retval.columns ();
+      while (order > 0 && k > 0)
+        {
+          idx_vector col1 (':'), col2 (':'), sl1 (1, k), sl2 (0, k-1);
+          retval = retval.index (col1, sl1, 0) - retval.index (col2, sl2, 0);
+          assert (retval.columns () == k-1);
+          order--;
+          k--;
+        }
+    }
+  else
+    {
+      octave_idx_type k = retval.rows ();
+      while (order > 0 && k > 0)
+        {
+          idx_vector col1 (':'), col2 (':'), sl1 (1, k), sl2 (0, k-1);
+          retval = retval.index (sl1, col1, 0) - retval.index (sl2, col2, 0);
+          assert (retval.rows () == k-1);
+          order--;
+          k--;
+        }
+    }
+
+  return retval;
+}
+
+static octave_value
+do_diff (const octave_value& array, octave_idx_type order,
+         int dim = -1)
+{
+  octave_value retval;
+
+  const dim_vector& dv = array.dims ();
+  if (dim == -1)
+    {
+      dim = array.dims ().first_non_singleton ();
+
+      // Bother Matlab. This behavior is really wicked.
+      if (dv(dim) <= order)
+        {
+          if (dv(dim) == 1)
+            retval = array.resize (dim_vector (0, 0));
+          else
+            {
+              retval = array;
+              while (order > 0)
+                {
+                  if (dim == dv.length ())
+                    {
+                      retval = do_diff (array, order, dim - 1);
+                      order = 0;
+                    }
+                  else if (dv(dim) == 1)
+                    dim++;
+                  else
+                    {
+                      retval = do_diff (array, dv(dim) - 1, dim);
+                      order -= dv(dim) - 1;
+                      dim++;
+                    }
+                }
+            }
+
+          return retval;
+        }
+    }
+
+  if (array.is_integer_type ())
+    {
+      if (array.is_int8_type ())
+        retval = array.int8_array_value ().diff (order, dim);
+      else if (array.is_int16_type ())
+        retval = array.int16_array_value ().diff (order, dim);
+      else if (array.is_int32_type ())
+        retval = array.int32_array_value ().diff (order, dim);
+      else if (array.is_int64_type ())
+        retval = array.int64_array_value ().diff (order, dim);
+      else if (array.is_uint8_type ())
+        retval = array.uint8_array_value ().diff (order, dim);
+      else if (array.is_uint16_type ())
+        retval = array.uint16_array_value ().diff (order, dim);
+      else if (array.is_uint32_type ())
+        retval = array.uint32_array_value ().diff (order, dim);
+      else if (array.is_uint64_type ())
+        retval = array.uint64_array_value ().diff (order, dim);
+      else
+        panic_impossible ();
+    }
+  else if (array.is_sparse_type ())
+    {
+      if (array.is_complex_type ())
+        retval = do_sparse_diff (array.sparse_complex_matrix_value (), order, dim);
+      else
+        retval = do_sparse_diff (array.sparse_matrix_value (), order, dim);
+    }
+  else if (array.is_single_type ())
+    {
+      if (array.is_complex_type ())
+        retval = array.float_complex_array_value ().diff (order, dim);
+      else
+        retval = array.float_array_value ().diff (order, dim);
+    }
+  else
+    {
+      if (array.is_complex_type ())
+        retval = array.complex_array_value ().diff (order, dim);
+      else
+        retval = array.array_value ().diff (order, dim);
+    }
+
+  return retval;
+}
+
+DEFUN (diff, args, ,
+  "-*- texinfo -*-\n\
+@deftypefn {Built-in Function} {} diff (@var{x}, @var{k}, @var{dim})\n\
+If @var{x} is a vector of length @var{n}, @code{diff (@var{x})} is the\n\
+vector of first differences\n\
+@tex\n\
+ $x_2 - x_1, \\ldots{}, x_n - x_{n-1}$.\n\
+@end tex\n\
+@ifnottex\n\
+ @var{x}(2) - @var{x}(1), @dots{}, @var{x}(n) - @var{x}(n-1).\n\
+@end ifnottex\n\
+\n\
+If @var{x} is a matrix, @code{diff (@var{x})} is the matrix of column\n\
+differences along the first non-singleton dimension.\n\
+\n\
+The second argument is optional.  If supplied, @code{diff (@var{x},\n\
+@var{k})}, where @var{k} is a non-negative integer, returns the\n\
+@var{k}-th differences.  It is possible that @var{k} is larger than\n\
+then first non-singleton dimension of the matrix.  In this case,\n\
+@code{diff} continues to take the differences along the next\n\
+non-singleton dimension.\n\
+\n\
+The dimension along which to take the difference can be explicitly\n\
+stated with the optional variable @var{dim}.  In this case the \n\
+@var{k}-th order differences are calculated along this dimension.\n\
+In the case where @var{k} exceeds @code{size (@var{x}, @var{dim})}\n\
+then an empty matrix is returned.\n\
+@end deftypefn")
+{
+  int nargin = args.length ();
+  octave_value retval;
+
+  if (nargin < 1 || nargin > 3)
+    print_usage ();
+  else if (! args(0).is_numeric_type ())
+    error ("diff: X must be numeric");
+
+  if (! error_state)
+    {
+      int dim = -1;
+      octave_idx_type order = 1;
+      if (nargin > 1)
+        {
+          if (args(1).is_scalar_type ())
+            order = args(1).idx_type_value (true, false);
+          else if (! args(1).is_zero_by_zero ())
+            error ("order must be a scalar or []");
+          if (! error_state && order < 0)
+            error ("order must be non-negative");
+        }
+
+      if (nargin > 2)
+        {
+          dim = args(2).int_value (true, false);
+          if (! error_state && (dim < 1 || dim > args(0).ndims ()))
+            error ("needs a valid dimension");
+          else
+            dim -= 1;
+        }
+
+      if (! error_state)
+        retval = do_diff (args(0), order, dim);
+    }
+
+  return retval;
+}
+
+/*
+
+%!assert (diff ([1, 2, 3, 4]), [1, 1, 1])
+%!assert (diff ([1, 3, 7, 19], 2), [2, 8])
+%!assert (diff ([1, 2; 5, 4; 8, 7; 9, 6; 3, 1]), [4, 2; 3, 3; 1, -1; -6, -5])
+%!assert (diff ([1, 2; 5, 4; 8, 7; 9, 6; 3, 1], 3), [-1, -5; -5, 0])
+%!assert (isempty (diff (1)));
+
+%!error diff ([1, 2; 3, 4], -1);
+
+%!error diff ("foo");
+
+%!error diff ();
+
+%!error diff (1, 2, 3, 4);
+
+*/
+
 /*
 ;;; Local Variables: ***
 ;;; mode: C++ ***