src/data.cc

/*

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.

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/>.

*/

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include "systime.h"

#ifdef HAVE_SYS_TYPES_H
#include <sys/types.h>
#endif

#ifdef HAVE_SYS_RESOURCE_H
#include <sys/resource.h>
#endif

#include <cfloat>

#include <string>

#include "lo-ieee.h"
#include "lo-math.h"
#include "str-vec.h"
#include "quit.h"
#include "mx-base.h"

#include "Cell.h"
#include "defun.h"
#include "error.h"
#include "gripes.h"
#include "oct-map.h"
#include "oct-obj.h"
#include "ov.h"
#include "ov-float.h"
#include "ov-complex.h"
#include "ov-flt-complex.h"
#include "ov-cx-mat.h"
#include "ov-flt-cx-mat.h"
#include "ov-cx-sparse.h"
#include "parse.h"
#include "pt-mat.h"
#include "utils.h"
#include "variables.h"
#include "pager.h"
#include "xnorm.h"

#if ! defined (HAVE_HYPOTF) && defined (HAVE__HYPOTF)
#define hypotf _hypotf
#define HAVE_HYPOTF 1
#endif

#define ANY_ALL(FCN) \
 \
  octave_value retval; \
 \
  int nargin = args.length (); \
 \
  if (nargin == 1 || nargin == 2) \
    { \
      int dim = (nargin == 1 ? -1 : args(1).int_value (true) - 1); \
 \
      if (! error_state) \
        { \
          if (dim >= -1) \
	    retval = args(0).FCN (dim); \
          else \
	    error (#FCN ": invalid dimension argument = %d", dim + 1); \
        } \
      else \
        error (#FCN ": expecting dimension argument to be an integer"); \
    } \
  else \
    print_usage (); \
 \
  return retval

DEFUN (all, args, ,
  "-*- texinfo -*-\n\
@deftypefn {Built-in Function} {} all (@var{x}, @var{dim})\n\
The function @code{all} behaves like the function @code{any}, except\n\
that it returns true only if all the elements of a vector, or all the\n\
elements along dimension @var{dim} of a matrix, are nonzero.\n\
@end deftypefn")
{
  ANY_ALL (all);
}

/*

%!test
%! x = ones (3);
%! x(1,1) = 0;
%! assert((all (all (rand (3) + 1) == [1, 1, 1]) == 1
%! && all (all (x) == [0, 1, 1]) == 1
%! && all (x, 1) == [0, 1, 1]
%! && all (x, 2) == [0; 1; 1]));

%!test
%! x = ones (3, 'single');
%! x(1,1) = 0;
%! assert((all (all (single (rand (3) + 1)) == [1, 1, 1]) == 1
%! && all (all (x) == [0, 1, 1]) == 1
%! && all (x, 1) == [0, 1, 1]
%! && all (x, 2) == [0; 1; 1]));

%!error <Invalid call to all.*> all ();
%!error <Invalid call to all.*> all (1, 2, 3);

 */

DEFUN (any, args, ,
  "-*- texinfo -*-\n\
@deftypefn {Built-in Function} {} any (@var{x}, @var{dim})\n\
For a vector argument, return 1 if any element of the vector is\n\
nonzero.\n\
\n\
For a matrix argument, return a row vector of ones and\n\
zeros with each element indicating whether any of the elements of the\n\
corresponding column of the matrix are nonzero.  For example,\n\
\n\
@example\n\
@group\n\
any (eye (2, 4))\n\
     @result{} [ 1, 1, 0, 0 ]\n\
@end group\n\
@end example\n\
\n\
If the optional argument @var{dim} is supplied, work along dimension\n\
@var{dim}.  For example,\n\
\n\
@example\n\
@group\n\
any (eye (2, 4), 2)\n\
     @result{} [ 1; 1 ]\n\
@end group\n\
@end example\n\
@end deftypefn")
{
  ANY_ALL (any);
}

/*

%!test
%! x = zeros (3);
%! x(3,3) = 1;
%! assert((all (any (x) == [0, 0, 1]) == 1
%! && all (any (ones (3)) == [1, 1, 1]) == 1
%! && any (x, 1) == [0, 0, 1]
%! && any (x, 2) == [0; 0; 1]));

%!test
%! x = zeros (3,'single');
%! x(3,3) = 1;
%! assert((all (any (x) == [0, 0, 1]) == 1
%! && all (any (ones (3, 'single')) == [1, 1, 1]) == 1
%! && any (x, 1) == [0, 0, 1]
%! && any (x, 2) == [0; 0; 1]));

%!error <Invalid call to any.*> any ();
%!error <Invalid call to any.*> any (1, 2, 3);

 */

// These mapping functions may also be useful in other places, eh?

typedef double (*d_dd_fcn) (double, double);
typedef float (*f_ff_fcn) (float, float);

static NDArray
map_d_m (d_dd_fcn f, double x, const NDArray& y)
{
  NDArray retval (y.dims ());
  double *r_data = retval.fortran_vec ();

  const double *y_data = y.data ();

  octave_idx_type nel = y.numel ();

  for (octave_idx_type i = 0; i < nel; i++)
    {
      OCTAVE_QUIT;
      r_data[i] = f (x, y_data[i]);
    }

  return retval;
}

static FloatNDArray
map_f_fm (f_ff_fcn f, float x, const FloatNDArray& y)
{
  FloatNDArray retval (y.dims ());
  float *r_data = retval.fortran_vec ();

  const float *y_data = y.data ();

  octave_idx_type nel = y.numel ();

  for (octave_idx_type i = 0; i < nel; i++)
    {
      OCTAVE_QUIT;
      r_data[i] = f (x, y_data[i]);
    }

  return retval;
}

static NDArray
map_m_d (d_dd_fcn f, const NDArray& x, double y)
{
  NDArray retval (x.dims ());
  double *r_data = retval.fortran_vec ();

  const double *x_data = x.data ();

  octave_idx_type nel = x.numel ();

  for (octave_idx_type i = 0; i < nel; i++)
    {
      OCTAVE_QUIT;
      r_data[i] = f (x_data[i], y);
    }

  return retval;
}

static FloatNDArray
map_fm_f (f_ff_fcn f, const FloatNDArray& x, float y)
{
  FloatNDArray retval (x.dims ());
  float *r_data = retval.fortran_vec ();

  const float *x_data = x.data ();

  octave_idx_type nel = x.numel ();

  for (octave_idx_type i = 0; i < nel; i++)
    {
      OCTAVE_QUIT;
      r_data[i] = f (x_data[i], y);
    }

  return retval;
}

static NDArray
map_m_m (d_dd_fcn f, const NDArray& x, const NDArray& y)
{
  assert (x.dims () == y.dims ());

  NDArray retval (x.dims ());
  double *r_data = retval.fortran_vec ();

  const double *x_data = x.data ();
  const double *y_data = y.data ();

  octave_idx_type nel = x.numel ();

  for (octave_idx_type i = 0; i < nel; i++)
    {
      OCTAVE_QUIT;
      r_data[i] = f (x_data[i], y_data[i]);
    }

  return retval;
}

static FloatNDArray
map_fm_fm (f_ff_fcn f, const FloatNDArray& x, const FloatNDArray& y)
{
  assert (x.dims () == y.dims ());

  FloatNDArray retval (x.dims ());
  float *r_data = retval.fortran_vec ();

  const float *x_data = x.data ();
  const float *y_data = y.data ();

  octave_idx_type nel = x.numel ();

  for (octave_idx_type i = 0; i < nel; i++)
    {
      OCTAVE_QUIT;
      r_data[i] = f (x_data[i], y_data[i]);
    }

  return retval;
}

static SparseMatrix
map_d_s (d_dd_fcn f, double x, const SparseMatrix& y)
{
  octave_idx_type nr = y.rows ();
  octave_idx_type nc = y.columns ();
  SparseMatrix retval;
  double f_zero = f (x, 0.);

  if (f_zero != 0)
    {
      retval = SparseMatrix (nr, nc, f_zero);

      for (octave_idx_type j = 0; j < nc; j++)
	for (octave_idx_type i = y.cidx (j); i < y.cidx (j+1); i++)
	  {
	    OCTAVE_QUIT;
	    retval.data (y.ridx(i) + j * nr) = f (x, y.data (i));
	  } 

      retval.maybe_compress (true);
    }
  else
    {
      octave_idx_type nz = y.nnz ();
      retval = SparseMatrix (nr, nc, nz);
      octave_idx_type ii = 0;
      retval.cidx (ii) = 0;

      for (octave_idx_type j = 0; j < nc; j++)
	{
	  for (octave_idx_type i = y.cidx (j); i < y.cidx (j+1); i++)
	    {
	      OCTAVE_QUIT;
	      double val = f (x, y.data (i));

	      if (val != 0.0)
		{
		  retval.data (ii) = val;
		  retval.ridx (ii++) = y.ridx (i);
		}
	    } 
	  retval.cidx (j + 1) = ii;
	}

      retval.maybe_compress (false);
    }

  return retval;
}

static SparseMatrix
map_s_d (d_dd_fcn f, const SparseMatrix& x, double y)
{
  octave_idx_type nr = x.rows ();
  octave_idx_type nc = x.columns ();
  SparseMatrix retval;
  double f_zero = f (0., y);

  if (f_zero != 0)
    {
      retval = SparseMatrix (nr, nc, f_zero);

      for (octave_idx_type j = 0; j < nc; j++)
	for (octave_idx_type i = x.cidx (j); i < x.cidx (j+1); i++)
	  {
	    OCTAVE_QUIT;
	    retval.data (x.ridx(i) + j * nr) = f (x.data (i), y);
	  } 

      retval.maybe_compress (true);
    }
  else
    {
      octave_idx_type nz = x.nnz ();
      retval = SparseMatrix (nr, nc, nz);
      octave_idx_type ii = 0;
      retval.cidx (ii) = 0;

      for (octave_idx_type j = 0; j < nc; j++)
	{
	  for (octave_idx_type i = x.cidx (j); i < x.cidx (j+1); i++)
	    {
	      OCTAVE_QUIT;
	      double val = f (x.data (i), y);

	      if (val != 0.0)
		{
		  retval.data (ii) = val;
		  retval.ridx (ii++) = x.ridx (i);
		}
	    } 
	  retval.cidx (j + 1) = ii;
	}

      retval.maybe_compress (false);
    }

  return retval;
}

static SparseMatrix
map_s_s (d_dd_fcn f, const SparseMatrix& x, const SparseMatrix& y)
{
  octave_idx_type nr = x.rows ();
  octave_idx_type nc = x.columns ();

  octave_idx_type y_nr = y.rows ();
  octave_idx_type y_nc = y.columns ();

  assert (nr == y_nr && nc == y_nc);

  SparseMatrix retval;
  double f_zero = f (0., 0.);

  if (f_zero != 0)
    {
      retval = SparseMatrix (nr, nc, f_zero);
      octave_idx_type k1 = 0, k2 = 0;

      for (octave_idx_type j = 0; j < nc; j++)
	{
	  while (k1 < x.cidx(j+1) && k2 < y.cidx(j+1))
	    {
	      OCTAVE_QUIT;
	      if (k1 >= x.cidx(j+1))
		{
		  retval.data (y.ridx(k2) + j * nr) = f (0.0, y.data (k2));
		  k2++;
		}
	      else if (k2 >= y.cidx(j+1))
		{
		  retval.data (x.ridx(k1) + j * nr) = f (x.data (k1), 0.0);
		  k1++;
		}
	      else
		{
		  octave_idx_type rx = x.ridx(k1);
		  octave_idx_type ry = y.ridx(k2);

		  if (rx < ry)
		    {
		      retval.data (rx + j * nr) = f (x.data (k1), 0.0);
		      k1++;
		    }
		  else if (rx > ry)
		    {
		      retval.data (ry + j * nr) = f (0.0, y.data (k2));
		      k2++;
		    }
		  else
		    {
		      retval.data (ry + j * nr) = f (x.data (k1), y.data (k2));
		      k1++;
		      k2++;
		    }
		}
	    }
	}

      retval.maybe_compress (true);
    }
  else
    {
      octave_idx_type nz = x.nnz () + y.nnz ();
      retval = SparseMatrix (nr, nc, nz);
      octave_idx_type ii = 0;
      retval.cidx (ii) = 0;
      octave_idx_type k1 = 0, k2 = 0;

      for (octave_idx_type j = 0; j < nc; j++)
	{
	  while (k1 < x.cidx(j+1) && k2 < y.cidx(j+1))
	    {
	      OCTAVE_QUIT;
	      double val;
	      octave_idx_type r;
	      if (k1 >= x.cidx(j+1))
		{
		  r = y.ridx (k2);
		  val = f (0.0, y.data (k2++));
		}
	      else if (k2 >= y.cidx(j+1))
		{
		  r = x.ridx (k1);
		  val = f (x.data (k1++), 0.0);
		}
	      else
		{
		  octave_idx_type rx = x.ridx(k1);
		  octave_idx_type ry = y.ridx(k2);

		  if (rx < ry)
		    {
		      r = x.ridx (k1);
		      val = f (x.data (k1++), 0.0);
		    }
		  else if (rx > ry)
		    {
		      r = y.ridx (k2);
		      val = f (0.0, y.data (k2++));
		    }
		  else
		    {
		      r = y.ridx (k2);
		      val = f (x.data (k1++), y.data (k2++));
		    }
		}
	      if (val != 0.0)
		{
		  retval.data (ii) = val;
		  retval.ridx (ii++) = r;
		}
	    }
	  retval.cidx (j + 1) = ii;
	}

      retval.maybe_compress (false);
    }

  return retval;
}

DEFUN (atan2, args, ,
  "-*- texinfo -*-\n\
@deftypefn {Mapping Function} {} atan2 (@var{y}, @var{x})\n\
Compute atan (@var{y} / @var{x}) for corresponding elements of @var{y}\n\
and @var{x}.  Signal an error if @var{y} and @var{x} do not match in size\n\
and orientation.\n\
@end deftypefn")
{
  octave_value retval;

  int nargin = args.length ();

  if (nargin == 2 && args(0).is_defined () && args(1).is_defined ())
    {
      if (args(0).is_integer_type () || args(1).is_integer_type ())
	error ("atan2: not defined for integer types");
      else
	{
	  octave_value arg_y = args(0);
	  octave_value arg_x = args(1);

	  dim_vector y_dims = arg_y.dims ();
	  dim_vector x_dims = arg_x.dims ();

	  bool y_is_scalar = y_dims.all_ones ();
	  bool x_is_scalar = x_dims.all_ones ();

	  bool is_float = arg_y.is_single_type () || arg_x.is_single_type ();

	  if (y_is_scalar && x_is_scalar)
	    {
	      if (is_float)
		{
		  float y = arg_y.float_value ();

		  if (! error_state)
		    {
		      float x = arg_x.float_value ();

		      if (! error_state)
			retval = atan2f (y, x);
		    }
		}
	      else
		{
		  double y = arg_y.double_value ();

		  if (! error_state)
		    {
		      double x = arg_x.double_value ();

		      if (! error_state)
			retval = atan2 (y, x);
		    }
		}
	    }
	  else if (y_is_scalar)
	    {
	      if (is_float)
		{
		  float y = arg_y.float_value ();

		  if (! error_state)
		    {
		      // Even if x is sparse return a full matrix here
		      FloatNDArray x = arg_x.float_array_value ();

		      if (! error_state)
			retval = map_f_fm (atan2f, y, x);
		    }
		}
	      else
		{
		  double y = arg_y.double_value ();

		  if (! error_state)
		    {
		      // Even if x is sparse return a full matrix here
		      NDArray x = arg_x.array_value ();

		      if (! error_state)
			retval = map_d_m (atan2, y, x);
		    }
		}
	    }
	  else if (x_is_scalar)
	    {
	      if (arg_y.is_sparse_type ())
		{
		  SparseMatrix y = arg_y.sparse_matrix_value ();

		  if (! error_state)
		    {
		      double x = arg_x.double_value ();
		      
		      if (! error_state)
			retval = map_s_d (atan2, y, x);
		    }
		}
	      else if (is_float)
		{
		  FloatNDArray y = arg_y.float_array_value ();
		  
		  if (! error_state)
		    {
		      float x = arg_x.float_value ();

		      if (! error_state)
			retval = map_fm_f (atan2f, y, x);
		    }
		}
	      else
		{
		  NDArray y = arg_y.array_value ();

		  if (! error_state)
		    {
		      double x = arg_x.double_value ();

		      if (! error_state)
			retval = map_m_d (atan2, y, x);
		    }
		}
	    }
	  else if (y_dims == x_dims)
	    {
	      // Even if y is sparse return a full matrix here
	      if (arg_x.is_sparse_type ())
		{
		  SparseMatrix y = arg_y.sparse_matrix_value ();

		  if (! error_state)
		    {
		      SparseMatrix x = arg_x.sparse_matrix_value ();

		      if (! error_state)
			retval = map_s_s (atan2, y, x);
		    }
		}
	      else if (is_float)
		{
		  FloatNDArray y = arg_y.array_value ();

		  if (! error_state)
		    {
		      FloatNDArray x = arg_x.array_value ();

		      if (! error_state)
			retval = map_fm_fm (atan2f, y, x);
		    }
		}
	      else
		{
		  NDArray y = arg_y.array_value ();

		  if (! error_state)
		    {
		      NDArray x = arg_x.array_value ();

		      if (! error_state)
			retval = map_m_m (atan2, y, x);
		    }
		}
	    }
	  else
	    error ("atan2: nonconformant matrices");
	}
    }
  else
    print_usage ();

  return retval;
}

/*
%!assert (size (atan2 (zeros (0, 2), zeros (0, 2))), [0, 2])
%!assert (size (atan2 (rand (2, 3, 4), zeros (2, 3, 4))), [2, 3, 4])
%!assert (size (atan2 (rand (2, 3, 4), 1)), [2, 3, 4])
%!assert (size (atan2 (1, rand (2, 3, 4))), [2, 3, 4])
%!assert (size (atan2 (1, 2)), [1, 1])

%!test
%! rt2 = sqrt (2);
%! rt3 = sqrt (3);
%! v = [0, pi/6, pi/4, pi/3, -pi/3, -pi/4, -pi/6, 0];
%! y = [0, rt3, 1, rt3, -rt3, -1, -rt3, 0];
%! x = [1, 3, 1, 1, 1, 1, 3, 1];
%! assert(atan2 (y, x), v, sqrt (eps));

%!test
%! rt2 = sqrt (2);
%! rt3 = sqrt (3);
%! v = single([0, pi/6, pi/4, pi/3, -pi/3, -pi/4, -pi/6, 0]);
%! y = single([0, rt3, 1, rt3, -rt3, -1, -rt3, 0]);
%! x = single([1, 3, 1, 1, 1, 1, 3, 1]);
%! assert(atan2 (y, x), v, sqrt (eps('single')));

%!error <Invalid call to atan2.*> atan2 ();
%!error <Invalid call to atan2.*> atan2 (1, 2, 3);

*/



DEFUN (hypot, args, ,
  "-*- texinfo -*-\n\
@deftypefn {Built-in Function} {} hypot (@var{x}, @var{y})\n\
Compute the element-by-element square root of the sum of the squares of\n\
@var{x} and @var{y}.  This is equivalent to\n\
@code{sqrt (@var{x}.^2 + @var{y}.^2)}, but calculated in a manner that\n\
avoids overflows for large values of @var{x} or @var{y}.\n\
@end deftypefn")
{
  octave_value retval;

  int nargin = args.length ();

  if (nargin == 2 && args(0).is_defined () && args(1).is_defined ())
    {
      if (args(0).is_integer_type () || args(1).is_integer_type ())
	error ("hypot: not defined for integer types");
      else
	{
	  octave_value arg_x = args(0);
	  octave_value arg_y = args(1);

	  dim_vector x_dims = arg_x.dims ();
	  dim_vector y_dims = arg_y.dims ();

	  bool x_is_scalar = x_dims.all_ones ();
	  bool y_is_scalar = y_dims.all_ones ();

	  bool is_float = arg_y.is_single_type () || arg_x.is_single_type ();

	  if (y_is_scalar && x_is_scalar)
	    {
	      if (is_float)
		{
		  float x;
		  if (arg_x.is_complex_type ())
		    x = abs (arg_x.float_complex_value ());
		  else
		    x = arg_x.float_value ();

		  if (! error_state)
		    {
		      float y;
		      if (arg_y.is_complex_type ())
			y = abs (arg_y.float_complex_value ());
		      else
			y = arg_y.float_value ();

		      if (! error_state)
			retval = hypotf (x, y);
		    }
		}
	      else
		{
		  double x;
		  if (arg_x.is_complex_type ())
		    x = abs (arg_x.complex_value ());
		  else
		    x = arg_x.double_value ();

		  if (! error_state)
		    {
		      double y;
		      if (arg_y.is_complex_type ())
			y = abs (arg_y.complex_value ());
		      else
			y = arg_y.double_value ();

		      if (! error_state)
			retval = hypot (x, y);
		    }
		}
	    }
	  else if (y_is_scalar)
	    {
	      if (is_float)
		{
		  FloatNDArray x;
		  if (arg_x.is_complex_type ())
		    x = arg_x.float_complex_array_value ().abs ();
		  else
		    x = arg_x.float_array_value ();

		  if (! error_state)
		    {
		      float y;
		      if (arg_y.is_complex_type ())
			y = abs (arg_y.float_complex_value ());
		      else
			y = arg_y.float_value ();

		      if (! error_state)
			retval = map_fm_f (hypotf, x, y);
		    }
		}
	      else
		{
		  NDArray x;
		  if (arg_x.is_complex_type ())
		    x = arg_x.complex_array_value ().abs ();
		  else
		    x = arg_x.array_value ();

		  if (! error_state)
		    {
		      double y;
		      if (arg_y.is_complex_type ())
			y = abs (arg_y.complex_value ());
		      else
			y = arg_y.double_value ();

		      if (! error_state)
			retval = map_m_d (hypot, x, y);
		    }
		}
	    }
	  else if (x_is_scalar)
	    {
	      if (is_float)
		{
		  float x;
		  if (arg_x.is_complex_type ())
		    x = abs (arg_x.float_complex_value ());
		  else
		    x = arg_x.float_value ();

		  if (! error_state)
		    {
		      FloatNDArray y;
		      if (arg_y.is_complex_type ())
			y = arg_y.float_complex_array_value ().abs ();
		      else
			y = arg_y.float_array_value ();

		      if (! error_state)
			retval = map_f_fm (hypotf, x, y);
		    }
		}
	      else
		{
		  double x;
		  if (arg_x.is_complex_type ())
		    x = abs (arg_x.complex_value ());
		  else
		    x = arg_x.double_value ();

		  if (! error_state)
		    {
		      NDArray y;
		      if (arg_y.is_complex_type ())
			y = arg_y.complex_array_value ().abs ();
		      else
			y = arg_y.array_value ();

		      if (! error_state)
			retval = map_d_m (hypot, x, y);
		    }
		}
	    }
	  else if (y_dims == x_dims)
	    {
	      if (arg_x.is_sparse_type () && arg_y.is_sparse_type ())
		{
		  SparseMatrix x;
		  if (arg_x.is_complex_type ())
		    x = arg_x.sparse_complex_matrix_value ().abs ();
		  else
		    x = arg_x.sparse_matrix_value ();

		  if (! error_state)
		    {
		      SparseMatrix y;
		      if (arg_y.is_complex_type ())
			y = arg_y.sparse_complex_matrix_value ().abs ();
		      else
			y = arg_y.sparse_matrix_value ();

		      if (! error_state)
			retval = map_s_s (hypot, x, y);
		    }
		}
	      else if (is_float)
		{
		  FloatNDArray x;
		  if (arg_x.is_complex_type ())
		    x = arg_x.float_complex_array_value ().abs ();
		  else
		    x = arg_x.float_array_value ();

		  if (! error_state)
		    {
		      FloatNDArray y;
		      if (arg_y.is_complex_type ())
			y = arg_y.float_complex_array_value ().abs ();
		      else
			y = arg_y.float_array_value ();

		      if (! error_state)
			retval = map_fm_fm (hypotf, x, y);
		    }
		}
	      else
		{
		  NDArray x;
		  if (arg_x.is_complex_type ())
		    x = arg_x.complex_array_value ().abs ();
		  else
		    x = arg_x.array_value ();

		  if (! error_state)
		    {
		      NDArray y;
		      if (arg_y.is_complex_type ())
			y = arg_y.complex_array_value ().abs ();
		      else
			y = arg_y.array_value ();

		      if (! error_state)
			retval = map_m_m (hypot, x, y);
		    }
		}
	    }
	  else
	    error ("hypot: nonconformant matrices");
	}
    }
  else
    print_usage ();

  return retval;
}

/*
%!assert (size (hypot (zeros (0, 2), zeros (0, 2))), [0, 2])
%!assert (size (hypot (rand (2, 3, 4), zeros (2, 3, 4))), [2, 3, 4])
%!assert (size (hypot (rand (2, 3, 4), 1)), [2, 3, 4])
%!assert (size (hypot (1, rand (2, 3, 4))), [2, 3, 4])
%!assert (size (hypot (1, 2)), [1, 1])
%!assert (hypot (1:10, 1:10), sqrt(2) * [1:10], 16*eps)
%!assert (hypot (single(1:10), single(1:10)), single(sqrt(2) * [1:10]));
*/

template<typename T, typename ET>
void 
map_2_xlog2 (const Array<T>& x, Array<T>& f, Array<ET>& e)
{
  f = Array<T>(x.dims ());
  e = Array<ET>(x.dims ());
  for (octave_idx_type i = 0; i < x.numel (); i++)
    {
      int exp;
      f.xelem (i) = xlog2 (x(i), exp);
      e.xelem (i) = exp;
    }
}

DEFUN (log2, args, nargout,
  "-*- texinfo -*-\n\
@deftypefn {Mapping Function} {} log2 (@var{x})\n\
@deftypefnx {Mapping Function} {[@var{f}, @var{e}] =} log2 (@var{x})\n\
Compute the base-2 logarithm of each element of @var{x}.\n\
\n\
If called with two output arguments, split @var{x} into\n\
binary mantissa and exponent so that\n\
@tex\n\
${1 \\over 2} \\le \\left| f \\right| < 1$\n\
@end tex\n\
@ifnottex\n\
@code{1/2 <= abs(f) < 1}\n\
@end ifnottex\n\
and @var{e} is an integer.  If\n\
@tex\n\
$x = 0$, $f = e = 0$.\n\
@end tex\n\
@ifnottex\n\
@code{x = 0}, @code{f = e = 0}.\n\
@end ifnottex\n\
@seealso{pow2, log, log10, exp}\n\
@end deftypefn")
{
  octave_value_list retval;

  if (args.length () == 1)
    {
      if (nargout < 2)
        retval(0) = args(0).log2 ();
      else if (args(0).is_single_type ())
	{
	  if (args(0).is_real_type ())
	    {
	      FloatNDArray f;
	      FloatNDArray x = args(0).float_array_value ();
	      // FIXME -- should E be an int value?
	      FloatMatrix e;
	      map_2_xlog2 (x, f, e);
	      retval (1) = e;
	      retval (0) = f;
	    }
	  else if (args(0).is_complex_type ())
	    {
	      FloatComplexNDArray f;
	      FloatComplexNDArray x = args(0).float_complex_array_value ();
	      // FIXME -- should E be an int value?
	      FloatNDArray e;
	      map_2_xlog2 (x, f, e);
	      retval (1) = e;
	      retval (0) = f;
	    }
	}
      else if (args(0).is_real_type ())
        {
          NDArray f;
          NDArray x = args(0).array_value ();
          // FIXME -- should E be an int value?
          Matrix e;
          map_2_xlog2 (x, f, e);
          retval (1) = e;
          retval (0) = f;
        }
      else if (args(0).is_complex_type ())
        {
          ComplexNDArray f;
          ComplexNDArray x = args(0).complex_array_value ();
          // FIXME -- should E be an int value?
          NDArray e;
          map_2_xlog2 (x, f, e);
          retval (1) = e;
          retval (0) = f;
        }
      else
        gripe_wrong_type_arg ("log2", args(0));
    }
  else
    print_usage ();

  return retval;
}

/*
%!assert(log2 ([1/4, 1/2, 1, 2, 4]), [-2, -1, 0, 1, 2]);
%!assert(log2(Inf), Inf);
%!assert(isnan(log2(NaN)));
%!assert(log2(4*i), 2 + log2(1*i));
%!assert(log2(complex(0,Inf)), Inf + log2(i));

%!test
%! [f, e] = log2 ([0,-1; 2,-4; Inf,-Inf]);
%! assert (f, [0,-0.5; 0.5,-0.5; Inf,-Inf]);
%! assert (e(1:2,:), [0,1;2,3])

%!test
%! [f, e] = log2 (complex (zeros (3, 2), [0,-1; 2,-4; Inf,-Inf]));
%! assert (f, complex (zeros (3, 2), [0,-0.5; 0.5,-0.5; Inf,-Inf]));
%! assert (e(1:2,:), [0,1; 2,3]);
*/

DEFUN (fmod, args, ,
  "-*- texinfo -*-\n\
@deftypefn {Mapping Function} {} fmod (@var{x}, @var{y})\n\
Compute the floating point remainder of dividing @var{x} by @var{y}\n\
using the C library function @code{fmod}.  The result has the same\n\
sign as @var{x}.  If @var{y} is zero, the result is implementation-dependent.\n\
@seealso{mod, rem}\n\
@end deftypefn")
{
  octave_value retval;

  int nargin = args.length ();

  if (nargin == 2 && args(0).is_defined () && args(1).is_defined ())
    {
      octave_value arg_x = args(0);
      octave_value arg_y = args(1);

      dim_vector y_dims = arg_y.dims ();
      dim_vector x_dims = arg_x.dims ();

      bool y_is_scalar = y_dims.all_ones ();
      bool x_is_scalar = x_dims.all_ones ();

      bool is_float = arg_y.is_single_type () || arg_x.is_single_type ();

      if (y_is_scalar && x_is_scalar)
	{
	  if (is_float)
	    {
	      float y = arg_y.float_value ();

	      if (! error_state)
		{
		  float x = arg_x.float_value ();

		  if (! error_state)
		    retval = fmod (x, y);
		}
	    }
	  else
	    {
	      double y = arg_y.double_value ();

	      if (! error_state)
		{
		  double x = arg_x.double_value ();

		  if (! error_state)
		    retval = fmod (x, y);
		}
	    }
	}
      else if (y_is_scalar)
	{
	  if (is_float)
	    {
	      float y = arg_y.float_value ();

	      if (! error_state)
		{
		  FloatNDArray x = arg_x.float_array_value ();

		  if (! error_state)
		    retval = map_fm_f (fmodf, x, y);
		}
	    }
	  else
	    {
	      double y = arg_y.double_value ();

	      if (! error_state)
		{
		  if (arg_x.is_sparse_type ())
		    {
		      SparseMatrix x = arg_x.sparse_matrix_value ();

		      if (! error_state)
			retval = map_s_d (fmod, x, y);
		    }
		  else
		    {
		      NDArray x = arg_x.array_value ();

		      if (! error_state)
			retval = map_m_d (fmod, x, y);
		    }
		}
	    }
	}
      else if (x_is_scalar)
	{
	  if (arg_y.is_sparse_type ())
	    {
	      SparseMatrix y = arg_y.sparse_matrix_value ();

	      if (! error_state)
		{
		  double x = arg_x.double_value ();

		  if (! error_state)
		    retval = map_d_s (fmod, x, y);
		}
	    }
	  else if (is_float)
	    {
	      FloatNDArray y = arg_y.float_array_value ();

	      if (! error_state)
		{
		  float x = arg_x.float_value ();

		  if (! error_state)
		    retval = map_f_fm (fmodf, x, y);
		}
	    }
	  else
	    {
	      NDArray y = arg_y.array_value ();

	      if (! error_state)
		{
		  double x = arg_x.double_value ();

		  if (! error_state)
		    retval = map_d_m (fmod, x, y);
		}
	    }
	}
      else if (y_dims == x_dims)
	{
	  if (arg_y.is_sparse_type () || arg_x.is_sparse_type ())
	    {
	      SparseMatrix y = arg_y.sparse_matrix_value ();

	      if (! error_state)
		{
		  SparseMatrix x = arg_x.sparse_matrix_value ();

		  if (! error_state)
		    retval = map_s_s (fmod, x, y);
		}
	    }
	  else if (is_float)
	    {
	      FloatNDArray y = arg_y.float_array_value ();

	      if (! error_state)
		{
		  FloatNDArray x = arg_x.float_array_value ();

		  if (! error_state)
		    retval = map_fm_fm (fmodf, x, y);
		}
	    }
	  else
	    {
	      NDArray y = arg_y.array_value ();

	      if (! error_state)
		{
		  NDArray x = arg_x.array_value ();

		  if (! error_state)
		    retval = map_m_m (fmod, x, y);
		}
	    }
	}
      else
	error ("fmod: nonconformant matrices");
    }
  else
    print_usage ();

  return retval;
}

/*
%!assert (size (fmod (zeros (0, 2), zeros (0, 2))), [0, 2])
%!assert (size (fmod (rand (2, 3, 4), zeros (2, 3, 4))), [2, 3, 4])
%!assert (size (fmod (rand (2, 3, 4), 1)), [2, 3, 4])
%!assert (size (fmod (1, rand (2, 3, 4))), [2, 3, 4])
%!assert (size (fmod (1, 2)), [1, 1])
*/

// FIXME Need to convert the reduction functions of this file for single precision

#define NATIVE_REDUCTION_1(FCN, TYPE, DIM) \
  (arg.is_ ## TYPE ## _type ()) \
    { \
      TYPE ## NDArray tmp = arg. TYPE ##_array_value (); \
      \
      if (! error_state) \
        { \
          octave_ ## TYPE::clear_conv_flags (); \
          retval = tmp.FCN (DIM); \
          if (octave_ ## TYPE::get_trunc_flag ()) \
            { \
              gripe_native_integer_math_truncated (#FCN, \
                                                   octave_ ## TYPE::type_name ()); \
              octave_ ## TYPE::clear_conv_flags (); \
            } \
        } \
    }

#define NATIVE_REDUCTION(FCN, BOOL_FCN) \
 \
  octave_value retval; \
 \
  int nargin = args.length (); \
 \
  bool isnative = false; \
  bool isdouble = false; \
  \
  if (nargin > 1 && args(nargin - 1).is_string ()) \
    { \
      std::string str = args(nargin - 1).string_value (); \
      \
      if (! error_state) \
	{ \
	  if (str == "native") \
	    isnative = true; \
	  else if (str == "double") \
            isdouble = true; \
          else \
	    error ("sum: unrecognized string argument"); \
          nargin --; \
	} \
    } \
  \
  if (nargin == 1 || nargin == 2) \
    { \
      octave_value arg = args(0); \
 \
      int dim = (nargin == 1 ? -1 : args(1).int_value (true) - 1); \
 \
      if (! error_state) \
	{ \
	  if (dim >= -1) \
	    { \
	      if (arg.is_sparse_type ()) \
		{ \
		  if (arg.is_real_type ()) \
		    { \
		      SparseMatrix tmp = arg.sparse_matrix_value (); \
		      \
		      if (! error_state) \
			retval = tmp.FCN (dim); \
		    } \
		  else \
		    { \
		      SparseComplexMatrix tmp = arg.sparse_complex_matrix_value (); \
                      \
		      if (! error_state) \
			retval = tmp.FCN (dim); \
		    } \
		} \
	      else \
		{ \
		  if (isnative)	\
		    { \
		      if NATIVE_REDUCTION_1 (FCN, uint8, dim) \
		      else if NATIVE_REDUCTION_1 (FCN, uint16, dim) \
                      else if NATIVE_REDUCTION_1 (FCN, uint32, dim) \
                      else if NATIVE_REDUCTION_1 (FCN, uint64, dim) \
                      else if NATIVE_REDUCTION_1 (FCN, int8, dim) \
                      else if NATIVE_REDUCTION_1 (FCN, int16, dim) \
                      else if NATIVE_REDUCTION_1 (FCN, int32, dim) \
                      else if NATIVE_REDUCTION_1 (FCN, int64, dim) \
                      else if (arg.is_bool_type ()) \
                        { \
                          boolNDArray tmp = arg.bool_array_value (); \
                          if (! error_state) \
                            retval = boolNDArray (tmp.BOOL_FCN (dim)); \
                        } \
                      else if (arg.is_char_matrix ()) \
                        { \
			  error (#FCN, ": invalid char type"); \
			} \
		      else if (!isdouble && arg.is_single_type ()) \
                        { \
	                  if (arg.is_complex_type ()) \
		            { \
		              FloatComplexNDArray tmp = \
				arg.float_complex_array_value (); \
                              \
		              if (! error_state) \
		                retval = tmp.FCN (dim); \
		            } \
	                  else if (arg.is_real_type ()) \
		            { \
		              FloatNDArray tmp = arg.float_array_value (); \
                              \
		              if (! error_state) \
		                retval = tmp.FCN (dim); \
		            } \
			} \
	              else if (arg.is_complex_type ()) \
		        { \
		          ComplexNDArray tmp = arg.complex_array_value (); \
                          \
		          if (! error_state) \
		            retval = tmp.FCN (dim); \
		        } \
	              else if (arg.is_real_type ()) \
		        { \
		          NDArray tmp = arg.array_value (); \
                          \
		          if (! error_state) \
		            retval = tmp.FCN (dim); \
		        } \
                      else \
		        { \
		          gripe_wrong_type_arg (#FCN, arg); \
		          return retval; \
		        } \
                    } \
                  else if (arg.is_bool_type ()) \
                    { \
                      boolNDArray tmp = arg.bool_array_value (); \
                      if (! error_state) \
                        retval = tmp.FCN (dim); \
                    } \
		  else if (!isdouble && arg.is_single_type ()) \
		    { \
	              if (arg.is_real_type ()) \
		        { \
		          FloatNDArray tmp = arg.float_array_value (); \
                          \
		          if (! error_state) \
		            retval = tmp.FCN (dim); \
		        } \
	              else if (arg.is_complex_type ()) \
		        { \
		          FloatComplexNDArray tmp = \
			    arg.float_complex_array_value (); \
                          \
		          if (! error_state) \
		            retval = tmp.FCN (dim); \
		        } \
		    } \
	          else if (arg.is_real_type ()) \
		    { \
		      NDArray tmp = arg.array_value (); \
                      \
		      if (! error_state) \
		        retval = tmp.FCN (dim); \
		    } \
	          else if (arg.is_complex_type ()) \
		    { \
		      ComplexNDArray tmp = arg.complex_array_value (); \
                      \
		      if (! error_state) \
		        retval = tmp.FCN (dim); \
		    } \
	          else \
		    { \
		      gripe_wrong_type_arg (#FCN, arg); \
		      return retval; \
		    } \
		} \
	    } \
	  else \
	    error (#FCN ": invalid dimension argument = %d", dim + 1); \
	} \
      \
    } \
  else \
    print_usage (); \
 \
  return retval

#define DATA_REDUCTION(FCN) \
 \
  octave_value retval; \
 \
  int nargin = args.length (); \
 \
  if (nargin == 1 || nargin == 2) \
    { \
      octave_value arg = args(0); \
 \
      int dim = (nargin == 1 ? -1 : args(1).int_value (true) - 1); \
 \
      if (! error_state) \
	{ \
	  if (dim >= -1) \
	    { \
	      if (arg.is_real_type ()) \
		{ \
		  if (arg.is_sparse_type ()) \
		    { \
		      SparseMatrix tmp = arg.sparse_matrix_value (); \
 \
		      if (! error_state) \
			retval = tmp.FCN (dim); \
		    } \
		  else if (arg.is_single_type ()) \
		    { \
		      FloatNDArray tmp = arg.float_array_value (); \
 \
		      if (! error_state) \
			retval = tmp.FCN (dim); \
		    } \
		  else \
		    { \
		      NDArray tmp = arg.array_value (); \
 \
		      if (! error_state) \
			retval = tmp.FCN (dim); \
		    } \
		} \
	      else if (arg.is_complex_type ()) \
		{ \
		  if (arg.is_sparse_type ()) \
		    { \
		      SparseComplexMatrix tmp = arg.sparse_complex_matrix_value (); \
 \
		      if (! error_state) \
			retval = tmp.FCN (dim); \
		    } \
		  else if (arg.is_single_type ()) \
		    { \
		      FloatComplexNDArray tmp = arg.float_complex_array_value (); \
 \
		      if (! error_state) \
			retval = tmp.FCN (dim); \
		    } \
		  else \
		    { \
		      ComplexNDArray tmp = arg.complex_array_value (); \
 \
		      if (! error_state) \
			retval = tmp.FCN (dim); \
		    } \
		} \
	      else \
		{ \
		  gripe_wrong_type_arg (#FCN, arg); \
		  return retval; \
		} \
	    } \
	  else \
	    error (#FCN ": invalid dimension argument = %d", dim + 1); \
	} \
    } \
  else \
    print_usage (); \
 \
  return retval

DEFUN (cumprod, args, ,
  "-*- texinfo -*-\n\
@deftypefn  {Built-in Function} {} cumprod (@var{x})\n\
@deftypefnx {Built-in Function} {} cumprod (@var{x}, @var{dim})\n\
Cumulative product of elements along dimension @var{dim}.  If\n\
@var{dim} is omitted, it defaults to the first non-singleton dimension.\n\
\n\
@seealso{prod, cumsum}\n\
@end deftypefn")
{
  DATA_REDUCTION (cumprod);
}

/*

%!assert (cumprod ([1, 2, 3]), [1, 2, 6]);
%!assert (cumprod ([-1; -2; -3]), [-1; 2; -6]);
%!assert (cumprod ([i, 2+i, -3+2i, 4]), [i, -1+2i, -1-8i, -4-32i]);
%!assert (cumprod ([1, 2, 3; i, 2i, 3i; 1+i, 2+2i, 3+3i]), [1, 2, 3; i, 4i, 9i; -1+i, -8+8i, -27+27i]);

%!assert (cumprod (single([1, 2, 3])), single([1, 2, 6]));
%!assert (cumprod (single([-1; -2; -3])), single([-1; 2; -6]));
%!assert (cumprod (single([i, 2+i, -3+2i, 4])), single([i, -1+2i, -1-8i, -4-32i]));
%!assert (cumprod (single([1, 2, 3; i, 2i, 3i; 1+i, 2+2i, 3+3i])), single([1, 2, 3; i, 4i, 9i; -1+i, -8+8i, -27+27i]));

%!error <Invalid call to cumprod.*> cumprod ();

%!assert (cumprod ([2, 3; 4, 5], 1), [2, 3; 8, 15]);
%!assert (cumprod ([2, 3; 4, 5], 2), [2, 6; 4, 20]);

%!assert (cumprod (single([2, 3; 4, 5]), 1), single([2, 3; 8, 15]));
%!assert (cumprod (single([2, 3; 4, 5]), 2), single([2, 6; 4, 20]));

 */

DEFUN (cumsum, args, ,
  "-*- texinfo -*-\n\
@deftypefn  {Built-in Function} {} cumsum (@var{x})\n\
@deftypefnx {Built-in Function} {} cumsum (@var{x}, @var{dim})\n\
@deftypefnx {Built-in Function} {} cumsum (@dots{}, 'native')\n\
Cumulative sum of elements along dimension @var{dim}.  If @var{dim}\n\
is omitted, it defaults to the first non-singleton dimension.\n\
\n\
The \"native\" argument implies the summation is performed in native type.\n\
 See @code{sum} for a complete description and example of the use of\n\
\"native\".\n\
@seealso{sum, cumprod}\n\
@end deftypefn")
{
  octave_value retval;

  int nargin = args.length ();

  bool isnative = false;
  bool isdouble = false;

  if (nargin > 1 && args(nargin - 1).is_string ())
    {
      std::string str = args(nargin - 1).string_value ();

      if (! error_state)
	{
	  if (str == "native")
	    isnative = true;
	  else if (str == "double")
            isdouble = true;
          else
	    error ("sum: unrecognized string argument");
          nargin --;
	}
    }

  if (error_state)
    return retval;

  if (nargin == 1 || nargin == 2)
    {
      octave_value arg = args(0);

      int dim = -1;
      if (nargin == 2)
        {
          dim = args(1).int_value () - 1;
          if (dim < 0)
	    error ("cumsum: invalid dimension argument = %d", dim + 1);
        }

      if (! error_state)
	{
          switch (arg.builtin_type ())
            {
            case btyp_double:
              if (arg.is_sparse_type ())
                retval = arg.sparse_matrix_value ().cumsum (dim);
              else
                retval = arg.array_value ().cumsum (dim);
              break;
            case btyp_complex:
              if (arg.is_sparse_type ())
                retval = arg.sparse_complex_matrix_value ().cumsum (dim);
              else
                retval = arg.complex_array_value ().cumsum (dim);
              break;
            case btyp_float:
              if (isdouble)
                retval = arg.array_value ().cumsum (dim);
              else
                retval = arg.float_array_value ().cumsum (dim);
              break;
            case btyp_float_complex:
              if (isdouble)
                retval = arg.complex_array_value ().cumsum (dim);
              else
                retval = arg.float_complex_array_value ().cumsum (dim);
              break;

#define MAKE_INT_BRANCH(X) \
            case btyp_ ## X: \
              if (isnative) \
                retval = arg.X ## _array_value ().cumsum (dim); \
              else \
                retval = arg.array_value ().cumsum (dim); \
              break
            MAKE_INT_BRANCH (int8);
            MAKE_INT_BRANCH (int16);
            MAKE_INT_BRANCH (int32);
            MAKE_INT_BRANCH (int64);
            MAKE_INT_BRANCH (uint8);
            MAKE_INT_BRANCH (uint16);
            MAKE_INT_BRANCH (uint32);
            MAKE_INT_BRANCH (uint64);
#undef MAKE_INT_BRANCH

            case btyp_bool:
              if (arg.is_sparse_type ())
                {
                  SparseMatrix cs = arg.sparse_matrix_value ().cumsum (dim);
                  if (isnative)
                    retval = cs != 0.0;
                  else
                    retval = cs;
                }
              else
                {
                  NDArray cs = arg.bool_array_value ().cumsum (dim);
                  if (isnative)
                    retval = cs != 0.0;
                  else
                    retval = cs;
                }
              break;

            default:
              gripe_wrong_type_arg ("cumsum", arg);
            }
	}
    }
  else
    print_usage ();

  return retval;
}

/*

%!assert (cumsum ([1, 2, 3]), [1, 3, 6]);
%!assert (cumsum ([-1; -2; -3]), [-1; -3; -6]);
%!assert (cumsum ([i, 2+i, -3+2i, 4]), [i, 2+2i, -1+4i, 3+4i]);
%!assert (cumsum ([1, 2, 3; i, 2i, 3i; 1+i, 2+2i, 3+3i]), [1, 2, 3; 1+i, 2+2i, 3+3i; 2+2i, 4+4i, 6+6i]);

%!assert (cumsum (single([1, 2, 3])), single([1, 3, 6]));
%!assert (cumsum (single([-1; -2; -3])), single([-1; -3; -6]));
%!assert (cumsum (single([i, 2+i, -3+2i, 4])), single([i, 2+2i, -1+4i, 3+4i]));
%!assert (cumsum (single([1, 2, 3; i, 2i, 3i; 1+i, 2+2i, 3+3i])), single([1, 2, 3; 1+i, 2+2i, 3+3i; 2+2i, 4+4i, 6+6i]));

%!error <Invalid call to cumsum.*> cumsum ();

%!assert (cumsum ([1, 2; 3, 4], 1), [1, 2; 4, 6]);
%!assert (cumsum ([1, 2; 3, 4], 2), [1, 3; 3, 7]);

%!assert (cumsum (single([1, 2; 3, 4]), 1), single([1, 2; 4, 6]));
%!assert (cumsum (single([1, 2; 3, 4]), 2), single([1, 3; 3, 7]));

 */

DEFUN (diag, args, ,
  "-*- texinfo -*-\n\
@deftypefn {Built-in Function} {} diag (@var{v}, @var{k})\n\
Return a diagonal matrix with vector @var{v} on diagonal @var{k}.  The\n\
second argument is optional.  If it is positive, the vector is placed on\n\
the @var{k}-th super-diagonal.  If it is negative, it is placed on the\n\
@var{-k}-th sub-diagonal.  The default value of @var{k} is 0, and the\n\
vector is placed on the main diagonal.  For example,\n\
\n\
@example\n\
@group\n\
diag ([1, 2, 3], 1)\n\
     @result{}  0  1  0  0\n\
         0  0  2  0\n\
         0  0  0  3\n\
         0  0  0  0\n\
@end group\n\
@end example\n\
\n\
@noindent\n\
Given a matrix argument, instead of a vector, @code{diag} extracts the\n\
@var{k}-th diagonal of the matrix.\n\
@end deftypefn")
{
  octave_value retval;

  int nargin = args.length ();

  if (nargin == 1 && args(0).is_defined ())
    retval = args(0).diag();
  else if (nargin == 2 && args(0).is_defined () && args(1).is_defined ())
    {
      octave_idx_type k = args(1).int_value ();

      if (error_state)
	error ("diag: invalid second argument");      
      else
	retval = args(0).diag(k);
    }
  else if (nargin == 3)
    {
      octave_value arg0 = args(0);
      if (arg0.ndims () == 2 && (args(0).rows () == 1 || args(0).columns () == 1))
        {
          octave_idx_type m = args(1).int_value (), n = args(2).int_value ();
          if (! error_state)
            retval = arg0.diag ().resize (dim_vector (m, n));
          else
            error ("diag: invalid dimensions");
        }
      else
        error ("diag: first argument must be a vector");
    }
  else
    print_usage ();

  return retval;
}

/*

%!assert(full (diag ([1; 2; 3])), [1, 0, 0; 0, 2, 0; 0, 0, 3]);
%!assert(diag ([1; 2; 3], 1), [0, 1, 0, 0; 0, 0, 2, 0; 0, 0, 0, 3; 0, 0, 0, 0]);
%!assert(diag ([1; 2; 3], 2), [0, 0, 1, 0, 0; 0, 0, 0, 2, 0; 0, 0, 0, 0, 3; 0, 0, 0, 0, 0; 0, 0, 0, 0, 0]);
%!assert(diag ([1; 2; 3],-1), [0, 0, 0, 0; 1, 0, 0, 0; 0, 2, 0, 0; 0, 0, 3, 0]);
%!assert(diag ([1; 2; 3],-2), [0, 0, 0, 0, 0; 0, 0, 0, 0, 0; 1, 0, 0, 0, 0; 0, 2, 0, 0, 0; 0, 0, 3, 0, 0]);

%!assert(diag ([1, 0, 0; 0, 2, 0; 0, 0, 3]), [1; 2; 3]);
%!assert(diag ([0, 1, 0, 0; 0, 0, 2, 0; 0, 0, 0, 3; 0, 0, 0, 0], 1), [1; 2; 3]);
%!assert(diag ([0, 0, 0, 0; 1, 0, 0, 0; 0, 2, 0, 0; 0, 0, 3, 0], -1), [1; 2; 3]);

%!assert(full (diag (single([1; 2; 3]))), single([1, 0, 0; 0, 2, 0; 0, 0, 3]));
%!assert(diag (single([1; 2; 3]), 1), single([0, 1, 0, 0; 0, 0, 2, 0; 0, 0, 0, 3; 0, 0, 0, 0]));
%!assert(diag (single([1; 2; 3]), 2), single([0, 0, 1, 0, 0; 0, 0, 0, 2, 0; 0, 0, 0, 0, 3; 0, 0, 0, 0, 0; 0, 0, 0, 0, 0]));
%!assert(diag (single([1; 2; 3]),-1), single([0, 0, 0, 0; 1, 0, 0, 0; 0, 2, 0, 0; 0, 0, 3, 0]));
%!assert(diag (single([1; 2; 3]),-2), single([0, 0, 0, 0, 0; 0, 0, 0, 0, 0; 1, 0, 0, 0, 0; 0, 2, 0, 0, 0; 0, 0, 3, 0, 0]));

%!assert(diag (single([1, 0, 0; 0, 2, 0; 0, 0, 3])), single([1; 2; 3]));
%!assert(diag (single([0, 1, 0, 0; 0, 0, 2, 0; 0, 0, 0, 3; 0, 0, 0, 0]), 1), single([1; 2; 3]));
%!assert(diag (single([0, 0, 0, 0; 1, 0, 0, 0; 0, 2, 0, 0; 0, 0, 3, 0]), -1), single([1; 2; 3]));

%!assert(diag (int8([1; 2; 3])), int8([1, 0, 0; 0, 2, 0; 0, 0, 3]));
%!assert(diag (int8([1; 2; 3]), 1), int8([0, 1, 0, 0; 0, 0, 2, 0; 0, 0, 0, 3; 0, 0, 0, 0]));
%!assert(diag (int8([1; 2; 3]), 2), int8([0, 0, 1, 0, 0; 0, 0, 0, 2, 0; 0, 0, 0, 0, 3; 0, 0, 0, 0, 0; 0, 0, 0, 0, 0]));
%!assert(diag (int8([1; 2; 3]),-1), int8([0, 0, 0, 0; 1, 0, 0, 0; 0, 2, 0, 0; 0, 0, 3, 0]));
%!assert(diag (int8([1; 2; 3]),-2), int8([0, 0, 0, 0, 0; 0, 0, 0, 0, 0; 1, 0, 0, 0, 0; 0, 2, 0, 0, 0; 0, 0, 3, 0, 0]));

%!assert(diag (int8([1, 0, 0; 0, 2, 0; 0, 0, 3])), int8([1; 2; 3]));
%!assert(diag (int8([0, 1, 0, 0; 0, 0, 2, 0; 0, 0, 0, 3; 0, 0, 0, 0]), 1), int8([1; 2; 3]));
%!assert(diag (int8([0, 0, 0, 0; 1, 0, 0, 0; 0, 2, 0, 0; 0, 0, 3, 0]), -1), int8([1; 2; 3]));

%!error <Invalid call to diag.*> diag ();

 */

DEFUN (prod, args, ,
  "-*- texinfo -*-\n\
@deftypefn  {Built-in Function} {} prod (@var{x})\n\
@deftypefnx {Built-in Function} {} prod (@var{x}, @var{dim})\n\
Product of elements along dimension @var{dim}.  If @var{dim} is\n\
omitted, it defaults to the first non-singleton dimension.\n\
@seealso{cumprod, sum}\n\
@end deftypefn")
{
  DATA_REDUCTION (prod);
}

/*

%!assert (prod ([1, 2, 3]), 6);
%!assert (prod ([-1; -2; -3]), -6);
%!assert (prod ([i, 2+i, -3+2i, 4]), -4 - 32i);
%!assert (prod ([1, 2, 3; i, 2i, 3i; 1+i, 2+2i, 3+3i]), [-1+i, -8+8i, -27+27i]);

%!assert (prod (single([1, 2, 3])), single(6));
%!assert (prod (single([-1; -2; -3])), single(-6));
%!assert (prod (single([i, 2+i, -3+2i, 4])), single(-4 - 32i));
%!assert (prod (single([1, 2, 3; i, 2i, 3i; 1+i, 2+2i, 3+3i])), single([-1+i, -8+8i, -27+27i]));

%!error <Invalid call to prod.*> prod ();

%!assert (prod ([1, 2; 3, 4], 1), [3, 8]);
%!assert (prod ([1, 2; 3, 4], 2), [2; 12]);
%!assert (prod (zeros (1, 0)), 1);
%!assert (prod (zeros (1, 0), 1), zeros (1, 0));
%!assert (prod (zeros (1, 0), 2), 1);
%!assert (prod (zeros (0, 1)), 1);
%!assert (prod (zeros (0, 1), 1), 1);
%!assert (prod (zeros (0, 1), 2), zeros (0, 1));
%!assert (prod (zeros (2, 0)), zeros (1, 0));
%!assert (prod (zeros (2, 0), 1), zeros (1, 0));
%!assert (prod (zeros (2, 0), 2), [1; 1]);
%!assert (prod (zeros (0, 2)), [1, 1]);
%!assert (prod (zeros (0, 2), 1), [1, 1]);
%!assert (prod (zeros (0, 2), 2), zeros(0, 1));

%!assert (prod (single([1, 2; 3, 4]), 1), single([3, 8]));
%!assert (prod (single([1, 2; 3, 4]), 2), single([2; 12]));
%!assert (prod (zeros (1, 0, 'single')), single(1));
%!assert (prod (zeros (1, 0, 'single'), 1), zeros (1, 0, 'single'));
%!assert (prod (zeros (1, 0, 'single'), 2), single(1));
%!assert (prod (zeros (0, 1, 'single')), single(1));
%!assert (prod (zeros (0, 1, 'single'), 1), single(1));
%!assert (prod (zeros (0, 1, 'single'), 2), zeros (0, 1, 'single'));
%!assert (prod (zeros (2, 0, 'single')), zeros (1, 0, 'single'));
%!assert (prod (zeros (2, 0, 'single'), 1), zeros (1, 0, 'single'));
%!assert (prod (zeros (2, 0, 'single'), 2), single([1; 1]));
%!assert (prod (zeros (0, 2, 'single')), single([1, 1]));
%!assert (prod (zeros (0, 2, 'single'), 1), single([1, 1]));
%!assert (prod (zeros (0, 2, 'single'), 2), zeros(0, 1, 'single'));

 */

#define SINGLE_TYPE_CONCAT(TYPE, EXTRACTOR) \
  do \
    { \
      int dv_len = dv.length (); \
      Array<octave_idx_type> ra_idx (dv_len > 1 ? dv_len : 2, 0); \
      \
      for (int j = 1; j < n_args; j++) \
	{ \
	  OCTAVE_QUIT; \
	  \
	  TYPE ra = args(j).EXTRACTOR ();	\
	  \
	  if (! error_state) \
	    { \
	      result.insert (ra, ra_idx); \
	      \
	      if (error_state) \
	        return retval; \
	      \
	      dim_vector dv_tmp = args (j).dims (); \
	      \
	      if (dim >= dv_len) \
	        { \
		  if (j > 1) \
		    error ("%s: indexing error", fname.c_str ()); \
		  break; \
		} \
	      else \
		ra_idx (dim) += (dim < dv_tmp.length () ? dv_tmp (dim) : 1); \
	    } \
	} \
    } \
 while (0)

#define DO_SINGLE_TYPE_CONCAT(TYPE, EXTRACTOR) \
  do \
    { \
      TYPE result (dv); \
      \
      SINGLE_TYPE_CONCAT(TYPE, EXTRACTOR); \
      \
      retval = result; \
    } \
 while (0)

static octave_value
do_cat (const octave_value_list& args, std::string fname)
{
  octave_value retval;

  int n_args = args.length (); 

  if (n_args == 1)
    retval = Matrix ();
  else if (n_args == 2)
    retval = args(1);
  else if (n_args > 2)
    {
      octave_idx_type dim = args(0).int_value () - 1;

      if (error_state)
	{
	  error ("cat: expecting first argument to be a integer");
	  return retval;
	}
  
      if (dim >= 0)
	{
 	  
 	  dim_vector  dv = args(1).dims ();
	  std::string result_type = args(1).class_name ();
	  
	  bool all_sq_strings_p = args(1).is_sq_string ();
	  bool all_dq_strings_p = args(1).is_dq_string ();
	  bool all_real_p = args(1).is_real_type ();
	  bool any_sparse_p = args(1).is_sparse_type();

 	  for (int i = 2; i < args.length (); i++)
  	    {
 	      // add_dims constructs a dimension vector which holds the
	      // dimensions of the final array after concatenation.

	      if (! dv.concat (args(i).dims (), dim))
		{
		  // Dimensions do not match. 
		  error ("cat: dimension mismatch");
		  return retval;
		}
	      
	      result_type = 
		get_concat_class (result_type, args(i).class_name ());

	      if (all_sq_strings_p && ! args(i).is_sq_string ())
		all_sq_strings_p = false;
	      if (all_dq_strings_p && ! args(i).is_dq_string ())
		all_dq_strings_p = false;
	      if (all_real_p && ! args(i).is_real_type ())
		all_real_p = false;
	      if (!any_sparse_p && args(i).is_sparse_type ())
		any_sparse_p = true;
	    }

	  if (result_type == "double")
	    {
	      if (any_sparse_p)
		{	    
		  if (all_real_p)
		    DO_SINGLE_TYPE_CONCAT (SparseMatrix, sparse_matrix_value);
		  else
		    DO_SINGLE_TYPE_CONCAT (SparseComplexMatrix, sparse_complex_matrix_value);
		}
	      else
		{
		  if (all_real_p)
		    DO_SINGLE_TYPE_CONCAT (NDArray, array_value);
		  else
		    DO_SINGLE_TYPE_CONCAT (ComplexNDArray, complex_array_value);
		}
	    }
	  else if (result_type == "single")
	    {
	      if (all_real_p)
		DO_SINGLE_TYPE_CONCAT (FloatNDArray, float_array_value);
	      else
		DO_SINGLE_TYPE_CONCAT (FloatComplexNDArray, 
				       float_complex_array_value);
	    }
	  else if (result_type == "char")
	    {
	      char type = all_dq_strings_p ? '"' : '\'';

	      maybe_warn_string_concat (all_dq_strings_p, all_sq_strings_p);

	      charNDArray result (dv, Vstring_fill_char);

	      SINGLE_TYPE_CONCAT (charNDArray, char_array_value);

	      retval = octave_value (result, type);
	    }
	  else if (result_type == "logical")
	    {
	      if (any_sparse_p)
		DO_SINGLE_TYPE_CONCAT (SparseBoolMatrix, sparse_bool_matrix_value);
	      else
		DO_SINGLE_TYPE_CONCAT (boolNDArray, bool_array_value);
	    }
	  else if (result_type == "int8")
	    DO_SINGLE_TYPE_CONCAT (int8NDArray, int8_array_value);
	  else if (result_type == "int16")
	    DO_SINGLE_TYPE_CONCAT (int16NDArray, int16_array_value);
	  else if (result_type == "int32")
	    DO_SINGLE_TYPE_CONCAT (int32NDArray, int32_array_value);
	  else if (result_type == "int64")
	    DO_SINGLE_TYPE_CONCAT (int64NDArray, int64_array_value);
	  else if (result_type == "uint8")
	    DO_SINGLE_TYPE_CONCAT (uint8NDArray, uint8_array_value);
	  else if (result_type == "uint16")
	    DO_SINGLE_TYPE_CONCAT (uint16NDArray, uint16_array_value);
	  else if (result_type == "uint32")
	    DO_SINGLE_TYPE_CONCAT (uint32NDArray, uint32_array_value);
	  else if (result_type == "uint64")
	    DO_SINGLE_TYPE_CONCAT (uint64NDArray, uint64_array_value);
	  else
	    {
	      // The lines below might seem crazy, since we take a copy
	      // of the first argument, resize it to be empty and then resize
	      // it to be full. This is done since it means that there is no
	      // recopying of data, as would happen if we used a single resize.
	      // It should be noted that resize operation is also significantly 
	      // slower than the do_cat_op function, so it makes sense to have
	      // an empty matrix and copy all data.
	      //
	      // We might also start with a empty octave_value using
	      //   tmp = octave_value_typeinfo::lookup_type 
	      //                                (args(1).type_name());
	      // and then directly resize. However, for some types there might
	      // be some additional setup needed, and so this should be avoided.

	      octave_value tmp = args (1);
	      tmp = tmp.resize (dim_vector (0,0)).resize (dv);

	      if (error_state)
		return retval;

	      int dv_len = dv.length ();
	      Array<octave_idx_type> ra_idx (dv_len, 0);

	      for (int j = 1; j < n_args; j++)
		{
		  // Can't fast return here to skip empty matrices as something
		  // like cat(1,[],single([])) must return an empty matrix of
		  // the right type.
		  tmp = do_cat_op (tmp, args (j), ra_idx);

		  if (error_state)
		    return retval;

		  dim_vector dv_tmp = args (j).dims ();

		  if (dim >= dv_len)
		    {
		      if (j > 1)
			error ("%s: indexing error", fname.c_str ());
		      break;
		    }
		  else
		    ra_idx (dim) += (dim < dv_tmp.length () ? 
				     dv_tmp (dim) : 1);
		}
	      retval = tmp;
	    }

	  if (! error_state)
	    {
	      // Reshape, chopping trailing singleton dimensions
	      dv.chop_trailing_singletons ();
	      retval = retval.reshape (dv);
	    }
	}
      else
	error ("%s: invalid dimension argument", fname.c_str ());
    }
  else
    print_usage ();

  return retval;
}

DEFUN (horzcat, args, ,
       "-*- texinfo -*-\n\
@deftypefn {Built-in Function} {} horzcat (@var{array1}, @var{array2}, @dots{}, @var{arrayN})\n\
Return the horizontal concatenation of N-d array objects, @var{array1},\n\
@var{array2}, @dots{}, @var{arrayN} along dimension 2.\n\
@seealso{cat, vertcat}\n\
@end deftypefn")
{
  octave_value_list args_tmp = args;
  
  int dim = 2;
  
  octave_value d (dim);
  
  args_tmp.prepend (d);
  
  return do_cat (args_tmp, "horzcat");
}

DEFUN (vertcat, args, ,
  "-*- texinfo -*-\n\
@deftypefn {Built-in Function} {} vertcat (@var{array1}, @var{array2}, @dots{}, @var{arrayN})\n\
Return the vertical concatenation of N-d array objects, @var{array1},\n\
@var{array2}, @dots{}, @var{arrayN} along dimension 1.\n\
@seealso{cat, horzcat}\n\
@end deftypefn")
{
  octave_value_list args_tmp = args;
  
  int dim = 1;
  
  octave_value d (dim);
  
  args_tmp.prepend (d);
  
  return do_cat (args_tmp, "vertcat");
}

DEFUN (cat, args, ,
  "-*- texinfo -*-\n\
@deftypefn {Built-in Function} {} cat (@var{dim}, @var{array1}, @var{array2}, @dots{}, @var{arrayN})\n\
Return the concatenation of N-d array objects, @var{array1},\n\
@var{array2}, @dots{}, @var{arrayN} along dimension @var{dim}.\n\
\n\
@example\n\
@group\n\
A = ones (2, 2);\n\
B = zeros (2, 2);\n\
cat (2, A, B)\n\
@result{} ans =\n\
\n\
     1 1 0 0\n\
     1 1 0 0\n\
@end group\n\
@end example\n\
\n\
Alternatively, we can concatenate @var{A} and @var{B} along the\n\
second dimension the following way:\n\
\n\
@example\n\
@group\n\
[A, B].\n\
@end group\n\
@end example\n\
\n\
@var{dim} can be larger than the dimensions of the N-d array objects\n\
and the result will thus have @var{dim} dimensions as the\n\
following example shows:\n\
@example\n\
@group\n\
cat (4, ones(2, 2), zeros (2, 2))\n\
@result{} ans =\n\
\n\
   ans(:,:,1,1) =\n\
\n\
     1 1\n\
     1 1\n\
\n\
   ans(:,:,1,2) =\n\
     0 0\n\
     0 0\n\
@end group\n\
@end example\n\
@seealso{horzcat, vertcat}\n\
@end deftypefn")
{
  return do_cat (args, "cat");
}

/*

%!function ret = testcat (t1, t2, tr, cmplx)
%! assert (cat (1, cast ([], t1), cast([], t2)), cast ([], tr));
%!
%! assert (cat (1, cast (1, t1), cast (2, t2)), cast ([1; 2], tr));
%! assert (cat (1, cast (1, t1), cast ([2; 3], t2)), cast ([1; 2; 3], tr));
%! assert (cat (1, cast ([1; 2], t1), cast (3, t2)), cast ([1; 2; 3], tr));
%! assert (cat (1, cast ([1; 2], t1), cast ([3; 4], t2)), cast ([1; 2; 3; 4], tr));
%! assert (cat (2, cast (1, t1), cast (2, t2)), cast ([1, 2], tr));
%! assert (cat (2, cast (1, t1), cast ([2, 3], t2)), cast ([1, 2, 3], tr));
%! assert (cat (2, cast ([1, 2], t1), cast (3, t2)), cast ([1, 2, 3], tr));
%! assert (cat (2, cast ([1, 2], t1), cast ([3, 4], t2)), cast ([1, 2, 3, 4], tr));
%! 
%! assert ([cast(1, t1); cast(2, t2)], cast ([1; 2], tr));
%! assert ([cast(1, t1); cast([2; 3], t2)], cast ([1; 2; 3], tr));
%! assert ([cast([1; 2], t1); cast(3, t2)], cast ([1; 2; 3], tr));
%! assert ([cast([1; 2], t1); cast([3; 4], t2)], cast ([1; 2; 3; 4], tr));
%! assert ([cast(1, t1), cast(2, t2)], cast ([1, 2], tr));
%! assert ([cast(1, t1), cast([2, 3], t2)], cast ([1, 2, 3], tr));
%! assert ([cast([1, 2], t1), cast(3, t2)], cast ([1, 2, 3], tr));
%! assert ([cast([1, 2], t1), cast([3, 4], t2)], cast ([1, 2, 3, 4], tr));
%!
%! if (nargin == 3 || cmplx)
%!   assert (cat (1, cast (1i, t1), cast (2, t2)), cast ([1i; 2], tr));
%!   assert (cat (1, cast (1i, t1), cast ([2; 3], t2)), cast ([1i; 2; 3], tr));
%!   assert (cat (1, cast ([1i; 2], t1), cast (3, t2)), cast ([1i; 2; 3], tr));
%!   assert (cat (1, cast ([1i; 2], t1), cast ([3; 4], t2)), cast ([1i; 2; 3; 4], tr));
%!   assert (cat (2, cast (1i, t1), cast (2, t2)), cast ([1i, 2], tr));
%!   assert (cat (2, cast (1i, t1), cast ([2, 3], t2)), cast ([1i, 2, 3], tr));
%!   assert (cat (2, cast ([1i, 2], t1), cast (3, t2)), cast ([1i, 2, 3], tr));
%!   assert (cat (2, cast ([1i, 2], t1), cast ([3, 4], t2)), cast ([1i, 2, 3, 4], tr));
%! 
%!   assert ([cast(1i, t1); cast(2, t2)], cast ([1i; 2], tr));
%!   assert ([cast(1i, t1); cast([2; 3], t2)], cast ([1i; 2; 3], tr));
%!   assert ([cast([1i; 2], t1); cast(3, t2)], cast ([1i; 2; 3], tr));
%!   assert ([cast([1i; 2], t1); cast([3; 4], t2)], cast ([1i; 2; 3; 4], tr));
%!   assert ([cast(1i, t1), cast(2, t2)], cast ([1i, 2], tr));
%!   assert ([cast(1i, t1), cast([2, 3], t2)], cast ([1i, 2, 3], tr));
%!   assert ([cast([1i, 2], t1), cast(3, t2)], cast ([1i, 2, 3], tr));
%!   assert ([cast([1i, 2], t1), cast([3, 4], t2)], cast ([1i, 2, 3, 4], tr));
%!
%!   assert (cat (1, cast (1, t1), cast (2i, t2)), cast ([1; 2i], tr));
%!   assert (cat (1, cast (1, t1), cast ([2i; 3], t2)), cast ([1; 2i; 3], tr));
%!   assert (cat (1, cast ([1; 2], t1), cast (3i, t2)), cast ([1; 2; 3i], tr));
%!   assert (cat (1, cast ([1; 2], t1), cast ([3i; 4], t2)), cast ([1; 2; 3i; 4], tr));
%!   assert (cat (2, cast (1, t1), cast (2i, t2)), cast ([1, 2i], tr));
%!   assert (cat (2, cast (1, t1), cast ([2i, 3], t2)), cast ([1, 2i, 3], tr));
%!   assert (cat (2, cast ([1, 2], t1), cast (3i, t2)), cast ([1, 2, 3i], tr));
%!   assert (cat (2, cast ([1, 2], t1), cast ([3i, 4], t2)), cast ([1, 2, 3i, 4], tr));
%! 
%!   assert ([cast(1, t1); cast(2i, t2)], cast ([1; 2i], tr));
%!   assert ([cast(1, t1); cast([2i; 3], t2)], cast ([1; 2i; 3], tr));
%!   assert ([cast([1; 2], t1); cast(3i, t2)], cast ([1; 2; 3i], tr));
%!   assert ([cast([1; 2], t1); cast([3i; 4], t2)], cast ([1; 2; 3i; 4], tr));
%!   assert ([cast(1, t1), cast(2i, t2)], cast ([1, 2i], tr));
%!   assert ([cast(1, t1), cast([2i, 3], t2)], cast ([1, 2i, 3], tr));
%!   assert ([cast([1, 2], t1), cast(3i, t2)], cast ([1, 2, 3i], tr));
%!   assert ([cast([1, 2], t1), cast([3i, 4], t2)], cast ([1, 2, 3i, 4], tr));
%!
%!   assert (cat (1, cast (1i, t1), cast (2i, t2)), cast ([1i; 2i], tr));
%!   assert (cat (1, cast (1i, t1), cast ([2i; 3], t2)), cast ([1i; 2i; 3], tr));
%!   assert (cat (1, cast ([1i; 2], t1), cast (3i, t2)), cast ([1i; 2; 3i], tr));
%!   assert (cat (1, cast ([1i; 2], t1), cast ([3i; 4], t2)), cast ([1i; 2; 3i; 4], tr));
%!   assert (cat (2, cast (1i, t1), cast (2i, t2)), cast ([1i, 2i], tr));
%!   assert (cat (2, cast (1i, t1), cast ([2i, 3], t2)), cast ([1i, 2i, 3], tr));
%!   assert (cat (2, cast ([1i, 2], t1), cast (3i, t2)), cast ([1i, 2, 3i], tr));
%!   assert (cat (2, cast ([1i, 2], t1), cast ([3i, 4], t2)), cast ([1i, 2, 3i, 4], tr));
%! 
%!   assert ([cast(1i, t1); cast(2i, t2)], cast ([1i; 2i], tr));
%!   assert ([cast(1i, t1); cast([2i; 3], t2)], cast ([1i; 2i; 3], tr));
%!   assert ([cast([1i; 2], t1); cast(3i, t2)], cast ([1i; 2; 3i], tr));
%!   assert ([cast([1i; 2], t1); cast([3i; 4], t2)], cast ([1i; 2; 3i; 4], tr));
%!   assert ([cast(1i, t1), cast(2i, t2)], cast ([1i, 2i], tr));
%!   assert ([cast(1i, t1), cast([2i, 3], t2)], cast ([1i, 2i, 3], tr));
%!   assert ([cast([1i, 2], t1), cast(3i, t2)], cast ([1i, 2, 3i], tr));
%!   assert ([cast([1i, 2], t1), cast([3i, 4], t2)], cast ([1i, 2, 3i, 4], tr));
%! endif
%! ret = true;

%!assert (testcat('double', 'double', 'double'));
%!assert (testcat('single', 'double', 'single'));
%!assert (testcat('double', 'single', 'single'));
%!assert (testcat('single', 'single', 'single'));

%!assert (testcat('double', 'int8', 'int8', false));
%!assert (testcat('int8', 'double', 'int8', false));
%!assert (testcat('single', 'int8', 'int8', false));
%!assert (testcat('int8', 'single', 'int8', false));
%!assert (testcat('int8', 'int8', 'int8', false));
%!assert (testcat('double', 'int16', 'int16', false));
%!assert (testcat('int16', 'double', 'int16', false));
%!assert (testcat('single', 'int16', 'int16', false));
%!assert (testcat('int16', 'single', 'int16', false));
%!assert (testcat('int16', 'int16', 'int16', false));
%!assert (testcat('double', 'int32', 'int32', false));
%!assert (testcat('int32', 'double', 'int32', false));
%!assert (testcat('single', 'int32', 'int32', false));
%!assert (testcat('int32', 'single', 'int32', false));
%!assert (testcat('int32', 'int32', 'int32', false));
%!assert (testcat('double', 'int64', 'int64', false));
%!assert (testcat('int64', 'double', 'int64', false));
%!assert (testcat('single', 'int64', 'int64', false));
%!assert (testcat('int64', 'single', 'int64', false));
%!assert (testcat('int64', 'int64', 'int64', false));

%!assert (testcat('double', 'uint8', 'uint8', false));
%!assert (testcat('uint8', 'double', 'uint8', false));
%!assert (testcat('single', 'uint8', 'uint8', false));
%!assert (testcat('uint8', 'single', 'uint8', false));
%!assert (testcat('uint8', 'uint8', 'uint8', false));
%!assert (testcat('double', 'uint16', 'uint16', false));
%!assert (testcat('uint16', 'double', 'uint16', false));
%!assert (testcat('single', 'uint16', 'uint16', false));
%!assert (testcat('uint16', 'single', 'uint16', false));
%!assert (testcat('uint16', 'uint16', 'uint16', false));
%!assert (testcat('double', 'uint32', 'uint32', false));
%!assert (testcat('uint32', 'double', 'uint32', false));
%!assert (testcat('single', 'uint32', 'uint32', false));
%!assert (testcat('uint32', 'single', 'uint32', false));
%!assert (testcat('uint32', 'uint32', 'uint32', false));
%!assert (testcat('double', 'uint64', 'uint64', false));
%!assert (testcat('uint64', 'double', 'uint64', false));
%!assert (testcat('single', 'uint64', 'uint64', false));
%!assert (testcat('uint64', 'single', 'uint64', false));
%!assert (testcat('uint64', 'uint64', 'uint64', false));

*/

static octave_value
do_permute (const octave_value_list& args, bool inv)
{
  octave_value retval;

  if (args.length () == 2 && args(1).length () >= args(1).ndims ())
    {
      Array<int> vec = args(1).int_vector_value ();

      // FIXME -- maybe we should create an idx_vector object
      // here and pass that to permute?

      int n = vec.length ();

      for (int i = 0; i < n; i++)
	vec(i)--;

      octave_value ret = args(0).permute (vec, inv);

      if (! error_state)
	retval = ret;
    }
  else
    print_usage ();

  return retval;
}

DEFUN (permute, args, ,
  "-*- texinfo -*-\n\
@deftypefn {Built-in Function} {} permute (@var{a}, @var{perm})\n\
Return the generalized transpose for an N-d array object @var{a}.\n\
The permutation vector @var{perm} must contain the elements\n\
@code{1:ndims(a)} (in any order, but each element must appear just once).\n\
@seealso{ipermute}\n\
@end deftypefn")
{
  return do_permute (args, false);
}

DEFUN (ipermute, args, ,
  "-*- texinfo -*-\n\
@deftypefn {Built-in Function} {} ipermute (@var{a}, @var{iperm})\n\
The inverse of the @code{permute} function.  The expression\n\
\n\
@example\n\
ipermute (permute (a, perm), perm)\n\
@end example\n\
returns the original array @var{a}.\n\
@seealso{permute}\n\
@end deftypefn")
{
  return do_permute (args, true);
}

DEFUN (length, args, ,
  "-*- texinfo -*-\n\
@deftypefn {Built-in Function} {} length (@var{a})\n\
Return the `length' of the object @var{a}.  For matrix objects, the\n\
length is the number of rows or columns, whichever is greater (this\n\
odd definition is used for compatibility with @sc{matlab}).\n\
@end deftypefn")
{
  octave_value retval;

  if (args.length () == 1)
    retval = args(0).length ();
  else
    print_usage ();

  return retval;
}

DEFUN (ndims, args, ,
  "-*- texinfo -*-\n\
@deftypefn {Built-in Function} {} ndims (@var{a})\n\
Returns the number of dimensions of array @var{a}.\n\
For any array, the result will always be larger than or equal to 2.\n\
Trailing singleton dimensions are not counted.\n\
@end deftypefn")
{
  octave_value retval;

  if (args.length () == 1)
    retval = args(0).ndims ();
  else
    print_usage ();

  return retval;
}

DEFUN (numel, args, ,
  "-*- texinfo -*-\n\
@deftypefn {Built-in Function} {} numel (@var{a})\n\
@deftypefnx {Built-in Function} {} numel (@var{a}, @var{idx1}, @var{idx2}, @dots{})\n\
Returns the number of elements in the object @var{a}.\n\
Optionally, if indices @var{idx1}, @var{idx2}, @dots{} are supplied,\n\
return the number of elements that would result from the indexing\n\
@example\n\
  @var{a}(@var{idx1}, @var{idx2}, @dots{})\n\
@end example\n\
This method is also called when an object appears as lvalue with cs-list\n\
indexing, i.e., @code{object@{@dots{}@}} or @code{object(@dots{}).field}.\n\
@seealso{size}\n\
@end deftypefn")
{
  octave_value retval;
  octave_idx_type nargin = args.length ();

  if (nargin == 1)
    retval = args(0).numel ();
  else if (nargin > 1)
    {
      // Don't use numel (const octave_value_list&) here as that corresponds to
      // an overloaded call, not to builtin!
      retval = dims_to_numel (args(0).dims (), args.slice (1, nargin-1));
    }
  else
    print_usage ();

  return retval;
}

DEFUN (size, args, nargout,
  "-*- texinfo -*-\n\
@deftypefn {Built-in Function} {} size (@var{a}, @var{n})\n\
Return the number rows and columns of @var{a}.\n\
\n\
With one input argument and one output argument, the result is returned\n\
in a row vector.  If there are multiple output arguments, the number of\n\
rows is assigned to the first, and the number of columns to the second,\n\
etc.  For example,\n\
\n\
@example\n\
@group\n\
size ([1, 2; 3, 4; 5, 6])\n\
     @result{} [ 3, 2 ]\n\
\n\
[nr, nc] = size ([1, 2; 3, 4; 5, 6])\n\
     @result{} nr = 3\n\
     @result{} nc = 2\n\
@end group\n\
@end example\n\
\n\
If given a second argument, @code{size} will return the size of the\n\
corresponding dimension.  For example\n\
\n\
@example\n\
@group\n\
size ([1, 2; 3, 4; 5, 6], 2)\n\
     @result{} 2\n\
@end group\n\
@end example\n\
\n\
@noindent\n\
returns the number of columns in the given matrix.\n\
@seealso{numel}\n\
@end deftypefn")
{
  octave_value_list retval;

  int nargin = args.length ();

  if (nargin == 1)
    {
      const dim_vector dimensions = args(0).dims ();

      if (nargout > 1)
	{
          const dim_vector rdims = dimensions.redim (nargout);
          retval.resize (nargout);
          for (int i = 0; i < nargout; i++)
            retval(i) = rdims(i);
	}
      else
	{
          int ndims = dimensions.length ();

          NoAlias<Matrix> m (1, ndims);

	  for (int i = 0; i < ndims; i++)
	    m(i) = dimensions(i);

	  retval(0) = m;
	}
    }
  else if (nargin == 2 && nargout < 2)
    {
      octave_idx_type nd = args(1).int_value (true);

      if (error_state)
	error ("size: expecting scalar as second argument");
      else
	{
	  const dim_vector dv = args(0).dims ();

	  if (nd > 0)
	    {
	      if (nd <= dv.length ())
		retval(0) = dv(nd-1);
	      else 
		retval(0) = 1;
	    }
	  else
	    error ("size: requested dimension (= %d) out of range", nd);
	}
    }
  else
    print_usage ();

  return retval;
}

DEFUN (size_equal, args, ,
   "-*- texinfo -*-\n\
@deftypefn {Built-in Function} {} size_equal (@var{a}, @var{b}, @dots{})\n\
Return true if the dimensions of all arguments agree.\n\
Trailing singleton dimensions are ignored.\n\
Called with a single argument, size_equal returns true.\n\
@seealso{size, numel}\n\
@end deftypefn")
{
  octave_value retval;

  int nargin = args.length ();

  if (nargin >= 1)
    {
      retval = true;

      dim_vector a_dims = args(0).dims ();

      for (int i = 1; i < nargin; ++i)
        {
          dim_vector b_dims = args(i).dims ();

          if (a_dims != b_dims)
	    {
	      retval = false;
	      break;
	    }
        }
    }
  else
    print_usage ();

  return retval;
}

DEFUN (nnz, args, ,
   "-*- texinfo -*-\n\
@deftypefn {Built-in Function} {@var{scalar} =} nnz (@var{a})\n\
Returns the number of non zero elements in @var{a}.\n\
@seealso{sparse}\n\
@end deftypefn")
{
  octave_value retval;

  if (args.length () == 1)
    retval = args(0).nnz ();
  else
    print_usage ();

  return retval;
}

DEFUN (nzmax, args, ,
   "-*- texinfo -*-\n\
@deftypefn {Built-in Function} {@var{scalar} =} nzmax (@var{SM})\n\
Return the amount of storage allocated to the sparse matrix @var{SM}.\n\
Note that Octave tends to crop unused memory at the first opportunity\n\
for sparse objects.  There are some cases of user created sparse objects\n\
where the value returned by @dfn{nzmax} will not be the same as @dfn{nnz},\n\
but in general they will give the same result.\n\
@seealso{sparse, spalloc}\n\
@end deftypefn")
{
  octave_value retval;

  if (args.length() == 1)
    retval = args(0).nzmax ();
  else
    print_usage ();

  return retval;
}

DEFUN (rows, args, ,
  "-*- texinfo -*-\n\
@deftypefn {Built-in Function} {} rows (@var{a})\n\
Return the number of rows of @var{a}.\n\
@seealso{size, numel, columns, length, isscalar, isvector, ismatrix}\n\
@end deftypefn")
{
  octave_value retval;

  if (args.length () == 1)
    retval = args(0).rows ();
  else
    print_usage ();

  return retval;
}

DEFUN (columns, args, ,
  "-*- texinfo -*-\n\
@deftypefn {Built-in Function} {} columns (@var{a})\n\
Return the number of columns of @var{a}.\n\
@seealso{size, numel, rows, length, isscalar, isvector, ismatrix}\n\
@end deftypefn")
{
  octave_value retval;

  if (args.length () == 1)
    retval = args(0).columns ();
  else
    print_usage ();

  return retval;
}

DEFUN (sum, args, ,
  "-*- texinfo -*-\n\
@deftypefn  {Built-in Function} {} sum (@var{x})\n\
@deftypefnx {Built-in Function} {} sum (@var{x}, @var{dim})\n\
@deftypefnx {Built-in Function} {} sum (@dots{}, 'native')\n\
@deftypefnx {Built-in Function} {} sum (@dots{}, 'double')\n\
@deftypefnx {Built-in Function} {} sum (@dots{}, 'extra')\n\
Sum of elements along dimension @var{dim}.  If @var{dim} is\n\
omitted, it defaults to the first non-singleton dimension.\n\
\n\
If the optional argument 'native' is given, then the sum is performed\n\
in the same type as the original argument, rather than in the default\n\
double type.  For example\n\
\n\
@example\n\
@group\n\
sum ([true, true])\n\
  @result{} 2\n\
sum ([true, true], 'native')\n\
  @result{} true\n\
@end group\n\
@end example\n\
  \n\
On the contrary, if 'double' is given, the sum is performed in double precision\n\
even for single precision inputs.\n\
\n\
For double precision inputs, 'extra' indicates that a more accurate algorithm\n\
than straightforward summation is to be used.  For single precision inputs, 'extra' is\n\
the same as 'double'.  Otherwise, 'extra' has no effect.\n\
@seealso{cumsum, sumsq, prod}\n\
@end deftypefn")
{
  octave_value retval;

  int nargin = args.length ();

  bool isnative = false;
  bool isdouble = false;
  bool isextra = false;

  if (nargin > 1 && args(nargin - 1).is_string ())
    {
      std::string str = args(nargin - 1).string_value ();

      if (! error_state)
	{
	  if (str == "native")
	    isnative = true;
	  else if (str == "double")
            isdouble = true;
          else if (str == "extra")
            isextra = true;
          else
	    error ("sum: unrecognized string argument");
          nargin --;
	}
    }

  if (error_state)
    return retval;

  if (nargin == 1 || nargin == 2)
    {
      octave_value arg = args(0);

      int dim = -1;
      if (nargin == 2)
        {
          dim = args(1).int_value () - 1;
          if (dim < 0)
	    error ("sum: invalid dimension argument = %d", dim + 1);
        }

      if (! error_state)
	{
          switch (arg.builtin_type ())
            {
            case btyp_double:
              if (arg.is_sparse_type ())
                {
                  if (isextra)
                    warning ("sum: 'extra' not yet implemented for sparse matrices");
                  retval = arg.sparse_matrix_value ().sum (dim);
                }
              else if (isextra)
                retval = arg.array_value ().xsum (dim);
              else
                retval = arg.array_value ().sum (dim);
              break;
            case btyp_complex:
              if (arg.is_sparse_type ())
                {
                  if (isextra)
                    warning ("sum: 'extra' not yet implemented for sparse matrices");
                  retval = arg.sparse_complex_matrix_value ().sum (dim);
                }
              else if (isextra)
                retval = arg.complex_array_value ().xsum (dim);
              else
                retval = arg.complex_array_value ().sum (dim);
              break;
            case btyp_float:
              if (isdouble || isextra)
                retval = arg.float_array_value ().dsum (dim);
              else
                retval = arg.float_array_value ().sum (dim);
              break;
            case btyp_float_complex:
              if (isdouble || isextra)
                retval = arg.float_complex_array_value ().dsum (dim);
              else
                retval = arg.float_complex_array_value ().sum (dim);
              break;

#define MAKE_INT_BRANCH(X) \
            case btyp_ ## X: \
              if (isnative) \
                retval = arg.X ## _array_value ().sum (dim); \
              else \
                retval = arg.X ## _array_value ().dsum (dim); \
              break
            MAKE_INT_BRANCH (int8);
            MAKE_INT_BRANCH (int16);
            MAKE_INT_BRANCH (int32);
            MAKE_INT_BRANCH (int64);
            MAKE_INT_BRANCH (uint8);
            MAKE_INT_BRANCH (uint16);
            MAKE_INT_BRANCH (uint32);
            MAKE_INT_BRANCH (uint64);
#undef MAKE_INT_BRANCH

            case btyp_bool:
              if (arg.is_sparse_type ())
                {
                  if (isnative)
                    retval = arg.sparse_bool_matrix_value ().any (dim);
                  else
                    retval = arg.sparse_matrix_value ().sum (dim);
                }
              else if (isnative)
                retval = arg.bool_array_value ().any (dim);
              else
                retval = arg.bool_array_value ().sum (dim);
              break;

            default:
              gripe_wrong_type_arg ("sum", arg);
            }
	}
    }
  else
    print_usage ();

  return retval;
}

/*

%!assert (sum([true,true]), 2)
%!assert (sum([true,true],'native'), true)
%!assert (sum(int8([127,10,-20])), 117);
%!assert (sum(int8([127,10,-20]),'native'), int8(107));

%!assert(sum ([1, 2, 3]), 6)
%!assert(sum ([-1; -2; -3]), -6);
%!assert(sum ([i, 2+i, -3+2i, 4]), 3+4i);
%!assert(sum ([1, 2, 3; i, 2i, 3i; 1+i, 2+2i, 3+3i]), [2+2i, 4+4i, 6+6i]);

%!assert(sum (single([1, 2, 3])), single(6))
%!assert(sum (single([-1; -2; -3])), single(-6));
%!assert(sum (single([i, 2+i, -3+2i, 4])), single(3+4i));
%!assert(sum (single([1, 2, 3; i, 2i, 3i; 1+i, 2+2i, 3+3i])), single([2+2i, 4+4i, 6+6i]));

%!error <Invalid call to sum.*> sum ();

%!assert (sum ([1, 2; 3, 4], 1), [4, 6]);
%!assert (sum ([1, 2; 3, 4], 2), [3; 7]);
%!assert (sum (zeros (1, 0)), 0);
%!assert (sum (zeros (1, 0), 1), zeros(1, 0));
%!assert (sum (zeros (1, 0), 2), 0);
%!assert (sum (zeros (0, 1)), 0);
%!assert (sum (zeros (0, 1), 1), 0);
%!assert (sum (zeros (0, 1), 2), zeros(0, 1));
%!assert (sum (zeros (2, 0)),  zeros(1, 0));
%!assert (sum (zeros (2, 0), 1), zeros(1, 0));
%!assert (sum (zeros (2, 0), 2),  [0; 0]);
%!assert (sum (zeros (0, 2)), [0, 0]);
%!assert (sum (zeros (0, 2), 1), [0, 0]);
%!assert (sum (zeros (0, 2), 2), zeros(0, 1));
%!assert (sum (zeros (2, 2, 0, 3)), zeros(1, 2, 0, 3));
%!assert (sum (zeros (2, 2, 0, 3), 2), zeros(2, 1, 0, 3));
%!assert (sum (zeros (2, 2, 0, 3), 3), zeros(2, 2, 1, 3));
%!assert (sum (zeros (2, 2, 0, 3), 4), zeros(2, 2, 0));
%!assert (sum (zeros (2, 2, 0, 3), 7), zeros(2, 2, 0, 3));

%!assert (sum (single([1, 2; 3, 4]), 1), single([4, 6]));
%!assert (sum (single([1, 2; 3, 4]), 2), single([3; 7]));
%!assert (sum (zeros (1, 0, 'single')), single(0));
%!assert (sum (zeros (1, 0, 'single'), 1), zeros(1, 0, 'single'));
%!assert (sum (zeros (1, 0, 'single'), 2), single(0));
%!assert (sum (zeros (0, 1, 'single')), single(0));
%!assert (sum (zeros (0, 1, 'single'), 1), single(0));
%!assert (sum (zeros (0, 1, 'single'), 2), zeros(0, 1, 'single'));
%!assert (sum (zeros (2, 0, 'single')),  zeros(1, 0, 'single'));
%!assert (sum (zeros (2, 0, 'single'), 1), zeros(1, 0, 'single'));
%!assert (sum (zeros (2, 0, 'single'), 2),  single([0; 0]));
%!assert (sum (zeros (0, 2, 'single')), single([0, 0]));
%!assert (sum (zeros (0, 2, 'single'), 1), single([0, 0]));
%!assert (sum (zeros (0, 2, 'single'), 2), zeros(0, 1, 'single'));
%!assert (sum (zeros (2, 2, 0, 3, 'single')), zeros(1, 2, 0, 3, 'single'));
%!assert (sum (zeros (2, 2, 0, 3, 'single'), 2), zeros(2, 1, 0, 3, 'single'));
%!assert (sum (zeros (2, 2, 0, 3, 'single'), 3), zeros(2, 2, 1, 3, 'single'));
%!assert (sum (zeros (2, 2, 0, 3, 'single'), 4), zeros(2, 2, 0, 'single'));
%!assert (sum (zeros (2, 2, 0, 3, 'single'), 7), zeros(2, 2, 0, 3, 'single'));

*/

DEFUN (sumsq, args, ,
  "-*- texinfo -*-\n\
@deftypefn  {Built-in Function} {} sumsq (@var{x})\n\
@deftypefnx {Built-in Function} {} sumsq (@var{x}, @var{dim})\n\
Sum of squares of elements along dimension @var{dim}.  If @var{dim}\n\
is omitted, it defaults to the first non-singleton dimension.\n\
\n\
This function is conceptually equivalent to computing\n\
@example\n\
sum (x .* conj (x), dim)\n\
@end example\n\
but it uses less memory and avoids calling @code{conj} if @var{x} is real.\n\
@seealso{sum}\n\
@end deftypefn")
{
  DATA_REDUCTION (sumsq);
}

/*

%!assert(sumsq ([1, 2, 3]), 14)
%!assert(sumsq ([-1; -2; 4i]), 21);
%!assert(sumsq ([1, 2, 3; 2, 3, 4; 4i, 6i, 2]), [21, 49, 29]);

%!assert(sumsq (single([1, 2, 3])), single(14))
%!assert(sumsq (single([-1; -2; 4i])), single(21));
%!assert(sumsq (single([1, 2, 3; 2, 3, 4; 4i, 6i, 2])), single([21, 49, 29]));

%!error <Invalid call to sumsq.*> sumsq ();

%!assert (sumsq ([1, 2; 3, 4], 1), [10, 20]);
%!assert (sumsq ([1, 2; 3, 4], 2), [5; 25]);

%!assert (sumsq (single([1, 2; 3, 4]), 1), single([10, 20]));
%!assert (sumsq (single([1, 2; 3, 4]), 2), single([5; 25]));

 */

DEFUN (islogical, args, ,
  "-*- texinfo -*-\n\
@deftypefn {Built-in Function} {} islogical (@var{x})\n\
Return true if @var{x} is a logical object.\n\
@end deftypefn")
{
  octave_value retval;

  if (args.length () == 1)
    retval = args(0).is_bool_type ();
  else
    print_usage ();

  return retval;
}

DEFALIAS (isbool, islogical);

/*

%!assert (islogical(true), true)
%!assert (islogical(false), true)
%!assert (islogical([true, false]), true)
%!assert (islogical(1), false)
%!assert (islogical(1i), false)
%!assert (islogical([1,1]), false)
%!assert (islogical(single(1)), false)
%!assert (islogical(single(1i)), false)
%!assert (islogical(single([1,1])), false)

 */

DEFUN (isinteger, args, ,
  "-*- texinfo -*-\n\
@deftypefn {Built-in Function} {} isinteger (@var{x})\n\
Return true if @var{x} is an integer object (int8, uint8, int16, etc.).\n\
Note that @code{isinteger (14)} is false because numeric constants in\n\
Octave are double precision floating point values.\n\
@seealso{isreal, isnumeric, class, isa}\n\
@end deftypefn")
{
  octave_value retval;

  if (args.length () == 1)
    retval = args(0).is_integer_type ();
  else
    print_usage ();

  return retval;
}

DEFUN (iscomplex, args, ,
  "-*- texinfo -*-\n\
@deftypefn {Built-in Function} {} iscomplex (@var{x})\n\
Return true if @var{x} is a complex-valued numeric object.\n\
@end deftypefn")
{
  octave_value retval;

  if (args.length () == 1)
    retval = args(0).is_complex_type ();
  else
    print_usage ();

  return retval;
}

DEFUN (isfloat, args, ,
  "-*- texinfo -*-\n\
@deftypefn {Built-in Function} {} isfloat (@var{x})\n\
Return true if @var{x} is a floating-point numeric object.\n\
@end deftypefn")
{
  octave_value retval;

  if (args.length () == 1)
    retval = args(0).is_float_type ();
  else
    print_usage ();

  return retval;
}

// FIXME -- perhaps this should be implemented with an
// octave_value member function?

DEFUN (complex, args, ,
  "-*- texinfo -*-\n\
@deftypefn  {Built-in Function} {} complex (@var{x})\n\
@deftypefnx {Built-in Function} {} complex (@var{re}, @var{im})\n\
Return a complex result from real arguments.  With 1 real argument @var{x},\n\
return the complex result @code{@var{x} + 0i}.  With 2 real arguments,\n\
return the complex result @code{@var{re} + @var{im}}.  @code{complex} can\n\
often be more convenient than expressions such as @code{a + i*b}.\n\
For example:\n\
\n\
@example\n\
@group\n\
complex ([1, 2], [3, 4])\n\
@result{}\n\
   1 + 3i   2 + 4i\n\
@end group\n\
@end example\n\
@seealso{real, imag, iscomplex}\n\
@end deftypefn")
{
  octave_value retval;

  int nargin = args.length ();

  if (nargin == 1)
    {
      octave_value arg = args(0);

      if (arg.is_complex_type ())
	retval = arg;
      else
	{
	  if (arg.is_sparse_type ())
	    {
	      SparseComplexMatrix val = arg.sparse_complex_matrix_value ();

	      if (! error_state)
		retval = octave_value (new octave_sparse_complex_matrix (val));
	    }
	  else if (arg.is_single_type ())
	    {
	      if (arg.numel () == 1)
		{
		  FloatComplex val = arg.float_complex_value ();

		  if (! error_state)
		    retval = octave_value (new octave_float_complex (val));
		}
	      else
		{
		  FloatComplexNDArray val = arg.float_complex_array_value ();

		  if (! error_state)
		    retval = octave_value (new octave_float_complex_matrix (val));
		}
	    }
	  else
	    {
	      if (arg.numel () == 1)
		{
		  Complex val = arg.complex_value ();

		  if (! error_state)
		    retval = octave_value (new octave_complex (val));
		}
	      else
		{
		  ComplexNDArray val = arg.complex_array_value ();

		  if (! error_state)
		    retval = octave_value (new octave_complex_matrix (val));
		}
	    }

	  if (error_state)
	    error ("complex: invalid conversion");
	}
    }
  else if (nargin == 2)
    {
      octave_value re = args(0);
      octave_value im = args(1);

      if (re.is_sparse_type () && im.is_sparse_type ())
	{
	  const SparseMatrix re_val = re.sparse_matrix_value ();
	  const SparseMatrix im_val = im.sparse_matrix_value ();

	  if (!error_state)
	    {
	      if (re.numel () == 1)
		{
		  SparseComplexMatrix result;
		  if (re_val.nnz () == 0)
		    result = Complex(0, 1) * SparseComplexMatrix (im_val);
		  else
		    {
		      result = SparseComplexMatrix (im_val.dims (), re_val (0));
		      octave_idx_type nr = im_val.rows ();
		      octave_idx_type nc = im_val.cols ();

		      for (octave_idx_type j = 0; j < nc; j++)
			{
			  octave_idx_type off = j * nr;
			  for (octave_idx_type i = im_val.cidx(j); 
			       i < im_val.cidx(j + 1); i++)
			    result.data (im_val.ridx(i) + off) =  
			      result.data (im_val.ridx(i) + off) + 
			      Complex (0, im_val.data (i));
			}
		    }
		  retval = octave_value (new octave_sparse_complex_matrix (result));
		}
	      else if (im.numel () == 1)
		{
		  SparseComplexMatrix result;
		  if (im_val.nnz () == 0)
		    result = SparseComplexMatrix (re_val);
		  else
		    {
		      result = SparseComplexMatrix (re_val.rows(), re_val.cols(), Complex(0, im_val (0)));
		      octave_idx_type nr = re_val.rows ();
		      octave_idx_type nc = re_val.cols ();

		      for (octave_idx_type j = 0; j < nc; j++)
			{
			  octave_idx_type off = j * nr;
			  for (octave_idx_type i = re_val.cidx(j); 
			       i < re_val.cidx(j + 1); i++)
			    result.data (re_val.ridx(i) + off) =  
			      result.data (re_val.ridx(i) + off) + 
			      re_val.data (i);
			}
		    }
		  retval = octave_value (new octave_sparse_complex_matrix (result));
		}
	      else
		{
		  if (re_val.dims () == im_val.dims ())
		    {
		      SparseComplexMatrix result = SparseComplexMatrix(re_val) 
			+ Complex(0, 1) * SparseComplexMatrix (im_val);
		      retval = octave_value (new octave_sparse_complex_matrix (result));
		    }
		  else
		    error ("complex: dimension mismatch");
		}
	    }
	}
      else if (re.is_single_type () || im.is_single_type ())
	{
	  if (re.numel () == 1)
	    {
	      float re_val = re.float_value ();

	      if (im.numel () == 1)
		{
		  float im_val = im.double_value ();

		  if (! error_state)
		    retval = octave_value (new octave_float_complex (FloatComplex (re_val, im_val)));
		}
	      else
		{
		  const FloatNDArray im_val = im.float_array_value ();

		  if (! error_state)
		    {
		      FloatComplexNDArray result (im_val.dims (), FloatComplex ());

		      for (octave_idx_type i = 0; i < im_val.numel (); i++)
			result.xelem (i) = FloatComplex (re_val, im_val(i));

		      retval = octave_value (new octave_float_complex_matrix (result));
		    }
		}
	    }
	  else
	    {
	      const FloatNDArray re_val = re.float_array_value ();

	      if (im.numel () == 1)
		{
		  float im_val = im.float_value ();

		  if (! error_state)
		    {
		      FloatComplexNDArray result (re_val.dims (), FloatComplex ());

		      for (octave_idx_type i = 0; i < re_val.numel (); i++)
			result.xelem (i) = FloatComplex (re_val(i), im_val);

		      retval = octave_value (new octave_float_complex_matrix (result));
		    }
		}
	      else
		{
		  const FloatNDArray im_val = im.float_array_value ();

		  if (! error_state)
		    {
		      if (re_val.dims () == im_val.dims ())
			{
			  FloatComplexNDArray result (re_val.dims (), FloatComplex ());

			  for (octave_idx_type i = 0; i < re_val.numel (); i++)
			    result.xelem (i) = FloatComplex (re_val(i), im_val(i));

			  retval = octave_value (new octave_float_complex_matrix (result));
			}
		      else
			error ("complex: dimension mismatch");
		    }
		}
	    }
	}
      else if (re.numel () == 1)
	{
	  double re_val = re.double_value ();

	  if (im.numel () == 1)
	    {
	      double im_val = im.double_value ();

	      if (! error_state)
		retval = octave_value (new octave_complex (Complex (re_val, im_val)));
	    }
	  else
	    {
	      const NDArray im_val = im.array_value ();

	      if (! error_state)
		{
		  ComplexNDArray result (im_val.dims (), Complex ());

		  for (octave_idx_type i = 0; i < im_val.numel (); i++)
		    result.xelem (i) = Complex (re_val, im_val(i));

		  retval = octave_value (new octave_complex_matrix (result));
		}
	    }
	}
      else
	{
	  const NDArray re_val = re.array_value ();

	  if (im.numel () == 1)
	    {
	      double im_val = im.double_value ();

	      if (! error_state)
		{
		  ComplexNDArray result (re_val.dims (), Complex ());

		  for (octave_idx_type i = 0; i < re_val.numel (); i++)
		    result.xelem (i) = Complex (re_val(i), im_val);

		  retval = octave_value (new octave_complex_matrix (result));
		}
	    }
	  else
	    {
	      const NDArray im_val = im.array_value ();

	      if (! error_state)
		{
		  if (re_val.dims () == im_val.dims ())
		    {
		      ComplexNDArray result (re_val.dims (), Complex ());

		      for (octave_idx_type i = 0; i < re_val.numel (); i++)
			result.xelem (i) = Complex (re_val(i), im_val(i));

		      retval = octave_value (new octave_complex_matrix (result));
		    }
		  else
		    error ("complex: dimension mismatch");
		}
	    }
	}

      if (error_state)
	error ("complex: invalid conversion");
    }
  else
    print_usage ();

  return retval;
}

DEFUN (isreal, args, ,
  "-*- texinfo -*-\n\
@deftypefn {Built-in Function} {} isreal (@var{x})\n\
Return true if @var{x} is a real-valued numeric object.\n\
@end deftypefn")
{
  octave_value retval;

  if (args.length () == 1)
    retval = args(0).is_real_type ();
  else
    print_usage ();

  return retval;
}

DEFUN (isempty, args, ,
  "-*- texinfo -*-\n\
@deftypefn {Built-in Function} {} isempty (@var{a})\n\
Return 1 if @var{a} is an empty matrix (either the number of rows, or\n\
the number of columns, or both are zero).  Otherwise, return 0.\n\
@end deftypefn")
{
  octave_value retval = false;

  if (args.length () == 1)
    retval = args(0).is_empty ();
  else
    print_usage ();

  return retval;
}

DEFUN (isnumeric, args, ,
  "-*- texinfo -*-\n\
@deftypefn {Built-in Function} {} isnumeric (@var{x})\n\
Return nonzero if @var{x} is a numeric object.\n\
@end deftypefn")
{
  octave_value retval;

  if (args.length () == 1)
    retval = args(0).is_numeric_type ();
  else
    print_usage ();

  return retval;
}

DEFUN (islist, args, ,
  "-*- texinfo -*-\n\
@deftypefn {Built-in Function} {} islist (@var{x})\n\
Return nonzero if @var{x} is a list.\n\
@end deftypefn")
{
  octave_value retval;

  if (args.length () == 1)
    retval = args(0).is_list ();
  else
    print_usage ();

  return retval;
}

DEFUN (ismatrix, args, ,
  "-*- texinfo -*-\n\
@deftypefn {Built-in Function} {} ismatrix (@var{a})\n\
Return 1 if @var{a} is a matrix.  Otherwise, return 0.\n\
@end deftypefn")
{
  octave_value retval = false;

  if (args.length () == 1)
    {
      octave_value arg = args(0);

      if (arg.is_scalar_type () || arg.is_range ())
	retval = true;
      else if (arg.is_matrix_type ())
	retval = (arg.rows () >= 1 && arg.columns () >= 1);
    }
  else
    print_usage ();

  return retval;
}

/*