Mercurial > octave-nkf

changeset 6869:f9c893831e68

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2007-09-06 16:38:44 by dbateman]
author dbateman
date Thu, 06 Sep 2007 16:38:44 +0000
parents 975fcdfb0d2d
children e2e5425905ac
files scripts/general/celldisp.m scripts/miscellaneous/gzip.m scripts/miscellaneous/swapbytes.m src/DLD-FUNCTIONS/bsxfun.cc src/DLD-FUNCTIONS/typecast.cc
diffstat 5 files changed, 924 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/general/celldisp.m	Thu Sep 06 16:38:44 2007 +0000
@@ -0,0 +1,63 @@
+## Copyright (C) 2007  David Bateman
+##
+## 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 2, 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, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} celldisp (@var{c}, @var{name})
+## Recursively display the contents of a cell array. By default the values
+## are displayed with the name of the variable @var{c}. However, this name
+## can be replaced with the variable @var{name}.
+## @seealso{disp}
+## @end deftypefn
+
+## This is ugly, but seems to be what matlab does..
+
+function celldisp (c, name)
+  if (nargin < 1 || nargin > 2)
+    print_usage ();
+  endif
+
+  if (! iscell (c))
+    error ("celldisp: argument must be a cell array");
+  endif
+
+  if (nargin == 1)
+    name = inputname (1);
+  endif
+
+  for i = 1: numel (c)
+    if (iscell (c{i}))
+      celldisp (c{i}, sprintf ("%s{%s}", name, indices (size (c), i)));
+    else
+      disp (sprintf ("%s{%s} = \n", name, indices (size (c), i)));
+      disp (c{i});
+      disp ("");
+    endif
+  endfor
+endfunction
+
+function s = indices (dv, i)
+  if (sum (dv != 1) > 1)
+    c = cell (size (dv));
+    [c{:}] = ind2sub (dv, i);
+    s = sprintf("%i,", c{:});
+    s(end) = [];
+  else
+    s = sprintf("%i", i);
+  endif
+endfunction
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/miscellaneous/gzip.m	Thu Sep 06 16:38:44 2007 +0000
@@ -0,0 +1,110 @@
+## Copyright (C) 2007  David Bateman
+##
+## 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 2, 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, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{entries} =} gzip (@var{files})
+## @deftypefnx {Function File} {@var{entries} =} gzip (@var{files}, @var{outdir})
+## Compress the list of files and/or directories specified in @var{files}.
+## Each file is compressed separately and a new file with a '.gz' extension
+## is create. The original file is not touch. If @var{rootdir} is defined 
+## the compressed versions of the files are placed in this directory.
+## @seealso{gunzip, zip, tar}
+## @end deftypefn
+
+function entries = gzip (files, outdir)
+
+  if (nargin == 1 || nargin == 2)
+
+    if (nargin == 2 && ! exist (outdir, "dir"))
+      error ("gzip: output directory does not exist");
+    endif
+
+    if (ischar (files))
+      files = cellstr (files);
+    endif
+
+    if (nargin == 1)
+      outdir = tmpnam ();
+      mkdir (outdir);
+    endif
+
+    cwd = pwd();
+    unwind_protect
+      if (iscellstr (files))
+	files = glob (files);
+
+	## Ignore any file with a .gz extension
+	files (cellfun (@(x) strcmp (x(end-2:end), ".gz"), files)) = [];
+    
+	copyfile (files, outdir);
+	[d, f] = myfileparts(files);
+	cd (outdir);
+
+	cmd = sprintf ("gzip -r %s", sprintf (" %s", f{:}));
+
+	[status, output] = system (cmd);
+
+	if (status == 0)
+
+	  if (nargin == 2)
+	    gzfiles = cellfun(@(x) fullfile (outdir, sprintf ("%s.gz", x)), ...
+			      f, "UniformOutput", false);
+	  else
+	    movefile (cellfun(@(x) sprintf ("%s.gz", x), f, ...
+			      "UniformOutput", false), cwd);
+	    gzfiles = cellfun(@(x) sprintf ("%s.gz", x), ...
+			      files, "UniformOutput", false);
+	  endif
+
+	  if (nargout > 0)
+            entries = gzfiles;
+	  endif
+	else
+	  error ("gzip: failed with exit status = %d", status);
+	endif
+    
+      else
+	error ("gzip: expecting all arguments to be character strings");
+      endif
+    unwind_protect_cleanup
+      cd(cwd);
+      if (nargin == 1)
+	crr = confirm_recursive_rmdir ();
+	unwind_protect
+	  confirm_recursive_rmdir (false);
+	  rmdir (outdir, "s");
+	unwind_protect_cleanup
+	  confirm_recursive_rmdir (crr);
+	end_unwind_protect
+      endif
+    end_unwind_protect
+  else
+    print_usage ();
+  endif
+
+endfunction
+
+function [d, f] = myfileparts (x)
+  [d, f, ext] = cellfun (@(x) fileparts (x), x, "UniformOutput", false);
+  f = cellfun (@(x, y) sprintf ("%s%s", x, y), f, ext, ...
+	       "UniformOutput", false); 
+  idx = cellfun (@(x) isdir (x), x);
+  d(idx) = "";
+  f(idx) = x(idx);
+endfunction
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/miscellaneous/swapbytes.m	Thu Sep 06 16:38:44 2007 +0000
@@ -0,0 +1,57 @@
+## Copyright (C) 2007  David Bateman
+##
+## 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 2, 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, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} swapbytes (@var{x})
+## Swaps the byte order on values, converting from little endian to big 
+## endian and visa-versa. For example
+##
+## @example
+## @group
+## swapbytes (uint16 (1:4))
+## @result{} [   256   512   768  1024]
+## @end group
+## @end example
+##
+## @seealso{typecast, cast}
+## @end deftypefn
+
+function y = swapbytes (x)
+  if (nargin != 1)
+    print_usage ();
+  endif
+
+  clx = class (x);
+  if (strcmp (clx, "int8") || strcmp (clx, "uint8") || isempty (x))
+    y = x;
+  else
+    if (strcmp (clx, "int16") || strcmp (clx, "uint16"))
+      nb = 2;
+    elseif (strcmp (clx, "int32") || strcmp (clx, "uint32"))
+      nb = 4;
+    elseif (strcmp (clx, "int64") || strcmp (clx, "uint64") ||
+	    strcmp (clx, "double"))
+      nb = 8;
+    else
+      error ("swapbytes: invalid class of object");
+    endif
+    y = reshape (typecast (reshape (typecast (x(:), "uint8"), nb, numel (x))
+			   ([nb : -1 : 1], :) (:), clx), size(x));
+  endif
+endfunction
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/DLD-FUNCTIONS/bsxfun.cc	Thu Sep 06 16:38:44 2007 +0000
@@ -0,0 +1,483 @@
+/*
+
+Copyright (C) 2007 David Bateman
+
+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 2, 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, write to the Free
+Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+02110-1301, USA.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <string>
+#include <vector>
+#include <list>
+
+#include "lo-mappers.h"
+
+#include "oct-map.h"
+#include "defun-dld.h"
+#include "parse.h"
+#include "variables.h"
+#include "ov-colon.h"
+#include "unwind-prot.h"
+
+static bool
+maybe_update_column (octave_value& Ac, const octave_value& A, 
+		     const dim_vector& dva, const dim_vector& dvc,
+		     octave_idx_type i, octave_value_list &idx)
+{
+  octave_idx_type nd = dva.length ();
+
+  if (i == 0)
+    {
+      idx(0) = octave_value (':');
+      for (octave_idx_type j = 1; j < nd; j++)
+	{
+	  if (dva (j) == 1)
+	    idx (j) = octave_value (1);
+	  else
+	    idx (j) = octave_value ((i % dvc(j)) + 1);
+
+	  i = i / dvc (j);
+	}
+
+      Ac = A;
+      Ac = Ac.single_subsref ("(", idx);
+      return true;
+    }
+  else
+    {
+      bool is_changed = false;
+      octave_idx_type k = i;
+      octave_idx_type k1 = i - 1;
+      for (octave_idx_type j = 1; j < nd; j++)
+	{
+	  if (dva(j) != 1 && k % dvc (j) != k1 % dvc (j))
+	    {
+	      idx (j) = octave_value ((k % dvc(j)) + 1);
+	      is_changed = true;
+	    }
+
+	  k = k / dvc (j);
+	  k1 = k1 / dvc (j);
+	}
+
+      if (is_changed)
+	{
+	  Ac = A;
+	  Ac = Ac.single_subsref ("(", idx);
+	  return true;
+	}
+      else
+	return false;
+    }
+}
+
+static void
+update_index (octave_value_list& idx, const dim_vector& dv, octave_idx_type i)
+{
+  octave_idx_type nd = dv.length ();
+
+  if (i == 0)
+    {
+      for (octave_idx_type j = nd - 1; j > 0; j--)
+	idx(j) = octave_value (static_cast<double>(1));
+      idx(0) = octave_value (':');
+    }
+  else
+    {
+      for (octave_idx_type j = 1; j < nd; j++)
+	{
+	  idx (j) = octave_value (i % dv (j) + 1);
+	  i = i / dv (j);
+	}
+    }
+}
+
+static void
+update_index (Array<int>& idx, const dim_vector& dv, octave_idx_type i)
+{
+  octave_idx_type nd = dv.length ();
+
+  idx(0) = 0;
+  for (octave_idx_type j = 1; j < nd; j++)
+    {
+      idx (j) = i % dv (j);
+      i = i / dv (j);
+    }
+}
+
+DEFUN_DLD (bsxfun, args, nargout,
+  " -*- texinfo -*-\n\
+@deftypefn {Lodable Function} {} bsxfun (@var{f}, @var{a}, @var{b})\n\
+Applies a binary function @var{f} element-wise to two matrix arguments\n\
+@var{a} and @var{b}. The function @var{f} must be capable of accepting\n\
+two column vector arguments of equal length, or one column vector\n\
+argument and a scalar.\n\
+\n\
+The dimensions of @var{a} and @var{b} must be equal or singleton. The\n\
+singleton dimensions a the matirces will be expanded to the same\n\
+dimensioanlity as the other matrix.\n\
+\n\
+@seealso{arrayfun, cellfun}\n\
+@end deftypefn")
+{
+  int nargin = args.length ();
+  octave_value_list retval;
+
+  if (nargin != 3)
+    print_usage ();
+  else
+    {
+      octave_function *func = 0;
+      std::string name;
+      std::string fcn_name;
+
+      if (args(0).is_function_handle () || args(0).is_inline_function ())
+	func = args(0).function_value ();
+      else if (args(0).is_string ())
+	{
+	  name = args(0).string_value ();
+	  fcn_name = unique_symbol_name ("__bsxfun_fcn_");
+	  std::string fname = "function y = ";
+	  fname.append (fcn_name);
+	  fname.append ("(x) y = ");
+	  func = extract_function (args(0), "bsxfun", fcn_name, fname,
+				   "; endfunction");
+	}
+      else
+	  error ("bsxfun: first argument must be a string or function handle");
+
+      if (! error_state)
+	{
+	  const octave_value A = args (1);
+	  dim_vector dva = A.dims ();
+	  octave_idx_type nda = dva.length ();
+	  const octave_value B = args (2);
+	  dim_vector dvb = B.dims ();
+	  octave_idx_type ndb = dvb.length ();
+	  octave_idx_type nd = nda;
+      
+	  if (nda > ndb)
+	      dvb.resize (nda, 1);
+	  else if (nda < ndb)
+	    {
+	      dva.resize (ndb, 1);
+	      nd = ndb;
+	    }
+
+	  for (octave_idx_type i = 0; i < nd; i++)
+	    if (dva (i) != dvb (i) && dva (i) != 1 && dvb (i) != 1)
+	      {
+		error ("bsxfun: dimensions don't match");
+		break;
+	      }
+
+	  if (!error_state)
+	    {
+	      // Find the size of the output
+	      dim_vector dvc;
+	      dvc.resize (nd);
+	  
+	      for (octave_idx_type i = 0; i < nd; i++)
+		dvc (i) = (dva (i) < 1  ? dva (i) : (dvb (i) < 1 ? dvb (i) :
+		      (dva (i) > dvb (i) ? dva (i) : dvb (i))));
+
+	      if (dva == dvb || dva.numel () == 1 || dvb.numel () == 1)
+		{
+		  octave_value_list inputs;
+		  inputs (0) = A;
+		  inputs (1) = B;
+		  retval = feval (func, inputs, 1);
+		}
+	      else if (dvc.numel () < 1)
+		{
+		  octave_value_list inputs;
+		  inputs (0) = A.resize (dvc);
+		  inputs (1) = B.resize (dvc);
+		  retval = feval (func, inputs, 1);	      
+		}
+	      else
+		{
+		  octave_idx_type ncount = 1;
+		  for (octave_idx_type i = 1; i < nd; i++)
+		    ncount *= dvc (i);
+
+#define BSXDEF(T) \
+		  T result_ ## T; \
+		  bool have_ ## T = false;
+
+		  BSXDEF(NDArray);
+		  BSXDEF(ComplexNDArray);
+		  BSXDEF(boolNDArray);
+		  BSXDEF(int8NDArray);
+		  BSXDEF(int16NDArray);
+		  BSXDEF(int32NDArray);
+		  BSXDEF(int64NDArray);
+		  BSXDEF(uint8NDArray);
+		  BSXDEF(uint16NDArray);
+		  BSXDEF(uint32NDArray);
+		  BSXDEF(uint64NDArray);
+
+		  octave_value Ac ;
+		  octave_value_list idxA;
+		  octave_value Bc;
+		  octave_value_list idxB;
+		  octave_value C;
+		  octave_value_list inputs;
+		  Array<int> ra_idx (dvc.length(), 0);
+
+
+		  for (octave_idx_type i = 0; i < ncount; i++)
+		    {
+		      if (maybe_update_column (Ac, A, dva, dvc, i, idxA))
+			inputs (0) = Ac;
+
+		      if (maybe_update_column (Bc, B, dvb, dvc, i, idxB))
+			inputs (1) = Bc;
+			
+		      octave_value_list tmp = feval (func, inputs, 1);
+
+		      if (error_state)
+			break;
+
+#define BSXINIT(T, CLS, EXTRACTOR) \
+ 		      (result_type == CLS) \
+			{ \
+			    have_ ## T = true; \
+			    result_ ## T = \
+				tmp (0). EXTRACTOR ## _array_value (); \
+			    result_ ## T .resize (dvc); \
+                        }
+
+		      if (i == 0)
+			{
+			  if (! tmp(0).is_sparse_type ())
+			    {
+			      std::string result_type = tmp(0).class_name ();
+			      if (result_type == "double")
+				{
+				  if (tmp(0).is_real_type ())
+				    {
+				      have_NDArray = true;
+				      result_NDArray = tmp(0).array_value ();
+				      result_NDArray.resize (dvc);
+				    }
+				  else
+				    {
+				      have_ComplexNDArray = true;
+				      result_ComplexNDArray = 
+					tmp(0).complex_array_value ();
+				      result_ComplexNDArray.resize (dvc);
+				    }
+				}
+			      else if BSXINIT(boolNDArray, "logical", bool)
+			      else if BSXINIT(int8NDArray, "int8", int8)
+			      else if BSXINIT(int16NDArray, "int16", int16)
+			      else if BSXINIT(int32NDArray, "int32", int32)
+			      else if BSXINIT(int64NDArray, "int64", int64)
+			      else if BSXINIT(uint8NDArray, "uint8", uint8)
+			      else if BSXINIT(uint16NDArray, "uint16", uint16)
+			      else if BSXINIT(uint32NDArray, "uint32", uint32)
+			      else if BSXINIT(uint64NDArray, "uint64", uint64)
+			      else
+				{
+				  C = tmp (0);
+				  C = C.resize (dvc);
+				}
+			    }
+			}
+		      else
+			{
+			  update_index (ra_idx, dvc, i);
+			  
+			  if (have_NDArray)
+			    {
+			      if (tmp(0).class_name () != "double")
+				{
+				  have_NDArray = false;
+				  C = result_NDArray;
+				  C = do_cat_op (C, tmp(0), ra_idx);
+				}
+			      else if (tmp(0).is_real_type ())
+				result_NDArray.insert (tmp(0).array_value(), 
+						       ra_idx);
+			      else
+				{
+				  result_ComplexNDArray = 
+				    ComplexNDArray (result_NDArray);
+				  result_ComplexNDArray.insert 
+				    (tmp(0).complex_array_value(), ra_idx);
+				  have_NDArray = false;
+				  have_ComplexNDArray = true;
+				}
+			    }
+
+#define BSXLOOP(T, CLS, EXTRACTOR) \
+			(have_ ## T) \
+			  { \
+			    if (tmp (0).class_name () != CLS) \
+			      { \
+				have_ ## T = false; \
+				C = result_ ## T; \
+				C = do_cat_op (C, tmp (0), ra_idx); \
+			      } \
+			    else \
+			      result_ ## T .insert \
+				(tmp(0). EXTRACTOR ## _array_value (), \
+				ra_idx); \
+			  }
+
+			  else if BSXLOOP(ComplexNDArray, "double", complex)
+			  else if BSXLOOP(boolNDArray, "logical", bool)
+			  else if BSXLOOP(int8NDArray, "int8", int8)
+			  else if BSXLOOP(int16NDArray, "int16", int16)
+			  else if BSXLOOP(int32NDArray, "int32", int32)
+			  else if BSXLOOP(int64NDArray, "int64", int64)
+			  else if BSXLOOP(uint8NDArray, "uint8", uint8)
+			  else if BSXLOOP(uint16NDArray, "uint16", uint16)
+			  else if BSXLOOP(uint32NDArray, "uint32", uint32)
+			  else if BSXLOOP(uint64NDArray, "uint64", uint64)
+			  else
+			    C = do_cat_op (C, tmp(0), ra_idx);
+			}
+		    }
+
+#define BSXEND(T) \
+		  (have_ ## T) \
+		    retval (0) = result_ ## T;
+
+		  if BSXEND(NDArray)
+		  else if BSXEND(ComplexNDArray)
+		  else if BSXEND(boolNDArray)
+		  else if BSXEND(int8NDArray)
+		  else if BSXEND(int16NDArray)
+		  else if BSXEND(int32NDArray)
+		  else if BSXEND(int64NDArray)
+		  else if BSXEND(uint8NDArray)
+		  else if BSXEND(uint16NDArray)
+		  else if BSXEND(uint32NDArray)
+		  else if BSXEND(uint64NDArray)
+		  else
+		    retval(0) = C;
+		}
+	    }
+	}
+
+      if (! fcn_name.empty ())
+	clear_function (fcn_name);
+    }	
+
+  return retval;
+}
+
+/*
+
+%!shared a, b, c, f
+%! a = randn (4, 4);
+%! b = mean (a, 1);
+%! c = mean (a, 2);
+%! f = @minus;
+%!error(bsxfun (f));
+%!error(bsxfun (f, a));
+%!error(bsxfun (a, b));
+%!error(bsxfun (a, b, c));
+%!error(bsxfun (f, a, b, c));
+%!error(bsxfun (f, ones(4, 0), ones(4, 4)))
+%!assert(bsxfun (f, ones(4, 0), ones(4, 1)), zeros(4, 0));
+%!assert(bsxfun (f, ones(1, 4), ones(4, 1)), zeros(4, 4));
+%!assert(bsxfun (f, a, b), a - repmat(b, 4, 1));
+%!assert(bsxfun (f, a, c), a - repmat(c, 1, 4));
+%!assert(bsxfun ("minus", ones(1, 4), ones(4, 1)), zeros(4, 4));
+
+%!shared a, b, c, f
+%! a = randn (4, 4);
+%! a(1) *= 1i;
+%! b = mean (a, 1);
+%! c = mean (a, 2);
+%! f = @minus;
+%!error(bsxfun (f));
+%!error(bsxfun (f, a));
+%!error(bsxfun (a, b));
+%!error(bsxfun (a, b, c));
+%!error(bsxfun (f, a, b, c));
+%!error(bsxfun (f, ones(4, 0), ones(4, 4)))
+%!assert(bsxfun (f, ones(4, 0), ones(4, 1)), zeros(4, 0));
+%!assert(bsxfun (f, ones(1, 4), ones(4, 1)), zeros(4, 4));
+%!assert(bsxfun (f, a, b), a - repmat(b, 4, 1));
+%!assert(bsxfun (f, a, c), a - repmat(c, 1, 4));
+%!assert(bsxfun ("minus", ones(1, 4), ones(4, 1)), zeros(4, 4));
+
+%!shared a, b, c, f
+%! a = randn (4, 4);
+%! a(end) *= 1i;
+%! b = mean (a, 1);
+%! c = mean (a, 2);
+%! f = @minus;
+%!error(bsxfun (f));
+%!error(bsxfun (f, a));
+%!error(bsxfun (a, b));
+%!error(bsxfun (a, b, c));
+%!error(bsxfun (f, a, b, c));
+%!error(bsxfun (f, ones(4, 0), ones(4, 4)))
+%!assert(bsxfun (f, ones(4, 0), ones(4, 1)), zeros(4, 0));
+%!assert(bsxfun (f, ones(1, 4), ones(4, 1)), zeros(4, 4));
+%!assert(bsxfun (f, a, b), a - repmat(b, 4, 1));
+%!assert(bsxfun (f, a, c), a - repmat(c, 1, 4));
+%!assert(bsxfun ("minus", ones(1, 4), ones(4, 1)), zeros(4, 4));
+
+%!shared a, b, c, f
+%! a = randn (4, 4);
+%! b = a (1, :);
+%! c = a (:, 1);
+%! f = @(x, y) x == y;
+%!error(bsxfun (f));
+%!error(bsxfun (f, a));
+%!error(bsxfun (a, b));
+%!error(bsxfun (a, b, c));
+%!error(bsxfun (f, a, b, c));
+%!error(bsxfun (f, ones(4, 0), ones(4, 4)))
+%!assert(bsxfun (f, ones(4, 0), ones(4, 1)), zeros(4, 0, "logical"));
+%!assert(bsxfun (f, ones(1, 4), ones(4, 1)), ones(4, 4, "logical"));
+%!assert(bsxfun (f, a, b), a == repmat(b, 4, 1));
+%!assert(bsxfun (f, a, c), a == repmat(c, 1, 4));
+
+%!shared a, b, c, d, f
+%! a = randn (4, 4, 4);
+%! b = mean (a, 1);
+%! c = mean (a, 2);
+%! d = mean (a, 3);
+%! f = @minus;
+%!error(bsxfun (f, ones([4, 0, 4]), ones([4, 4, 4])));
+%!assert(bsxfun (f, ones([4, 0, 4]), ones([4, 1, 4])), zeros([4, 0, 4]));
+%!assert(bsxfun (f, ones([4, 4, 0]), ones([4, 1, 1])), zeros([4, 4, 0]));
+%!assert(bsxfun (f, ones([1, 4, 4]), ones([4, 1, 4])), zeros([4, 4, 4]));
+%!assert(bsxfun (f, ones([4, 4, 1]), ones([4, 1, 4])), zeros([4, 4, 4]));
+%!assert(bsxfun (f, ones([4, 1, 4]), ones([1, 4, 4])), zeros([4, 4, 4]));
+%!assert(bsxfun (f, ones([4, 1, 4]), ones([1, 4, 1])), zeros([4, 4, 4]));
+%!assert(bsxfun (f, a, b), a - repmat(b, [4, 1, 1]));
+%!assert(bsxfun (f, a, c), a - repmat(c, [1, 4, 1]));
+%!assert(bsxfun (f, a, d), a - repmat(d, [1, 1, 4]));
+%!assert(bsxfun ("minus", ones([4, 0, 4]), ones([4, 1, 4])), zeros([4, 0, 4]));
+
+%% The below is a very hard case to treat
+%!assert(bsxfun (f, ones([4, 1, 4, 1]), ones([1, 4, 1, 4])), zeros([4, 4, 4, 4]));
+
+*/
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/DLD-FUNCTIONS/typecast.cc	Thu Sep 06 16:38:44 2007 +0000
@@ -0,0 +1,211 @@
+/*
+
+Copyright (C) 2007 David Bateman
+
+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 2, 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, write to the Free
+Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+02110-1301, USA.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <string>
+
+#include "oct.h"
+
+template <class LT, class RT>
+static void
+typecast (const Array<RT>& x, Array<LT>& y)
+{
+  octave_idx_type n = x.length ();
+  size_t ns = sizeof (RT);
+  size_t ms = sizeof (LT);
+
+  if (n * ns % ms != 0)
+      error ("typecast: incorrect number of input values to make output value");
+  else
+    {
+      octave_idx_type m = n * ns / ms;
+      dim_vector dv = x.dims ();
+      for (octave_idx_type i = 0; i < dv.length(); i++)
+	if (dv(i) == n)
+	  {
+	    dv(i) = m;
+	    break;
+	  }
+      y.resize (dv);
+      const unsigned char *xp = reinterpret_cast<const unsigned char *>
+	(x.fortran_vec ());
+      unsigned char *yp = reinterpret_cast<unsigned char *>(y.fortran_vec ());
+      for (octave_idx_type i = 0; 
+	   i < n * static_cast<octave_idx_type>(ns); i++)
+	*yp++ = *xp++;
+    }
+}
+
+template <class T>
+static octave_value
+typecast (const T& x, std::string type)
+{
+  octave_value retval;
+  if (type == "uint8")
+    {
+      uint8NDArray y;
+      typecast (x, y);
+      retval = octave_value (y);
+    }
+  else if (type == "uint16")
+    {
+      uint16NDArray y;
+      typecast (x, y);
+      retval = octave_value (y);
+    }
+  else if (type == "uint32")
+    {
+      uint32NDArray y;
+      typecast (x, y);
+      retval = octave_value (y);
+    }
+  else if (type == "uint64")
+    {
+      uint64NDArray y;
+      typecast (x, y);
+      retval = octave_value (y);
+    }
+  else if (type == "int8")
+    {
+      int8NDArray y;
+      typecast (x, y);
+      retval = octave_value (y);
+    }
+  else if (type == "int16")
+    {
+      int16NDArray y;
+      typecast (x, y);
+      retval = octave_value (y);
+    }
+  else if (type == "int32")
+    {
+      int32NDArray y;
+      typecast (x, y);
+      retval = octave_value (y);
+    }
+  else if (type == "int64")
+    {
+      int64NDArray y;
+      typecast (x, y);
+      retval = octave_value (y);
+    }
+  else
+    {
+      NDArray y;
+      typecast (x, y);
+      retval = octave_value (y);
+    }
+
+  return retval;
+}
+
+DEFUN_DLD (typecast, args, ,
+  " -*- texinfo -*-\n\
+@deftypefn {Loadable Function} {} typecast (@var{x}, @var{type})\n\
+Converts from one datatype to another without changing the underlying\n\
+data. The argument @var{type} defines the type of the return argument\n\
+and must be one of 'uint8', 'uint16', 'uint32', 'uint64', 'int8', 'int16',\n\
+'int32', 'int64', 'single' or 'double'.\n\
+\n\
+An example of the use of typecast on a little-endian machine is\n\
+\n\
+@example\n\
+@group\n\
+@var{x} = uint16 ([1, 65535]);\n\
+typecast (@var{x}, 'uint8')\n\
+@result{} [   0,   1, 255, 255]\n\
+@end group\n\
+@end example\n\
+@seealso{cast, swapbytes}\n\
+@end deftypefn")
+{
+  int nargin = args.length ();
+  octave_value retval;
+
+  if (nargin != 2)
+    print_usage ();
+  else
+    {
+      std::string type = args (1).string_value ();
+
+      if (! error_state)
+	{
+	  std::transform (type.begin (), type.end (), type.begin (), tolower);
+
+	  if (type == "single")
+	    error ("typecast: type 'single' is not supported");
+	  else if (type != "uint8" && type != "uint16" &&
+	      type != "uint32" && type != "uint64" &&
+	      type != "int8" && type != "int16" &&
+	      type != "int32" && type != "int64" &&
+	      type != "single" && type != "double")
+	    error ("typecast: unrecognized or invalid type");
+	  else if (args(0).is_sparse_type () || args(0).is_complex_type ())
+	    error ("typecast: sparse and complex types are invalid");
+	  else
+	    {
+	      dim_vector dv = args(0).dims ();
+	      bool seen = false;
+
+	      for (octave_idx_type i = 0; i < dv.length(); i++)
+		if (dv (i) != 1)
+		  {
+		    if (seen)
+		      {
+			error ("typecast: argument must be a vector");
+			break;
+		      }
+		    else
+		      seen = true;
+		  }
+
+	      if (!error_state)
+		{
+		  if (args(0).is_uint8_type ())
+		    retval = typecast (args(0).uint8_array_value (), type); 
+		  else if (args(0).is_uint16_type ())
+		    retval = typecast (args(0).uint16_array_value (), type); 
+		  else if (args(0).is_uint32_type ())
+		    retval = typecast (args(0).uint32_array_value (), type); 
+		  else if (args(0).is_uint64_type ())
+		    retval = typecast (args(0).uint64_array_value (), type); 
+		  else if (args(0).is_int8_type ())
+		    retval = typecast (args(0).int8_array_value (), type); 
+		  else if (args(0).is_int16_type ())
+		    retval = typecast (args(0).int16_array_value (), type); 
+		  else if (args(0).is_int32_type ())
+		    retval = typecast (args(0).int32_array_value (), type); 
+		  else if (args(0).is_int64_type ())
+		    retval = typecast (args(0).int64_array_value (), type); 
+		  else
+		    retval = typecast (args(0).array_value (), type);
+		}
+	    }
+	}
+    }
+
+  return retval;
+}