GNU Octave  4.4.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
data.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1994-2018 John W. Eaton
4 Copyright (C) 2009 Jaroslav Hajek
5 Copyright (C) 2009-2010 VZLU Prague
6 Copyright (C) 2012 Carlo de Falco
7 
8 This file is part of Octave.
9 
10 Octave is free software: you can redistribute it and/or modify it
11 under the terms of the GNU General Public License as published by
12 the Free Software Foundation, either version 3 of the License, or
13 (at your option) any later version.
14 
15 Octave is distributed in the hope that it will be useful, but
16 WITHOUT ANY WARRANTY; without even the implied warranty of
17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 GNU General Public License for more details.
19 
20 You should have received a copy of the GNU General Public License
21 along with Octave; see the file COPYING. If not, see
22 <https://www.gnu.org/licenses/>.
23 
24 */
25 
26 #if defined (HAVE_CONFIG_H)
27 # include "config.h"
28 #endif
29 
30 #include <cmath>
31 #include <cstddef>
32 #include <cstdint>
33 #include <ctime>
34 
35 #include <algorithm>
36 #include <limits>
37 #include <string>
38 
39 #include "lo-ieee.h"
40 #include "mx-base.h"
41 #include "oct-base64.h"
42 #include "oct-binmap.h"
43 #include "oct-time.h"
44 #include "quit.h"
45 
46 #include "Cell.h"
47 #include "data.h"
48 #include "defun.h"
49 #include "error.h"
50 #include "errwarn.h"
51 #include "interpreter-private.h"
52 #include "oct-map.h"
53 #include "ov-class.h"
54 #include "ov-complex.h"
55 #include "ov-cx-mat.h"
56 #include "ov-cx-sparse.h"
57 #include "ov-float.h"
58 #include "ov-flt-complex.h"
59 #include "ov-flt-cx-mat.h"
60 #include "ov.h"
61 #include "ovl.h"
62 #include "pager.h"
63 #include "parse.h"
64 #include "pt-mat.h"
65 #include "utils.h"
66 #include "variables.h"
67 #include "xnorm.h"
68 
69 static void
70 index_error (const char *fmt, const std::string& idx, const std::string& msg)
71 {
72  error (fmt, idx.c_str (), msg.c_str ());
73 }
74 
75 DEFUN (all, args, ,
76  doc: /* -*- texinfo -*-
77 @deftypefn {} {} all (@var{x})
78 @deftypefnx {} {} all (@var{x}, @var{dim})
79 For a vector argument, return true (logical 1) if all elements of the vector
80 are nonzero.
81 
82 For a matrix argument, return a row vector of logical ones and
83 zeros with each element indicating whether all of the elements of the
84 corresponding column of the matrix are nonzero. For example:
85 
86 @example
87 @group
88 all ([2, 3; 1, 0])
89  @result{} [ 1, 0 ]
90 @end group
91 @end example
92 
93 If the optional argument @var{dim} is supplied, work along dimension
94 @var{dim}.
95 @seealso{any}
96 @end deftypefn */)
97 {
98  int nargin = args.length ();
99 
101  print_usage ();
102 
103  int dim = (nargin == 1 ? -1
104  : args(1).int_value ("all: DIM must be an integer")-1);
105 
106  if (dim < -1)
107  error ("all: invalid dimension argument = %d", dim + 1);
108 
109  return ovl (args(0).all (dim));
110 }
111 
112 /*
113 %!test
114 %! x = ones (3);
115 %! x(1,1) = 0;
116 %! assert (all (all (rand (3) + 1) == [1, 1, 1]) == 1);
117 %! assert (all (all (x) == [0, 1, 1]) == 1);
118 %! assert (all (x, 1) == [0, 1, 1]);
119 %! assert (all (x, 2) == [0; 1; 1]);
120 
121 %!test
122 %! x = ones (3, "single");
123 %! x(1,1) = 0;
124 %! assert (all (all (single (rand (3) + 1)) == [1, 1, 1]) == 1);
125 %! assert (all (all (x) == [0, 1, 1]) == 1);
126 %! assert (all (x, 1) == [0, 1, 1]);
127 %! assert (all (x, 2) == [0; 1; 1]);
128 
129 %!error all ()
130 %!error all (1, 2, 3)
131 */
132 
133 DEFUN (any, args, ,
134  doc: /* -*- texinfo -*-
135 @deftypefn {} {} any (@var{x})
136 @deftypefnx {} {} any (@var{x}, @var{dim})
137 For a vector argument, return true (logical 1) if any element of the vector
138 is nonzero.
139 
140 For a matrix argument, return a row vector of logical ones and
141 zeros with each element indicating whether any of the elements of the
142 corresponding column of the matrix are nonzero. For example:
143 
144 @example
145 @group
146 any (eye (2, 4))
147  @result{} [ 1, 1, 0, 0 ]
148 @end group
149 @end example
150 
151 If the optional argument @var{dim} is supplied, work along dimension
152 @var{dim}. For example:
153 
154 @example
155 @group
156 any (eye (2, 4), 2)
157  @result{} [ 1; 1 ]
158 @end group
159 @end example
160 @seealso{all}
161 @end deftypefn */)
162 {
163  int nargin = args.length ();
164 
166  print_usage ();
167 
168  int dim = (nargin == 1 ? -1
169  : args(1).int_value ("any: DIM must be an integer")-1);
170 
171  if (dim < -1)
172  error ("any: invalid dimension argument = %d", dim + 1);
173 
174  return ovl (args(0).any (dim));
175 }
176 
177 /*
178 %!test
179 %! x = zeros (3);
180 %! x(3,3) = 1;
181 %! assert (all (any (x) == [0, 0, 1]) == 1);
182 %! assert (all (any (ones (3)) == [1, 1, 1]) == 1);
183 %! assert (any (x, 1) == [0, 0, 1]);
184 %! assert (any (x, 2) == [0; 0; 1]);
185 
186 %!test
187 %! x = zeros (3, "single");
188 %! x(3,3) = 1;
189 %! assert (all (any (x) == [0, 0, 1]) == 1);
190 %! assert (all (any (ones (3, "single")) == [1, 1, 1]) == 1);
191 %! assert (any (x, 1) == [0, 0, 1]);
192 %! assert (any (x, 2) == [0; 0; 1]);
193 
194 %!error any ()
195 %!error any (1, 2, 3)
196 */
197 
198 // These mapping functions may also be useful in other places, eh?
199 
200 DEFUN (atan2, args, ,
201  doc: /* -*- texinfo -*-
202 @deftypefn {} {} atan2 (@var{y}, @var{x})
203 Compute atan (@var{y} / @var{x}) for corresponding elements of @var{y}
204 and @var{x}.
205 
206 @var{y} and @var{x} must match in size and orientation. The signs of
207 elements of @var{y} and @var{x} are used to determine the quadrants of each
208 resulting value.
209 
210 This function is equivalent to @code{arg (complex (@var{x}, @var{y}))}.
211 @seealso{tan, tand, tanh, atanh}
212 @end deftypefn */)
213 {
214  if (args.length () != 2)
215  print_usage ();
216 
218 
219  if (! args(0).isnumeric ())
220  err_wrong_type_arg ("atan2", args(0));
221 
222  if (! args(1).isnumeric ())
223  err_wrong_type_arg ("atan2", args(1));
224 
225  if (args(0).iscomplex () || args(1).iscomplex ())
226  error ("atan2: not defined for complex numbers");
227 
228  if (args(0).is_single_type () || args(1).is_single_type ())
229  {
230  if (args(0).is_scalar_type () && args(1).is_scalar_type ())
231  retval = atan2f (args(0).float_value (), args(1).float_value ());
232  else
233  {
234  FloatNDArray a0 = args(0).float_array_value ();
235  FloatNDArray a1 = args(1).float_array_value ();
236  retval = binmap<float> (a0, a1, std::atan2, "atan2");
237  }
238  }
239  else
240  {
241  if (args(0).is_scalar_type () && args(1).is_scalar_type ())
242  retval = atan2 (args(0).scalar_value (), args(1).scalar_value ());
243  else if (args(0).issparse ())
244  {
245  SparseMatrix m0 = args(0).sparse_matrix_value ();
246  SparseMatrix m1 = args(1).sparse_matrix_value ();
247  retval = binmap<double> (m0, m1, std::atan2, "atan2");
248  }
249  else
250  {
251  NDArray a0 = args(0).array_value ();
252  NDArray a1 = args(1).array_value ();
253  retval = binmap<double> (a0, a1, std::atan2, "atan2");
254  }
255  }
256 
257  return retval;
258 }
259 
260 /*
261 %!assert (size (atan2 (zeros (0, 2), zeros (0, 2))), [0, 2])
262 %!assert (size (atan2 (rand (2, 3, 4), zeros (2, 3, 4))), [2, 3, 4])
263 %!assert (size (atan2 (rand (2, 3, 4), 1)), [2, 3, 4])
264 %!assert (size (atan2 (1, rand (2, 3, 4))), [2, 3, 4])
265 %!assert (size (atan2 (1, 2)), [1, 1])
266 
267 %!test
268 %! rt2 = sqrt (2);
269 %! rt3 = sqrt (3);
270 %! v = [0, pi/6, pi/4, pi/3, -pi/3, -pi/4, -pi/6, 0];
271 %! y = [0, rt3, 1, rt3, -rt3, -1, -rt3, 0];
272 %! x = [1, 3, 1, 1, 1, 1, 3, 1];
273 %! assert (atan2 (y, x), v, sqrt (eps));
274 
275 %!test
276 %! rt2 = sqrt (2);
277 %! rt3 = sqrt (3);
278 %! v = single ([0, pi/6, pi/4, pi/3, -pi/3, -pi/4, -pi/6, 0]);
279 %! y = single ([0, rt3, 1, rt3, -rt3, -1, -rt3, 0]);
280 %! x = single ([1, 3, 1, 1, 1, 1, 3, 1]);
281 %! assert (atan2 (y, x), v, sqrt (eps ("single")));
282 
283 ## Test sparse implementations
284 %!shared xs
285 %! xs = sparse (0:3);
286 %!test
287 %! y = atan2 (1, xs);
288 %! assert (issparse (y), false);
289 %! assert (nnz (y), 4);
290 %! assert (y, atan2 (1, 0:3));
291 %!test
292 %! y = atan2 (0, xs);
293 %! assert (issparse (y), false);
294 %! assert (nnz (y), 0);
295 %! assert (y, zeros (1,4));
296 %!test
297 %! y = atan2 (xs, 1);
298 %! assert (issparse (y));
299 %! assert (nnz (y), 3);
300 %! assert (y, sparse (atan2 (0:3, 1)));
301 %!test
302 %! y = atan2 (xs, 0);
303 %! assert (issparse (y));
304 %! assert (nnz (y), 3);
305 %! assert (y, sparse (atan2 (0:3, 0)));
306 %!test
307 %! y = atan2 (xs, sparse (ones (1, 4)));
308 %! assert (issparse (y));
309 %! assert (nnz (y), 3);
310 %! assert (y, sparse (atan2 (0:3, ones (1,4))));
311 %!test
312 %! y = atan2 (xs, sparse (zeros (1,4)));
313 %! assert (issparse (y));
314 %! assert (nnz (y), 3);
315 %! assert (y, sparse (atan2 (0:3, zeros (1,4))));
316 
317 %!error atan2 ()
318 %!error atan2 (1, 2, 3)
319 */
320 
321 static octave_value
322 do_hypot (const octave_value& x, const octave_value& y)
323 {
325 
326  octave_value arg0 = x;
327  octave_value arg1 = y;
328  if (! arg0.isnumeric ())
329  err_wrong_type_arg ("hypot", arg0);
330  if (! arg1.isnumeric ())
331  err_wrong_type_arg ("hypot", arg1);
332 
333  if (arg0.iscomplex ())
334  arg0 = arg0.abs ();
335  if (arg1.iscomplex ())
336  arg1 = arg1.abs ();
337 
338  if (arg0.is_single_type () || arg1.is_single_type ())
339  {
340  if (arg0.is_scalar_type () && arg1.is_scalar_type ())
341  retval = hypotf (arg0.float_value (), arg1.float_value ());
342  else
343  {
344  FloatNDArray a0 = arg0.float_array_value ();
345  FloatNDArray a1 = arg1.float_array_value ();
346  retval = binmap<float> (a0, a1, std::hypot, "hypot");
347  }
348  }
349  else
350  {
351  if (arg0.is_scalar_type () && arg1.is_scalar_type ())
352  retval = hypot (arg0.scalar_value (), arg1.scalar_value ());
353  else if (arg0.issparse () || arg1.issparse ())
354  {
355  SparseMatrix m0 = arg0.sparse_matrix_value ();
356  SparseMatrix m1 = arg1.sparse_matrix_value ();
357  retval = binmap<double> (m0, m1, std::hypot, "hypot");
358  }
359  else
360  {
361  NDArray a0 = arg0.array_value ();
362  NDArray a1 = arg1.array_value ();
363  retval = binmap<double> (a0, a1, std::hypot, "hypot");
364  }
365  }
366 
367  return retval;
368 }
369 
370 DEFUN (hypot, args, ,
371  doc: /* -*- texinfo -*-
372 @deftypefn {} {} hypot (@var{x}, @var{y})
373 @deftypefnx {} {} hypot (@var{x}, @var{y}, @var{z}, @dots{})
374 Compute the element-by-element square root of the sum of the squares of
375 @var{x} and @var{y}.
376 
377 This is equivalent to
378 @code{sqrt (@var{x}.^2 + @var{y}.^2)}, but is calculated in a manner that
379 avoids overflows for large values of @var{x} or @var{y}.
380 
381 @code{hypot} can also be called with more than 2 arguments; in this case,
382 the arguments are accumulated from left to right:
383 
384 @example
385 @group
386 hypot (hypot (@var{x}, @var{y}), @var{z})
387 hypot (hypot (hypot (@var{x}, @var{y}), @var{z}), @var{w}), etc.
388 @end group
389 @end example
390 @end deftypefn */)
391 {
392  int nargin = args.length ();
393 
394  if (nargin < 2)
395  print_usage ();
396 
398 
399  if (nargin == 2)
400  retval = do_hypot (args(0), args(1));
401  else
402  {
403  retval = args(0);
404 
405  for (int i = 1; i < nargin; i++)
406  retval = do_hypot (retval, args(i));
407  }
408 
409  return retval;
410 }
411 
412 /*
413 %!assert (size (hypot (zeros (0, 2), zeros (0, 2))), [0, 2])
414 %!assert (size (hypot (rand (2, 3, 4), zeros (2, 3, 4))), [2, 3, 4])
415 %!assert (size (hypot (rand (2, 3, 4), 1)), [2, 3, 4])
416 %!assert (size (hypot (1, rand (2, 3, 4))), [2, 3, 4])
417 %!assert (size (hypot (1, 2)), [1, 1])
418 %!assert (hypot (1:10, 1:10), sqrt (2) * [1:10], 16*eps)
419 %!assert (hypot (single (1:10), single (1:10)), single (sqrt (2) * [1:10]))
420 
421 ## Test sparse implementations
422 %!shared xs
423 %! xs = sparse (0:3);
424 %!test
425 %! y = hypot (1, xs);
426 %! assert (nnz (y), 4);
427 %! assert (y, sparse (hypot (1, 0:3)));
428 %!test
429 %! y = hypot (0, xs);
430 %! assert (nnz (y), 3);
431 %! assert (y, xs);
432 %!test
433 %! y = hypot (xs, 1);
434 %! assert (nnz (y), 4);
435 %! assert (y, sparse (hypot (0:3, 1)));
436 %!test
437 %! y = hypot (xs, 0);
438 %! assert (nnz (y), 3);
439 %! assert (y, xs);
440 %!test
441 %! y = hypot (sparse ([0 0]), sparse ([0 1]));
442 %! assert (nnz (y), 1);
443 %! assert (y, sparse ([0 1]));
444 %!test
445 %! y = hypot (sparse ([0 1]), sparse ([0 0]));
446 %! assert (nnz (y), 1);
447 %! assert (y, sparse ([0 1]));
448 
449 */
450 
451 template <typename T, typename ET>
452 void
453 map_2_xlog2 (const Array<T>& x, Array<T>& f, Array<ET>& e)
454 {
455  f = Array<T>(x.dims ());
456  e = Array<ET>(x.dims ());
457  for (octave_idx_type i = 0; i < x.numel (); i++)
458  {
459  int exp;
460  f.xelem (i) = octave::math::log2 (x(i), exp);
461  e.xelem (i) = exp;
462  }
463 }
464 
465 DEFUN (log2, args, nargout,
466  doc: /* -*- texinfo -*-
467 @deftypefn {} {} log2 (@var{x})
468 @deftypefnx {} {[@var{f}, @var{e}] =} log2 (@var{x})
469 Compute the base-2 logarithm of each element of @var{x}.
470 
471 If called with two output arguments, split @var{x} into
472 binary mantissa and exponent so that
473 @tex
474 ${1 \over 2} \le \left| f \right| < 1$
475 @end tex
476 @ifnottex
477 @w{@code{1/2 <= abs(f) < 1}}
478 @end ifnottex
479 and @var{e} is an integer. If
480 @tex
481 $x = 0$, $f = e = 0$.
482 @end tex
483 @ifnottex
484 @w{@code{x = 0}}, @w{@code{f = e = 0}}.
485 @end ifnottex
486 @seealso{pow2, log, log10, exp}
487 @end deftypefn */)
488 {
489  if (args.length () != 1)
490  print_usage ();
491 
493 
494  if (nargout < 2)
495  retval = ovl (args(0).log2 ());
496  else if (args(0).is_single_type ())
497  {
498  if (args(0).isreal ())
499  {
500  FloatNDArray f;
501  FloatNDArray x = args(0).float_array_value ();
502  // FIXME: should E be an int value?
503  FloatMatrix e;
504  map_2_xlog2 (x, f, e);
505  retval = ovl (f, e);
506  }
507  else if (args(0).iscomplex ())
508  {
510  FloatComplexNDArray x = args(0).float_complex_array_value ();
511  // FIXME: should E be an int value?
512  FloatNDArray e;
513  map_2_xlog2 (x, f, e);
514  retval = ovl (f, e);
515  }
516  }
517  else if (args(0).isreal ())
518  {
519  NDArray f;
520  NDArray x = args(0).array_value ();
521  // FIXME: should E be an int value?
522  Matrix e;
523  map_2_xlog2 (x, f, e);
524  retval = ovl (f, e);
525  }
526  else if (args(0).iscomplex ())
527  {
529  ComplexNDArray x = args(0).complex_array_value ();
530  // FIXME: should E be an int value?
531  NDArray e;
532  map_2_xlog2 (x, f, e);
533  retval = ovl (f, e);
534  }
535  else
536  err_wrong_type_arg ("log2", args(0));
537 
538  return retval;
539 }
540 
541 /*
542 %!assert (log2 ([1/4, 1/2, 1, 2, 4]), [-2, -1, 0, 1, 2])
543 %!assert (log2 (Inf), Inf)
544 %!assert (isnan (log2 (NaN)))
545 %!assert (log2 (4*i), 2 + log2 (1*i))
546 %!assert (log2 (complex (0,Inf)), Inf + log2 (i))
547 
548 %!test
549 %! [f, e] = log2 ([0,-1; 2,-4; Inf,-Inf]);
550 %! assert (f, [0,-0.5; 0.5,-0.5; Inf,-Inf]);
551 %! assert (e(1:2,:), [0,1;2,3]);
552 
553 %!test
554 %! [f, e] = log2 (complex (zeros (3, 2), [0,-1; 2,-4; Inf,-Inf]));
555 %! assert (f, complex (zeros (3, 2), [0,-0.5; 0.5,-0.5; Inf,-Inf]));
556 %! assert (e(1:2,:), [0,1; 2,3]);
557 
558 %!assert <*42583> (all (log2 (pow2 (-1074:1023)) == -1074:1023))
559 */
560 
561 DEFUN (rem, args, ,
562  doc: /* -*- texinfo -*-
563 @deftypefn {} {} rem (@var{x}, @var{y})
564 Return the remainder of the division @code{@var{x} / @var{y}}.
565 
566 The remainder is computed using the expression
567 
568 @example
569 x - y .* fix (x ./ y)
570 @end example
571 
572 An error message is printed if the dimensions of the arguments do not agree,
573 or if either argument is complex.
574 
575 Programming Notes: Floating point numbers within a few eps of an integer
576 will be rounded to an integer before computation for compatibility with
577 @sc{matlab}.
578 
579 By convention,
580 
581 @example
582 @group
583 rem (@var{x}, 0) = NaN if @var{x} is a floating point variable
584 rem (@var{x}, 0) = 0 if @var{x} is an integer variable
585 rem (@var{x}, @var{y}) returns a value with the signbit from @var{x}
586 @end group
587 @end example
588 
589 For the opposite conventions see the @code{mod} function. In general,
590 @code{rem} is best when computing the remainder after division of two
591 @emph{positive} numbers. For negative numbers, or when the values are
592 periodic, @code{mod} is a better choice.
593 @seealso{mod}
594 @end deftypefn */)
595 {
596  if (args.length () != 2)
597  print_usage ();
598 
600 
601  if (! args(0).isnumeric ())
602  err_wrong_type_arg ("rem", args(0));
603 
604  if (! args(1).isnumeric ())
605  err_wrong_type_arg ("rem", args(1));
606 
607  if (args(0).iscomplex () || args(1).iscomplex ())
608  error ("rem: not defined for complex numbers");
609 
610  if (args(0).isinteger () || args(1).isinteger ())
611  {
612  builtin_type_t btyp0 = args(0).builtin_type ();
613  builtin_type_t btyp1 = args(1).builtin_type ();
614  if (btyp0 == btyp_double || btyp0 == btyp_float)
615  btyp0 = btyp1;
616  if (btyp1 == btyp_double || btyp1 == btyp_float)
617  btyp1 = btyp0;
618 
619  if (btyp0 != btyp1)
620  error ("rem: cannot combine %s and %d",
621  args(0).class_name ().c_str (),
622  args(1).class_name ().c_str ());
623 
624  switch (btyp0)
625  {
626 #define MAKE_INT_BRANCH(X) \
627  case btyp_ ## X: \
628  { \
629  X##NDArray a0 = args(0).X##_array_value (); \
630  X##NDArray a1 = args(1).X##_array_value (); \
631  retval = binmap<octave_##X,octave_##X,octave_##X> (a0, a1, rem, "rem"); \
632  } \
633  break
634 
635  MAKE_INT_BRANCH (int8);
636  MAKE_INT_BRANCH (int16);
637  MAKE_INT_BRANCH (int32);
638  MAKE_INT_BRANCH (int64);
639  MAKE_INT_BRANCH (uint8);
640  MAKE_INT_BRANCH (uint16);
641  MAKE_INT_BRANCH (uint32);
642  MAKE_INT_BRANCH (uint64);
643 
644 #undef MAKE_INT_BRANCH
645 
646  default:
647  panic_impossible ();
648  }
649  }
650  else if (args(0).is_single_type () || args(1).is_single_type ())
651  {
652  if (args(0).is_scalar_type () && args(1).is_scalar_type ())
653  retval = octave::math::rem (args(0).float_value (), args(1).float_value ());
654  else
655  {
656  FloatNDArray a0 = args(0).float_array_value ();
657  FloatNDArray a1 = args(1).float_array_value ();
658  retval = binmap<float> (a0, a1, octave::math::rem<float>, "rem");
659  }
660  }
661  else
662  {
663  if (args(0).is_scalar_type () && args(1).is_scalar_type ())
664  retval = octave::math::rem (args(0).scalar_value (), args(1).scalar_value ());
665  else if (args(0).issparse () || args(1).issparse ())
666  {
667  SparseMatrix m0 = args(0).sparse_matrix_value ();
668  SparseMatrix m1 = args(1).sparse_matrix_value ();
669  retval = binmap<double> (m0, m1, octave::math::rem<double>, "rem");
670  }
671  else
672  {
673  NDArray a0 = args(0).array_value ();
674  NDArray a1 = args(1).array_value ();
675  retval = binmap<double> (a0, a1, octave::math::rem<double>, "rem");
676  }
677  }
678 
679  return retval;
680 }
681 
682 /*
683 %!assert (size (rem (zeros (0, 2), zeros (0, 2))), [0, 2])
684 %!assert (size (rem (rand (2, 3, 4), zeros (2, 3, 4))), [2, 3, 4])
685 %!assert (size (rem (rand (2, 3, 4), 1)), [2, 3, 4])
686 %!assert (size (rem (1, rand (2, 3, 4))), [2, 3, 4])
687 %!assert (size (rem (1, 2)), [1, 1])
688 
689 %!assert (rem ([1, 2, 3; -1, -2, -3], 2), [1, 0, 1; -1, 0, -1])
690 %!assert (rem ([1, 2, 3; -1, -2, -3], 2 * ones (2, 3)),[1, 0, 1; -1, 0, -1])
691 %!assert (rem ([0, 1, 2], [0, 0, 1]), [NaN, NaN, 0])
692 %!assert (rem (uint8 ([1, 2, 3; -1, -2, -3]), uint8 (2)), uint8 ([1, 0, 1; -1, 0, -1]))
693 %!assert (uint8 (rem ([1, 2, 3; -1, -2, -3], 2 * ones (2, 3))),uint8 ([1, 0, 1; -1, 0, -1]))
694 %!assert (rem (uint8 ([0, 1, 2]), [0, 0, 1]), uint8 ([0, 0, 0]))
695 
696 ## Test sparse implementations
697 %!shared xs
698 %! xs = sparse (0:3);
699 %!test
700 %! y = rem (11, xs);
701 %! assert (isnan (y(1)));
702 %! assert (y, sparse (rem (11, 0:3)));
703 %!test
704 %! y = rem (0, xs);
705 %! assert (nnz (y), 1);
706 %! assert (y, sparse ([NaN 0 0 0]));
707 %!test
708 %! y = rem (xs, 2);
709 %! assert (nnz (y), 2);
710 %! assert (y, sparse (rem (0:3, 2)));
711 %!test
712 %! y = rem (xs, 1);
713 %! assert (nnz (y), 0);
714 %! assert (y, sparse (rem (0:3, 1)));
715 %!test
716 %! y = rem (sparse ([11 11 11 11]), xs);
717 %! assert (nnz (y), 3);
718 %! assert (y, sparse (rem (11, 0:3)));
719 %!test
720 %! y = rem (sparse ([0 0 0 0]), xs);
721 %! assert (nnz (y), 1);
722 %! assert (y, sparse ([NaN 0 0 0]));
723 
724 %!assert <*45587> (signbit (rem (-0, 1)))
725 %!assert <*45587> (! signbit (rem (0, 1)))
726 
727 %!assert <*42627> (rem (0.94, 0.01), 0.0)
728 
729 %!error rem (uint (8), int8 (5))
730 %!error rem (uint8 ([1, 2]), uint8 ([3, 4, 5]))
731 %!error rem ()
732 %!error rem (1, 2, 3)
733 %!error rem ([1, 2], [3, 4, 5])
734 %!error rem (i, 1)
735 */
736 
737 DEFUN (mod, args, ,
738  doc: /* -*- texinfo -*-
739 @deftypefn {} {} mod (@var{x}, @var{y})
740 Compute the modulo of @var{x} and @var{y}.
741 
742 Conceptually this is given by
743 
744 @example
745 x - y .* floor (x ./ y)
746 @end example
747 
748 @noindent
749 and is written such that the correct modulus is returned for integer types.
750 This function handles negative values correctly. That is,
751 @w{@code{mod (-1, 3)}} is 2, not -1, as @w{@code{rem (-1, 3)}} returns.
752 
753 An error results if the dimensions of the arguments do not agree, or if
754 either of the arguments is complex.
755 
756 Programming Notes: Floating point numbers within a few eps of an integer
757 will be rounded to an integer before computation for compatibility with
758 @sc{matlab}.
759 
760 By convention,
761 
762 @example
763 @group
764 mod (@var{x}, 0) = @var{x}
765 mod (@var{x}, @var{y}) returns a value with the signbit from @var{y}
766 @end group
767 @end example
768 
769 For the opposite conventions see the @code{rem} function. In general,
770 @code{mod} is a better choice than @code{rem} when any of the inputs are
771 negative numbers or when the values are periodic.
772 @seealso{rem}
773 @end deftypefn */)
774 {
775  if (args.length () != 2)
776  print_usage ();
777 
779 
780  if (! args(0).isnumeric ())
781  err_wrong_type_arg ("mod", args(0));
782 
783  if (! args(1).isnumeric ())
784  err_wrong_type_arg ("mod", args(1));
785 
786  if (args(0).iscomplex () || args(1).iscomplex ())
787  error ("mod: not defined for complex numbers");
788 
789  if (args(0).isinteger () || args(1).isinteger ())
790  {
791  builtin_type_t btyp0 = args(0).builtin_type ();
792  builtin_type_t btyp1 = args(1).builtin_type ();
793  if (btyp0 == btyp_double || btyp0 == btyp_float)
794  btyp0 = btyp1;
795  if (btyp1 == btyp_double || btyp1 == btyp_float)
796  btyp1 = btyp0;
797 
798  if (btyp0 != btyp1)
799  error ("mod: cannot combine %s and %d",
800  args(0).class_name ().c_str (),
801  args(1).class_name ().c_str ());
802 
803  switch (btyp0)
804  {
805 #define MAKE_INT_BRANCH(X) \
806  case btyp_ ## X: \
807  { \
808  X##NDArray a0 = args(0).X##_array_value (); \
809  X##NDArray a1 = args(1).X##_array_value (); \
810  retval = binmap<octave_##X,octave_##X,octave_##X> (a0, a1, mod, "mod"); \
811  } \
812  break
813 
814  MAKE_INT_BRANCH (int8);
815  MAKE_INT_BRANCH (int16);
816  MAKE_INT_BRANCH (int32);
817  MAKE_INT_BRANCH (int64);
818  MAKE_INT_BRANCH (uint8);
819  MAKE_INT_BRANCH (uint16);
820  MAKE_INT_BRANCH (uint32);
821  MAKE_INT_BRANCH (uint64);
822 
823 #undef MAKE_INT_BRANCH
824 
825  default:
826  panic_impossible ();
827  }
828  }
829  else if (args(0).is_single_type () || args(1).is_single_type ())
830  {
831  if (args(0).is_scalar_type () && args(1).is_scalar_type ())
832  retval = octave::math::mod (args(0).float_value (), args(1).float_value ());
833  else
834  {
835  FloatNDArray a0 = args(0).float_array_value ();
836  FloatNDArray a1 = args(1).float_array_value ();
837  retval = binmap<float> (a0, a1, octave::math::mod<float>, "mod");
838  }
839  }
840  else
841  {
842  if (args(0).is_scalar_type () && args(1).is_scalar_type ())
843  retval = octave::math::mod (args(0).scalar_value (), args(1).scalar_value ());
844  else if (args(0).issparse () || args(1).issparse ())
845  {
846  SparseMatrix m0 = args(0).sparse_matrix_value ();
847  SparseMatrix m1 = args(1).sparse_matrix_value ();
848  retval = binmap<double> (m0, m1, octave::math::mod<double>, "mod");
849  }
850  else
851  {
852  NDArray a0 = args(0).array_value ();
853  NDArray a1 = args(1).array_value ();
854  retval = binmap<double> (a0, a1, octave::math::mod<double>, "mod");
855  }
856  }
857 
858  return retval;
859 }
860 
861 /*
862 ## empty input test
863 %!assert (isempty (mod ([], [])))
864 
865 ## x mod y, y != 0 tests
866 %!assert (mod (5, 3), 2)
867 %!assert (mod (-5, 3), 1)
868 %!assert (mod (0, 3), 0)
869 %!assert (mod ([-5, 5, 0], [3, 3, 3]), [1, 2, 0])
870 %!assert (mod ([-5; 5; 0], [3; 3; 3]), [1; 2; 0])
871 %!assert (mod ([-5, 5; 0, 3], [3, 3 ; 3, 1]), [1, 2 ; 0, 0])
872 
873 ## x mod 0 tests
874 %!assert (mod (5, 0), 5)
875 %!assert (mod (-5, 0), -5)
876 %!assert (mod ([-5, 5, 0], [3, 0, 3]), [1, 5, 0])
877 %!assert (mod ([-5; 5; 0], [3; 0; 3]), [1; 5; 0])
878 %!assert (mod ([-5, 5; 0, 3], [3, 0 ; 3, 1]), [1, 5 ; 0, 0])
879 %!assert (mod ([-5, 5; 0, 3], [0, 0 ; 0, 0]), [-5, 5; 0, 3])
880 
881 ## mixed scalar/matrix tests
882 %!assert (mod ([-5, 5; 0, 3], 0), [-5, 5; 0, 3])
883 %!assert (mod ([-5, 5; 0, 3], 3), [1, 2; 0, 0])
884 %!assert (mod (-5, [0,0; 0,0]), [-5, -5; -5, -5])
885 %!assert (mod (-5, [3,0; 3,1]), [1, -5; 1, 0])
886 %!assert (mod (-5, [3,2; 3,1]), [1, 1; 1, 0])
887 
888 ## integer types
889 %!assert (mod (uint8 (5), uint8 (4)), uint8 (1))
890 %!assert (mod (uint8 ([1:5]), uint8 (4)), uint8 ([1,2,3,0,1]))
891 %!assert (mod (uint8 ([1:5]), uint8 (0)), uint8 ([1:5]))
892 %!error (mod (uint8 (5), int8 (4)))
893 
894 ## mixed integer/real types
895 %!assert (mod (uint8 (5), 4), uint8 (1))
896 %!assert (mod (5, uint8 (4)), uint8 (1))
897 %!assert (mod (uint8 ([1:5]), 4), uint8 ([1,2,3,0,1]))
898 
899 ## non-integer real numbers
900 %!assert (mod (2.1, 0.1), 0)
901 %!assert (mod (2.1, 0.2), 0.1, eps)
902 
903 %!assert <*45587> (signbit (mod (-0, 0)))
904 %!assert <*45587> (! signbit (mod (0, -0)))
905 
906 %!assert <*42627> (mod (0.94, 0.01), 0.0)
907 */
908 
909 #define DATA_REDUCTION(FCN) \
910  \
911  int nargin = args.length (); \
912  \
913  if (nargin < 1 || nargin > 2) \
914  print_usage (); \
915  \
916  octave_value retval; \
917  \
918  octave_value arg = args(0); \
919  \
920  int dim = (nargin == 1 ? -1 : args(1).int_value (true) - 1); \
921  \
922  if (dim < -1) \
923  error (#FCN ": invalid dimension argument = %d", dim + 1); \
924  \
925  if (arg.isreal ()) \
926  { \
927  if (arg.issparse ()) \
928  { \
929  SparseMatrix tmp = arg.sparse_matrix_value (); \
930  \
931  retval = tmp.FCN (dim); \
932  } \
933  else if (arg.is_single_type ()) \
934  { \
935  FloatNDArray tmp = arg.float_array_value (); \
936  \
937  retval = tmp.FCN (dim); \
938  } \
939  else \
940  { \
941  NDArray tmp = arg.array_value (); \
942  \
943  retval = tmp.FCN (dim); \
944  } \
945  } \
946  else if (arg.iscomplex ()) \
947  { \
948  if (arg.issparse ()) \
949  { \
950  SparseComplexMatrix tmp = arg.sparse_complex_matrix_value (); \
951  \
952  retval = tmp.FCN (dim); \
953  } \
954  else if (arg.is_single_type ()) \
955  { \
956  FloatComplexNDArray tmp \
957  = arg.float_complex_array_value (); \
958  \
959  retval = tmp.FCN (dim); \
960  } \
961  else \
962  { \
963  ComplexNDArray tmp = arg.complex_array_value (); \
964  \
965  retval = tmp.FCN (dim); \
966  } \
967  } \
968  else \
969  err_wrong_type_arg (#FCN, arg); \
970  \
971  return retval
972 
973 DEFUN (cumprod, args, ,
974  doc: /* -*- texinfo -*-
975 @deftypefn {} {} cumprod (@var{x})
976 @deftypefnx {} {} cumprod (@var{x}, @var{dim})
977 Cumulative product of elements along dimension @var{dim}.
978 
979 If @var{dim} is omitted, it defaults to the first non-singleton dimension.
980 For example:
981 
982 @example
983 @group
984 cumprod ([1, 2; 3, 4; 5, 6])
985  @result{} 1 2
986  3 8
987  15 48
988 @end group
989 @end example
990 @seealso{prod, cumsum}
991 @end deftypefn */)
992 {
993  DATA_REDUCTION (cumprod);
994 }
995 
996 /*
997 %!assert (cumprod ([1, 2, 3]), [1, 2, 6])
998 %!assert (cumprod ([-1; -2; -3]), [-1; 2; -6])
999 %!assert (cumprod ([i, 2+i, -3+2i, 4]), [i, -1+2i, -1-8i, -4-32i])
1000 %!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])
1001 
1002 %!assert (cumprod (single ([1, 2, 3])), single ([1, 2, 6]))
1003 %!assert (cumprod (single ([-1; -2; -3])), single ([-1; 2; -6]))
1004 %!assert (cumprod (single ([i, 2+i, -3+2i, 4])), single ([i, -1+2i, -1-8i, -4-32i]))
1005 %!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]))
1006 
1007 %!assert (cumprod ([2, 3; 4, 5], 1), [2, 3; 8, 15])
1008 %!assert (cumprod ([2, 3; 4, 5], 2), [2, 6; 4, 20])
1009 
1010 %!assert (cumprod (single ([2, 3; 4, 5]), 1), single ([2, 3; 8, 15]))
1011 %!assert (cumprod (single ([2, 3; 4, 5]), 2), single ([2, 6; 4, 20]))
1012 
1013 %!error cumprod ()
1014 */
1015 
1016 DEFUN (cumsum, args, ,
1017  doc: /* -*- texinfo -*-
1018 @deftypefn {} {} cumsum (@var{x})
1019 @deftypefnx {} {} cumsum (@var{x}, @var{dim})
1020 @deftypefnx {} {} cumsum (@dots{}, "native")
1021 @deftypefnx {} {} cumsum (@dots{}, "double")
1022 Cumulative sum of elements along dimension @var{dim}.
1023 
1024 If @var{dim} is omitted, it defaults to the first non-singleton dimension.
1025 For example:
1026 
1027 @example
1028 @group
1029 cumsum ([1, 2; 3, 4; 5, 6])
1030  @result{} 1 2
1031  4 6
1032  9 12
1033 @end group
1034 @end example
1035 
1036 See @code{sum} for an explanation of the optional parameters @qcode{"native"}
1037 and @qcode{"double"}.
1038 @seealso{sum, cumprod}
1039 @end deftypefn */)
1040 {
1041  int nargin = args.length ();
1042 
1043  bool isnative = false;
1044  bool isdouble = false;
1045 
1046  if (nargin > 1 && args(nargin - 1).is_string ())
1047  {
1048  std::string str = args(nargin - 1).string_value ();
1049 
1050  if (str == "native")
1051  isnative = true;
1052  else if (str == "double")
1053  isdouble = true;
1054  else
1055  error ("cumsum: unrecognized string argument");
1056 
1057  nargin--;
1058  }
1059 
1060  if (nargin < 1 || nargin > 2)
1061  print_usage ();
1062 
1063  int dim = -1;
1064  if (nargin == 2)
1065  {
1066  dim = args(1).int_value () - 1;
1067  if (dim < 0)
1068  error ("cumsum: invalid dimension argument = %d", dim + 1);
1069  }
1070 
1072  octave_value arg = args(0);
1073 
1074  switch (arg.builtin_type ())
1075  {
1076  case btyp_double:
1077  if (arg.issparse ())
1078  retval = arg.sparse_matrix_value ().cumsum (dim);
1079  else
1080  retval = arg.array_value ().cumsum (dim);
1081  break;
1082  case btyp_complex:
1083  if (arg.issparse ())
1085  else
1086  retval = arg.complex_array_value ().cumsum (dim);
1087  break;
1088  case btyp_float:
1089  if (isdouble)
1090  retval = arg.array_value ().cumsum (dim);
1091  else
1092  retval = arg.float_array_value ().cumsum (dim);
1093  break;
1094  case btyp_float_complex:
1095  if (isdouble)
1096  retval = arg.complex_array_value ().cumsum (dim);
1097  else
1099  break;
1100 
1101 #define MAKE_INT_BRANCH(X) \
1102  case btyp_ ## X: \
1103  if (isnative) \
1104  retval = arg.X ## _array_value ().cumsum (dim); \
1105  else \
1106  retval = arg.array_value ().cumsum (dim); \
1107  break;
1108 
1109  MAKE_INT_BRANCH (int8);
1110  MAKE_INT_BRANCH (int16);
1111  MAKE_INT_BRANCH (int32);
1112  MAKE_INT_BRANCH (int64);
1113  MAKE_INT_BRANCH (uint8);
1114  MAKE_INT_BRANCH (uint16);
1115  MAKE_INT_BRANCH (uint32);
1116  MAKE_INT_BRANCH (uint64);
1117 
1118 #undef MAKE_INT_BRANCH
1119 
1120  case btyp_bool:
1121  if (arg.issparse ())
1122  {
1124  if (isnative)
1125  retval = cs != 0.0;
1126  else
1127  retval = cs;
1128  }
1129  else
1130  {
1131  NDArray cs = arg.array_value ().cumsum (dim);
1132  if (isnative)
1133  retval = cs != 0.0;
1134  else
1135  retval = cs;
1136  }
1137  break;
1138 
1139  default:
1140  err_wrong_type_arg ("cumsum", arg);
1141  }
1142 
1143  return retval;
1144 }
1145 
1146 /*
1147 %!assert (cumsum ([1, 2, 3]), [1, 3, 6])
1148 %!assert (cumsum ([-1; -2; -3]), [-1; -3; -6])
1149 %!assert (cumsum ([i, 2+i, -3+2i, 4]), [i, 2+2i, -1+4i, 3+4i])
1150 %!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])
1151 
1152 %!assert (cumsum (single ([1, 2, 3])), single ([1, 3, 6]))
1153 %!assert (cumsum (single ([-1; -2; -3])), single ([-1; -3; -6]))
1154 %!assert (cumsum (single ([i, 2+i, -3+2i, 4])), single ([i, 2+2i, -1+4i, 3+4i]))
1155 %!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]))
1156 
1157 %!assert (cumsum ([1, 2; 3, 4], 1), [1, 2; 4, 6])
1158 %!assert (cumsum ([1, 2; 3, 4], 2), [1, 3; 3, 7])
1159 
1160 %!assert (cumsum (single ([1, 2; 3, 4]), 1), single ([1, 2; 4, 6]))
1161 %!assert (cumsum (single ([1, 2; 3, 4]), 2), single ([1, 3; 3, 7]))
1162 
1163 %!error cumsum ()
1164 */
1165 
1166 DEFUN (diag, args, ,
1167  doc: /* -*- texinfo -*-
1168 @deftypefn {} {@var{M} =} diag (@var{v})
1169 @deftypefnx {} {@var{M} =} diag (@var{v}, @var{k})
1170 @deftypefnx {} {@var{M} =} diag (@var{v}, @var{m}, @var{n})
1171 @deftypefnx {} {@var{v} =} diag (@var{M})
1172 @deftypefnx {} {@var{v} =} diag (@var{M}, @var{k})
1173 Return a diagonal matrix with vector @var{v} on diagonal @var{k}.
1174 
1175 The second argument is optional. If it is positive, the vector is placed on
1176 the @var{k}-th superdiagonal. If it is negative, it is placed on the
1177 @var{-k}-th subdiagonal. The default value of @var{k} is 0, and the vector
1178 is placed on the main diagonal. For example:
1179 
1180 @example
1181 @group
1182 diag ([1, 2, 3], 1)
1183  @result{} 0 1 0 0
1184  0 0 2 0
1185  0 0 0 3
1186  0 0 0 0
1187 @end group
1188 @end example
1189 
1190 @noindent
1191 The 3-input form returns a diagonal matrix with vector @var{v} on the main
1192 diagonal and the resulting matrix being of size @var{m} rows x @var{n}
1193 columns.
1194 
1195 Given a matrix argument, instead of a vector, @code{diag} extracts the
1196 @var{k}-th diagonal of the matrix.
1197 @end deftypefn */)
1198 {
1199  int nargin = args.length ();
1200 
1201  if (nargin < 1 || nargin > 3)
1202  print_usage ();
1203 
1205 
1206  if (nargin == 1)
1207  retval = args(0).diag ();
1208  else if (nargin == 2)
1209  {
1210  octave_idx_type k = args(1).xidx_type_value ("diag: invalid argument K");
1211 
1212  retval = args(0).diag (k);
1213  }
1214  else
1215  {
1216  octave_value arg0 = args(0);
1217 
1218  if (arg0.ndims () != 2 || (arg0.rows () != 1 && arg0.columns () != 1))
1219  error ("diag: V must be a vector");
1220 
1221  octave_idx_type m = args(1).xidx_type_value ("diag: invalid dimensions");
1222  octave_idx_type n = args(2).xidx_type_value ("diag: invalid dimensions");
1223 
1224  retval = arg0.diag (m, n);
1225  }
1226 
1227  return retval;
1228 }
1229 
1230 /*
1231 
1232 %!assert (full (diag ([1; 2; 3])), [1, 0, 0; 0, 2, 0; 0, 0, 3])
1233 %!assert (diag ([1; 2; 3], 1), [0, 1, 0, 0; 0, 0, 2, 0; 0, 0, 0, 3; 0, 0, 0, 0])
1234 %!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])
1235 %!assert (diag ([1; 2; 3],-1), [0, 0, 0, 0; 1, 0, 0, 0; 0, 2, 0, 0; 0, 0, 3, 0])
1236 %!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])
1237 
1238 %!assert (diag ([1, 0, 0; 0, 2, 0; 0, 0, 3]), [1; 2; 3])
1239 %!assert (diag ([0, 1, 0, 0; 0, 0, 2, 0; 0, 0, 0, 3; 0, 0, 0, 0], 1), [1; 2; 3])
1240 %!assert (diag ([0, 0, 0, 0; 1, 0, 0, 0; 0, 2, 0, 0; 0, 0, 3, 0], -1), [1; 2; 3])
1241 %!assert (diag (ones (1, 0), 2), zeros (2))
1242 %!assert (diag (1:3, 4, 2), [1, 0; 0, 2; 0, 0; 0, 0])
1243 
1244 %!assert (full (diag (single ([1; 2; 3]))), single ([1, 0, 0; 0, 2, 0; 0, 0, 3]))
1245 %!assert (diag (single ([1; 2; 3]), 1), single ([0, 1, 0, 0; 0, 0, 2, 0; 0, 0, 0, 3; 0, 0, 0, 0]))
1246 %!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]))
1247 %!assert (diag (single ([1; 2; 3]),-1), single ([0, 0, 0, 0; 1, 0, 0, 0; 0, 2, 0, 0; 0, 0, 3, 0]))
1248 %!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]))
1249 
1250 %!assert (diag (single ([1, 0, 0; 0, 2, 0; 0, 0, 3])), single ([1; 2; 3]))
1251 %!assert (diag (single ([0, 1, 0, 0; 0, 0, 2, 0; 0, 0, 0, 3; 0, 0, 0, 0]), 1), single ([1; 2; 3]))
1252 %!assert (diag (single ([0, 0, 0, 0; 1, 0, 0, 0; 0, 2, 0, 0; 0, 0, 3, 0]), -1), single ([1; 2; 3]))
1253 
1254 %!assert (diag (int8 ([1; 2; 3])), int8 ([1, 0, 0; 0, 2, 0; 0, 0, 3]))
1255 %!assert (diag (int8 ([1; 2; 3]), 1), int8 ([0, 1, 0, 0; 0, 0, 2, 0; 0, 0, 0, 3; 0, 0, 0, 0]))
1256 %!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]))
1257 %!assert (diag (int8 ([1; 2; 3]),-1), int8 ([0, 0, 0, 0; 1, 0, 0, 0; 0, 2, 0, 0; 0, 0, 3, 0]))
1258 %!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]))
1259 
1260 %!assert (diag (int8 ([1, 0, 0; 0, 2, 0; 0, 0, 3])), int8 ([1; 2; 3]))
1261 %!assert (diag (int8 ([0, 1, 0, 0; 0, 0, 2, 0; 0, 0, 0, 3; 0, 0, 0, 0]), 1), int8 ([1; 2; 3]))
1262 %!assert (diag (int8 ([0, 0, 0, 0; 1, 0, 0, 0; 0, 2, 0, 0; 0, 0, 3, 0]), -1), int8 ([1; 2; 3]))
1263 
1264 %!assert <*37411> (diag (diag ([5, 2, 3])(:,1)), diag([5 0 0 ]))
1265 %!assert <*37411> (diag (diag ([5, 2, 3])(:,1), 2), [0 0 5 0 0; zeros(4, 5)])
1266 %!assert <*37411> (diag (diag ([5, 2, 3])(:,1), -2), [[0 0 5 0 0]', zeros(5, 4)])
1267 
1268 ## Test non-square size
1269 %!assert (diag ([1,2,3], 6, 3), [1 0 0; 0 2 0; 0 0 3; 0 0 0; 0 0 0; 0 0 0])
1270 %!assert (diag (1, 2, 3), [1,0,0; 0,0,0])
1271 %!assert (diag ({1}, 2, 3), {1,[],[]; [],[],[]})
1272 %!assert (diag ({1,2}, 3, 4), {1,[],[],[]; [],2,[],[]; [],[],[],[]})
1273 
1274 ## Test out-of-range diagonals
1275 %!assert (diag (ones (3,3), 4), zeros (0, 1))
1276 %!assert (diag (cell (3,3), 4), cell (0, 1))
1277 %!assert (diag (sparse (ones (3,3)), 4), sparse (zeros (0, 1)))
1278 
1279 ## Test input validation
1280 %!error <Invalid call to diag> diag ()
1281 %!error <Invalid call to diag> diag (1,2,3,4)
1282 %!error diag (ones (2), 3, 3)
1283 %!error diag (1:3, -4, 3)
1284 
1285 %!assert (diag (1, 3, 3), diag ([1, 0, 0]))
1286 %!assert (diag (i, 3, 3), diag ([i, 0, 0]))
1287 %!assert (diag (single (1), 3, 3), diag ([single(1), 0, 0]))
1288 %!assert (diag (single (i), 3, 3), diag ([single(i), 0, 0]))
1289 %!assert (diag ([1, 2], 3, 3), diag ([1, 2, 0]))
1290 %!assert (diag ([1, 2]*i, 3, 3), diag ([1, 2, 0]*i))
1291 %!assert (diag (single ([1, 2]), 3, 3), diag (single ([1, 2, 0])))
1292 %!assert (diag (single ([1, 2]*i), 3, 3), diag (single ([1, 2, 0]*i)))
1293 */
1294 
1295 DEFUN (prod, args, ,
1296  doc: /* -*- texinfo -*-
1297 @deftypefn {} {} prod (@var{x})
1298 @deftypefnx {} {} prod (@var{x}, @var{dim})
1299 @deftypefnx {} {} prod (@dots{}, "native")
1300 @deftypefnx {} {} prod (@dots{}, "double")
1301 Product of elements along dimension @var{dim}.
1302 
1303 If @var{dim} is omitted, it defaults to the first non-singleton dimension.
1304 
1305 The optional @qcode{"type"} input determines the class of the variable
1306 used for calculations. If the argument @qcode{"native"} is given, then
1307 the operation is performed in the same type as the original argument, rather
1308 than the default double type.
1309 
1310 For example:
1311 
1312 @example
1313 @group
1314 prod ([true, true])
1315  @result{} 1
1316 prod ([true, true], "native")
1317  @result{} true
1318 @end group
1319 @end example
1320 
1321 On the contrary, if @qcode{"double"} is given, the operation is performed
1322 in double precision even for single precision inputs.
1323 @seealso{cumprod, sum}
1324 @end deftypefn */)
1325 {
1326  int nargin = args.length ();
1327 
1328  bool isnative = false;
1329  bool isdouble = false;
1330 
1331  if (nargin > 1 && args(nargin - 1).is_string ())
1332  {
1333  std::string str = args(nargin - 1).string_value ();
1334 
1335  if (str == "native")
1336  isnative = true;
1337  else if (str == "double")
1338  isdouble = true;
1339  else
1340  error ("prod: unrecognized type argument '%s'", str.c_str ());
1341 
1342  nargin--;
1343  }
1344 
1345  if (nargin < 1 || nargin > 2)
1346  print_usage ();
1347 
1349 
1350  octave_value arg = args(0);
1351 
1352  int dim = -1;
1353  if (nargin == 2)
1354  {
1355  dim = args(1).int_value () - 1;
1356  if (dim < 0)
1357  error ("prod: invalid dimension DIM = %d", dim + 1);
1358  }
1359 
1360  switch (arg.builtin_type ())
1361  {
1362  case btyp_double:
1363  if (arg.issparse ())
1364  retval = arg.sparse_matrix_value ().prod (dim);
1365  else
1366  retval = arg.array_value ().prod (dim);
1367  break;
1368  case btyp_complex:
1369  if (arg.issparse ())
1371  else
1372  retval = arg.complex_array_value ().prod (dim);
1373  break;
1374  case btyp_float:
1375  if (isdouble)
1376  retval = arg.float_array_value ().dprod (dim);
1377  else
1378  retval = arg.float_array_value ().prod (dim);
1379  break;
1380  case btyp_float_complex:
1381  if (isdouble)
1383  else
1385  break;
1386 
1387 #define MAKE_INT_BRANCH(X) \
1388  case btyp_ ## X: \
1389  if (isnative) \
1390  retval = arg.X ## _array_value ().prod (dim); \
1391  else \
1392  retval = arg.array_value ().prod (dim); \
1393  break;
1394 
1395  MAKE_INT_BRANCH (int8);
1396  MAKE_INT_BRANCH (int16);
1397  MAKE_INT_BRANCH (int32);
1398  MAKE_INT_BRANCH (int64);
1399  MAKE_INT_BRANCH (uint8);
1400  MAKE_INT_BRANCH (uint16);
1401  MAKE_INT_BRANCH (uint32);
1402  MAKE_INT_BRANCH (uint64);
1403 
1404 #undef MAKE_INT_BRANCH
1405 
1406  // GAGME: Accursed Matlab compatibility...
1407  case btyp_char:
1408  retval = arg.array_value (true).prod (dim);
1409  break;
1410 
1411  case btyp_bool:
1412  if (arg.issparse ())
1413  {
1414  if (isnative)
1416  else
1417  retval = arg.sparse_matrix_value ().prod (dim);
1418  }
1419  else if (isnative)
1420  retval = arg.bool_array_value ().all (dim);
1421  else
1422  retval = NDArray (arg.bool_array_value ().all (dim));
1423  break;
1424 
1425  default:
1426  err_wrong_type_arg ("prod", arg);
1427  }
1428 
1429  return retval;
1430 }
1431 
1432 /*
1433 %!assert (prod ([1, 2, 3]), 6)
1434 %!assert (prod ([-1; -2; -3]), -6)
1435 %!assert (prod ([i, 2+i, -3+2i, 4]), -4 - 32i)
1436 %!assert (prod ([1, 2, 3; i, 2i, 3i; 1+i, 2+2i, 3+3i]), [-1+i, -8+8i, -27+27i])
1437 
1438 %!assert (prod (single ([1, 2, 3])), single (6))
1439 %!assert (prod (single ([-1; -2; -3])), single (-6))
1440 %!assert (prod (single ([i, 2+i, -3+2i, 4])), single (-4 - 32i))
1441 %!assert (prod (single ([1, 2, 3; i, 2i, 3i; 1+i, 2+2i, 3+3i])), single ([-1+i, -8+8i, -27+27i]))
1442 
1443 ## Test sparse
1444 %!assert (prod (sparse ([1, 2, 3])), sparse (6))
1445 %!assert (prod (sparse ([-1; -2; -3])), sparse (-6))
1446 ## Commented out until bug #42290 is fixed
1447 #%!assert (prod (sparse ([i, 2+i, -3+2i, 4])), sparse (-4 - 32i))
1448 #%!assert (prod (sparse ([1, 2, 3; i, 2i, 3i; 1+i, 2+2i, 3+3i])), sparse ([-1+i, -8+8i, -27+27i]))
1449 
1450 %!assert (prod ([1, 2; 3, 4], 1), [3, 8])
1451 %!assert (prod ([1, 2; 3, 4], 2), [2; 12])
1452 %!assert (prod (zeros (1, 0)), 1)
1453 %!assert (prod (zeros (1, 0), 1), zeros (1, 0))
1454 %!assert (prod (zeros (1, 0), 2), 1)
1455 %!assert (prod (zeros (0, 1)), 1)
1456 %!assert (prod (zeros (0, 1), 1), 1)
1457 %!assert (prod (zeros (0, 1), 2), zeros (0, 1))
1458 %!assert (prod (zeros (2, 0)), zeros (1, 0))
1459 %!assert (prod (zeros (2, 0), 1), zeros (1, 0))
1460 %!assert (prod (zeros (2, 0), 2), [1; 1])
1461 %!assert (prod (zeros (0, 2)), [1, 1])
1462 %!assert (prod (zeros (0, 2), 1), [1, 1])
1463 %!assert (prod (zeros (0, 2), 2), zeros (0, 1))
1464 
1465 %!assert (prod (single ([1, 2; 3, 4]), 1), single ([3, 8]))
1466 %!assert (prod (single ([1, 2; 3, 4]), 2), single ([2; 12]))
1467 %!assert (prod (zeros (1, 0, "single")), single (1))
1468 %!assert (prod (zeros (1, 0, "single"), 1), zeros (1, 0, "single"))
1469 %!assert (prod (zeros (1, 0, "single"), 2), single (1))
1470 %!assert (prod (zeros (0, 1, "single")), single (1))
1471 %!assert (prod (zeros (0, 1, "single"), 1), single (1))
1472 %!assert (prod (zeros (0, 1, "single"), 2), zeros (0, 1, "single"))
1473 %!assert (prod (zeros (2, 0, "single")), zeros (1, 0, "single"))
1474 %!assert (prod (zeros (2, 0, "single"), 1), zeros (1, 0, "single"))
1475 %!assert (prod (zeros (2, 0, "single"), 2), single ([1; 1]))
1476 %!assert (prod (zeros (0, 2, "single")), single ([1, 1]))
1477 %!assert (prod (zeros (0, 2, "single"), 1), single ([1, 1]))
1478 %!assert (prod (zeros (0, 2, "single"), 2), zeros (0, 1, "single"))
1479 
1480 ## Test "double" type argument
1481 %!assert (prod (single ([1, 2, 3]), "double"), 6)
1482 %!assert (prod (single ([-1; -2; -3]), "double"), -6)
1483 %!assert (prod (single ([i, 2+i, -3+2i, 4]), "double"), -4 - 32i)
1484 %!assert (prod (single ([1, 2, 3; i, 2i, 3i; 1+i, 2+2i, 3+3i]), "double"), [-1+i, -8+8i, -27+27i])
1485 
1486 ## Test "native" type argument
1487 %!assert (prod (uint8 ([1, 2, 3]), "native"), uint8 (6))
1488 %!assert (prod (uint8 ([-1; -2; -3]), "native"), uint8 (0))
1489 %!assert (prod (int8 ([1, 2, 3]), "native"), int8 (6))
1490 %!assert (prod (int8 ([-1; -2; -3]), "native"), int8 (-6))
1491 %!assert (prod ([true false; true true], "native"), [true false])
1492 %!assert (prod ([true false; true true], 2, "native"), [false; true])
1493 
1494 ## Test input validation
1495 %!error prod ()
1496 %!error prod (1,2,3)
1497 %!error <unrecognized type argument 'foobar'> prod (1, "foobar")
1498 */
1499 
1500 static bool
1501 all_scalar_1x1 (const octave_value_list& args)
1502 {
1503  int n_args = args.length ();
1504  for (int i = 0; i < n_args; i++)
1505  if (args(i).numel () != 1)
1506  return false;
1507 
1508  return true;
1509 }
1510 
1511 template <typename TYPE, typename T>
1512 static void
1514  const octave_value_list& args,
1515  int dim)
1516 {
1517  int n_args = args.length ();
1520  && all_scalar_1x1 (args))
1521  {
1522  // Optimize all scalars case.
1523  dim_vector dv (1, 1);
1524  if (dim == -1 || dim == -2)
1525  dim = -dim - 1;
1526  else if (dim >= 2)
1527  dv.resize (dim+1, 1);
1528  dv(dim) = n_args;
1529 
1530  result.clear (dv);
1531 
1532  for (int j = 0; j < n_args; j++)
1533  {
1534  octave_quit ();
1535 
1536  result(j) = octave_value_extract<T> (args(j));
1537  }
1538  }
1539  else
1540  {
1541  OCTAVE_LOCAL_BUFFER (Array<T>, array_list, n_args);
1542 
1543  for (int j = 0; j < n_args; j++)
1544  {
1545  octave_quit ();
1546 
1547  array_list[j] = octave_value_extract<TYPE> (args(j));
1548  }
1549 
1550  result = Array<T>::cat (dim, n_args, array_list);
1551  }
1552 }
1553 
1554 template <typename TYPE, typename T>
1555 static void
1557  const octave_value_list& args,
1558  int dim)
1559 {
1560  int n_args = args.length ();
1561  OCTAVE_LOCAL_BUFFER (Sparse<T>, sparse_list, n_args);
1562 
1563  for (int j = 0; j < n_args; j++)
1564  {
1565  octave_quit ();
1566 
1567  sparse_list[j] = octave_value_extract<TYPE> (args(j));
1568  }
1569 
1570  result = Sparse<T>::cat (dim, n_args, sparse_list);
1571 }
1572 
1573 // Dispatcher.
1574 template <typename TYPE>
1575 static TYPE
1576 do_single_type_concat (const octave_value_list& args, int dim)
1577 {
1578  TYPE result;
1579 
1580  single_type_concat<TYPE, typename TYPE::element_type> (result, args, dim);
1581 
1582  return result;
1583 }
1584 
1585 template <typename MAP>
1586 static void
1587 single_type_concat_map (octave_map& result,
1588  const octave_value_list& args,
1589  int dim)
1590 {
1591  int n_args = args.length ();
1592  OCTAVE_LOCAL_BUFFER (MAP, map_list, n_args);
1593 
1594  for (int j = 0; j < n_args; j++)
1595  {
1596  octave_quit ();
1597 
1598  map_list[j] = octave_value_extract<MAP> (args(j));
1599  }
1600 
1601  result = octave_map::cat (dim, n_args, map_list);
1602 }
1603 
1604 static octave_map
1605 do_single_type_concat_map (const octave_value_list& args,
1606  int dim)
1607 {
1609  if (all_scalar_1x1 (args)) // optimize all scalars case.
1610  single_type_concat_map<octave_scalar_map> (result, args, dim);
1611  else
1612  single_type_concat_map<octave_map> (result, args, dim);
1613 
1614  return result;
1615 }
1616 
1617 static octave_value
1618 attempt_type_conversion (const octave_value& ov, std::string dtype)
1619 {
1621 
1622  // First try to find function in the class of OV that can convert to
1623  // the dispatch type dtype. It will have the name of the dispatch
1624  // type.
1625 
1626  std::string cname = ov.class_name ();
1627 
1628  octave::symbol_table& symtab =
1629  octave::__get_symbol_table__ ("attempt_type_conversion");
1630 
1631  octave_value fcn = symtab.find_method (dtype, cname);
1632 
1633  if (fcn.is_defined ())
1634  {
1636 
1637  try
1638  {
1639  result = octave::feval (fcn, ovl (ov), 1);
1640  }
1641  catch (octave::execution_exception& e)
1642  {
1643  error (e, "conversion from %s to %s failed", dtype.c_str (),
1644  cname.c_str ());
1645  }
1646 
1647  if (result.empty ())
1648  error ("conversion from %s to %s failed", dtype.c_str (),
1649  cname.c_str ());
1650 
1651  retval = result(0);
1652  }
1653  else
1654  {
1655  // No conversion function available. Try the constructor for the
1656  // dispatch type.
1657 
1658  fcn = symtab.find_method (dtype, dtype);
1659 
1660  if (! fcn.is_defined ())
1661  error ("no constructor for %s!", dtype.c_str ());
1662 
1664 
1665  try
1666  {
1667  result = octave::feval (fcn, ovl (ov), 1);
1668  }
1669  catch (octave::execution_exception& e)
1670  {
1671  error (e, "%s constructor failed for %s argument", dtype.c_str (),
1672  cname.c_str ());
1673  }
1674 
1675  if (result.empty ())
1676  error ("%s constructor failed for %s argument", dtype.c_str (),
1677  cname.c_str ());
1678 
1679  retval = result(0);
1680  }
1681 
1682  return retval;
1683 }
1684 
1686 do_class_concat (const octave_value_list& ovl, std::string cattype, int dim)
1687 {
1689 
1690  // Get dominant type for list
1691 
1693 
1694  octave::symbol_table& symtab = octave::__get_symbol_table__ ("do_class_concat");
1695 
1696  octave_value fcn = symtab.find_method (cattype, dtype);
1697 
1698  if (fcn.is_defined ())
1699  {
1700  // Have method for dominant type. Call it and let it handle conversions.
1701 
1702  octave_value_list tmp2;
1703 
1704  try
1705  {
1706  tmp2 = octave::feval (fcn, ovl, 1);
1707  }
1708  catch (octave::execution_exception& e)
1709  {
1710  error (e, "%s/%s method failed", dtype.c_str (), cattype.c_str ());
1711  }
1712 
1713  if (tmp2.empty ())
1714  error ("%s/%s method did not return a value", dtype.c_str (),
1715  cattype.c_str ());
1716 
1717  retval = tmp2(0);
1718  }
1719  else
1720  {
1721  // No method for dominant type, so attempt type conversions for
1722  // all elements that are not of the dominant type, then do the
1723  // default operation for octave_class values.
1724 
1725  octave_idx_type j = 0;
1726  octave_idx_type len = ovl.length ();
1728  for (octave_idx_type k = 0; k < len; k++)
1729  {
1730  octave_value elt = ovl(k);
1731 
1732  std::string t1_type = elt.class_name ();
1733 
1734  if (t1_type == dtype)
1735  tmp(j++) = elt;
1736  else if (elt.isobject () || ! elt.isempty ())
1737  tmp(j++) = attempt_type_conversion (elt, dtype);
1738  }
1739 
1740  tmp.resize (j);
1741 
1742  octave_map m = do_single_type_concat_map (tmp, dim);
1743 
1744  std::string cname = tmp(0).class_name ();
1745  std::list<std::string> parents = tmp(0).parent_class_name_list ();
1746 
1747  retval = octave_value (new octave_class (m, cname, parents));
1748  }
1749 
1750  return retval;
1751 }
1752 
1753 static octave_value
1754 do_cat (const octave_value_list& xargs, int dim, std::string fname)
1755 {
1757 
1758  // We may need to convert elements of the list to cells, so make a copy.
1759  // This should be efficient, it is done mostly by incrementing reference
1760  // counts.
1761  octave_value_list args = xargs;
1762 
1763  int n_args = args.length ();
1764 
1765  if (n_args == 0)
1766  retval = Matrix ();
1767  else if (n_args == 1)
1768  retval = args(0);
1769  else if (n_args > 1)
1770  {
1772 
1773  bool all_strings_p = true;
1774  bool all_sq_strings_p = true;
1775  bool all_dq_strings_p = true;
1776  bool all_real_p = true;
1777  bool all_cmplx_p = true;
1778  bool any_sparse_p = false;
1779  bool any_cell_p = false;
1780  bool any_class_p = false;
1781 
1782  bool first_elem_is_struct = false;
1783 
1784  for (int i = 0; i < n_args; i++)
1785  {
1786  if (i == 0)
1787  {
1788  result_type = args(i).class_name ();
1789 
1790  first_elem_is_struct = args(i).isstruct ();
1791  }
1792  else
1793  result_type = octave::get_concat_class (result_type, args(i).class_name ());
1794 
1795  if (all_strings_p && ! args(i).is_string ())
1796  all_strings_p = false;
1797  if (all_sq_strings_p && ! args(i).is_sq_string ())
1798  all_sq_strings_p = false;
1799  if (all_dq_strings_p && ! args(i).is_dq_string ())
1800  all_dq_strings_p = false;
1801  if (all_real_p && ! args(i).isreal ())
1802  all_real_p = false;
1803  if (all_cmplx_p && ! (args(i).iscomplex ()
1804  || args(i).isreal ()))
1805  all_cmplx_p = false;
1806  if (! any_sparse_p && args(i).issparse ())
1807  any_sparse_p = true;
1808  if (! any_cell_p && args(i).iscell ())
1809  any_cell_p = true;
1810  if (! any_class_p && args(i).isobject ())
1811  any_class_p = true;
1812  }
1813 
1814  if (any_cell_p && ! any_class_p && ! first_elem_is_struct)
1815  {
1816  int j = 0;
1817  for (int i = 0; i < n_args; i++)
1818  {
1819  if (args(i).iscell ())
1820  args(j++) = args(i);
1821  else
1822  {
1823  if (args(i).isempty ())
1824  continue; // Delete empty non-cell arg
1825  else
1826  args(j++) = Cell (args(i));
1827  }
1828  }
1829  n_args = j;
1830  args.resize (n_args);
1831  }
1832 
1833  if (any_class_p)
1834  {
1835  retval = do_class_concat (args, fname, dim);
1836  }
1837  else if (result_type == "double")
1838  {
1839  if (any_sparse_p)
1840  {
1841  if (all_real_p)
1842  retval = do_single_type_concat<SparseMatrix> (args, dim);
1843  else
1844  retval = do_single_type_concat<SparseComplexMatrix> (args, dim);
1845  }
1846  else
1847  {
1848  if (all_real_p)
1849  retval = do_single_type_concat<NDArray> (args, dim);
1850  else
1851  retval = do_single_type_concat<ComplexNDArray> (args, dim);
1852  }
1853  }
1854  else if (result_type == "single")
1855  {
1856  if (all_real_p)
1857  retval = do_single_type_concat<FloatNDArray> (args, dim);
1858  else
1859  retval = do_single_type_concat<FloatComplexNDArray> (args, dim);
1860  }
1861  else if (result_type == "char")
1862  {
1863  char type = (all_dq_strings_p ? '"' : '\'');
1864 
1865  if (! all_strings_p)
1866  warn_implicit_conversion ("Octave:num-to-str",
1867  "numeric", result_type);
1868  else
1869  octave::maybe_warn_string_concat (all_dq_strings_p, all_sq_strings_p);
1870 
1871  charNDArray result = do_single_type_concat<charNDArray> (args, dim);
1872 
1874  }
1875  else if (result_type == "logical")
1876  {
1877  if (any_sparse_p)
1878  retval = do_single_type_concat<SparseBoolMatrix> (args, dim);
1879  else
1880  retval = do_single_type_concat<boolNDArray> (args, dim);
1881  }
1882  else if (result_type == "int8")
1883  retval = do_single_type_concat<int8NDArray> (args, dim);
1884  else if (result_type == "int16")
1885  retval = do_single_type_concat<int16NDArray> (args, dim);
1886  else if (result_type == "int32")
1887  retval = do_single_type_concat<int32NDArray> (args, dim);
1888  else if (result_type == "int64")
1889  retval = do_single_type_concat<int64NDArray> (args, dim);
1890  else if (result_type == "uint8")
1891  retval = do_single_type_concat<uint8NDArray> (args, dim);
1892  else if (result_type == "uint16")
1893  retval = do_single_type_concat<uint16NDArray> (args, dim);
1894  else if (result_type == "uint32")
1895  retval = do_single_type_concat<uint32NDArray> (args, dim);
1896  else if (result_type == "uint64")
1897  retval = do_single_type_concat<uint64NDArray> (args, dim);
1898  else if (result_type == "cell")
1899  retval = do_single_type_concat<Cell> (args, dim);
1900  else if (result_type == "struct")
1901  retval = do_single_type_concat_map (args, dim);
1902  else
1903  {
1904  dim_vector dv = args(0).dims ();
1905 
1906  // Default concatenation.
1907  bool (dim_vector::*concat_rule) (const dim_vector&, int)
1908  = &dim_vector::concat;
1909 
1910  if (dim == -1 || dim == -2)
1911  {
1912  concat_rule = &dim_vector::hvcat;
1913  dim = -dim - 1;
1914  }
1915 
1916  for (int i = 1; i < args.length (); i++)
1917  {
1918  if (! (dv.*concat_rule) (args(i).dims (), dim))
1919  error ("cat: dimension mismatch");
1920  }
1921 
1922  // The lines below might seem crazy, since we take a copy
1923  // of the first argument, resize it to be empty and then resize
1924  // it to be full. This is done since it means that there is no
1925  // recopying of data, as would happen if we used a single resize.
1926  // It should be noted that resize operation is also significantly
1927  // slower than the do_cat_op function, so it makes sense to have
1928  // an empty matrix and copy all data.
1929  //
1930  // We might also start with a empty octave_value using
1931  //
1932  // tmp = octave::type_info::lookup_type (args(1).type_name());
1933  //
1934  // and then directly resize. However, for some types there might
1935  // be some additional setup needed, and so this should be avoided.
1936 
1937  octave_value tmp = args(0);
1938  tmp = tmp.resize (dim_vector (0,0)).resize (dv);
1939 
1940  int dv_len = dv.ndims ();
1941  Array<octave_idx_type> ra_idx (dim_vector (dv_len, 1), 0);
1942 
1943  for (int j = 0; j < n_args; j++)
1944  {
1945  // Can't fast return here to skip empty matrices as something
1946  // like cat (1,[],single ([])) must return an empty matrix of
1947  // the right type.
1948  tmp = do_cat_op (tmp, args(j), ra_idx);
1949 
1950  dim_vector dv_tmp = args(j).dims ();
1951 
1952  if (dim >= dv_len)
1953  {
1954  if (j > 1)
1955  error ("%s: indexing error", fname.c_str ());
1956 
1957  break;
1958  }
1959  else
1960  ra_idx(dim) += (dim < dv_tmp.ndims () ? dv_tmp(dim) : 1);
1961  }
1962  retval = tmp;
1963  }
1964  }
1965  else
1966  print_usage ();
1967 
1968  return retval;
1969 }
1970 
1971 DEFUN (horzcat, args, ,
1972  doc: /* -*- texinfo -*-
1973 @deftypefn {} {} horzcat (@var{array1}, @var{array2}, @dots{}, @var{arrayN})
1974 Return the horizontal concatenation of N-D array objects, @var{array1},
1975 @var{array2}, @dots{}, @var{arrayN} along dimension 2.
1976 
1977 Arrays may also be concatenated horizontally using the syntax for creating
1978 new matrices. For example:
1979 
1980 @example
1981 @var{hcat} = [ @var{array1}, @var{array2}, @dots{} ]
1982 @end example
1983 @seealso{cat, vertcat}
1984 @end deftypefn */)
1985 {
1986  return do_cat (args, -2, "horzcat");
1987 }
1988 
1989 /*
1990 ## Test concatenation with all zero matrices
1991 %!test
1992 %! warning ("off", "Octave:num-to-str", "local");
1993 %! assert (horzcat ("", 65*ones (1,10)), "AAAAAAAAAA");
1994 %! assert (horzcat (65*ones (1,10), ""), "AAAAAAAAAA");
1995 
1996 %!assert (class (horzcat (int64 (1), int64 (1))), "int64")
1997 %!assert (class (horzcat (int64 (1), int32 (1))), "int64")
1998 %!assert (class (horzcat (int64 (1), int16 (1))), "int64")
1999 %!assert (class (horzcat (int64 (1), int8 (1))), "int64")
2000 %!assert (class (horzcat (int64 (1), uint64 (1))), "int64")
2001 %!assert (class (horzcat (int64 (1), uint32 (1))), "int64")
2002 %!assert (class (horzcat (int64 (1), uint16 (1))), "int64")
2003 %!assert (class (horzcat (int64 (1), uint8 (1))), "int64")
2004 %!assert (class (horzcat (int64 (1), single (1))), "int64")
2005 %!assert (class (horzcat (int64 (1), double (1))), "int64")
2006 %!assert (class (horzcat (int64 (1), cell (1))), "cell")
2007 %!assert (class (horzcat (int64 (1), true)), "int64")
2008 %!test
2009 %! warning ("off", "Octave:num-to-str", "local");
2010 %! assert (class (horzcat (int64 (1), "a")), "char");
2011 
2012 %!assert (class (horzcat (int32 (1), int64 (1))), "int32")
2013 %!assert (class (horzcat (int32 (1), int32 (1))), "int32")
2014 %!assert (class (horzcat (int32 (1), int16 (1))), "int32")
2015 %!assert (class (horzcat (int32 (1), int8 (1))), "int32")
2016 %!assert (class (horzcat (int32 (1), uint64 (1))), "int32")
2017 %!assert (class (horzcat (int32 (1), uint32 (1))), "int32")
2018 %!assert (class (horzcat (int32 (1), uint16 (1))), "int32")
2019 %!assert (class (horzcat (int32 (1), uint8 (1))), "int32")
2020 %!assert (class (horzcat (int32 (1), single (1))), "int32")
2021 %!assert (class (horzcat (int32 (1), double (1))), "int32")
2022 %!assert (class (horzcat (int32 (1), cell (1))), "cell")
2023 %!assert (class (horzcat (int32 (1), true)), "int32")
2024 %!test
2025 %! warning ("off", "Octave:num-to-str", "local");
2026 %! assert (class (horzcat (int32 (1), "a")), "char");
2027 
2028 %!assert (class (horzcat (int16 (1), int64 (1))), "int16")
2029 %!assert (class (horzcat (int16 (1), int32 (1))), "int16")
2030 %!assert (class (horzcat (int16 (1), int16 (1))), "int16")
2031 %!assert (class (horzcat (int16 (1), int8 (1))), "int16")
2032 %!assert (class (horzcat (int16 (1), uint64 (1))), "int16")
2033 %!assert (class (horzcat (int16 (1), uint32 (1))), "int16")
2034 %!assert (class (horzcat (int16 (1), uint16 (1))), "int16")
2035 %!assert (class (horzcat (int16 (1), uint8 (1))), "int16")
2036 %!assert (class (horzcat (int16 (1), single (1))), "int16")
2037 %!assert (class (horzcat (int16 (1), double (1))), "int16")
2038 %!assert (class (horzcat (int16 (1), cell (1))), "cell")
2039 %!assert (class (horzcat (int16 (1), true)), "int16")
2040 %!test
2041 %! warning ("off", "Octave:num-to-str", "local");
2042 %! assert (class (horzcat (int16 (1), "a")), "char");
2043 
2044 %!assert (class (horzcat (int8 (1), int64 (1))), "int8")
2045 %!assert (class (horzcat (int8 (1), int32 (1))), "int8")
2046 %!assert (class (horzcat (int8 (1), int16 (1))), "int8")
2047 %!assert (class (horzcat (int8 (1), int8 (1))), "int8")
2048 %!assert (class (horzcat (int8 (1), uint64 (1))), "int8")
2049 %!assert (class (horzcat (int8 (1), uint32 (1))), "int8")
2050 %!assert (class (horzcat (int8 (1), uint16 (1))), "int8")
2051 %!assert (class (horzcat (int8 (1), uint8 (1))), "int8")
2052 %!assert (class (horzcat (int8 (1), single (1))), "int8")
2053 %!assert (class (horzcat (int8 (1), double (1))), "int8")
2054 %!assert (class (horzcat (int8 (1), cell (1))), "cell")
2055 %!assert (class (horzcat (int8 (1), true)), "int8")
2056 %!test
2057 %! warning ("off", "Octave:num-to-str", "local");
2058 %! assert (class (horzcat (int8 (1), "a")), "char");
2059 
2060 %!assert (class (horzcat (uint64 (1), int64 (1))), "uint64")
2061 %!assert (class (horzcat (uint64 (1), int32 (1))), "uint64")
2062 %!assert (class (horzcat (uint64 (1), int16 (1))), "uint64")
2063 %!assert (class (horzcat (uint64 (1), int8 (1))), "uint64")
2064 %!assert (class (horzcat (uint64 (1), uint64 (1))), "uint64")
2065 %!assert (class (horzcat (uint64 (1), uint32 (1))), "uint64")
2066 %!assert (class (horzcat (uint64 (1), uint16 (1))), "uint64")
2067 %!assert (class (horzcat (uint64 (1), uint8 (1))), "uint64")
2068 %!assert (class (horzcat (uint64 (1), single (1))), "uint64")
2069 %!assert (class (horzcat (uint64 (1), double (1))), "uint64")
2070 %!assert (class (horzcat (uint64 (1), cell (1))), "cell")
2071 %!assert (class (horzcat (uint64 (1), true)), "uint64")
2072 %!test
2073 %! warning ("off", "Octave:num-to-str", "local");
2074 %! assert (class (horzcat (uint64 (1), "a")), "char");
2075 
2076 %!assert (class (horzcat (uint32 (1), int64 (1))), "uint32")
2077 %!assert (class (horzcat (uint32 (1), int32 (1))), "uint32")
2078 %!assert (class (horzcat (uint32 (1), int16 (1))), "uint32")
2079 %!assert (class (horzcat (uint32 (1), int8 (1))), "uint32")
2080 %!assert (class (horzcat (uint32 (1), uint64 (1))), "uint32")
2081 %!assert (class (horzcat (uint32 (1), uint32 (1))), "uint32")
2082 %!assert (class (horzcat (uint32 (1), uint16 (1))), "uint32")
2083 %!assert (class (horzcat (uint32 (1), uint8 (1))), "uint32")
2084 %!assert (class (horzcat (uint32 (1), single (1))), "uint32")
2085 %!assert (class (horzcat (uint32 (1), double (1))), "uint32")
2086 %!assert (class (horzcat (uint32 (1), cell (1))), "cell")
2087 %!assert (class (horzcat (uint32 (1), true)), "uint32")
2088 %!test
2089 %! warning ("off", "Octave:num-to-str", "local");
2090 %! assert (class (horzcat (uint32 (1), "a")), "char");
2091 
2092 %!assert (class (horzcat (uint16 (1), int64 (1))), "uint16")
2093 %!assert (class (horzcat (uint16 (1), int32 (1))), "uint16")
2094 %!assert (class (horzcat (uint16 (1), int16 (1))), "uint16")
2095 %!assert (class (horzcat (uint16 (1), int8 (1))), "uint16")
2096 %!assert (class (horzcat (uint16 (1), uint64 (1))), "uint16")
2097 %!assert (class (horzcat (uint16 (1), uint32 (1))), "uint16")
2098 %!assert (class (horzcat (uint16 (1), uint16 (1))), "uint16")
2099 %!assert (class (horzcat (uint16 (1), uint8 (1))), "uint16")
2100 %!assert (class (horzcat (uint16 (1), single (1))), "uint16")
2101 %!assert (class (horzcat (uint16 (1), double (1))), "uint16")
2102 %!assert (class (horzcat (uint16 (1), cell (1))), "cell")
2103 %!assert (class (horzcat (uint16 (1), true)), "uint16")
2104 %!test
2105 %! warning ("off", "Octave:num-to-str", "local");
2106 %! assert (class (horzcat (uint16 (1), "a")), "char");
2107 
2108 %!assert (class (horzcat (uint8 (1), int64 (1))), "uint8")
2109 %!assert (class (horzcat (uint8 (1), int32 (1))), "uint8")
2110 %!assert (class (horzcat (uint8 (1), int16 (1))), "uint8")
2111 %!assert (class (horzcat (uint8 (1), int8 (1))), "uint8")
2112 %!assert (class (horzcat (uint8 (1), uint64 (1))), "uint8")
2113 %!assert (class (horzcat (uint8 (1), uint32 (1))), "uint8")
2114 %!assert (class (horzcat (uint8 (1), uint16 (1))), "uint8")
2115 %!assert (class (horzcat (uint8 (1), uint8 (1))), "uint8")
2116 %!assert (class (horzcat (uint8 (1), single (1))), "uint8")
2117 %!assert (class (horzcat (uint8 (1), double (1))), "uint8")
2118 %!assert (class (horzcat (uint8 (1), cell (1))), "cell")
2119 %!assert (class (horzcat (uint8 (1), true)), "uint8")
2120 %!test
2121 %! warning ("off", "Octave:num-to-str", "local");
2122 %! assert (class (horzcat (uint8 (1), "a")), "char");
2123 
2124 %!assert (class (horzcat (single (1), int64 (1))), "int64")
2125 %!assert (class (horzcat (single (1), int32 (1))), "int32")
2126 %!assert (class (horzcat (single (1), int16 (1))), "int16")
2127 %!assert (class (horzcat (single (1), int8 (1))), "int8")
2128 %!assert (class (horzcat (single (1), uint64 (1))), "uint64")
2129 %!assert (class (horzcat (single (1), uint32 (1))), "uint32")
2130 %!assert (class (horzcat (single (1), uint16 (1))), "uint16")
2131 %!assert (class (horzcat (single (1), uint8 (1))), "uint8")
2132 %!assert (class (horzcat (single (1), single (1))), "single")
2133 %!assert (class (horzcat (single (1), double (1))), "single")
2134 %!assert (class (horzcat (single (1), cell (1))), "cell")
2135 %!assert (class (horzcat (single (1), true)), "single")
2136 %!test
2137 %! warning ("off", "Octave:num-to-str", "local");
2138 %! assert (class (horzcat (single (1), "a")), "char");
2139 
2140 %!assert (class (horzcat (double (1), int64 (1))), "int64")
2141 %!assert (class (horzcat (double (1), int32 (1))), "int32")
2142 %!assert (class (horzcat (double (1), int16 (1))), "int16")
2143 %!assert (class (horzcat (double (1), int8 (1))), "int8")
2144 %!assert (class (horzcat (double (1), uint64 (1))), "uint64")
2145 %!assert (class (horzcat (double (1), uint32 (1))), "uint32")
2146 %!assert (class (horzcat (double (1), uint16 (1))), "uint16")
2147 %!assert (class (horzcat (double (1), uint8 (1))), "uint8")
2148 %!assert (class (horzcat (double (1), single (1))), "single")
2149 %!assert (class (horzcat (double (1), double (1))), "double")
2150 %!assert (class (horzcat (double (1), cell (1))), "cell")
2151 %!assert (class (horzcat (double (1), true)), "double")
2152 %!test
2153 %! warning ("off", "Octave:num-to-str", "local");
2154 %! assert (class (horzcat (double (1), "a")), "char");
2155 
2156 %!assert (class (horzcat (cell (1), int64 (1))), "cell")
2157 %!assert (class (horzcat (cell (1), int32 (1))), "cell")
2158 %!assert (class (horzcat (cell (1), int16 (1))), "cell")
2159 %!assert (class (horzcat (cell (1), int8 (1))), "cell")
2160 %!assert (class (horzcat (cell (1), uint64 (1))), "cell")
2161 %!assert (class (horzcat (cell (1), uint32 (1))), "cell")
2162 %!assert (class (horzcat (cell (1), uint16 (1))), "cell")
2163 %!assert (class (horzcat (cell (1), uint8 (1))), "cell")
2164 %!assert (class (horzcat (cell (1), single (1))), "cell")
2165 %!assert (class (horzcat (cell (1), double (1))), "cell")
2166 %!assert (class (horzcat (cell (1), cell (1))), "cell")
2167 %!assert (class (horzcat (cell (1), true)), "cell")
2168 %!assert (class (horzcat (cell (1), "a")), "cell")
2169 
2170 %!assert (class (horzcat (true, int64 (1))), "int64")
2171 %!assert (class (horzcat (true, int32 (1))), "int32")
2172 %!assert (class (horzcat (true, int16 (1))), "int16")
2173 %!assert (class (horzcat (true, int8 (1))), "int8")
2174 %!assert (class (horzcat (true, uint64 (1))), "uint64")
2175 %!assert (class (horzcat (true, uint32 (1))), "uint32")
2176 %!assert (class (horzcat (true, uint16 (1))), "uint16")
2177 %!assert (class (horzcat (true, uint8 (1))), "uint8")
2178 %!assert (class (horzcat (true, single (1))), "single")
2179 %!assert (class (horzcat (true, double (1))), "double")
2180 %!assert (class (horzcat (true, cell (1))), "cell")
2181 %!assert (class (horzcat (true, true)), "logical")
2182 %!test
2183 %! warning ("off", "Octave:num-to-str", "local");
2184 %! assert (class (horzcat (true, "a")), "char");
2185 
2186 %!test
2187 %! warning ("off", "Octave:num-to-str", "local");
2188 %! assert (class (horzcat ("a", int64 (1))), "char");
2189 %! assert (class (horzcat ("a", int32 (1))), "char");
2190 %! assert (class (horzcat ("a", int16 (1))), "char");
2191 %! assert (class (horzcat ("a", int8 (1))), "char");
2192 %! assert (class (horzcat ("a", int64 (1))), "char");
2193 %! assert (class (horzcat ("a", int32 (1))), "char");
2194 %! assert (class (horzcat ("a", int16 (1))), "char");
2195 %! assert (class (horzcat ("a", int8 (1))), "char");
2196 %! assert (class (horzcat ("a", single (1))), "char");
2197 %! assert (class (horzcat ("a", double (1))), "char");
2198 %! assert (class (horzcat ("a", cell (1))), "cell");
2199 %! assert (class (horzcat ("a", true)), "char");
2200 %! assert (class (horzcat ("a", "a")), "char");
2201 
2202 %!assert (class (horzcat (cell (1), struct ("foo", "bar"))), "cell")
2203 
2204 %!error horzcat (struct ("foo", "bar"), cell (1))
2205 
2206 %!test <*39041> assert (class (horzcat (cell(0), struct())), "cell")
2207 %!test <51086> assert (class (horzcat (struct(), cell(0))), "struct")
2208 */
2209 
2210 DEFUN (vertcat, args, ,
2211  doc: /* -*- texinfo -*-
2212 @deftypefn {} {} vertcat (@var{array1}, @var{array2}, @dots{}, @var{arrayN})
2213 Return the vertical concatenation of N-D array objects, @var{array1},
2214 @var{array2}, @dots{}, @var{arrayN} along dimension 1.
2215 
2216 Arrays may also be concatenated vertically using the syntax for creating
2217 new matrices. For example:
2218 
2219 @example
2220 @var{vcat} = [ @var{array1}; @var{array2}; @dots{} ]
2221 @end example
2222 @seealso{cat, horzcat}
2223 @end deftypefn */)
2224 {
2225  return do_cat (args, -1, "vertcat");
2226 }
2227 
2228 /*
2229 %!test
2230 %! c = {"foo"; "bar"; "bazoloa"};
2231 %! assert (vertcat (c, "a", "bc", "def"), {"foo"; "bar"; "bazoloa"; "a"; "bc"; "def"});
2232 */
2233 
2234 DEFUN (cat, args, ,
2235  doc: /* -*- texinfo -*-
2236 @deftypefn {} {} cat (@var{dim}, @var{array1}, @var{array2}, @dots{}, @var{arrayN})
2237 Return the concatenation of N-D array objects, @var{array1},
2238 @var{array2}, @dots{}, @var{arrayN} along dimension @var{dim}.
2239 
2240 @example
2241 @group
2242 A = ones (2, 2);
2243 B = zeros (2, 2);
2244 cat (2, A, B)
2245  @result{} 1 1 0 0
2246  1 1 0 0
2247 @end group
2248 @end example
2249 
2250 Alternatively, we can concatenate @var{A} and @var{B} along the
2251 second dimension in the following way:
2252 
2253 @example
2254 @group
2255 [A, B]
2256 @end group
2257 @end example
2258 
2259 @var{dim} can be larger than the dimensions of the N-D array objects
2260 and the result will thus have @var{dim} dimensions as the
2261 following example shows:
2262 
2263 @example
2264 @group
2265 cat (4, ones (2, 2), zeros (2, 2))
2266  @result{} ans(:,:,1,1) =
2267 
2268  1 1
2269  1 1
2270 
2271  ans(:,:,1,2) =
2272 
2273  0 0
2274  0 0
2275 @end group
2276 @end example
2277 @seealso{horzcat, vertcat}
2278 @end deftypefn */)
2279 {
2280  if (args.length () == 0)
2281  print_usage ();
2282 
2283  int dim = args(0).xint_value ("cat: DIM must be an integer") - 1;
2284 
2285  if (dim < 0)
2286  error ("cat: DIM must be a valid dimension");
2287 
2288  return ovl (do_cat (args.slice (1, args.length () - 1), dim, "cat"));
2289 }
2290 
2291 /*
2292 %!function ret = __testcat (t1, t2, tr, cmplx)
2293 %! assert (cat (1, cast ([], t1), cast ([], t2)), cast ([], tr));
2294 %!
2295 %! assert (cat (1, cast (1, t1), cast (2, t2)), cast ([1; 2], tr));
2296 %! assert (cat (1, cast (1, t1), cast ([2; 3], t2)), cast ([1; 2; 3], tr));
2297 %! assert (cat (1, cast ([1; 2], t1), cast (3, t2)), cast ([1; 2; 3], tr));
2298 %! assert (cat (1, cast ([1; 2], t1), cast ([3; 4], t2)), cast ([1; 2; 3; 4], tr));
2299 %! assert (cat (2, cast (1, t1), cast (2, t2)), cast ([1, 2], tr));
2300 %! assert (cat (2, cast (1, t1), cast ([2, 3], t2)), cast ([1, 2, 3], tr));
2301 %! assert (cat (2, cast ([1, 2], t1), cast (3, t2)), cast ([1, 2, 3], tr));
2302 %! assert (cat (2, cast ([1, 2], t1), cast ([3, 4], t2)), cast ([1, 2, 3, 4], tr));
2303 %!
2304 %! assert ([cast(1, t1); cast(2, t2)], cast ([1; 2], tr));
2305 %! assert ([cast(1, t1); cast([2; 3], t2)], cast ([1; 2; 3], tr));
2306 %! assert ([cast([1; 2], t1); cast(3, t2)], cast ([1; 2; 3], tr));
2307 %! assert ([cast([1; 2], t1); cast([3; 4], t2)], cast ([1; 2; 3; 4], tr));
2308 %! assert ([cast(1, t1), cast(2, t2)], cast ([1, 2], tr));
2309 %! assert ([cast(1, t1), cast([2, 3], t2)], cast ([1, 2, 3], tr));
2310 %! assert ([cast([1, 2], t1), cast(3, t2)], cast ([1, 2, 3], tr));
2311 %! assert ([cast([1, 2], t1), cast([3, 4], t2)], cast ([1, 2, 3, 4], tr));
2312 %!
2313 %! if (nargin == 3 || cmplx)
2314 %! assert (cat (1, cast (1i, t1), cast (2, t2)), cast ([1i; 2], tr));
2315 %! assert (cat (1, cast (1i, t1), cast ([2; 3], t2)), cast ([1i; 2; 3], tr));
2316 %! assert (cat (1, cast ([1i; 2], t1), cast (3, t2)), cast ([1i; 2; 3], tr));
2317 %! assert (cat (1, cast ([1i; 2], t1), cast ([3; 4], t2)), cast ([1i; 2; 3; 4], tr));
2318 %! assert (cat (2, cast (1i, t1), cast (2, t2)), cast ([1i, 2], tr));
2319 %! assert (cat (2, cast (1i, t1), cast ([2, 3], t2)), cast ([1i, 2, 3], tr));
2320 %! assert (cat (2, cast ([1i, 2], t1), cast (3, t2)), cast ([1i, 2, 3], tr));
2321 %! assert (cat (2, cast ([1i, 2], t1), cast ([3, 4], t2)), cast ([1i, 2, 3, 4], tr));
2322 %!
2323 %! assert ([cast(1i, t1); cast(2, t2)], cast ([1i; 2], tr));
2324 %! assert ([cast(1i, t1); cast([2; 3], t2)], cast ([1i; 2; 3], tr));
2325 %! assert ([cast([1i; 2], t1); cast(3, t2)], cast ([1i; 2; 3], tr));
2326 %! assert ([cast([1i; 2], t1); cast([3; 4], t2)], cast ([1i; 2; 3; 4], tr));
2327 %! assert ([cast(1i, t1), cast(2, t2)], cast ([1i, 2], tr));
2328 %! assert ([cast(1i, t1), cast([2, 3], t2)], cast ([1i, 2, 3], tr));
2329 %! assert ([cast([1i, 2], t1), cast(3, t2)], cast ([1i, 2, 3], tr));
2330 %! assert ([cast([1i, 2], t1), cast([3, 4], t2)], cast ([1i, 2, 3, 4], tr));
2331 %!
2332 %! assert (cat (1, cast (1, t1), cast (2i, t2)), cast ([1; 2i], tr));
2333 %! assert (cat (1, cast (1, t1), cast ([2i; 3], t2)), cast ([1; 2i; 3], tr));
2334 %! assert (cat (1, cast ([1; 2], t1), cast (3i, t2)), cast ([1; 2; 3i], tr));
2335 %! assert (cat (1, cast ([1; 2], t1), cast ([3i; 4], t2)), cast ([1; 2; 3i; 4], tr));
2336 %! assert (cat (2, cast (1, t1), cast (2i, t2)), cast ([1, 2i], tr));
2337 %! assert (cat (2, cast (1, t1), cast ([2i, 3], t2)), cast ([1, 2i, 3], tr));
2338 %! assert (cat (2, cast ([1, 2], t1), cast (3i, t2)), cast ([1, 2, 3i], tr));
2339 %! assert (cat (2, cast ([1, 2], t1), cast ([3i, 4], t2)), cast ([1, 2, 3i, 4], tr));
2340 %!
2341 %! assert ([cast(1, t1); cast(2i, t2)], cast ([1; 2i], tr));
2342 %! assert ([cast(1, t1); cast([2i; 3], t2)], cast ([1; 2i; 3], tr));
2343 %! assert ([cast([1; 2], t1); cast(3i, t2)], cast ([1; 2; 3i], tr));
2344 %! assert ([cast([1; 2], t1); cast([3i; 4], t2)], cast ([1; 2; 3i; 4], tr));
2345 %! assert ([cast(1, t1), cast(2i, t2)], cast ([1, 2i], tr));
2346 %! assert ([cast(1, t1), cast([2i, 3], t2)], cast ([1, 2i, 3], tr));
2347 %! assert ([cast([1, 2], t1), cast(3i, t2)], cast ([1, 2, 3i], tr));
2348 %! assert ([cast([1, 2], t1), cast([3i, 4], t2)], cast ([1, 2, 3i, 4], tr));
2349 %!
2350 %! assert (cat (1, cast (1i, t1), cast (2i, t2)), cast ([1i; 2i], tr));
2351 %! assert (cat (1, cast (1i, t1), cast ([2i; 3], t2)), cast ([1i; 2i; 3], tr));
2352 %! assert (cat (1, cast ([1i; 2], t1), cast (3i, t2)), cast ([1i; 2; 3i], tr));
2353 %! assert (cat (1, cast ([1i; 2], t1), cast ([3i; 4], t2)), cast ([1i; 2; 3i; 4], tr));
2354 %! assert (cat (2, cast (1i, t1), cast (2i, t2)), cast ([1i, 2i], tr));
2355 %! assert (cat (2, cast (1i, t1), cast ([2i, 3], t2)), cast ([1i, 2i, 3], tr));
2356 %! assert (cat (2, cast ([1i, 2], t1), cast (3i, t2)), cast ([1i, 2, 3i], tr));
2357 %! assert (cat (2, cast ([1i, 2], t1), cast ([3i, 4], t2)), cast ([1i, 2, 3i, 4], tr));
2358 %!
2359 %! assert ([cast(1i, t1); cast(2i, t2)], cast ([1i; 2i], tr));
2360 %! assert ([cast(1i, t1); cast([2i; 3], t2)], cast ([1i; 2i; 3], tr));
2361 %! assert ([cast([1i; 2], t1); cast(3i, t2)], cast ([1i; 2; 3i], tr));
2362 %! assert ([cast([1i; 2], t1); cast([3i; 4], t2)], cast ([1i; 2; 3i; 4], tr));
2363 %! assert ([cast(1i, t1), cast(2i, t2)], cast ([1i, 2i], tr));
2364 %! assert ([cast(1i, t1), cast([2i, 3], t2)], cast ([1i, 2i, 3], tr));
2365 %! assert ([cast([1i, 2], t1), cast(3i, t2)], cast ([1i, 2, 3i], tr));
2366 %! assert ([cast([1i, 2], t1), cast([3i, 4], t2)], cast ([1i, 2, 3i, 4], tr));
2367 %! endif
2368 %! ret = true;
2369 %!endfunction
2370 
2371 %!assert (__testcat ("double", "double", "double"))
2372 %!assert (__testcat ("single", "double", "single"))
2373 %!assert (__testcat ("double", "single", "single"))
2374 %!assert (__testcat ("single", "single", "single"))
2375 
2376 %!assert (__testcat ("double", "int8", "int8", false))
2377 %!assert (__testcat ("int8", "double", "int8", false))
2378 %!assert (__testcat ("single", "int8", "int8", false))
2379 %!assert (__testcat ("int8", "single", "int8", false))
2380 %!assert (__testcat ("int8", "int8", "int8", false))
2381 %!assert (__testcat ("double", "int16", "int16", false))
2382 %!assert (__testcat ("int16", "double", "int16", false))
2383 %!assert (__testcat ("single", "int16", "int16", false))
2384 %!assert (__testcat ("int16", "single", "int16", false))
2385 %!assert (__testcat ("int16", "int16", "int16", false))
2386 %!assert (__testcat ("double", "int32", "int32", false))
2387 %!assert (__testcat ("int32", "double", "int32", false))
2388 %!assert (__testcat ("single", "int32", "int32", false))
2389 %!assert (__testcat ("int32", "single", "int32", false))
2390 %!assert (__testcat ("int32", "int32", "int32", false))
2391 %!assert (__testcat ("double", "int64", "int64", false))
2392 %!assert (__testcat ("int64", "double", "int64", false))
2393 %!assert (__testcat ("single", "int64", "int64", false))
2394 %!assert (__testcat ("int64", "single", "int64", false))
2395 %!assert (__testcat ("int64", "int64", "int64", false))
2396 
2397 %!assert (__testcat ("double", "uint8", "uint8", false))
2398 %!assert (__testcat ("uint8", "double", "uint8", false))
2399 %!assert (__testcat ("single", "uint8", "uint8", false))
2400 %!assert (__testcat ("uint8", "single", "uint8", false))
2401 %!assert (__testcat ("uint8", "uint8", "uint8", false))
2402 %!assert (__testcat ("double", "uint16", "uint16", false))
2403 %!assert (__testcat ("uint16", "double", "uint16", false))
2404 %!assert (__testcat ("single", "uint16", "uint16", false))
2405 %!assert (__testcat ("uint16", "single", "uint16", false))
2406 %!assert (__testcat ("uint16", "uint16", "uint16", false))
2407 %!assert (__testcat ("double", "uint32", "uint32", false))
2408 %!assert (__testcat ("uint32", "double", "uint32", false))
2409 %!assert (__testcat ("single", "uint32", "uint32", false))
2410 %!assert (__testcat ("uint32", "single", "uint32", false))
2411 %!assert (__testcat ("uint32", "uint32", "uint32", false))
2412 %!assert (__testcat ("double", "uint64", "uint64", false))
2413 %!assert (__testcat ("uint64", "double", "uint64", false))
2414 %!assert (__testcat ("single", "uint64", "uint64", false))
2415 %!assert (__testcat ("uint64", "single", "uint64", false))
2416 %!assert (__testcat ("uint64", "uint64", "uint64", false))
2417 
2418 %!assert (cat (3, [], [1,2;3,4]), [1,2;3,4])
2419 %!assert (cat (3, [1,2;3,4], []), [1,2;3,4])
2420 %!assert (cat (3, [], [1,2;3,4], []), [1,2;3,4])
2421 %!assert (cat (3, [], [], []), zeros (0, 0, 3))
2422 
2423 %!assert (cat (3, [], [], 1, 2), cat (3, 1, 2))
2424 %!assert (cat (3, [], [], [1,2;3,4]), [1,2;3,4])
2425 %!assert (cat (4, [], [], [1,2;3,4]), [1,2;3,4])
2426 
2427 %!assert ([zeros(3,2,2); ones(1,2,2)], repmat ([0;0;0;1],[1,2,2]))
2428 %!assert ([zeros(3,2,2); ones(1,2,2)], vertcat (zeros (3,2,2), ones (1,2,2)))
2429 
2430 %!test <*49759>
2431 %! A = [];
2432 %! B = {1; 2};
2433 %! assert (cat (1, A, B), {1; 2});
2434 %! assert (cat (2, A, B), {1; 2});
2435 
2436 %!error <dimension mismatch> cat (3, cat (3, [], []), [1,2;3,4])
2437 %!error <dimension mismatch> cat (3, zeros (0, 0, 2), [1,2;3,4])
2438 */
2439 
2440 static octave_value
2441 do_permute (const octave_value_list& args, bool inv)
2442 {
2443  if (args.length () != 2 || args(1).length () < args(1).ndims ())
2444  print_usage ();
2445 
2446  Array<int> vec = args(1).int_vector_value ();
2447 
2448  // FIXME: maybe we should create an idx_vector object here
2449  // and pass that to permute?
2450  int n = vec.numel ();
2451  for (int i = 0; i < n; i++)
2452  vec(i)--;
2453 
2454  return octave_value (args(0).permute (vec, inv));
2455 }
2456 
2457 DEFUN (permute, args, ,
2458  doc: /* -*- texinfo -*-
2459 @deftypefn {} {} permute (@var{A}, @var{perm})
2460 Return the generalized transpose for an N-D array object @var{A}.
2461 
2462 The permutation vector @var{perm} must contain the elements
2463 @w{@code{1:ndims (A)}} (in any order, but each element must appear only
2464 once). The @var{N}th dimension of @var{A} gets remapped to dimension
2465 @code{@var{PERM}(@var{N})}. For example:
2466 
2467 @example
2468 @group
2469 @var{x} = zeros ([2, 3, 5, 7]);
2470 size (@var{x})
2471  @result{} 2 3 5 7
2472 
2473 size (permute (@var{x}, [2, 1, 3, 4]))
2474  @result{} 3 2 5 7
2475 
2476 size (permute (@var{x}, [1, 3, 4, 2]))
2477  @result{} 2 5 7 3
2478 
2479 ## The identity permutation
2480 size (permute (@var{x}, [1, 2, 3, 4]))
2481  @result{} 2 3 5 7
2482 @end group
2483 @end example
2484 @seealso{ipermute}
2485 @end deftypefn */)
2486 {
2487  return do_permute (args, false);
2488 }
2489 
2490 DEFUN (ipermute, args, ,
2491  doc: /* -*- texinfo -*-
2492 @deftypefn {} {} ipermute (@var{A}, @var{iperm})
2493 The inverse of the @code{permute} function.
2494 
2495 The expression
2496 
2497 @example
2498 ipermute (permute (A, perm), perm)
2499 @end example
2500 
2501 @noindent
2502 returns the original array @var{A}.
2503 @seealso{permute}
2504 @end deftypefn */)
2505 {
2506  return do_permute (args, true);
2507 }
2508 
2509 DEFUN (length, args, ,
2510  doc: /* -*- texinfo -*-
2511 @deftypefn {} {} length (@var{a})
2512 Return the length of the object @var{a}.
2513 
2514 The length is 0 for empty objects, 1 for scalars, and the number of elements
2515 for vectors. For matrix or N-dimensional objects, the length is the number
2516 of elements along the largest dimension
2517 (equivalent to @w{@code{max (size (@var{a}))}}).
2518 @seealso{numel, size}
2519 @end deftypefn */)
2520 {
2521  if (args.length () != 1)
2522  print_usage ();
2523 
2524  return ovl (args(0).length ());
2525 }
2526 
2527 DEFUN (ndims, args, ,
2528  doc: /* -*- texinfo -*-
2529 @deftypefn {} {} ndims (@var{a})
2530 Return the number of dimensions of @var{a}.
2531 
2532 For any array, the result will always be greater than or equal to 2.
2533 Trailing singleton dimensions are not counted.
2534 
2535 @example
2536 @group
2537 ndims (ones (4, 1, 2, 1))
2538  @result{} 3
2539 @end group
2540 @end example
2541 @seealso{size}
2542 @end deftypefn */)
2543 {
2544  if (args.length () != 1)
2545  print_usage ();
2546 
2547  return ovl (args(0).ndims ());
2548 }
2549 
2550 DEFUN (numel, args, ,
2551  doc: /* -*- texinfo -*-
2552 @deftypefn {} {} numel (@var{a})
2553 @deftypefnx {} {} numel (@var{a}, @var{idx1}, @var{idx2}, @dots{})
2554 Return the number of elements in the object @var{a}.
2555 
2556 Optionally, if indices @var{idx1}, @var{idx2}, @dots{} are supplied,
2557 return the number of elements that would result from the indexing
2558 
2559 @example
2560 @var{a}(@var{idx1}, @var{idx2}, @dots{})
2561 @end example
2562 
2563 Note that the indices do not have to be scalar numbers. For example,
2564 
2565 @example
2566 @group
2567 @var{a} = 1;
2568 @var{b} = ones (2, 3);
2569 numel (@var{a}, @var{b})
2570 @end group
2571 @end example
2572 
2573 @noindent
2574 will return 6, as this is the number of ways to index with @var{b}.
2575 Or the index could be the string @qcode{":"} which represents the colon
2576 operator. For example,
2577 
2578 @example
2579 @group
2580 @var{a} = ones (5, 3);
2581 numel (@var{a}, 2, ":")
2582 @end group
2583 @end example
2584 
2585 @noindent
2586 will return 3 as the second row has three column entries.
2587 
2588 This method is also called when an object appears as lvalue with cs-list
2589 indexing, i.e., @code{object@{@dots{}@}} or @code{object(@dots{}).field}.
2590 @seealso{size, length, ndims}
2591 @end deftypefn */)
2592 {
2593  int nargin = args.length ();
2594 
2595  if (nargin == 0)
2596  print_usage ();
2597 
2599 
2600  if (nargin == 1)
2601  retval = args(0).numel ();
2602  else if (nargin > 1)
2603  {
2604  // Don't use numel (const octave_value_list&) here as that corresponds to
2605  // an overloaded call, not to builtin!
2606  retval = dims_to_numel (args(0).dims (), args.slice (1, nargin-1));
2607  }
2608 
2609  return retval;
2610 }
2611 
2612 DEFUN (size, args, nargout,
2613  doc: /* -*- texinfo -*-
2614 @deftypefn {} {@var{sz} =} size (@var{a})
2615 @deftypefnx {} {@var{dim_sz} =} size (@var{a}, @var{dim})
2616 @deftypefnx {} {[@var{rows}, @var{cols}, @dots{}, @var{dim_N_sz}] =} size (@dots{})
2617 Return a row vector with the size (number of elements) of each dimension for
2618 the object @var{a}.
2619 
2620 When given a second argument, @var{dim}, return the size of the corresponding
2621 dimension.
2622 
2623 With a single output argument, @code{size} returns a row vector. When called
2624 with multiple output arguments, @code{size} returns the size of dimension N
2625 in the Nth argument. The number of rows, dimension 1, is returned in the
2626 first argument, the number of columns, dimension 2, is returned in the
2627 second argument, etc. If there are more dimensions in @var{a} than there are
2628 output arguments, @code{size} returns the total number of elements in the
2629 remaining dimensions in the final output argument.
2630 
2631 Example 1: single row vector output
2632 
2633 @example
2634 @group
2635 size ([1, 2; 3, 4; 5, 6])
2636  @result{} [ 3, 2 ]
2637 @end group
2638 @end example
2639 
2640 Example 2: number of elements in 2nd dimension (columns)
2641 
2642 @example
2643 @group
2644 size ([1, 2; 3, 4; 5, 6], 2)
2645  @result{} 2
2646 @end group
2647 @end example
2648 
2649 Example 3: number of output arguments == number of dimensions
2650 
2651 @example
2652 @group
2653 [nr, nc] = size ([1, 2; 3, 4; 5, 6])
2654  @result{} nr = 3
2655  @result{} nc = 2
2656 @end group
2657 @end example
2658 
2659 Example 4: number of output arguments < number of dimensions
2660 
2661 @example
2662 @group
2663 [nr, remainder] = size (ones (2, 3, 4, 5))
2664  @result{} nr = 2
2665  @result{} remainder = 60
2666 @end group
2667 @end example
2668 
2669 @seealso{numel, ndims, length, rows, columns, size_equal, common_size}
2670 @end deftypefn */)
2671 {
2673 
2674  int nargin = args.length ();
2675 
2676  if (nargin == 1)
2677  {
2678  const dim_vector dimensions = args(0).dims ();
2679 
2680  if (nargout > 1)
2681  {
2682  const dim_vector rdims = dimensions.redim (nargout);
2683  retval.resize (nargout);
2684  for (int i = 0; i < nargout; i++)
2685  retval(i) = rdims(i);
2686  }
2687  else
2688  {
2689  int ndims = dimensions.ndims ();
2690 
2691  NoAlias<Matrix> m (1, ndims);
2692 
2693  for (int i = 0; i < ndims; i++)
2694  m(i) = dimensions(i);
2695 
2696  retval(0) = m;
2697  }
2698  }
2699  else if (nargin == 2 && nargout < 2)
2700  {
2701  if (! args(1).is_real_scalar ())
2702  error ("size: DIM must be a positive integer");
2703 
2704  octave_idx_type nd = args(1).idx_type_value ();
2705 
2706  const dim_vector dv = args(0).dims ();
2707 
2708  if (nd < 1)
2709  error ("size: requested dimension DIM (= %d) out of range", nd);
2710 
2711  if (nd <= dv.ndims ())
2712  retval(0) = dv(nd-1);
2713  else
2714  retval(0) = 1;
2715  }
2716  else
2717  print_usage ();
2718 
2719  return retval;
2720 }
2721 
2722 DEFUN (size_equal, args, ,
2723  doc: /* -*- texinfo -*-
2724 @deftypefn {} {} size_equal (@var{a}, @var{b}, @dots{})
2725 Return true if the dimensions of all arguments agree.
2726 
2727 Trailing singleton dimensions are ignored. When called with a single argument,
2728 or no argument, @code{size_equal} returns true.
2729 @seealso{size, numel, ndims, common_size}
2730 @end deftypefn */)
2731 {
2732  int nargin = args.length ();
2733 
2734  if (nargin >= 1)
2735  {
2736  dim_vector a_dims = args(0).dims ();
2737 
2738  for (int i = 1; i < nargin; ++i)
2739  {
2740  dim_vector b_dims = args(i).dims ();
2741 
2742  if (a_dims != b_dims)
2743  return ovl (false);
2744  }
2745  }
2746 
2747  return ovl (true);
2748 }
2749 
2750 DEFUN (nnz, args, ,
2751  doc: /* -*- texinfo -*-
2752 @deftypefn {} {@var{n} =} nnz (@var{a})
2753 Return the number of nonzero elements in @var{a}.
2754 @seealso{nzmax, nonzeros, find}
2755 @end deftypefn */)
2756 {
2757  if (args.length () != 1)
2758  print_usage ();
2759 
2760  return ovl (args(0).nnz ());
2761 }
2762 
2763 /*
2764 %!assert (nnz (1:5), 5)
2765 %!assert (nnz (-5:-1), 5)
2766 %!assert (nnz (0:5), 5)
2767 %!assert (nnz (-5:0), 5)
2768 %!assert (nnz (-5:5), 10)
2769 %!assert (nnz (-2:1:2), 4)
2770 %!assert (nnz (-2+eps(2):1:2), 5)
2771 %!assert (nnz (-2-eps(2):1:2), 5)
2772 %!assert (nnz (-2:1+eps(1):2), 5)
2773 %!assert (nnz (-2:1-eps(1):2), 5)
2774 %!assert (nnz ([1:5] * 0), 0)
2775 %!assert (nnz ([-5:-1] * 0), 0)
2776 %!assert (nnz ([-1:1] * 0), 0)
2777 */
2778 
2779 DEFUN (nzmax, args, ,
2780  doc: /* -*- texinfo -*-
2781 @deftypefn {} {@var{n} =} nzmax (@var{SM})
2782 Return the amount of storage allocated to the sparse matrix @var{SM}.
2783 
2784 Note that Octave tends to crop unused memory at the first opportunity
2785 for sparse objects. Thus, in general the value of @code{nzmax} will be the
2786 same as @code{nnz} except for some cases of user-created sparse objects.
2787 @seealso{nnz, spalloc, sparse}
2788 @end deftypefn */)
2789 {
2790  if (args.length () != 1)
2791  print_usage ();
2792 
2793  return ovl (args(0).nzmax ());
2794 }
2795 
2796 DEFUN (rows, args, ,
2797  doc: /* -*- texinfo -*-
2798 @deftypefn {} {} rows (@var{a})
2799 Return the number of rows of @var{a}.
2800 @seealso{columns, size, length, numel, isscalar, isvector, ismatrix}
2801 @end deftypefn */)
2802 {
2803  if (args.length () != 1)
2804  print_usage ();
2805 
2806  return ovl (args(0).rows ());
2807 }
2808 
2809 /*
2810 %!assert (rows (ones (2,5)), 2)
2811 %!assert (rows (ones (5,2)), 5)
2812 %!assert (rows (ones (5,4,3,2)), 5)
2813 %!assert (rows (ones (3,4,5,2)), 3)
2814 
2815 %!assert (rows (cell (2,5)), 2)
2816 %!assert (rows (cell (5,2)), 5)
2817 %!assert (rows (cell (5,4,3,2)), 5)
2818 %!assert (rows (cell (3,4,5,2)), 3)
2819 
2820 %!test
2821 %! x(2,5,3).a = 1;
2822 %! assert (rows (x), 2);
2823 %! y(5,4,3).b = 2;
2824 %! assert (rows (y), 5);
2825 
2826 %!assert (rows ("Hello World"), 1)
2827 
2828 %!assert (rows ([]), 0)
2829 %!assert (rows (zeros (2,0)), 2)
2830 
2831 ## Test input validation
2832 %!error rows ()
2833 %!error rows (1,2)
2834 */
2835 
2836 DEFUN (columns, args, ,
2837  doc: /* -*- texinfo -*-
2838 @deftypefn {} {} columns (@var{a})
2839 Return the number of columns of @var{a}.
2840 @seealso{rows, size, length, numel, isscalar, isvector, ismatrix}
2841 @end deftypefn */)
2842 {
2843  if (args.length () != 1)
2844  print_usage ();
2845 
2846  return ovl (args(0).columns ());
2847 }
2848 
2849 DEFUN (sum, args, ,
2850  doc: /* -*- texinfo -*-
2851 @deftypefn {} {} sum (@var{x})
2852 @deftypefnx {} {} sum (@var{x}, @var{dim})
2853 @deftypefnx {} {} sum (@dots{}, "native")
2854 @deftypefnx {} {} sum (@dots{}, "double")
2855 @deftypefnx {} {} sum (@dots{}, "extra")
2856 Sum of elements along dimension @var{dim}.
2857 
2858 If @var{dim} is omitted, it defaults to the first non-singleton dimension.
2859 
2860 The optional @qcode{"type"} input determines the class of the variable
2861 used for calculations. By default, operations on floating point inputs (double
2862 or single) are performed in their native data type, while operations on
2863 integer, logical, and character data types are performed using doubles. If the
2864 argument @qcode{"native"} is given, then the operation is performed in the same
2865 type as the original argument.
2866 
2867 For example:
2868 
2869 @example
2870 @group
2871 sum ([true, true])
2872  @result{} 2
2873 sum ([true, true], "native")
2874  @result{} true
2875 @end group
2876 @end example
2877 
2878 If @qcode{"double"} is given the sum is performed in double precision even for
2879 single precision inputs.
2880 
2881 For double precision inputs, the @qcode{"extra"} option will use a more
2882 accurate algorithm than straightforward summation. For single precision
2883 inputs, @qcode{"extra"} is the same as @qcode{"double"}. For all other data
2884 type @qcode{"extra"} has no effect.
2885 @seealso{cumsum, sumsq, prod}
2886 @end deftypefn */)
2887 {
2888  int nargin = args.length ();
2889 
2890  bool isnative = false;
2891  bool isdouble = false;
2892  bool isextra = false;
2893 
2894  if (nargin > 1 && args(nargin - 1).is_string ())
2895  {
2896  std::string str = args(nargin - 1).string_value ();
2897 
2898  if (str == "native")
2899  isnative = true;
2900  else if (str == "double")
2901  isdouble = true;
2902  else if (str == "extra")
2903  isextra = true;
2904  else
2905  error ("sum: unrecognized type argument '%s'", str.c_str ());
2906 
2907  nargin--;
2908  }
2909 
2910  if (nargin < 1 || nargin > 2)
2911  print_usage ();
2912 
2913  int dim = -1;
2914  if (nargin == 2)
2915  {
2916  dim = args(1).int_value () - 1;
2917  if (dim < 0)
2918  error ("sum: invalid dimension DIM = %d", dim + 1);
2919  }
2920 
2922  octave_value arg = args(0);
2923 
2924  switch (arg.builtin_type ())
2925  {
2926  case btyp_double:
2927  if (arg.issparse ())
2928  {
2929  if (isextra)
2930  warning ("sum: 'extra' not yet implemented for sparse matrices");
2931  retval = arg.sparse_matrix_value ().sum (dim);
2932  }
2933  else if (isextra)
2934  retval = arg.array_value ().xsum (dim);
2935  else
2936  retval = arg.array_value ().sum (dim);
2937  break;
2938 
2939  case btyp_complex:
2940  if (arg.issparse ())
2941  {
2942  if (isextra)
2943  warning ("sum: 'extra' not yet implemented for sparse matrices");
2945  }
2946  else if (isextra)
2947  retval = arg.complex_array_value ().xsum (dim);
2948  else
2949  retval = arg.complex_array_value ().sum (dim);
2950  break;
2951 
2952  case btyp_float:
2953  if (isdouble || isextra)
2954  retval = arg.float_array_value ().dsum (dim);
2955  else
2956  retval = arg.float_array_value ().sum (dim);
2957  break;
2958 
2959  case btyp_float_complex:
2960  if (isdouble || isextra)
2962  else
2964  break;
2965 
2966 #define MAKE_INT_BRANCH(X) \
2967  case btyp_ ## X: \
2968  if (isnative) \
2969  retval = arg.X ## _array_value ().sum (dim); \
2970  else \
2971  retval = arg.X ## _array_value ().dsum (dim); \
2972  break;
2973 
2974  MAKE_INT_BRANCH (int8);
2975  MAKE_INT_BRANCH (int16);
2976  MAKE_INT_BRANCH (int32);
2977  MAKE_INT_BRANCH (int64);
2978  MAKE_INT_BRANCH (uint8);
2979  MAKE_INT_BRANCH (uint16);
2980  MAKE_INT_BRANCH (uint32);
2981  MAKE_INT_BRANCH (uint64);
2982 
2983 #undef MAKE_INT_BRANCH
2984 
2985  // GAGME: Accursed Matlab compatibility...
2986  case btyp_char:
2987  if (isextra)
2988  retval = arg.array_value (true).xsum (dim);
2989  else
2990  retval = arg.array_value (true).sum (dim);
2991  break;
2992 
2993  case btyp_bool:
2994  if (arg.issparse ())
2995  {
2996  if (isnative)
2998  else
3000  }
3001  else if (isnative)
3002  retval = arg.bool_array_value ().any (dim);
3003  else
3004  retval = arg.array_value ().sum (dim);
3005  break;
3006 
3007  default:
3008  err_wrong_type_arg ("sum", arg);
3009  }
3010 
3011  return retval;
3012 }
3013 
3014 /*
3015 %!assert (sum ([1, 2, 3]), 6)
3016 %!assert (sum ([-1; -2; -3]), -6)
3017 %!assert (sum ([i, 2+i, -3+2i, 4]), 3+4i)
3018 %!assert (sum ([1, 2, 3; i, 2i, 3i; 1+i, 2+2i, 3+3i]), [2+2i, 4+4i, 6+6i])
3019 
3020 %!assert (sum (single ([1, 2, 3])), single (6))
3021 %!assert (sum (single ([-1; -2; -3])), single (-6))
3022 %!assert (sum (single ([i, 2+i, -3+2i, 4])), single (3+4i))
3023 %!assert (sum (single ([1, 2, 3; i, 2i, 3i; 1+i, 2+2i, 3+3i])), single ([2+2i, 4+4i, 6+6i]))
3024 
3025 %!assert (sum ([1, 2; 3, 4], 1), [4, 6])
3026 %!assert (sum ([1, 2; 3, 4], 2), [3; 7])
3027 %!assert (sum (zeros (1, 0)), 0)
3028 %!assert (sum (zeros (1, 0), 1), zeros (1, 0))
3029 %!assert (sum (zeros (1, 0), 2), 0)
3030 %!assert (sum (zeros (0, 1)), 0)
3031 %!assert (sum (zeros (0, 1), 1), 0)
3032 %!assert (sum (zeros (0, 1), 2), zeros (0, 1))
3033 %!assert (sum (zeros (2, 0)), zeros (1, 0))
3034 %!assert (sum (zeros (2, 0), 1), zeros (1, 0))
3035 %!assert (sum (zeros (2, 0), 2), [0; 0])
3036 %!assert (sum (zeros (0, 2)), [0, 0])
3037 %!assert (sum (zeros (0, 2), 1), [0, 0])
3038 %!assert (sum (zeros (0, 2), 2), zeros (0, 1))
3039 %!assert (sum (zeros (2, 2, 0, 3)), zeros (1, 2, 0, 3))
3040 %!assert (sum (zeros (2, 2, 0, 3), 2), zeros (2, 1, 0, 3))
3041 %!assert (sum (zeros (2, 2, 0, 3), 3), zeros (2, 2, 1, 3))
3042 %!assert (sum (zeros (2, 2, 0, 3), 4), zeros (2, 2, 0))
3043 %!assert (sum (zeros (2, 2, 0, 3), 7), zeros (2, 2, 0, 3))
3044 
3045 %!assert (sum (single ([1, 2; 3, 4]), 1), single ([4, 6]))
3046 %!assert (sum (single ([1, 2; 3, 4]), 2), single ([3; 7]))
3047 %!assert (sum (zeros (1, 0, "single")), single (0))
3048 %!assert (sum (zeros (1, 0, "single"), 1), zeros (1, 0, "single"))
3049 %!assert (sum (zeros (1, 0, "single"), 2), single (0))
3050 %!assert (sum (zeros (0, 1, "single")), single (0))
3051 %!assert (sum (zeros (0, 1, "single"), 1), single (0))
3052 %!assert (sum (zeros (0, 1, "single"), 2), zeros (0, 1, "single"))
3053 %!assert (sum (zeros (2, 0, "single")), zeros (1, 0, "single"))
3054 %!assert (sum (zeros (2, 0, "single"), 1), zeros (1, 0, "single"))
3055 %!assert (sum (zeros (2, 0, "single"), 2), single ([0; 0]))
3056 %!assert (sum (zeros (0, 2, "single")), single ([0, 0]))
3057 %!assert (sum (zeros (0, 2, "single"), 1), single ([0, 0]))
3058 %!assert (sum (zeros (0, 2, "single"), 2), zeros (0, 1, "single"))
3059 %!assert (sum (zeros (2, 2, 0, 3, "single")), zeros (1, 2, 0, 3, "single"))
3060 %!assert (sum (zeros (2, 2, 0, 3, "single"), 2), zeros (2, 1, 0, 3, "single"))
3061 %!assert (sum (zeros (2, 2, 0, 3, "single"), 3), zeros (2, 2, 1, 3, "single"))
3062 %!assert (sum (zeros (2, 2, 0, 3, "single"), 4), zeros (2, 2, 0, "single"))
3063 %!assert (sum (zeros (2, 2, 0, 3, "single"), 7), zeros (2, 2, 0, 3, "single"))
3064 
3065 ## Test "native"
3066 %!assert (sum ([true,true]), 2)
3067 %!assert (sum ([true,true], "native"), true)
3068 %!assert (sum (int8 ([127,10,-20])), 117)
3069 %!assert (sum (int8 ([127,10,-20]), "native"), int8 (107))
3070 
3071 ;-)
3072 %!assert (sum ("Octave") + "8", sumsq (primes (17)))
3073 
3074 %!error sum ()
3075 %!error sum (1,2,3)
3076 %!error <unrecognized type argument 'foobar'> sum (1, "foobar")
3077 */
3078 
3079 DEFUN (sumsq, args, ,
3080  doc: /* -*- texinfo -*-
3081 @deftypefn {} {} sumsq (@var{x})
3082 @deftypefnx {} {} sumsq (@var{x}, @var{dim})
3083 Sum of squares of elements along dimension @var{dim}.
3084 
3085 If @var{dim} is omitted, it defaults to the first non-singleton dimension.
3086 
3087 This function is conceptually equivalent to computing
3088 
3089 @example
3090 sum (x .* conj (x), dim)
3091 @end example
3092 
3093 @noindent
3094 but it uses less memory and avoids calling @code{conj} if @var{x} is real.
3095 @seealso{sum, prod}
3096 @end deftypefn */)
3097 {
3098  DATA_REDUCTION (sumsq);
3099 }
3100 
3101 /*
3102 %!assert (sumsq ([1, 2, 3]), 14)
3103 %!assert (sumsq ([-1; -2; 4i]), 21)
3104 %!assert (sumsq ([1, 2, 3; 2, 3, 4; 4i, 6i, 2]), [21, 49, 29])
3105 
3106 %!assert (sumsq (single ([1, 2, 3])), single (14))
3107 %!assert (sumsq (single ([-1; -2; 4i])), single (21))
3108 %!assert (sumsq (single ([1, 2, 3; 2, 3, 4; 4i, 6i, 2])), single ([21, 49, 29]))
3109 
3110 %!assert (sumsq ([1, 2; 3, 4], 1), [10, 20])
3111 %!assert (sumsq ([1, 2; 3, 4], 2), [5; 25])
3112 
3113 %!assert (sumsq (single ([1, 2; 3, 4]), 1), single ([10, 20]))
3114 %!assert (sumsq (single ([1, 2; 3, 4]), 2), single ([5; 25]))
3115 
3116 %!error sumsq ()
3117 */
3118 
3119 DEFUN (islogical, args, ,
3120  doc: /* -*- texinfo -*-
3121 @deftypefn {} {} islogical (@var{x})
3122 @deftypefnx {} {} isbool (@var{x})
3123 Return true if @var{x} is a logical object.
3124 @seealso{ischar, isfloat, isinteger, isstring, isnumeric, isa}
3125 @end deftypefn */)
3126 {
3127  if (args.length () != 1)
3128  print_usage ();
3129 
3130  return ovl (args(0).islogical ());
3131 }
3132 
3133 DEFALIAS (isbool, islogical);
3134 
3135 /*
3136 %!assert (islogical (true), true)
3137 %!assert (islogical (false), true)
3138 %!assert (islogical ([true, false]), true)
3139 %!assert (islogical (1), false)
3140 %!assert (islogical (1i), false)
3141 %!assert (islogical ([1,1]), false)
3142 %!assert (islogical (single (1)), false)
3143 %!assert (islogical (single (1i)), false)
3144 %!assert (islogical (single ([1,1])), false)
3145 %!assert (islogical (sparse ([true, false])), true)
3146 %!assert (islogical (sparse ([1, 0])), false)
3147 */
3148 
3149 DEFUN (isinteger, args, ,
3150  doc: /* -*- texinfo -*-
3151 @deftypefn {} {} isinteger (@var{x})
3152 Return true if @var{x} is an integer object (int8, uint8, int16, etc.).
3153 
3154 Note that @w{@code{isinteger (14)}} is false because numeric constants in
3155 Octave are double precision floating point values.
3156 @seealso{isfloat, ischar, islogical, isstring, isnumeric, isa}
3157 @end deftypefn */)
3158 {
3159  if (args.length () != 1)
3160  print_usage ();
3161 
3162  return ovl (args(0).isinteger ());
3163 }
3164 
3165 /*
3166 %!assert (isinteger (int8 (16)))
3167 %!assert (isinteger (int16 (16)))
3168 %!assert (isinteger (int32 (16)))
3169 %!assert (isinteger (int64 (16)))
3170 
3171 %!assert (isinteger (uint8 (16)))
3172 %!assert (isinteger (uint16 (16)))
3173 %!assert (isinteger (uint32 (16)))
3174 %!assert (isinteger (uint64 (16)))
3175 
3176 %!assert (isinteger (intmax ("int8")))
3177 %!assert (isinteger (intmax ("int16")))
3178 %!assert (isinteger (intmax ("int32")))
3179 %!assert (isinteger (intmax ("int64")))
3180 
3181 %!assert (isinteger (intmax ("uint8")))
3182 %!assert (isinteger (intmax ("uint16")))
3183 %!assert (isinteger (intmax ("uint32")))
3184 %!assert (isinteger (intmax ("uint64")))
3185 
3186 %!assert (isinteger (intmin ("int8")))
3187 %!assert (isinteger (intmin ("int16")))
3188 %!assert (isinteger (intmin ("int32")))
3189 %!assert (isinteger (intmin ("int64")))
3190 
3191 %!assert (isinteger (intmin ("uint8")))
3192 %!assert (isinteger (intmin ("uint16")))
3193 %!assert (isinteger (intmin ("uint32")))
3194 %!assert (isinteger (intmin ("uint64")))
3195 
3196 %!assert (isinteger (uint8 ([1:10])))
3197 %!assert (isinteger (uint8 ([1:10; 1:10])))
3198 
3199 %!assert (! isinteger (16))
3200 %!assert (! isinteger ("parrot"))
3201 %!assert (! isinteger ([1, 2, 3]))
3202 
3203 %!error isinteger ()
3204 %!error isinteger ("multiple", "parameters")
3205 */
3206 
3207 DEFUN (iscomplex, args, ,
3208  doc: /* -*- texinfo -*-
3209 @deftypefn {} {} iscomplex (@var{x})
3210 Return true if @var{x} is a complex-valued numeric object.
3211 @seealso{isreal, isnumeric, ischar, isfloat, islogical, isstring, isa}
3212 @end deftypefn */)
3213 {
3214  if (args.length () != 1)
3215  print_usage ();
3216 
3217  return ovl (args(0).iscomplex ());
3218 }
3219 
3220 DEFUN (isfloat, args, ,
3221  doc: /* -*- texinfo -*-
3222 @deftypefn {} {} isfloat (@var{x})
3223 Return true if @var{x} is a floating-point numeric object.
3224 
3225 Objects of class double or single are floating-point objects.
3226 @seealso{isinteger, ischar, islogical, isnumeric, isstring, isa}
3227 @end deftypefn */)
3228 {
3229  if (args.length () != 1)
3230  print_usage ();
3231 
3232  return ovl (args(0).isfloat ());
3233 }
3234 
3235 // FIXME: perhaps this should be implemented with an
3236 // octave_value member function?
3237 
3238 DEFUN (complex, args, ,
3239  doc: /* -*- texinfo -*-
3240 @deftypefn {} {} complex (@var{x})
3241 @deftypefnx {} {} complex (@var{re}, @var{im})
3242 Return a complex value from real arguments.
3243 
3244 With 1 real argument @var{x}, return the complex result
3245 @w{@code{@var{x} + 0i}}.
3246 
3247 With 2 real arguments, return the complex result
3248 @w{@code{@var{re} + @var{im}i}}.
3249 @code{complex} can often be more convenient than expressions such as
3250 @w{@code{a + b*i}}.
3251 For example:
3252 
3253 @example
3254 @group
3255 complex ([1, 2], [3, 4])
3256  @result{} [ 1 + 3i 2 + 4i ]
3257 @end group
3258 @end example
3259 @seealso{real, imag, iscomplex, abs, arg}
3260 @end deftypefn */)
3261 {
3262  int nargin = args.length ();
3263 
3264  if (nargin < 1 || nargin > 2)
3265  print_usage ();
3266 
3268 
3269  if (nargin == 1)
3270  {
3271  octave_value arg = args(0);
3272 
3273  if (arg.iscomplex ())
3274  retval = arg;
3275  else
3276  {
3277  if (arg.issparse ())
3278  {
3279  SparseComplexMatrix val = arg.xsparse_complex_matrix_value ("complex: invalid conversion");
3280 
3282  }
3283  else if (arg.is_single_type ())
3284  {
3285  if (arg.numel () == 1)
3286  {
3287  FloatComplex val = arg.xfloat_complex_value ("complex: invalid conversion");
3288 
3290  }
3291  else
3292  {
3293  FloatComplexNDArray val = arg.xfloat_complex_array_value ("complex: invalid conversion");
3294 
3296  }
3297  }
3298  else
3299  {
3300  if (arg.numel () == 1)
3301  {
3302  Complex val = arg.xcomplex_value ("complex: invalid conversion");
3303 
3305  }
3306  else
3307  {
3308  ComplexNDArray val = arg.xcomplex_array_value ("complex: invalid conversion");
3309 
3311  }
3312  }
3313  }
3314  }
3315  else
3316  {
3317  octave_value re = args(0);
3318  octave_value im = args(1);
3319 
3320  if (re.issparse () && im.issparse ())
3321  {
3322  const SparseMatrix re_val = re.sparse_matrix_value ();
3323  const SparseMatrix im_val = im.sparse_matrix_value ();
3324 
3325  if (re.numel () == 1)
3326  {
3328  if (re_val.nnz () == 0)
3329  result = Complex (0, 1) * SparseComplexMatrix (im_val);
3330  else
3331  {
3332  octave_idx_type nr = im_val.rows ();
3333  octave_idx_type nc = im_val.cols ();
3334  result = SparseComplexMatrix (nr, nc, re_val(0));
3335 
3336  for (octave_idx_type j = 0; j < nc; j++)
3337  {
3338  octave_idx_type off = j * nr;
3339  for (octave_idx_type i = im_val.cidx (j);
3340  i < im_val.cidx (j + 1); i++)
3341  result.data (im_val.ridx (i) + off) +=
3342  Complex (0, im_val.data (i));
3343  }
3344  }
3346  }
3347  else if (im.numel () == 1)
3348  {
3350  if (im_val.nnz () == 0)
3351  result = SparseComplexMatrix (re_val);
3352  else
3353  {
3354  octave_idx_type nr = re_val.rows ();
3355  octave_idx_type nc = re_val.cols ();
3356  result = SparseComplexMatrix (nr, nc,
3357  Complex (0, im_val(0)));
3358 
3359  for (octave_idx_type j = 0; j < nc; j++)
3360  {
3361  octave_idx_type off = j * nr;
3362  for (octave_idx_type i = re_val.cidx (j);
3363  i < re_val.cidx (j + 1); i++)
3364  result.data (re_val.ridx (i) + off) +=
3365  re_val.data (i);
3366  }
3367  }
3369  }
3370  else
3371  {
3372  if (re_val.dims () != im_val.dims ())
3373  error ("complex: dimension mismatch");
3374 
3376  result = SparseComplexMatrix (re_val)
3377  + Complex (0, 1) * SparseComplexMatrix (im_val);
3379  }
3380  }
3381  else if (re.is_single_type () || im.is_single_type ())
3382  {
3383  if (re.numel () == 1)
3384  {
3385  float re_val = re.float_value ();
3386 
3387  if (im.numel () == 1)
3388  {
3389  float im_val = im.double_value ();
3390 
3392  (FloatComplex (re_val, im_val)));
3393  }
3394  else
3395  {
3396  const FloatNDArray im_val = im.float_array_value ();
3397 
3398  FloatComplexNDArray result (im_val.dims (),
3399  FloatComplex ());
3400 
3401  for (octave_idx_type i = 0; i < im_val.numel (); i++)
3402  result.xelem (i) = FloatComplex (re_val, im_val(i));
3403 
3405  (result));
3406  }
3407  }
3408  else
3409  {
3410  const FloatNDArray re_val = re.float_array_value ();
3411 
3412  if (im.numel () == 1)
3413  {
3414  float im_val = im.float_value ();
3415 
3416  FloatComplexNDArray result (re_val.dims (),
3417  FloatComplex ());
3418 
3419  for (octave_idx_type i = 0; i < re_val.numel (); i++)
3420  result.xelem (i) = FloatComplex (re_val(i), im_val);
3421 
3423  (result));
3424  }
3425  else
3426  {
3427  const FloatNDArray im_val = im.float_array_value ();
3428 
3429  if (re_val.dims () != im_val.dims ())
3430  error ("complex: dimension mismatch");
3431 
3432  FloatComplexNDArray result (re_val.dims (),
3433  FloatComplex ());
3434 
3435  for (octave_idx_type i = 0; i < re_val.numel (); i++)
3436  result.xelem (i) = FloatComplex (re_val(i),
3437  im_val(i));
3438 
3440  (result));
3441  }
3442  }
3443  }
3444  else if (re.numel () == 1)
3445  {
3446  double re_val = re.double_value ();
3447 
3448  if (im.numel () == 1)
3449  {
3450  double im_val = im.double_value ();
3451 
3453  (Complex (re_val, im_val)));
3454  }
3455  else
3456  {
3457  const NDArray im_val = im.array_value ();
3458 
3459  ComplexNDArray result (im_val.dims (), Complex ());
3460 
3461  for (octave_idx_type i = 0; i < im_val.numel (); i++)
3462  result.xelem (i) = Complex (re_val, im_val(i));
3463 
3465  }
3466  }
3467  else
3468  {
3469  const NDArray re_val = re.array_value ();
3470 
3471  if (im.numel () == 1)
3472  {
3473  double im_val = im.double_value ();
3474 
3475  ComplexNDArray result (re_val.dims (), Complex ());
3476 
3477  for (octave_idx_type i = 0; i < re_val.numel (); i++)
3478  result.xelem (i) = Complex (re_val(i), im_val);
3479 
3481  }
3482  else
3483  {
3484  const NDArray im_val = im.array_value ();
3485 
3486  if (re_val.dims () != im_val.dims ())
3487  error ("complex: dimension mismatch");
3488 
3489  ComplexNDArray result (re_val.dims (), Complex ());
3490 
3491  for (octave_idx_type i = 0; i < re_val.numel (); i++)
3492  result.xelem (i) = Complex (re_val(i), im_val(i));
3493 
3495  }
3496  }
3497  }
3498 
3499  return retval;
3500 }
3501 
3502 /*
3503 %!error <undefined> 1+Infj
3504 %!error <undefined> 1+Infi
3505 
3506 %!test <31974>
3507 %! assert (Inf + Inf*i, complex (Inf, Inf))
3508 %!
3509 %! assert (1 + Inf*i, complex (1, Inf))
3510 %! assert (1 + Inf*j, complex (1, Inf))
3511 %!
3512 %! ## whitespace should not affect parsing
3513 %! assert (1+Inf*i, complex (1, Inf))
3514 %! assert (1+Inf*j, complex (1, Inf))
3515 %!
3516 %! assert (NaN*j, complex (0, NaN))
3517 %!
3518 %! assert (Inf * 4j, complex (0, Inf))
3519 
3520 %!test <31974>
3521 %! x = Inf;
3522 %! assert (x * j, complex (0, Inf))
3523 %! j = complex (0, 1);
3524 %! assert (Inf * j, complex (0, Inf))
3525 
3526 %!test <31974>
3527 %! exp = complex (zeros (2, 2), Inf (2, 2));
3528 %! assert (Inf (2, 2) * j, exp)
3529 %! assert (Inf (2, 2) .* j, exp)
3530 %! assert (Inf * (ones (2, 2) * j), exp)
3531 %! assert (Inf (2, 2) .* (ones (2, 2) * j), exp)
3532 
3533 %!test <31974>
3534 %! assert ([Inf; 0] * [i, 0], complex ([NaN NaN; 0 0], [Inf NaN; 0 0]))
3535 %! assert ([Inf, 0] * [i; 0], complex (NaN, Inf))
3536 %! assert ([Inf, 0] .* [i, 0], complex ([0 0], [Inf 0]))
3537 
3538 %!test <31974>
3539 %! m = @(x, y) x * y;
3540 %! d = @(x, y) x / y;
3541 %! assert (m (Inf, i), complex (0, +Inf))
3542 %! assert (d (Inf, i), complex (0, -Inf))
3543 */
3544 
3545 DEFUN (isreal, args, ,
3546  doc: /* -*- texinfo -*-
3547 @deftypefn {} {} isreal (@var{x})
3548 Return true if @var{x} is a non-complex matrix or scalar.
3549 
3550 For compatibility with @sc{matlab}, this includes logical and character
3551 matrices.
3552 @seealso{iscomplex, isnumeric, isa}
3553 @end deftypefn */)
3554 {
3555  if (args.length () != 1)
3556  print_usage ();
3557 
3558  return ovl (args(0).isreal ());
3559 }
3560 
3561 DEFUN (isempty, args, ,
3562  doc: /* -*- texinfo -*-
3563 @deftypefn {} {} isempty (@var{a})
3564 Return true if @var{a} is an empty matrix (any one of its dimensions is
3565 zero).
3566 @seealso{isnull, isa}
3567 @end deftypefn */)
3568 {
3569  if (args.length () != 1)
3570  print_usage ();
3571 
3572  return ovl (args(0).isempty ());
3573 }
3574 
3575 /*
3576 ## Debian bug #706376
3577 %!assert (isempty (speye(2^16)), false)
3578 */
3579 
3580 DEFUN (isnumeric, args, ,
3581  doc: /* -*- texinfo -*-
3582 @deftypefn {} {} isnumeric (@var{x})
3583 Return true if @var{x} is a numeric object, i.e., an integer, real, or
3584 complex array.
3585 
3586 Logical and character arrays are not considered to be numeric.
3587 @seealso{isinteger, isfloat, isreal, iscomplex, ischar, islogical, isstring, iscell, isstruct, isa}
3588 @end deftypefn */)
3589 {
3590  if (args.length () != 1)
3591  print_usage ();
3592 
3593  return ovl (args(0).isnumeric ());
3594 }
3595 
3596 /*
3597 %!assert (isnumeric (1), true)
3598 %!assert (isnumeric (1i), true)
3599 %!assert (isnumeric ([1,1]), true)
3600 %!assert (isnumeric (single (1)), true)
3601 %!assert (isnumeric (single (1i)), true)
3602 %!assert (isnumeric (single ([1,1])), true)
3603 %!assert (isnumeric (int8 (1)), true)
3604 %!assert (isnumeric (uint8 ([1,1])), true)
3605 %!assert (isnumeric ("Hello World"), false)
3606 %!assert (isnumeric (true), false)
3607 %!assert (isnumeric (false), false)
3608 %!assert (isnumeric ([true, false]), false)
3609 %!assert (isnumeric (sparse ([true, false])), false)
3610 */
3611 
3612 DEFUN (isscalar, args, ,
3613  doc: /* -*- texinfo -*-
3614 @deftypefn {} {} isscalar (@var{x})
3615 Return true if @var{x} is a scalar.
3616 @seealso{isvector, ismatrix}
3617 @end deftypefn */)
3618 {
3619  if (args.length () != 1)
3620  print_usage ();
3621 
3622  return ovl (args(0).numel () == 1);
3623 }
3624 
3625 /*
3626 %!assert (isscalar (1))
3627 %!assert (isscalar ([1, 2]), false)
3628 %!assert (isscalar ([]), false)
3629 %!assert (isscalar ([1, 2; 3, 4]), false)
3630 
3631 %!assert (isscalar ("t"))
3632 %!assert (isscalar ("test"), false)
3633 %!assert (isscalar (["test"; "ing"]), false)
3634 
3635 %!test
3636 %! s.a = 1;
3637 %! assert (isscalar (s));
3638 
3639 ## Test input validation
3640 %!error isscalar ()
3641 %!error isscalar (1, 2)
3642 */
3643 
3644 DEFUN (isvector, args, ,
3645  doc: /* -*- texinfo -*-
3646 @deftypefn {} {} isvector (@var{x})
3647 Return true if @var{x} is a vector.
3648 
3649 A vector is a 2-D array where one of the dimensions is equal to 1. As a
3650 consequence a 1x1 array, or scalar, is also a vector.
3651 @seealso{isscalar, ismatrix, size, rows, columns, length}
3652 @end deftypefn */)
3653 {
3654  if (args.length () != 1)
3655  print_usage ();
3656 
3657  dim_vector sz = args(0).dims ();
3658 
3659  return ovl (sz.ndims () == 2 && (sz(0) == 1 || sz(1) == 1));
3660 }
3661 
3662 /*
3663 %!assert (isvector (1), true)
3664 %!assert (isvector ([1; 2; 3]), true)
3665 %!assert (isvector ([1, 2, 3]), true)
3666 %!assert (isvector ([]), false)
3667 %!assert (isvector ([1, 2; 3, 4]), false)
3668 
3669 %!assert (isvector ("t"), true)
3670 %!assert (isvector ("test"), true)
3671 %!assert (isvector (["test"; "ing"]), false)
3672 
3673 %!test
3674 %! s.a = 1;
3675 %! assert (isvector (s), true);
3676 
3677 ## Test input validation
3678 %!error isvector ()
3679 %!error isvector ([1, 2], 2)
3680 */
3681 
3682 DEFUN (isrow, args, ,
3683  doc: /* -*- texinfo -*-
3684 @deftypefn {} {} isrow (@var{x})
3685 Return true if @var{x} is a row vector 1xN with non-negative N.
3686 @seealso{iscolumn, isscalar, isvector, ismatrix}
3687 @end deftypefn */)
3688 {
3689  if (args.length () != 1)
3690  print_usage ();
3691 
3692  dim_vector sz = args(0).dims ();
3693 
3694  return ovl (sz.ndims () == 2 && sz(0) == 1);
3695 }
3696 
3697 /*
3698 %!assert (isrow ([1, 2, 3]))
3699 %!assert (isrow ([1; 2; 3]), false)
3700 %!assert (isrow (1))
3701 %!assert (isrow ([]), false)
3702 %!assert (isrow ([1, 2; 3, 4]), false)
3703 
3704 %!assert (isrow (ones (1, 0)), true)
3705 %!assert (isrow (ones (1, 1)), true)
3706 %!assert (isrow (ones (1, 2)), true)
3707 %!assert (isrow (ones (1, 1, 1)), true)
3708 %!assert (isrow (ones (1, 1, 1, 1)), true)
3709 
3710 %!assert (isrow (ones (0, 0)), false)
3711 %!assert (isrow (ones (1, 1, 0)), false)
3712 
3713 %!assert (isrow ("t"), true)
3714 %!assert (isrow ("test"), true)
3715 %!assert (isrow (["test"; "ing"]), false)
3716 
3717 %!test
3718 %! s.a = 1;
3719 %! assert (isrow (s), true);
3720 
3721 ## Test input validation
3722 %!error isrow ()
3723 %!error isrow ([1, 2], 2)
3724 */
3725 
3726 DEFUN (iscolumn, args, ,
3727  doc: /* -*- texinfo -*-
3728 @deftypefn {} {} iscolumn (@var{x})
3729 Return true if @var{x} is a column vector Nx1 with non-negative N.
3730 @seealso{isrow, isscalar, isvector, ismatrix}
3731 @end deftypefn */)
3732 {
3733  if (args.length () != 1)
3734  print_usage ();
3735 
3736  dim_vector sz = args(0).dims ();
3737 
3738  return ovl (sz.ndims () == 2 && sz(1) == 1);
3739 }
3740 
3741 /*
3742 %!assert (iscolumn ([1, 2, 3]), false)
3743 %!assert (iscolumn ([1; 2; 3]), true)
3744 %!assert (iscolumn (1), true)
3745 %!assert (iscolumn ([]), false)
3746 %!assert (iscolumn ([1, 2; 3, 4]), false)
3747 
3748 %!assert (iscolumn ("t"), true)
3749 %!assert (iscolumn ("test"), false)
3750 %!assert (iscolumn (["test"; "ing"]), false)
3751 
3752 %!assert (iscolumn (ones (0, 1)), true)
3753 %!assert (iscolumn (ones (1, 1)), true)
3754 %!assert (iscolumn (ones (2, 1)), true)
3755 %!assert (iscolumn (ones (1, 1, 1)), true)
3756 %!assert (iscolumn (ones (1, 1, 1, 1)), true)
3757 
3758 %!assert (iscolumn (ones (0, 0)), false)
3759 %!assert (iscolumn (ones (0, 1, 0)), false)
3760 
3761 %!test
3762 %! s.a = 1;
3763 %! assert (iscolumn (s));
3764 
3765 ## Test input validation
3766 %!error iscolumn ()
3767 %!error iscolumn ([1, 2], 2)
3768 */
3769 
3770 DEFUN (ismatrix, args, ,
3771  doc: /* -*- texinfo -*-
3772 @deftypefn {} {} ismatrix (@var{a})
3773 Return true if @var{a} is a 2-D array.
3774 @seealso{isscalar, isvector, iscell, isstruct, issparse, isa}
3775 @end deftypefn */)
3776 {
3777  if (args.length () != 1)
3778  print_usage ();
3779 
3780  dim_vector sz = args(0).dims ();
3781 
3782  return ovl (sz.ndims () == 2 && sz(0) >= 0 && sz(1) >= 0);
3783 }
3784 
3785 /*
3786 %!assert (ismatrix ([]), true)
3787 %!assert (ismatrix (1), true)
3788 %!assert (ismatrix ([1, 2, 3]), true)
3789 %!assert (ismatrix ([1, 2; 3, 4]), true)
3790 
3791 %!assert (ismatrix (zeros (0)), true)
3792 %!assert (ismatrix (zeros (0, 0)), true)
3793 %!assert (ismatrix (zeros (0, 0, 0)), false)
3794 %!assert (ismatrix (zeros (3, 2, 4)), false)
3795 
3796 %!assert (ismatrix (single ([])), true)
3797 %!assert (ismatrix (single (1)), true)
3798 %!assert (ismatrix (single ([1, 2, 3])), true)
3799 %!assert (ismatrix (single ([1, 2; 3, 4])), true)
3800 
3801 %!assert (ismatrix ("t"), true)
3802 %!assert (ismatrix ("test"), true)
3803 %!assert (ismatrix (["test"; "ing"]), true)
3804 
3805 %!test
3806 %! s.a = 1;
3807 %! assert (ismatrix (s), true);
3808 
3809 %!error ismatrix ()
3810 %!error ismatrix ([1, 2; 3, 4], 2)
3811 */
3812 
3813 DEFUN (issquare, args, ,
3814  doc: /* -*- texinfo -*-
3815 @deftypefn {} {} issquare (@var{x})
3816 Return true if @var{x} is a square matrix.
3817 @seealso{isscalar, isvector, ismatrix, size}
3818 @end deftypefn */)
3819 {
3820  if (args.length () != 1)
3821  print_usage ();
3822 
3823  dim_vector sz = args(0).dims ();
3824 
3825  return ovl (sz.ndims () == 2 && sz(0) == sz(1));
3826 }
3827 
3828 /*
3829 %!assert (issquare ([]))
3830 %!assert (issquare (1))
3831 %!assert (! issquare ([1, 2]))
3832 %!assert (issquare ([1, 2; 3, 4]))
3833 %!assert (! issquare ([1, 2; 3, 4; 5, 6]))
3834 %!assert (! issquare (ones (3,3,3)))
3835 %!assert (issquare ("t"))
3836 %!assert (! issquare ("test"))
3837 %!assert (issquare (["test"; "ing"; "1"; "2"]))
3838 %!test
3839 %! s.a = 1;
3840 %! assert (issquare (s));
3841 %!assert (issquare ({1, 2; 3, 4}))
3842 %!assert (sparse (([1, 2; 3, 4])))
3843 
3844 ## Test input validation
3845 %!error issquare ()
3846 %!error issquare ([1, 2; 3, 4], 2)
3847 */
3848 
3849 static octave_value
3850 fill_matrix (const octave_value_list& args, int val, const char *fcn)
3851 {
3853 
3854  int nargin = args.length ();
3855 
3857 
3858  dim_vector dims (1, 1);
3859 
3860  if (nargin > 0 && args(nargin-1).is_string ())
3861  {
3862  std::string nm = args(nargin-1).string_value ();
3863  nargin--;
3864 
3866  }
3867 
3868  switch (nargin)
3869  {
3870  case 0:
3871  break;
3872 
3873  case 1:
3874  get_dimensions (args(0), fcn, dims);
3875  break;
3876 
3877  default:
3878  {
3879  dims.resize (nargin);
3880 
3881  for (int i = 0; i < nargin; i++)
3882  dims(i) = (args(i).isempty ()
3883  ? 0 : args(i).xidx_type_value ("%s: dimension arguments must be scalar integers", fcn));
3884  }
3885  break;
3886  }
3887 
3888  dims.chop_trailing_singletons ();
3889 
3891 
3892  // FIXME: perhaps this should be made extensible by
3893  // using the class name to lookup a function to call to create
3894  // the new value.
3895 
3896  // Note that automatic narrowing will handle conversion from
3897  // NDArray to scalar.
3898 
3899  switch (dt)
3900  {
3902  retval = int8NDArray (dims, val);
3903  break;
3904 
3906  retval = uint8NDArray (dims, val);
3907  break;
3908 
3910  retval = int16NDArray (dims, val);
3911  break;
3912 
3915  break;
3916 
3918  retval = int32NDArray (dims, val);
3919  break;
3920 
3923  break;
3924 
3926  retval = int64NDArray (dims, val);
3927  break;
3928 
3931  break;
3932 
3934  retval = FloatNDArray (dims, val);
3935  break;
3936 
3938  {
3939  if (val == 1 && dims.ndims () == 2 && dims(0) == 1)
3940  retval = Range (1.0, 0.0, dims(1)); // packed form
3941  else
3942  retval = NDArray (dims, val);
3943  }
3944  break;
3945 
3947  retval = boolNDArray (dims, val);
3948  break;
3949 
3950  default:
3951  error ("%s: invalid class name", fcn);
3952  break;
3953  }
3954 
3955  return retval;
3956 }
3957 
3958 static octave_value
3959 fill_matrix (const octave_value_list& args, double val, float fval,
3960  const char *fcn)
3961 {
3963 
3964  int nargin = args.length ();
3965 
3967 
3968  dim_vector dims (1, 1);
3969 
3970  if (nargin > 0 && args(nargin-1).is_string ())
3971  {
3972  std::string nm = args(nargin-1).string_value ();
3973  nargin--;
3974 
3976  }
3977 
3978  switch (nargin)
3979  {
3980  case 0:
3981  break;
3982 
3983  case 1:
3984  get_dimensions (args(0), fcn, dims);
3985  break;
3986 
3987  default:
3988  {
3989  dims.resize (nargin);
3990 
3991  for (int i = 0; i < nargin; i++)
3992  dims(i) = (args(i).isempty ()
3993  ? 0 : args(i).xidx_type_value ("%s: dimension arguments must be scalar integers", fcn));
3994  }
3995  break;
3996  }
3997 
3998  dims.chop_trailing_singletons ();
3999 
4001 
4002  // Note that automatic narrowing will handle conversion from
4003  // NDArray to scalar.
4004 
4005  switch (dt)
4006  {
4008  retval = FloatNDArray (dims, fval);
4009  break;
4010 
4012  retval = NDArray (dims, val);
4013  break;
4014 
4015  default:
4016  error ("%s: invalid class name", fcn);
4017  break;
4018  }
4019 
4020  return retval;
4021 }
4022 
4023 static octave_value
4024 fill_matrix (const octave_value_list& args, double val, const char *fcn)
4025 {
4027 
4028  int nargin = args.length ();
4029 
4031 
4032  dim_vector dims (1, 1);
4033 
4034  if (nargin > 0 && args(nargin-1).is_string ())
4035  {
4036  std::string nm = args(nargin-1).string_value ();
4037  nargin--;
4038 
4040  }
4041 
4042  switch (nargin)
4043  {
4044  case 0:
4045  break;
4046 
4047  case 1:
4048  get_dimensions (args(0), fcn, dims);
4049  break;
4050 
4051  default:
4052  {
4053  dims.resize (nargin);
4054 
4055  for (int i = 0; i < nargin; i++)
4056  dims(i) = (args(i).isempty ()
4057  ? 0 : args(i).xidx_type_value ("%s: dimension arguments must be scalar integers", fcn));
4058  }
4059  break;
4060  }
4061 
4062  dims.chop_trailing_singletons ();
4063 
4065 
4066  // Note that automatic narrowing will handle conversion from
4067  // NDArray to scalar.
4068 
4069  switch (dt)
4070  {
4072  retval = FloatNDArray (dims, static_cast<float> (val));
4073  break;
4074 
4076  retval = NDArray (dims, val);
4077  break;
4078 
4079  default:
4080  error ("%s: invalid class name", fcn);
4081  break;
4082  }
4083 
4084  return retval;
4085 }
4086 
4087 static octave_value
4088 fill_matrix (const octave_value_list& args, const Complex& val,
4089  const char *fcn)
4090 {
4092 
4093  int nargin = args.length ();
4094 
4096 
4097  dim_vector dims (1, 1);
4098 
4099  if (nargin > 0 && args(nargin-1).is_string ())
4100  {
4101  std::string nm = args(nargin-1).string_value ();
4102  nargin--;
4103 
4105  }
4106 
4107  switch (nargin)
4108  {
4109  case 0:
4110  break;
4111 
4112  case 1:
4113  get_dimensions (args(0), fcn, dims);
4114  break;
4115 
4116  default:
4117  {
4118  dims.resize (nargin);
4119 
4120  for (int i = 0; i < nargin; i++)
4121  dims(i) = (args(i).isempty ()
4122  ? 0 : args(i).xidx_type_value ("%s: dimension arguments must be scalar integers", fcn));
4123  }
4124  break;
4125  }
4126 
4127  dims.chop_trailing_singletons ();
4128 
4130 
4131  // Note that automatic narrowing will handle conversion from
4132  // NDArray to scalar.
4133 
4134  switch (dt)
4135  {
4138  static_cast<FloatComplex> (val));
4139  break;
4140 
4143  break;
4144 
4145  default:
4146  error ("%s: invalid class name", fcn);
4147  break;
4148  }
4149 
4150  return retval;
4151 }
4152 
4153 static octave_value
4154 fill_matrix (const octave_value_list& args, bool val, const char *fcn)
4155 {
4157 
4158  int nargin = args.length ();
4159 
4160  dim_vector dims (1, 1);
4161 
4162  switch (nargin)
4163  {
4164  case 0:
4165  break;
4166 
4167  case 1:
4168  get_dimensions (args(0), fcn, dims);
4169  break;
4170 
4171  default:
4172  {
4173  dims.resize (nargin);
4174 
4175  for (int i = 0; i < nargin; i++)
4176  dims(i) = (args(i).isempty ()
4177  ? 0 : args(i).xidx_type_value ("%s: dimension arguments must be scalar integers", fcn));
4178  }
4179  break;
4180  }
4181 
4182  dims.chop_trailing_singletons ();
4183 
4185 
4186  // Note that automatic narrowing will handle conversion from
4187  // NDArray to scalar.
4188 
4189  retval = boolNDArray (dims, val);
4190 
4191  return retval;
4192 }
4193 
4194 DEFUN (ones, args, ,
4195  doc: /* -*- texinfo -*-
4196 @deftypefn {} {} ones (@var{n})
4197 @deftypefnx {} {} ones (@var{m}, @var{n})
4198 @deftypefnx {} {} ones (@var{m}, @var{n}, @var{k}, @dots{})
4199 @deftypefnx {} {} ones ([@var{m} @var{n} @dots{}])
4200 @deftypefnx {} {} ones (@dots{}, @var{class})
4201 Return a matrix or N-dimensional array whose elements are all 1.
4202 
4203 If invoked with a single scalar integer argument @var{n}, return a square
4204 @nospell{NxN} matrix.
4205 
4206 If invoked with two or more scalar integer arguments, or a vector of integer
4207 values, return an array with the given dimensions.
4208 
4209 To create a constant matrix whose values are all the same use an expression
4210 such as
4211 
4212 @example
4213 val_matrix = val * ones (m, n)
4214 @end example
4215 
4216 The optional argument @var{class} specifies the class of the return array
4217 and defaults to double. For example:
4218 
4219 @example
4220 val = ones (m,n, "uint8")
4221 @end example
4222 @seealso{zeros}
4223 @end deftypefn */)
4224 {
4225  return fill_matrix (args, 1, "ones");
4226 }
4227 
4228 /*
4229 %!assert (ones (3), [1, 1, 1; 1, 1, 1; 1, 1, 1])
4230 %!assert (ones (2, 3), [1, 1, 1; 1, 1, 1])
4231 %!assert (ones (3, 2), [1, 1; 1, 1; 1, 1])
4232 %!assert (size (ones (3, 4, 5)), [3, 4, 5])
4233 
4234 %!assert (ones (3, "single"), single ([1, 1, 1; 1, 1, 1; 1, 1, 1]))
4235 %!assert (ones (2, 3, "single"), single ([1, 1, 1; 1, 1, 1]))
4236 %!assert (ones (3, 2, "single"), single ([1, 1; 1, 1; 1, 1]))
4237 %!assert (size (ones (3, 4, 5, "single")), [3, 4, 5])
4238 
4239 %!assert (ones (3, "int8"), int8 ([1, 1, 1; 1, 1, 1; 1, 1, 1]))
4240 %!assert (ones (2, 3, "int8"), int8 ([1, 1, 1; 1, 1, 1]))
4241 %!assert (ones (3, 2, "int8"), int8 ([1, 1; 1, 1; 1, 1]))
4242 %!assert (size (ones (3, 4, 5, "int8")), [3, 4, 5])
4243 */
4244 
4245 /*
4246 ## Tests for bug #47298
4247 ## Matlab requires the size to be a row vector. In that logic, it supports
4248 ## n to be a 1x0 vector (returns 0x0) but not a 0x1 vector. Octave supports
4249 ## any vector and therefore must support 0x1, 1x0, and 0x0x1 (but not 0x1x1).
4250 %!test <*47298>
4251 %! funcs = {@zeros, @ones, @inf, @nan, @NA, @i, @pi, @e};
4252 %! for idx = 1:numel (funcs)
4253 %! func = funcs{idx};
4254 %! assert (func (zeros (1, 0)), zeros (0, 0));
4255 %! assert (func (zeros (0, 1)), zeros (0, 0));
4256 %! assert (func (zeros (0, 1, 1)), zeros (0, 0));
4257 %! fail ([func2str(func) " ([])"]);
4258 %! fail ([func2str(func) " (zeros (0, 0, 1))"]);
4259 %! endfor
4260 */
4261 
4262 DEFUN (zeros, args, ,
4263  doc: /* -*- texinfo -*-
4264 @deftypefn {} {} zeros (@var{n})
4265 @deftypefnx {} {} zeros (@var{m}, @var{n})
4266 @deftypefnx {} {} zeros (@var{m}, @var{n}, @var{k}, @dots{})
4267 @deftypefnx {} {} zeros ([@var{m} @var{n} @dots{}])
4268 @deftypefnx {} {} zeros (@dots{}, @var{class})
4269 Return a matrix or N-dimensional array whose elements are all 0.
4270 
4271 If invoked with a single scalar integer argument, return a square
4272 @nospell{NxN} matrix.
4273 
4274 If invoked with two or more scalar integer arguments, or a vector of integer
4275 values, return an array with the given dimensions.
4276 
4277 The optional argument @var{class} specifies the class of the return array
4278 and defaults to double. For example:
4279 
4280 @example
4281 val = zeros (m,n, "uint8")
4282 @end example
4283 @seealso{ones}
4284 @end deftypefn */)
4285 {
4286  return fill_matrix (args, 0, "zeros");
4287 }
4288 
4289 /*
4290 %!assert (zeros (3), [0, 0, 0; 0, 0, 0; 0, 0, 0])
4291 %!assert (zeros (2, 3), [0, 0, 0; 0, 0, 0])
4292 %!assert (zeros (3, 2), [0, 0; 0, 0; 0, 0])
4293 %!assert (size (zeros (3, 4, 5)), [3, 4, 5])
4294 
4295 %!assert (zeros (3, "single"), single ([0, 0, 0; 0, 0, 0; 0, 0, 0]))
4296 %!assert (zeros (2, 3, "single"), single ([0, 0, 0; 0, 0, 0]))
4297 %!assert (zeros (3, 2, "single"), single ([0, 0; 0, 0; 0, 0]))
4298 %!assert (size (zeros (3, 4, 5, "single")), [3, 4, 5])
4299 
4300 %!assert (zeros (3, "int8"), int8 ([0, 0, 0; 0, 0, 0; 0, 0, 0]))
4301 %!assert (zeros (2, 3, "int8"), int8 ([0, 0, 0; 0, 0, 0]))
4302 %!assert (zeros (3, 2, "int8"), int8 ([0, 0; 0, 0; 0, 0]))
4303 %!assert (size (zeros (3, 4, 5, "int8")), [3, 4, 5])
4304 */
4305 
4306 DEFUN (Inf, args, ,
4307  doc: /* -*- texinfo -*-
4308 @c List other form of function in documentation index
4309 @findex inf
4310 
4311 @deftypefn {} {} Inf
4312 @deftypefnx {} {} Inf (@var{n})
4313 @deftypefnx {} {} Inf (@var{n}, @var{m})
4314 @deftypefnx {} {} Inf (@var{n}, @var{m}, @var{k}, @dots{})
4315 @deftypefnx {} {} Inf (@dots{}, @var{class})
4316 Return a scalar, matrix or N-dimensional array whose elements are all equal
4317 to the IEEE representation for positive infinity.
4318 
4319 Infinity is produced when results are too large to be represented using the
4320 IEEE floating point format for numbers. Two common examples which produce
4321 infinity are division by zero and overflow.
4322 
4323 @example
4324 @group
4325 [ 1/0 e^800 ]
4326 @result{} Inf Inf
4327 @end group
4328 @end example
4329 
4330 When called with no arguments, return a scalar with the value @samp{Inf}.
4331 
4332 When called with a single argument, return a square matrix with the
4333 dimension specified.
4334 
4335 When called with more than one scalar argument the first two arguments are
4336 taken as the number of rows and columns and any further arguments specify
4337 additional matrix dimensions.
4338 
4339 The optional argument @var{class} specifies the return type and may be
4340 either @qcode{"double"} or @qcode{"single"}.
4341 @seealso{isinf, NaN}
4342 @end deftypefn */)
4343 {
4344  return fill_matrix (args, lo_ieee_inf_value (),
4345  lo_ieee_float_inf_value (), "Inf");
4346 }
4347 
4348 DEFALIAS (inf, Inf);
4349 
4350 /*
4351 %!assert (inf (3), [Inf, Inf, Inf; Inf, Inf, Inf; Inf, Inf, Inf])
4352 %!assert (inf (2, 3), [Inf, Inf, Inf; Inf, Inf, Inf])
4353 %!assert (inf (3, 2), [Inf, Inf; Inf, Inf; Inf, Inf])
4354 %!assert (size (inf (3, 4, 5)), [3, 4, 5])
4355 
4356 %!assert (inf (3, "single"), single ([Inf, Inf, Inf; Inf, Inf, Inf; Inf, Inf, Inf]))
4357 %!assert (inf (2, 3, "single"), single ([Inf, Inf, Inf; Inf, Inf, Inf]))
4358 %!assert (inf (3, 2, "single"), single ([Inf, Inf; Inf, Inf; Inf, Inf]))
4359 %!assert (size (inf (3, 4, 5, "single")), [3, 4, 5])
4360 
4361 %!error (inf (3, "int8"))
4362 %!error (inf (2, 3, "int8"))
4363 %!error (inf (3, 2, "int8"))
4364 %!error (inf (3, 4, 5, "int8"))
4365 */
4366 
4367 DEFUN (NaN, args, ,
4368  doc: /* -*- texinfo -*-
4369 @c List other form of function in documentation index
4370 @findex nan
4371 
4372 @deftypefn {} {} NaN
4373 @deftypefnx {} {} NaN (@var{n})
4374 @deftypefnx {} {} NaN (@var{n}, @var{m})
4375 @deftypefnx {} {} NaN (@var{n}, @var{m}, @var{k}, @dots{})
4376 @deftypefnx {} {} NaN (@dots{}, @var{class})
4377 Return a scalar, matrix, or N-dimensional array whose elements are all equal
4378 to the IEEE symbol NaN (Not a Number).
4379 
4380 NaN is the result of operations which do not produce a well defined
4381 numerical result. Common operations which produce a NaN are arithmetic
4382 with infinity
4383 @tex
4384 ($\infty - \infty$), zero divided by zero ($0/0$),
4385 @end tex
4386 @ifnottex
4387 (Inf - Inf), zero divided by zero (0/0),
4388 @end ifnottex
4389 and any operation involving another NaN value (5 + NaN).
4390 
4391 Note that NaN always compares not equal to NaN (NaN != NaN). This behavior
4392 is specified by the IEEE standard for floating point arithmetic. To find
4393 NaN values, use the @code{isnan} function.
4394 
4395 When called with no arguments, return a scalar with the value @samp{NaN}.
4396 
4397 When called with a single argument, return a square matrix with the
4398 dimension specified.
4399 
4400 When called with more than one scalar argument the first two arguments are
4401 taken as the number of rows and columns and any further arguments specify
4402 additional matrix dimensions.
4403 
4404 The optional argument @var{class} specifies the return type and may be
4405 either @qcode{"double"} or @qcode{"single"}.
4406 @seealso{isnan, Inf}
4407 @end deftypefn */)
4408 {
4409  return fill_matrix (args, lo_ieee_nan_value (),
4410  lo_ieee_float_nan_value (), "NaN");
4411 }
4412 
4413 DEFALIAS (nan, NaN);
4414 
4415 /*
4416 %!assert (NaN (3), [NaN, NaN, NaN; NaN, NaN, NaN; NaN, NaN, NaN])
4417 %!assert (NaN (2, 3), [NaN, NaN, NaN; NaN, NaN, NaN])
4418 %!assert (NaN (3, 2), [NaN, NaN; NaN, NaN; NaN, NaN])
4419 %!assert (size (NaN (3, 4, 5)), [3, 4, 5])
4420 
4421 %!assert (NaN (3, "single"), single ([NaN, NaN, NaN; NaN, NaN, NaN; NaN, NaN, NaN]))
4422 %!assert (NaN (2, 3, "single"), single ([NaN, NaN, NaN; NaN, NaN, NaN]))
4423 %!assert (NaN (3, 2, "single"), single ([NaN, NaN; NaN, NaN; NaN, NaN]))
4424 %!assert (size (NaN (3, 4, 5, "single")), [3, 4, 5])
4425 
4426 %!error (NaN (3, "int8"))
4427 %!error (NaN (2, 3, "int8"))
4428 %!error (NaN (3, 2, "int8"))
4429 %!error (NaN (3, 4, 5, "int8"))
4430 */
4431 
4432 DEFUN (e, args, ,
4433  doc: /* -*- texinfo -*-
4434 @deftypefn {} {} e
4435 @deftypefnx {} {} e (@var{n})
4436 @deftypefnx {} {} e (@var{n}, @var{m})
4437 @deftypefnx {} {} e (@var{n}, @var{m}, @var{k}, @dots{})
4438 @deftypefnx {} {} e (@dots{}, @var{class})
4439 Return a scalar, matrix, or N-dimensional array whose elements are all equal
4440 to the base of natural logarithms.
4441 
4442 The constant
4443 @tex
4444 $e$ satisfies the equation $\log (e) = 1$.
4445 @end tex
4446 @ifnottex
4447 @samp{e} satisfies the equation @code{log} (e) = 1.
4448 @end ifnottex
4449 
4450 When called with no arguments, return a scalar with the value @math{e}.
4451 
4452 When called with a single argument, return a square matrix with the
4453 dimension specified.
4454 
4455 When called with more than one scalar argument the first two arguments are
4456 taken as the number of rows and columns and any further arguments specify
4457 additional matrix dimensions.
4458 
4459 The optional argument @var{class} specifies the return type and may be
4460 either @qcode{"double"} or @qcode{"single"}.
4461 @seealso{log, exp, pi, I}
4462 @end deftypefn */)
4463 {
4464 #if defined (M_E)
4465  double e_val = M_E;
4466 #else
4467  double e_val = exp (1.0);
4468 #endif
4469 
4470  return fill_matrix (args, e_val, "e");
4471 }
4472 
4473 template <typename T>
4474 T
4475 eps (const T& x)
4476 {
4477  T epsval = x.abs ();
4478  typedef typename T::value_type P;
4479  for (octave_idx_type i = 0; i < x.numel (); i++)
4480  {
4481  P val = epsval.xelem (i);
4483  epsval(i) = octave::numeric_limits<P>::NaN ();
4484  else if (val < std::numeric_limits<P>::min ())
4485  epsval(i) = std::numeric_limits<P>::denorm_min ();
4486  else
4487  {
4488  int exponent;
4489  octave::math::frexp (val, &exponent);
4490  const P digits = std::numeric_limits<P>::digits;
4491  epsval(i) = std::pow (static_cast<P> (2.0),
4492  static_cast<P> (exponent - digits));
4493  }
4494  }
4495  return epsval;
4496 }
4497 
4498 DEFUN (eps, args, ,
4499  doc: /* -*- texinfo -*-
4500 @deftypefn {} {} eps
4501 @deftypefnx {} {} eps (@var{x})
4502 @deftypefnx {} {} eps (@var{n}, @var{m})
4503 @deftypefnx {} {} eps (@var{n}, @var{m}, @var{k}, @dots{})
4504 @deftypefnx {} {} eps (@dots{}, @var{class})
4505 Return a scalar, matrix or N-dimensional array whose elements are all eps,
4506 the machine precision.
4507 
4508 More precisely, @code{eps} is the relative spacing between any two adjacent
4509 numbers in the machine's floating point system. This number is obviously
4510 system dependent. On machines that support IEEE floating point arithmetic,
4511 @code{eps} is approximately
4512 @tex
4513 $2.2204\times10^{-16}$ for double precision and $1.1921\times10^{-7}$
4514 @end tex
4515 @ifnottex
4516 2.2204e-16 for double precision and 1.1921e-07
4517 @end ifnottex
4518 for single precision.
4519 
4520 When called with no arguments, return a scalar with the value
4521 @code{eps (1.0)}.
4522 
4523 Given a single argument @var{x}, return the distance between @var{x} and the
4524 next largest value.
4525 
4526 When called with more than one argument the first two arguments are taken as
4527 the number of rows and columns and any further arguments specify additional
4528 matrix dimensions. The optional argument @var{class} specifies the return
4529 type and may be either @qcode{"double"} or @qcode{"single"}.
4530 @seealso{realmax, realmin, intmax, flintmax}
4531 @end deftypefn */)
4532 {
4534 
4535  if (args.length () == 1 && ! args(0).is_string ())
4536  {
4537  octave_value arg0 = args(0);
4538  if (arg0.is_single_type ())
4539  {
4540  FloatNDArray epsval = eps (arg0.float_array_value ());
4541  retval = epsval;
4542  }
4543  else if (arg0.is_double_type ())
4544  {
4545  NDArray epsval = eps (arg0.array_value ());
4546  retval = epsval;
4547  }
4548  else
4549  error ("eps: X must be of a floating point type");
4550  }
4551  else
4552  retval = fill_matrix (args, std::numeric_limits<double>::epsilon (),
4553  std::numeric_limits<float>::epsilon (), "eps");
4554 
4555  return retval;
4556 }
4557 
4558 /*
4559 %!assert (eps (1/2), 2^(-53))
4560 %!assert (eps (1), 2^(-52))
4561 %!assert (eps (2), 2^(-51))
4562 %!assert (eps (realmax), 2^971)
4563 %!assert (eps (0), 2^(-1074))
4564 %!assert (eps (realmin/2), 2^(-1074))
4565 %!assert (eps (realmin/16), 2^(-1074))
4566 %!assert (eps (Inf), NaN)
4567 %!assert (eps (NaN), NaN)
4568 %!assert (eps ([1/2 1 2 realmax 0 realmin/2 realmin/16 Inf NaN]),
4569 %! [2^(-53) 2^(-52) 2^(-51) 2^971 2^(-1074) 2^(-1074) 2^(-1074) NaN NaN])
4570 %!assert (eps (single (1/2)), single (2^(-24)))
4571 %!assert (eps (single (1)), single (2^(-23)))
4572 %!assert (eps (single (2)), single (2^(-22)))
4573 %!assert (eps (realmax ("single")), single (2^104))
4574 %!assert (eps (single (0)), single (2^(-149)))
4575 %!assert (eps (realmin ("single")/2), single (2^(-149)))
4576 %!assert (eps (realmin ("single")/16), single (2^(-149)))
4577 %!assert (eps (single (Inf)), single (NaN))
4578 %!assert (eps (single (NaN)), single (NaN))
4579 %!assert (eps (single ([1/2 1 2 realmax("single") 0 realmin("single")/2 realmin("single")/16 Inf NaN])),
4580 %! single ([2^(-24) 2^(-23) 2^(-22) 2^104 2^(-149) 2^(-149) 2^(-149) NaN NaN]))
4581 %!error <X must be of a floating point type> eps (uint8 ([0 1 2]))
4582 */
4583 
4584 DEFUN (pi, args, ,
4585  doc: /* -*- texinfo -*-
4586 @deftypefn {} {} pi
4587 @deftypefnx {} {} pi (@var{n})
4588 @deftypefnx {} {} pi (@var{n}, @var{m})
4589 @deftypefnx {} {} pi (@var{n}, @var{m}, @var{k}, @dots{})
4590 @deftypefnx {} {} pi (@dots{}, @var{class})
4591 Return a scalar, matrix, or N-dimensional array whose elements are all equal
4592 to the ratio of the circumference of a circle to its
4593 @tex
4594 diameter($\pi$).
4595 @end tex
4596 @ifnottex
4597 diameter.
4598 @end ifnottex
4599 
4600 Internally, @code{pi} is computed as @samp{4.0 * atan (1.0)}.
4601 
4602 When called with no arguments, return a scalar with the value of
4603 @tex
4604 $\pi$.
4605 @end tex
4606 @ifnottex
4607 pi.
4608 @end ifnottex
4609 
4610 When called with a single argument, return a square matrix with the
4611 dimension specified.
4612 
4613 When called with more than one scalar argument the first two arguments are
4614 taken as the number of rows and columns and any further arguments specify
4615 additional matrix dimensions.
4616 
4617 The optional argument @var{class} specifies the return type and may be
4618 either @qcode{"double"} or @qcode{"single"}.
4619 @seealso{e, I}
4620 @end deftypefn */)
4621 {
4622 #if defined (M_PI)
4623  double pi_val = M_PI;
4624 #else
4625  double pi_val = 4.0 * atan (1.0);
4626 #endif
4627 
4628  return fill_matrix (args, pi_val, "pi");
4629 }
4630 
4631 DEFUN (realmax, args, ,
4632  doc: /* -*- texinfo -*-
4633 @deftypefn {} {} realmax
4634 @deftypefnx {} {} realmax (@var{n})
4635 @deftypefnx {} {} realmax (@var{n}, @var{m})
4636 @deftypefnx {} {} realmax (@var{n}, @var{m}, @var{k}, @dots{})
4637 @deftypefnx {} {} realmax (@dots{}, @var{class})
4638 Return a scalar, matrix, or N-dimensional array whose elements are all equal
4639 to the largest floating point number that is representable.
4640 
4641 The actual value is system dependent. On machines that support IEEE
4642 floating point arithmetic, @code{realmax} is approximately
4643 @tex
4644 $1.7977\times10^{308}$ for double precision and $3.4028\times10^{38}$
4645 @end tex
4646 @ifnottex
4647 1.7977e+308 for double precision and 3.4028e+38
4648 @end ifnottex
4649 for single precision.
4650 
4651 When called with no arguments, return a scalar with the value
4652 @code{realmax (@qcode{"double"})}.
4653 
4654 When called with a single argument, return a square matrix with the
4655 dimension specified.
4656 
4657 When called with more than one scalar argument the first two arguments are
4658 taken as the number of rows and columns and any further arguments specify
4659 additional matrix dimensions.
4660 
4661 The optional argument @var{class} specifies the return type and may be
4662 either @qcode{"double"} or @qcode{"single"}.
4663 @seealso{realmin, intmax, flintmax, eps}
4664 @end deftypefn */)
4665 {
4666  return fill_matrix (args, std::numeric_limits<double>::max (),
4667  std::numeric_limits<float>::max (), "realmax");
4668 }
4669 
4670 DEFUN (realmin, args, ,
4671  doc: /* -*- texinfo -*-
4672 @deftypefn {} {} realmin
4673 @deftypefnx {} {} realmin (@var{n})
4674 @deftypefnx {} {} realmin (@var{n}, @var{m})
4675 @deftypefnx {} {} realmin (@var{n}, @var{m}, @var{k}, @dots{})
4676 @deftypefnx {} {} realmin (@dots{}, @var{class})
4677 Return a scalar, matrix, or N-dimensional array whose elements are all equal
4678 to the smallest normalized floating point number that is representable.
4679 
4680 The actual value is system dependent. On machines that support
4681 IEEE floating point arithmetic, @code{realmin} is approximately
4682 @tex
4683 $2.2251\times10^{-308}$ for double precision and $1.1755\times10^{-38}$
4684 @end tex
4685 @ifnottex
4686 2.2251e-308 for double precision and 1.1755e-38
4687 @end ifnottex
4688 for single precision.
4689 
4690 When called with no arguments, return a scalar with the value
4691 @code{realmin (@qcode{"double"})}.
4692 
4693 When called with a single argument, return a square matrix with the
4694 dimension specified.
4695 
4696 When called with more than one scalar argument the first two arguments are
4697 taken as the number of rows and columns and any further arguments specify
4698 additional matrix dimensions.
4699 
4700 The optional argument @var{class} specifies the return type and may be
4701 either @qcode{"double"} or @qcode{"single"}.
4702 @seealso{realmax, intmin, eps}
4703 @end deftypefn */)
4704 {
4705  return fill_matrix (args, std::numeric_limits<double>::min (),
4706  std::numeric_limits<float>::min (), "realmin");
4707 }
4708 
4709 DEFUN (I, args, ,
4710  doc: /* -*- texinfo -*-
4711 @c List other forms of function in documentation index
4712 @findex i
4713 @findex j
4714 @findex J
4715 
4716 @deftypefn {} {} I
4717 @deftypefnx {} {} I (@var{n})
4718 @deftypefnx {} {} I (@var{n}, @var{m})
4719 @deftypefnx {} {} I (@var{n}, @var{m}, @var{k}, @dots{})
4720 @deftypefnx {} {} I (@dots{}, @var{class})
4721 Return a scalar, matrix, or N-dimensional array whose elements are all equal
4722 to the pure imaginary unit, defined as
4723 @tex
4724 $\sqrt{-1}$.
4725 @end tex
4726 @ifnottex
4727 @w{@code{sqrt (-1)}}.
4728 @end ifnottex
4729 
4730 I, and its equivalents i, j, and J, are functions so any of the names may
4731 be reused for other purposes (such as i for a counter variable).
4732 
4733 When called with no arguments, return a scalar with the value @math{i}.
4734 
4735 When called with a single argument, return a square matrix with the
4736 dimension specified.
4737 
4738 When called with more than one scalar argument the first two arguments are
4739 taken as the number of rows and columns and any further arguments specify
4740 additional matrix dimensions.
4741 
4742 The optional argument @var{class} specifies the return type and may be
4743 either @qcode{"double"} or @qcode{"single"}.
4744 @seealso{e, pi, log, exp}
4745 @end deftypefn */)
4746 {
4747  return fill_matrix (args, Complex (0.0, 1.0), "I");
4748 }
4749 
4750 DEFALIAS (i, I);
4751 DEFALIAS (J, I);
4752 DEFALIAS (j, I);
4753 
4754 DEFUN (NA, args, ,
4755  doc: /* -*- texinfo -*-
4756 @deftypefn {} {} NA
4757 @deftypefnx {} {} NA (@var{n})
4758 @deftypefnx {} {} NA (@var{n}, @var{m})
4759 @deftypefnx {} {} NA (@var{n}, @var{m}, @var{k}, @dots{})
4760 @deftypefnx {} {} NA (@dots{}, @var{class})
4761 Return a scalar, matrix, or N-dimensional array whose elements are all equal
4762 to the special constant used to designate missing values.
4763 
4764 Note that NA always compares not equal to NA (NA != NA).
4765 To find NA values, use the @code{isna} function.
4766 
4767 When called with no arguments, return a scalar with the value @samp{NA}.
4768 
4769 When called with a single argument, return a square matrix with the
4770 dimension specified.
4771 
4772 When called with more than one scalar argument the first two arguments are
4773 taken as the number of rows and columns and any further arguments specify
4774 additional matrix dimensions.
4775 
4776 The optional argument @var{class} specifies the return type and may be
4777 either @qcode{"double"} or @qcode{"single"}.
4778 @seealso{isna}
4779 @end deftypefn */)
4780 {
4781  return fill_matrix (args, lo_ieee_na_value (),
4782  lo_ieee_float_na_value (), "NA");
4783 }
4784 
4785 /*
4786 %!assert (single (NA ("double")), NA ("single"))
4787 %!assert (double (NA ("single")), NA ("double"))
4788 */
4789 
4790 DEFUN (false, args, ,
4791  doc: /* -*- texinfo -*-
4792 @deftypefn {} {} false (@var{x})
4793 @deftypefnx {} {} false (@var{n}, @var{m})
4794 @deftypefnx {} {} false (@var{n}, @var{m}, @var{k}, @dots{})
4795 Return a matrix or N-dimensional array whose elements are all logical 0.
4796 
4797 If invoked with a single scalar integer argument, return a square
4798 matrix of the specified size.
4799 
4800 If invoked with two or more scalar integer arguments, or a vector of integer
4801 values, return an array with given dimensions.
4802 @seealso{true}
4803 @end deftypefn */)
4804 {
4805  return fill_matrix (args, false, "false");
4806 }
4807 
4808 DEFUN (true, args, ,
4809  doc: /* -*- texinfo -*-
4810 @deftypefn {} {} true (@var{x})
4811 @deftypefnx {} {} true (@var{n}, @var{m})
4812 @deftypefnx {} {} true (@var{n}, @var{m}, @var{k}, @dots{})
4813 Return a matrix or N-dimensional array whose elements are all logical 1.
4814 
4815 If invoked with a single scalar integer argument, return a square
4816 matrix of the specified size.
4817 
4818 If invoked with two or more scalar integer arguments, or a vector of integer
4819 values, return an array with given dimensions.
4820 @seealso{false}
4821 @end deftypefn */)
4822 {
4823  return fill_matrix (args, true, "true");
4824 }
4825 
4826 template <typename MT>
4828 identity_matrix (int nr, int nc)
4829 {
4831 
4832  typename MT::element_type one (1);
4833 
4834  if (nr == 1 && nc == 1)
4835  retval = one;
4836  else
4837  {
4838  dim_vector dims (nr, nc);
4839 
4840  typename MT::element_type zero (0);
4841 
4842  MT m (dims, zero);
4843 
4844  if (nr > 0 && nc > 0)
4845  {
4846  int n = std::min (nr, nc);
4847 
4848  for (int i = 0; i < n; i++)
4849  m(i,i) = one;
4850  }
4851 
4852  retval = m;
4853  }
4854 
4855  return retval;
4856 }
4857 
4858 #define INSTANTIATE_EYE(T) \
4859  template octave_value identity_matrix<T> (int, int)
4860 
4872 
4873 static octave_value
4875 {
4877 
4878  // FIXME: perhaps this should be made extensible by using
4879  // the class name to lookup a function to call to create the new
4880  // value.
4881 
4882  switch (dt)
4883  {
4885  retval = identity_matrix<int8NDArray> (nr, nc);
4886  break;
4887 
4890  break;
4891 
4894  break;
4895 
4898  break;
4899 
4902  break;
4903 
4906  break;
4907 
4910  break;
4911 
4914  break;
4915 
4917  retval = FloatDiagMatrix (nr, nc, 1.0f);
4918  break;
4919 
4921  retval = DiagMatrix (nr, nc, 1.0);
4922  break;
4923 
4926  break;
4927 
4928  default:
4929  error ("eye: invalid class name");
4930  break;
4931  }
4932 
4933  return retval;
4934 }
4935 
4936 #undef INT_EYE_MATRIX
4937 
4938 DEFUN (eye, args, ,
4939  doc: /* -*- texinfo -*-
4940 @deftypefn {} {} eye (@var{n})
4941 @deftypefnx {} {} eye (@var{m}, @var{n})
4942 @deftypefnx {} {} eye ([@var{m} @var{n}])
4943 @deftypefnx {} {} eye (@dots{}, @var{class})
4944 Return an identity matrix.
4945 
4946 If invoked with a single scalar argument @var{n}, return a square
4947 @nospell{NxN} identity matrix.
4948 
4949 If supplied two scalar arguments (@var{m}, @var{n}), @code{eye} takes them
4950 to be the number of rows and columns. If given a vector with two elements,
4951 @code{eye} uses the values of the elements as the number of rows and
4952 columns, respectively. For example:
4953 
4954 @example
4955 @group
4956 eye (3)
4957  @result{} 1 0 0
4958  0 1 0
4959  0 0 1
4960 @end group
4961 @end example
4962 
4963 The following expressions all produce the same result:
4964 
4965 @example
4966 @group
4967 eye (2)
4968 @equiv{}
4969 eye (2, 2)
4970 @equiv{}
4971 eye (size ([1, 2; 3, 4]))
4972 @end group
4973 @end example
4974 
4975 The optional argument @var{class}, allows @code{eye} to return an array of
4976 the specified type, like
4977 
4978 @example
4979 val = zeros (n,m, "uint8")
4980 @end example
4981 
4982 Calling @code{eye} with no arguments is equivalent to calling it with an
4983 argument of 1. Any negative dimensions are treated as zero. These odd
4984 definitions are for compatibility with @sc{matlab}.
4985 @seealso{speye, ones, zeros}
4986 @end deftypefn */)
4987 {
4988  int nargin = args.length ();
4989 
4991 
4992  // Check for type information.
4993 
4994  if (nargin > 0 && args(nargin-1).is_string ())
4995  {
4996  std::string nm = args(nargin-1).string_value ();
4997  nargin--;
4998 
5000  }
5001 
5002  if (nargin > 2)
5003  print_usage ();
5004 
5006 
5007  if (nargin == 0)
5008  retval = identity_matrix (1, 1, dt);
5009  else if (nargin == 1)
5010  {
5011  octave_idx_type nr, nc;
5012  get_dimensions (args(0), "eye", nr, nc);
5013 
5014  retval = identity_matrix (nr, nc, dt);
5015  }
5016  else
5017  {
5018  octave_idx_type nr, nc;
5019  get_dimensions (args(0), args(1), "eye", nr, nc);
5020 
5021  retval = identity_matrix (nr, nc, dt);
5022  }
5023 
5024  return retval;
5025 }
5026 
5027 /*
5028 %!assert (full (eye (3)), [1, 0, 0; 0, 1, 0; 0, 0, 1])
5029 %!assert (full (eye (2, 3)), [1, 0, 0; 0, 1, 0])
5030 
5031 %!assert (full (eye (3,"single")), single ([1, 0, 0; 0, 1, 0; 0, 0, 1]))
5032 %!assert (full (eye (2, 3,"single")), single ([1, 0, 0; 0, 1, 0]))
5033 
5034 %!assert (eye (3, "int8"), int8 ([1, 0, 0; 0, 1, 0; 0, 0, 1]))
5035 %!assert (eye (2, 3, "int8"), int8 ([1, 0, 0; 0, 1, 0]))
5036 
5037 %!error eye (1, 2, 3)
5038 */
5039 
5040 template <typename MT>
5041 static octave_value
5042 do_linspace (const octave_value& base, const octave_value& limit,
5043  octave_idx_type n)
5044 {
5045  typedef typename MT::column_vector_type CVT;
5046  typedef typename MT::element_type T;
5047 
5049 
5050  if (base.is_scalar_type ())
5051  {
5052  T bs = octave_value_extract<T> (base);
5053  if (limit.is_scalar_type ())
5054  {
5055  T ls = octave_value_extract<T> (limit);
5056  retval = linspace (bs, ls, n);
5057  }
5058  else
5059  {
5060  CVT lv = octave_value_extract<CVT> (limit);
5061  CVT bv (lv.numel (), bs);
5062  retval = linspace (bv, lv, n);
5063  }
5064  }
5065  else
5066  {
5067  CVT bv = octave_value_extract<CVT> (base);
5068  if (limit.is_scalar_type ())
5069  {
5070  T ls = octave_value_extract<T> (limit);
5071  CVT lv (bv.numel (), ls);
5072  retval = linspace (bv, lv, n);
5073  }
5074  else
5075  {
5076  CVT lv = octave_value_extract<CVT> (limit);
5077  retval = linspace (bv, lv, n);
5078  }
5079  }
5080 
5081  return retval;
5082 }
5083 
5084 DEFUN (linspace, args, ,
5085  doc: /* -*- texinfo -*-
5086 @deftypefn {} {} linspace (@var{start}, @var{end})
5087 @deftypefnx {} {} linspace (@var{start}, @var{end}, @var{n})
5088 Return a row vector with @var{n} linearly spaced elements between @var{start}
5089 and @var{end}.
5090 
5091 If the number of elements is greater than one, then the endpoints @var{start}
5092 and @var{end} are always included in the range. If @var{start} is greater than
5093 @var{end}, the elements are stored in decreasing order. If the number of
5094 points is not specified, a value of 100 is used.
5095 
5096 The @code{linspace} function returns a row vector when both @var{start} and
5097 @var{end} are scalars. If one, or both, inputs are vectors, then
5098 @code{linspace} transforms them to column vectors and returns a matrix where
5099 each row is an independent sequence between
5100 @w{@code{@var{start}(@var{row_n}), @var{end}(@var{row_n})}}.
5101 
5102 For compatibility with @sc{matlab}, return the second argument (@var{end}) when
5103 only a single value (@var{n} = 1) is requested.
5104 @seealso{colon, logspace}
5105 @end deftypefn */)
5106 {
5107  int nargin = args.length ();
5108 
5109  if (nargin != 2 && nargin != 3)
5110  print_usage ();
5111 
5112  octave_idx_type npoints = 100;
5113  if (nargin == 3)
5114  {
5115  // Apparently undocumented Matlab. If the third arg is an empty
5116  // numeric value, the number of points defaults to 1.
5117  octave_value arg_3 = args(2);
5118 
5119  if (arg_3.isnumeric () && arg_3.isempty ())
5120  npoints = 1;
5121  else if (! arg_3.is_scalar_type ())
5122  error ("linspace: N must be a scalar");
5123  else
5124  // Even if third arg is not an integer, it must be cast to int
5125  npoints = arg_3.idx_type_value ();
5126  }
5127 
5128  octave_value arg_1 = args(0);
5129  octave_value arg_2 = args(1);
5130 
5131  dim_vector sz1 = arg_1.dims ();
5132  bool isvector1 = sz1.ndims () == 2 && (sz1(0) == 1 || sz1(1) == 1);
5133  dim_vector sz2 = arg_2.dims ();
5134  bool isvector2 = sz2.ndims () == 2 && (sz2(0) == 1 || sz2(1) == 1);
5135 
5136  if (! isvector1 || ! isvector2)
5137  error ("linspace: START, END must be scalars or vectors");
5138 
5140 
5141  if (arg_1.is_single_type () || arg_2.is_single_type ())
5142  {
5143  if (arg_1.iscomplex () || arg_2.iscomplex ())
5144  retval = do_linspace<FloatComplexMatrix> (arg_1, arg_2, npoints);
5145  else
5146  retval = do_linspace<FloatMatrix> (arg_1, arg_2, npoints);
5147  }
5148  else
5149  {
5150  if (arg_1.iscomplex () || arg_2.iscomplex ())
5151  retval = do_linspace<ComplexMatrix> (arg_1, arg_2, npoints);
5152  else
5153  retval = do_linspace<Matrix> (arg_1, arg_2, npoints);
5154  }
5155 
5156  return retval;
5157 }
5158 
5159 /*
5160 %!test
5161 %! x1 = linspace (1, 2);
5162 %! x2 = linspace (1, 2, 10);
5163 %! x3 = linspace (1, -2, 10);
5164 %! assert (size (x1) == [1, 100] && x1(1) == 1 && x1(100) == 2);
5165 %! assert (x1(2) - x1(1), (2 - 1)/ (100 - 1), eps);
5166 %! assert (size (x2) == [1, 10] && x2(1) == 1 && x2(10) == 2);
5167 %! assert (x2(2) - x2(1), (2 - 1)/ (10 - 1), eps);
5168 %! assert (size (x3) == [1, 10] && x3(1) == 1 && x3(10) == -2);
5169 %! assert (x3(2) - x3(1), (-2 - 1)/ (10 - 1), eps);
5170 
5171 ## Test complex values
5172 %!test
5173 %! exp = [1+0i, 2-1.25i, 3-2.5i, 4-3.75i, 5-5i];
5174 %! obs = linspace (1, 5-5i, 5);
5175 %! assert (obs, exp);
5176 
5177 ## Test support for vectors in START and END
5178 %!assert (linspace ([1 2 3], [7 8 9]),
5179 %! [linspace(1, 7); linspace(2, 8); linspace(3, 9)], 10*eps)
5180 %!assert (linspace ([1 2 3]', [7 8 9]'),
5181 %! [linspace(1, 7); linspace(2, 8); linspace(3, 9)], 10*eps)
5182 %!assert (linspace ([1 2 3], 9),
5183 %! [linspace(1, 9); linspace(2, 9); linspace(3, 9)], 10*eps)
5184 %!assert (linspace ([1 2 3]', 9),
5185 %! [linspace(1, 9); linspace(2, 9); linspace(3, 9)], 10*eps)
5186 %!assert (linspace (1, [7 8 9]),
5187 %! [linspace(1, 7); linspace(1, 8); linspace(1, 9)], 10*eps)
5188 %!assert (linspace (1, [7 8 9]'),
5189 %! [linspace(1, 7); linspace(1, 8); linspace(1, 9)], 10*eps)
5190 
5191 ## Test class of output
5192 %!assert (class (linspace (1, 2)), "double")
5193 %!assert (class (linspace (single (1), 2)), "single")
5194 %!assert (class (linspace (1, single (2))), "single")
5195 
5196 ## Test obscure Matlab compatibility options
5197 %!assert (linspace (0, 1, []), 1)
5198 %!assert (linspace (10, 20, 2), [10 20])
5199 %!assert (linspace (10, 20, 1), [20])
5200 %!assert (linspace (10, 20, 0), zeros (1, 0))
5201 %!assert (linspace (10, 20, -1), zeros (1, 0))
5202 %!assert (numel (linspace (0, 1, 2+eps)), 2)
5203 %!assert (numel (linspace (0, 1, 2-eps)), 1)
5204 %!assert (linspace (10, 20, 2.1), [10 20])
5205 %!assert (linspace (10, 20, 2.9), [10 20])
5206 %!assert (1 ./ linspace (-0, 0, 4), [-Inf, Inf, Inf, Inf])
5207 
5208 %!error linspace ()
5209 %!error linspace (1, 2, 3, 4)
5210 %!error <N must be a scalar> linspace (1, 2, [3, 4])
5211 %!error <START, END must be scalars or vectors> linspace (ones (2,2), 2, 3)
5212 %!error <START, END must be scalars or vectors> linspace (2, ones (2,2), 3)
5213 %!error <START, END must be scalars or vectors> linspace (1, [], 3)
5214 */
5215 
5216 // FIXME: should accept dimensions as separate args for N-D
5217 // arrays as well as 1-D and 2-D arrays.
5218 
5219 DEFUN (resize, args, ,
5220  doc: /* -*- texinfo -*-
5221 @deftypefn {} {} resize (@var{x}, @var{m})
5222 @deftypefnx {} {} resize (@var{x}, @var{m}, @var{n}, @dots{})
5223 @deftypefnx {} {} resize (@var{x}, [@var{m} @var{n} @dots{}])
5224 Resize @var{x} cutting off elements as necessary.
5225 
5226 In the result, element with certain indices is equal to the corresponding
5227 element of @var{x} if the indices are within the bounds of @var{x};
5228 otherwise, the element is set to zero.
5229 
5230 In other words, the statement
5231 
5232 @example
5233 y = resize (x, dv)
5234 @end example
5235 
5236 @noindent
5237 is equivalent to the following code:
5238 
5239 @example
5240 @group
5241 y = zeros (dv, class (x));
5242 sz = min (dv, size (x));
5243 for i = 1:length (sz)
5244  idx@{i@} = 1:sz(i);
5245 endfor
5246 y(idx@{:@}) = x(idx@{:@});
5247 @end group
5248 @end example
5249 
5250 @noindent
5251 but is performed more efficiently.
5252 
5253 If only @var{m} is supplied, and it is a scalar, the dimension of the
5254 result is @var{m}-by-@var{m}.
5255 If @var{m}, @var{n}, @dots{} are all scalars, then the dimensions of
5256 the result are @var{m}-by-@var{n}-by-@dots{}.
5257 If given a vector as input, then the
5258 dimensions of the result are given by the elements of that vector.
5259 
5260 An object can be resized to more dimensions than it has;
5261 in such case the missing dimensions are assumed to be 1.
5262 Resizing an object to fewer dimensions is not possible.
5263 @seealso{reshape, postpad, prepad, cat}
5264 @end deftypefn */)
5265 {
5266  int nargin = args.length ();
5267 
5268  if (nargin < 2)
5269  print_usage ();
5270 
5272 
5273  if (nargin == 2)
5274  {
5275  Array<double> vec = args(1).vector_value ();
5276  int ndim = vec.numel ();
5277  if (ndim == 1)
5278  {
5279  octave_idx_type m = static_cast<octave_idx_type> (vec(0));
5280  retval = args(0);
5281  retval = retval.resize (dim_vector (m, m), true);
5282  }
5283  else
5284  {
5285  dim_vector dv;
5286  dv.resize (ndim);
5287  for (int i = 0; i < ndim; i++)
5288  dv(i) = static_cast<octave_idx_type> (vec(i));
5289  retval = args(0);
5290  retval = retval.resize (dv, true);
5291  }
5292  }
5293  else
5294  {
5295  dim_vector dv;
5296  dv.resize (nargin - 1);
5297  for (octave_idx_type i = 1; i < nargin; i++)
5298  dv(i-1) = static_cast<octave_idx_type> (args(i).scalar_value ());
5299 
5300  retval = args(0);
5301  retval = retval.resize (dv, true);
5302  }
5303 
5304  return retval;
5305 }
5306 
5307 // FIXME: should use octave_idx_type for dimensions.
5308 
5309 DEFUN (reshape, args, ,
5310  doc: /* -*- texinfo -*-
5311 @deftypefn {} {} reshape (@var{A}, @var{m}, @var{n}, @dots{})
5312 @deftypefnx {} {} reshape (@var{A}, [@var{m} @var{n} @dots{}])
5313 @deftypefnx {} {} reshape (@var{A}, @dots{}, [], @dots{})
5314 @deftypefnx {} {} reshape (@var{A}, @var{size})
5315 Return a matrix with the specified dimensions (@var{m}, @var{n}, @dots{})
5316 whose elements are taken from the matrix @var{A}.
5317 
5318 The elements of the matrix are accessed in column-major order (like Fortran
5319 arrays are stored).
5320 
5321 The following code demonstrates reshaping a 1x4 row vector into a 2x2 square
5322 matrix.
5323 
5324 @example
5325 @group
5326 reshape ([1, 2, 3, 4], 2, 2)
5327  @result{} 1 3
5328  2 4
5329 @end group
5330 @end example
5331 
5332 @noindent
5333 Note that the total number of elements in the original matrix
5334 (@code{prod (size (@var{A}))}) must match the total number of elements
5335 in the new matrix (@code{prod ([@var{m} @var{n} @dots{}])}).
5336 
5337 A single dimension of the return matrix may be left unspecified and Octave
5338 will determine its size automatically. An empty matrix ([]) is used to flag
5339 the unspecified dimension.
5340 @seealso{resize, vec, postpad, cat, squeeze}
5341 @end deftypefn */)
5342 {
5343  int nargin = args.length ();
5344 
5345  if (nargin < 2)
5346  print_usage ();
5347 
5349 
5350  dim_vector new_dims;
5351 
5352  if (nargin == 2)
5353  {
5354  Array<octave_idx_type> new_size = args(1).octave_idx_type_vector_value ();
5355 
5356  if (new_size.numel () < 2)
5357  error ("reshape: SIZE must have 2 or more dimensions");
5358 
5359  new_dims = dim_vector::alloc (new_size.numel ());
5360 
5361  for (octave_idx_type i = 0; i < new_size.numel (); i++)
5362  {
5363  if (new_size(i) < 0)
5364  error ("reshape: SIZE must be non-negative");
5365 
5366  new_dims(i) = new_size(i);
5367  }
5368  }
5369  else
5370  {
5371  new_dims = dim_vector::alloc (nargin-1);
5372  int empty_dim = -1;
5373 
5374  for (int i = 1; i < nargin; i++)
5375  {
5376  if (args(i).isempty ())
5377  {
5378  if (empty_dim > 0)
5379  error ("reshape: only a single dimension can be unknown");
5380 
5381  empty_dim = i;
5382  new_dims(i-1) = 1;
5383  }
5384  else
5385  {
5386  new_dims(i-1) = args(i).idx_type_value ();
5387 
5388  if (new_dims(i-1) < 0)
5389  error ("reshape: SIZE must be non-negative");
5390  }
5391  }
5392 
5393  if (empty_dim > 0)
5394  {
5395  octave_idx_type nel = new_dims.numel ();
5396 
5397  if (nel == 0)
5398  new_dims(empty_dim-1) = 0;
5399  else
5400  {
5401  octave_idx_type a_nel = args(0).numel ();
5402  octave_idx_type size_empty_dim = a_nel / nel;
5403 
5404  if (a_nel != size_empty_dim * nel)
5405  error ("reshape: SIZE is not divisible by the product of known dimensions (= %d)",
5406  nel);
5407 
5408  new_dims(empty_dim-1) = size_empty_dim;
5409  }
5410  }
5411  }
5412 
5413  retval = args(0).reshape (new_dims);
5414 
5415  return retval;
5416 }
5417 
5418 /*
5419 %!assert (size (reshape (ones (4, 4), 2, 8)), [2, 8])
5420 %!assert (size (reshape (ones (4, 4), 8, 2)), [8, 2])
5421 %!assert (size (reshape (ones (15, 4), 1, 60)), [1, 60])
5422 %!assert (size (reshape (ones (15, 4), 60, 1)), [60, 1])
5423 
5424 %!assert (size (reshape (ones (4, 4, "single"), 2, 8)), [2, 8])
5425 %!assert (size (reshape (ones (4, 4, "single"), 8, 2)), [8, 2])
5426 %!assert (size (reshape (ones (15, 4, "single"), 1, 60)), [1, 60])
5427 %!assert (size (reshape (ones (15, 4, "single"), 60, 1)), [60, 1])
5428 
5429 %!test
5430 %! s.a = 1;
5431 %! fail ("reshape (s, 2, 3)", "can't reshape 1x1 array to 2x3 array");
5432 
5433 %!error reshape ()
5434 %!error reshape (1, 2, 3, 4)
5435 %!error <SIZE must have 2 or more dimensions> reshape (1:3, 3)
5436 %!error <SIZE must be non-negative> reshape (1:3, [3 -1])
5437 %!error <only a single dimension can be unknown> reshape (1:3, 1,[],[],3)
5438 %!error <SIZE must be non-negative> reshape (1:3, 3, -1)
5439 %!error <SIZE is not divisible> reshape (1:3, 3, [], 2)
5440 */
5441 
5442 DEFUN (vec, args, ,
5443  doc: /* -*- texinfo -*-
5444 @deftypefn {} {@var{v} =} vec (@var{x})
5445 @deftypefnx {} {@var{v} =} vec (@var{x}, @var{dim})
5446 Return the vector obtained by stacking the columns of the matrix @var{x}
5447 one above the other.
5448 
5449 Without @var{dim} this is equivalent to @code{@var{x}(:)}.
5450 
5451 If @var{dim} is supplied, the dimensions of @var{v} are set to @var{dim}
5452 with all elements along the last dimension. This is equivalent to
5453 @code{shiftdim (@var{x}(:), 1-@var{dim})}.
5454 @seealso{vech, resize, cat}
5455 @end deftypefn */)
5456 {
5457  int nargin = args.length ();
5458 
5459  if (nargin < 1 || nargin > 2)
5460  print_usage ();
5461 
5462  int dim = 1;
5463  if (nargin == 2)
5464  {
5465  dim = args(1).idx_type_value ();
5466 
5467  if (dim < 1)
5468  error ("vec: DIM must be greater than zero");
5469  }
5470 
5472  octave_value arg = args(0);
5473 
5474  octave_value retval = arg.single_subsref ("(", colon);
5475 
5476  if (dim > 1)
5477  {
5478  dim_vector new_dims = dim_vector::alloc (dim);
5479 
5480  for (int i = 0; i < dim-1; i++)
5481  new_dims(i) = 1;
5482 
5483  new_dims(dim-1) = retval.numel ();
5484 
5485  retval = retval.reshape (new_dims);
5486  }
5487 
5488  return retval;
5489 }
5490 
5491 /*
5492 %!assert (vec ([1, 2; 3, 4]), [1; 3; 2; 4])
5493 %!assert (vec ([1, 3, 2, 4]), [1; 3; 2; 4])
5494 %!assert (vec ([1, 2, 3, 4], 2), [1, 2, 3, 4])
5495 %!assert (vec ([1, 2; 3, 4]), vec ([1, 2; 3, 4], 1))
5496 %!assert (vec ([1, 2; 3, 4], 1), [1; 3; 2; 4])
5497 %!assert (vec ([1, 2; 3, 4], 2), [1, 3, 2, 4])
5498 %!assert (vec ([1, 3; 2, 4], 3), reshape ([1, 2, 3, 4], 1, 1, 4))
5499 %!assert (vec ([1, 3; 2, 4], 3), shiftdim (vec ([1, 3; 2, 4]), -2))
5500 
5501 %!error vec ()
5502 %!error vec (1, 2, 3)
5503 %!error vec ([1, 2; 3, 4], 0)
5504 */
5505 
5506 DEFUN (squeeze, args, ,
5507  doc: /* -*- texinfo -*-
5508 @deftypefn {} {} squeeze (@var{x})
5509 Remove singleton dimensions from @var{x} and return the result.
5510 
5511 Note that for compatibility with @sc{matlab}, all objects have
5512 a minimum of two dimensions and row vectors are left unchanged.
5513 @seealso{reshape}
5514 @end deftypefn */)
5515 {
5516  if (args.length () != 1)
5517  print_usage ();
5518 
5519  return ovl (args(0).squeeze ());
5520 }
5521 
5522 DEFUN (full, args, ,
5523  doc: /* -*- texinfo -*-
5524 @deftypefn {} {@var{FM} =} full (@var{SM})
5525 Return a full storage matrix from a sparse, diagonal, or permutation matrix,
5526 or a range.
5527 @seealso{sparse, issparse}
5528 @end deftypefn */)
5529 {
5530  if (args.length () != 1)
5531  print_usage ();
5532 
5533  return ovl (args(0).full_value ());
5534 }
5535 
5536 // Compute various norms of the vector X.
5537 
5538 DEFUN (norm, args, ,
5539  doc: /* -*- texinfo -*-
5540 @deftypefn {} {} norm (@var{A})
5541 @deftypefnx {} {} norm (@var{A}, @var{p})
5542 @deftypefnx {} {} norm (@var{A}, @var{p}, @var{opt})
5543 Compute the p-norm of the matrix @var{A}.
5544 
5545 If the second argument is not given, @w{@code{p = 2}} is used.
5546 
5547 If @var{A} is a matrix (or sparse matrix):
5548 
5549 @table @asis
5550 @item @var{p} = @code{1}
5551 1-norm, the largest column sum of the absolute values of @var{A}.
5552 
5553 @item @var{p} = @code{2}
5554 Largest singular value of @var{A}.
5555 
5556 @item @var{p} = @code{Inf} or @qcode{"inf"}
5557 @cindex infinity norm
5558 Infinity norm, the largest row sum of the absolute values of @var{A}.
5559 
5560 @item @var{p} = @qcode{"fro"}
5561 @cindex @nospell{Frobenius} norm
5562 @nospell{Frobenius} norm of @var{A},
5563 @code{sqrt (sum (diag (@var{A}' * @var{A})))}.
5564 
5565 @item other @var{p}, @code{@var{p} > 1}
5566 @cindex general p-norm
5567 maximum @code{norm (A*x, p)} such that @code{norm (x, p) == 1}
5568 @end table
5569 
5570 If @var{A} is a vector or a scalar:
5571 
5572 @table @asis
5573 @item @var{p} = @code{Inf} or @qcode{"inf"}
5574 @code{max (abs (@var{A}))}.
5575 
5576 @item @var{p} = @code{-Inf}
5577 @code{min (abs (@var{A}))}.
5578 
5579 @item @var{p} = @qcode{"fro"}
5580 @nospell{Frobenius} norm of @var{A}, @code{sqrt (sumsq (abs (A)))}.
5581 
5582 @item @var{p} = 0
5583 Hamming norm---the number of nonzero elements.
5584 
5585 @item other @var{p}, @code{@var{p} > 1}
5586 p-norm of @var{A}, @code{(sum (abs (@var{A}) .^ @var{p})) ^ (1/@var{p})}.
5587 
5588 @item other @var{p} @code{@var{p} < 1}
5589 the p-pseudonorm defined as above.
5590 @end table
5591 
5592 If @var{opt} is the value @qcode{"rows"}, treat each row as a vector and
5593 compute its norm. The result is returned as a column vector.
5594 Similarly, if @var{opt} is @qcode{"columns"} or @qcode{"cols"} then
5595 compute the norms of each column and return a row vector.
5596 @seealso{normest, normest1, vecnorm, cond, svd}
5597 @end deftypefn */)
5598 {
5599  int nargin = args.length ();
5600 
5601  if (nargin < 1 || nargin > 3)
5602  print_usage ();
5603 
5604  octave_value x_arg = args(0);
5605 
5606  if (x_arg.ndims () != 2)
5607  error ("norm: only valid for 2-D objects");
5608 
5609  enum {sfmatrix, sfcols, sfrows, sffrob, sfinf, sfneginf} strflag = sfmatrix;
5610  if (nargin > 1 && args(nargin-1).is_string ())
5611  {
5612  std::string str = args(nargin-1).string_value ();
5613  std::transform (str.begin (), str.end (), str.begin (), tolower);
5614  if (str == "cols" || str == "columns")
5615  strflag = sfcols;
5616  else if (str == "rows")
5617  strflag = sfrows;
5618  else if (str == "fro")
5619  strflag = sffrob;
5620  else if (str == "inf")
5621  strflag = sfinf;
5622  else if (str == "-inf")
5623  strflag = sfneginf;
5624  else
5625  error ("norm: unrecognized option: %s", str.c_str ());
5626 
5627  // we've handled the last parameter, so act as if it was removed
5628  nargin--;
5629  }
5630 
5631  octave_value p_arg = (nargin > 1) ? args(1) : octave_value (2);
5632 
5633  if (p_arg.isempty ())
5634  p_arg = octave_value (2);
5635  else if (p_arg.is_string ())
5636  {
5637  std::string str = p_arg.string_value ();
5638  std::transform (str.begin (), str.end (), str.begin (), tolower);
5639  if (strflag != sfcols && strflag != sfrows)
5640  error ("norm: invalid combination of options");
5641 
5642  if (str == "cols" || str == "columns" || str == "rows")
5643  error ("norm: invalid combination of options");
5644 
5645  if (str == "fro")
5646  p_arg = octave_value (2);
5647  else if (str == "inf")
5649  else if (str == "-inf")
5651  else
5652  error ("norm: unrecognized option: %s", str.c_str ());
5653  }
5654  else if (! p_arg.is_scalar_type ())
5655  err_wrong_type_arg ("norm", p_arg);
5656 
5658 
5659  switch (strflag)
5660  {
5661  case sfmatrix:
5662  retval = xnorm (x_arg, p_arg);
5663  break;
5664 
5665  case sfcols:
5666  retval = xcolnorms (x_arg, p_arg);
5667  break;
5668 
5669  case sfrows:
5670  retval = xrownorms (x_arg, p_arg);
5671  break;
5672 
5673  case sffrob:
5674  retval = xfrobnorm (x_arg);
5675  break;
5676 
5677  case sfinf:
5679  break;
5680 
5681  case sfneginf:
5683  break;
5684  }
5685 
5686  return retval;
5687 }
5688 
5689 /*
5690 %!shared x
5691 %! x = [1, -3, 4, 5, -7];
5692 %!assert (norm (x,0), 5)
5693 %!assert (norm (x,1), 20)
5694 %!assert (norm (x,2), 10)
5695 %!assert (norm (x,3), 8.24257059961711, -4*eps)
5696 %!assert (norm (x,Inf), 7)
5697 %!assert (norm (x,-Inf), 1)
5698 %!assert (norm (x,"inf"), 7)
5699 %!assert (norm (x,"-Inf"), 1)
5700 %!assert (norm (x,"fro"), 10, -eps)
5701 %!assert (norm (x), 10)
5702 %!assert (norm ([1e200, 1]), 1e200)
5703 %!assert (norm ([3+4i, 3-4i, sqrt(31)]), 9, -4*eps)
5704 %!shared m
5705 %! m = magic (4);
5706 %!assert (norm (m,1), 34)
5707 %!assert (norm (m,2), 34, -eps)
5708 %!assert (norm (m,3), 34, -sqrt (eps))
5709 %!assert (norm (m,Inf), 34)
5710 %!assert (norm (m,"inf"), 34)
5711 %!shared m2, flo, fhi
5712 %! m2 = [1,2;3,4];
5713 %! flo = 1e-300;
5714 %! fhi = 1e+300;
5715 %!assert (norm (flo*m2,"fro"), sqrt (30)*flo, -eps)
5716 %!assert (norm (fhi*m2,"fro"), sqrt (30)*fhi, -eps)
5717 
5718 %!shared x
5719 %! x = single ([1, -3, 4, 5, -7]);
5720 %!assert (norm (x,0), single (5))
5721 %!assert (norm (x,1), single (20))
5722 %!assert (norm (x,2), single (10))
5723 %!assert (norm (x,3), single (8.24257059961711), -4*eps ("single"))
5724 %!assert (norm (x,Inf), single (7))
5725 %!assert (norm (x,-Inf), single (1))
5726 %!assert (norm (x,"inf"), single (7))
5727 %!assert (norm (x,"-Inf"), single (1))
5728 %!assert (norm (x,"fro"), single (10), -eps ("single"))
5729 %!assert (norm (x), single (10))
5730 %!assert (norm (single ([1e38, 1])), single (1e38))
5731 %!assert (norm (single ([3+4i, 3-4i, sqrt(31)])), single (9), -4*eps ("single"))
5732 %!shared m
5733 %! m = single (magic (4));
5734 %!assert (norm (m,1), single (34))
5735 %!assert (norm (m,2), single (34), -eps ("single"))
5736 %!assert (norm (m,3), single (34), -sqrt (eps ("single")))
5737 %!assert (norm (m,Inf), single (34))
5738 %!assert (norm (m,"inf"), single (34))
5739 %!shared m2, flo, fhi
5740 %! m2 = single ([1,2;3,4]);
5741 %! flo = single (1e-300);
5742 %! fhi = single (1e+300);
5743 %!assert (norm (flo*m2,"fro"), single (sqrt (30)*flo), -eps ("single"))
5744 %!assert (norm (fhi*m2,"fro"), single (sqrt (30)*fhi), -eps ("single"))
5745 
5746 ## Hamming norm (p == 0)
5747 %!assert (norm ([1, 0, 0, 0, 1], 0), 2)
5748 
5749 %!shared q
5750 %! q = rand (1e3, 3);
5751 %!assert (norm (q, 3, "rows"), sum (q.^3, 2).^(1/3), sqrt (eps))
5752 %!assert (norm (q, "fro", "rows"), sum (q.^2, 2).^(1/2), sqrt (eps))
5753 %!assert (norm (q, "fro", "rows"), sqrt (sumsq (q, 2)), sqrt (eps))
5754 %!assert (norm (q, "fro", "cols"), sqrt (sumsq (q, 1)), sqrt (eps))
5755 %!assert (norm (q, 3, "cols"), sum (q.^3, 1).^(1/3), sqrt (eps))
5756 %!assert (norm (q, "inf", "rows"), norm (q, Inf, "rows"))
5757 %!assert (norm (q, "inf", "cols"), norm (q, Inf, "cols"))
5758 %!assert (norm (q, [], "rows"), norm (q, 2, "rows"))
5759 %!assert (norm (q, [], "cols"), norm (q, 2, "cols"))
5760 
5761 %!test <30631>
5762 %! ## Test for norm returning NaN on sparse matrix
5763 %! A = sparse (2,2);
5764 %! A(2,1) = 1;
5765 %! assert (norm (A), 1);
5766 
5767 ## Test input validation
5768 %!error norm ()
5769 %!error norm (1,2,3,4)
5770 %!error <unrecognized option> norm (1, "invalid")
5771 %!error <unrecognized option> norm (1, "rows", "invalid")
5772 %!error <unrecognized option> norm (1, "invalid", "rows")
5773 %!error <invalid combination of options> norm (1, "cols", "rows")
5774 %!error <invalid combination of options> norm (1, "rows", "rows")
5775 %!error <p must be .= 1> norm (ones (2,2), -Inf)
5776 */
5777 
5778 static octave_value
5779 unary_op_defun_body (octave_value::unary_op op,
5780  const octave_value_list& args)
5781 {
5782  if (args.length () != 1)
5783  print_usage ();
5784 
5785  return do_unary_op (op, args(0));
5786 }
5787 
5788 DEFUN (not, args, ,
5789  doc: /* -*- texinfo -*-
5790 @deftypefn {} {@var{z} =} not (@var{x})
5791 Return the logical NOT of @var{x}.
5792 
5793 This function is equivalent to the operator syntax @w{@code{! @var{x}}}.
5794 @seealso{and, or, xor}
5795 @end deftypefn */)
5796 {
5797  return unary_op_defun_body (octave_value::op_not, args);
5798 }
5799 
5800 DEFUN (uplus, args, ,
5801  doc: /* -*- texinfo -*-
5802 @deftypefn {} {} uplus (@var{x})
5803 This function and @w{@tcode{+ @var{x}}} are equivalent.
5804 @seealso{uminus, plus, minus}
5805 @end deftypefn */)
5806 {
5807  return unary_op_defun_body (octave_value::op_uplus, args);
5808 }
5809 
5810 DEFUN (uminus, args, ,
5811  doc: /* -*- texinfo -*-
5812 @deftypefn {} {} uminus (@var{x})
5813 This function and @w{@tcode{- @var{x}}} are equivalent.
5814 @seealso{uplus, minus}
5815 @end deftypefn */)
5816 {
5817  return unary_op_defun_body (octave_value::op_uminus, args);
5818 }
5819 
5820 DEFUN (transpose, args, ,
5821  doc: /* -*- texinfo -*-
5822 @deftypefn {} {} transpose (@var{x})
5823 Return the transpose of @var{x}.
5824 
5825 This function and @tcode{@var{x}.'} are equivalent.
5826 @seealso{ctranspose}
5827 @end deftypefn */)
5828 {
5829  return unary_op_defun_body (octave_value::op_transpose, args);
5830 }
5831 
5832 /*
5833 %!assert (2.', 2)
5834 %!assert (2i.', 2i)
5835 %!assert ([1:4].', [1;2;3;4])
5836 %!assert ([1;2;3;4].', [1:4])
5837 %!assert ([1,2;3,4].', [1,3;2,4])
5838 %!assert ([1,2i;3,4].', [1,3;2i,4])
5839 
5840 %!assert (transpose ([1,2;3,4]), [1,3;2,4])
5841 
5842 %!assert (single (2).', single (2))
5843 %!assert (single (2i).', single (2i))
5844 %!assert (single ([1:4]).', single ([1;2;3;4]))
5845 %!assert (single ([1;2;3;4]).', single ([1:4]))
5846 %!assert (single ([1,2;3,4]).', single ([1,3;2,4]))
5847 %!assert (single ([1,2i;3,4]).', single ([1,3;2i,4]))
5848 
5849 %!assert (transpose (single ([1,2;3,4])), single ([1,3;2,4]))
5850 */
5851 
5852 DEFUN (ctranspose, args, ,
5853  doc: /* -*- texinfo -*-
5854 @deftypefn {} {} ctranspose (@var{x})
5855 Return the complex conjugate transpose of @var{x}.
5856 
5857 This function and @tcode{@var{x}'} are equivalent.
5858 @seealso{transpose}
5859 @end deftypefn */)
5860 {
5861  return unary_op_defun_body (octave_value::op_hermitian, args);
5862 }
5863 
5864 /*
5865 %!assert (2', 2)
5866 %!assert (2i', -2i)
5867 %!assert ([1:4]', [1;2;3;4])
5868 %!assert ([1;2;3;4]', [1:4])
5869 %!assert ([1,2;3,4]', [1,3;2,4])
5870 %!assert ([1,2i;3,4]', [1,3;-2i,4])
5871 
5872 %!assert (ctranspose ([1,2i;3,4]), [1,3;-2i,4])
5873 
5874 %!assert (single (2)', single (2))
5875 %!assert (single (2i)', single (-2i))
5876 %!assert (single ([1:4])', single ([1;2;3;4]))
5877 %!assert (single ([1;2;3;4])', single ([1:4]))
5878 %!assert (single ([1,2;3,4])', single ([1,3;2,4]))
5879 %!assert (single ([1,2i;3,4])', single ([1,3;-2i,4]))
5880 
5881 %!assert (ctranspose (single ([1,2i;3,4])), single ([1,3;-2i,4]))
5882 */
5883 
5884 static octave_value
5885 binary_op_defun_body (octave_value::binary_op op,
5886  const octave_value_list& args)
5887 {
5888  if (args.length () != 2)
5889  print_usage ();
5890 
5891  return do_binary_op (op, args(0), args(1));
5892 }
5893 
5894 static octave_value
5895 binary_assoc_op_defun_body (octave_value::binary_op op,
5897  const octave_value_list& args)
5898 {
5899  int nargin = args.length ();
5900 
5901  if (nargin == 0)
5902  print_usage ();
5903 
5905 
5906  if (nargin == 1)
5907  retval = args(0);
5908  else if (nargin == 2)
5909  retval = do_binary_op (op, args(0), args(1));
5910  else
5911  {
5912  retval = do_binary_op (op, args(0), args(1));
5913 
5914  for (int i = 2; i < nargin; i++)
5915  retval.assign (aop, args(i));
5916  }
5917 
5918  return retval;
5919 }
5920 
5921 DEFUN (plus, args, ,
5922  doc: /* -*- texinfo -*-
5923 @deftypefn {} {} plus (@var{x}, @var{y})
5924 @deftypefnx {} {} plus (@var{x1}, @var{x2}, @dots{})
5925 This function and @w{@tcode{@var{x} + @var{y}}} are equivalent.
5926 
5927 If more arguments are given, the summation is applied
5928 cumulatively from left to right:
5929 
5930 @example
5931 (@dots{}((@var{x1} + @var{x2}) + @var{x3}) + @dots{})
5932 @end example
5933 
5934 At least one argument is required.
5935 @seealso{minus, uplus}
5936 @end deftypefn */)
5937 {
5938  return binary_assoc_op_defun_body (octave_value::op_add,
5939  octave_value::op_add_eq, args);
5940 }
5941 
5942 DEFUN (minus, args, ,
5943  doc: /* -*- texinfo -*-
5944 @deftypefn {} {} minus (@var{x}, @var{y})
5945 This function and @w{@tcode{@var{x} - @var{y}}} are equivalent.
5946 @seealso{plus, uminus}
5947 @end deftypefn */)
5948 {
5949  return binary_op_defun_body (octave_value::op_sub, args);
5950 }
5951 
5952 DEFUN (mtimes, args, ,
5953  doc: /* -*- texinfo -*-
5954 @deftypefn {} {} mtimes (@var{x}, @var{y})
5955 @deftypefnx {} {} mtimes (@var{x1}, @var{x2}, @dots{})
5956 Return the matrix multiplication product of inputs.
5957 
5958 This function and @w{@tcode{@var{x} * @var{y}}} are equivalent.
5959 If more arguments are given, the multiplication is applied
5960 cumulatively from left to right:
5961 
5962 @example
5963 (@dots{}((@var{x1} * @var{x2}) * @var{x3}) * @dots{})
5964 @end example
5965 
5966 At least one argument is required.
5967 @seealso{times, plus, minus, rdivide, mrdivide, mldivide, mpower}
5968 @end deftypefn */)
5969 {
5970  return binary_assoc_op_defun_body (octave_value::op_mul,
5971  octave_value::op_mul_eq, args);
5972 }
5973 
5974 DEFUN (mrdivide, args, ,
5975  doc: /* -*- texinfo -*-
5976 @deftypefn {} {} mrdivide (@var{x}, @var{y})
5977 Return the matrix right division of @var{x} and @var{y}.
5978 
5979 This function and @w{@tcode{@var{x} / @var{y}}} are equivalent.
5980 @seealso{mldivide, rdivide, plus, minus}
5981 @end deftypefn */)
5982 {
5983  return binary_op_defun_body (octave_value::op_div, args);
5984 }
5985 
5986 DEFUN (mpower, args, ,
5987  doc: /* -*- texinfo -*-
5988 @deftypefn {} {} mpower (@var{x}, @var{y})
5989 Return the matrix power operation of @var{x} raised to the @var{y} power.
5990 
5991 This function and @w{@tcode{@var{x} ^ @var{y}}} are equivalent.
5992 @seealso{power, mtimes, plus, minus}
5993 @end deftypefn */)
5994 {
5995  return binary_op_defun_body (octave_value::op_pow, args);
5996 }
5997 
5998 DEFUN (mldivide, args, ,
5999  doc: /* -*- texinfo -*-
6000 @deftypefn {} {} mldivide (@var{x}, @var{y})
6001 Return the matrix left division of @var{x} and @var{y}.
6002 
6003 This function and @w{@tcode{@var{x} @xbackslashchar{} @var{y}}} are equivalent.
6004 @seealso{mrdivide, ldivide, rdivide}
6005 @end deftypefn */)
6006 {
6007  return binary_op_defun_body (octave_value::op_ldiv, args);
6008 }
6009 
6010 DEFUN (lt, args, ,
6011  doc: /* -*- texinfo -*-
6012 @deftypefn {} {} lt (@var{x}, @var{y})
6013 This function is equivalent to @w{@code{@var{x} < @var{y}}}.
6014 @seealso{le, eq, ge, gt, ne}
6015 @end deftypefn */)
6016 {
6017  return binary_op_defun_body (octave_value::op_lt, args);
6018 }
6019 
6020 DEFUN (le, args, ,
6021  doc: /* -*- texinfo -*-
6022 @deftypefn {} {} le (@var{x}, @var{y})
6023 This function is equivalent to @w{@code{@var{x} <= @var{y}}}.
6024 @seealso{eq, ge, gt, ne, lt}
6025 @end deftypefn */)
6026 {
6027  return binary_op_defun_body (octave_value::op_le, args);
6028 }
6029 
6030 DEFUN (eq, args, ,
6031  doc: /* -*- texinfo -*-
6032 @deftypefn {} {} eq (@var{x}, @var{y})
6033 Return true if the two inputs are equal.
6034 
6035 This function is equivalent to @w{@code{@var{x} == @var{y}}}.
6036 @seealso{ne, isequal, le, ge, gt, ne, lt}
6037 @end deftypefn */)
6038 {
6039  return binary_op_defun_body (octave_value::op_eq, args);
6040 }
6041 
6042 DEFUN (ge, args, ,
6043  doc: /* -*- texinfo -*-
6044 @deftypefn {} {} ge (@var{x}, @var{y})
6045 This function is equivalent to @w{@code{@var{x} >= @var{y}}}.
6046 @seealso{le, eq, gt, ne, lt}
6047 @end deftypefn */)
6048 {
6049  return binary_op_defun_body (octave_value::op_ge, args);
6050 }
6051 
6052 DEFUN (gt, args, ,
6053  doc: /* -*- texinfo -*-
6054 @deftypefn {} {} gt (@var{x}, @var{y})
6055 This function is equivalent to @w{@code{@var{x} > @var{y}}}.
6056 @seealso{le, eq, ge, ne, lt}
6057 @end deftypefn */)
6058 {
6059  return binary_op_defun_body (octave_value::op_gt, args);
6060 }
6061 
6062 DEFUN (ne, args, ,
6063  doc: /* -*- texinfo -*-
6064 @deftypefn {} {} ne (@var{x}, @var{y})
6065 Return true if the two inputs are not equal.
6066 
6067 This function is equivalent to @w{@code{@var{x} != @var{y}}}.
6068 @seealso{eq, isequal, le, ge, lt}
6069 @end deftypefn */)
6070 {
6071  return binary_op_defun_body (octave_value::op_ne, args);
6072 }
6073 
6074 DEFUN (times, args, ,
6075  doc: /* -*- texinfo -*-
6076 @deftypefn {} {} times (@var{x}, @var{y})
6077 @deftypefnx {} {} times (@var{x1}, @var{x2}, @dots{})
6078 Return the element-by-element multiplication product of inputs.
6079 
6080 This function and @w{@tcode{@var{x} .* @var{y}}} are equivalent.
6081 If more arguments are given, the multiplication is applied
6082 cumulatively from left to right:
6083 
6084 @example
6085 (@dots{}((@var{x1} .* @var{x2}) .* @var{x3}) .* @dots{})
6086 @end example
6087 
6088 At least one argument is required.
6089 @seealso{mtimes, rdivide}
6090 @end deftypefn */)
6091 {
6092  return binary_assoc_op_defun_body (octave_value::op_el_mul,
6094 }
6095 
6096 DEFUN (rdivide, args, ,
6097  doc: /* -*- texinfo -*-
6098 @deftypefn {} {} rdivide (@var{x}, @var{y})
6099 Return the element-by-element right division of @var{x} and @var{y}.
6100 
6101 This function and @w{@tcode{@var{x} ./ @var{y}}} are equivalent.
6102 @seealso{ldivide, mrdivide, times, plus}
6103 @end deftypefn */)
6104 {
6105  return binary_op_defun_body (octave_value::op_el_div, args);
6106 }
6107 
6108 DEFUN (power, args, ,
6109  doc: /* -*- texinfo -*-
6110 @deftypefn {} {} power (@var{x}, @var{y})
6111 Return the element-by-element operation of @var{x} raised to the
6112 @var{y} power.
6113 
6114 This function and @w{@tcode{@var{x} .^ @var{y}}} are equivalent.
6115 
6116 If several complex results are possible, returns the one with smallest
6117 non-negative argument (angle). Use @code{realpow}, @code{realsqrt},
6118 @code{cbrt}, or @code{nthroot} if a real result is preferred.
6119 
6120 @seealso{mpower, realpow, realsqrt, cbrt, nthroot}
6121 @end deftypefn */)
6122 {
6123  return binary_op_defun_body (octave_value::op_el_pow, args);
6124 }
6125 
6126 DEFUN (ldivide, args, ,
6127  doc: /* -*- texinfo -*-
6128 @deftypefn {} {} ldivide (@var{x}, @var{y})
6129 Return the element-by-element left division of @var{x} and @var{y}.
6130 
6131 This function and @w{@tcode{@var{x} .@xbackslashchar{} @var{y}}} are
6132 equivalent.
6133 @seealso{rdivide, mldivide, times, plus}
6134 @end deftypefn */)
6135 {
6136  return binary_op_defun_body (octave_value::op_el_ldiv, args);
6137 }
6138 
6139 DEFUN (and, args, ,
6140  doc: /* -*- texinfo -*-
6141 @deftypefn {} {@var{z} =} and (@var{x}, @var{y})
6142 @deftypefnx {} {@var{z} =} and (@var{x1}, @var{x2}, @dots{})
6143 Return the logical AND of @var{x} and @var{y}.
6144 
6145 This function is equivalent to the operator syntax
6146 @w{@code{@var{x} & @var{y}}}. If more than two arguments are given, the
6147 logical AND is applied cumulatively from left to right:
6148 
6149 @example
6150 (@dots{}((@var{x1} & @var{x2}) & @var{x3}) & @dots{})
6151 @end example
6152 
6153 At least one argument is required.
6154 @seealso{or, not, xor}
6155 @end deftypefn */)
6156 {
6157  return binary_assoc_op_defun_body (octave_value::op_el_and,
6159 }
6160 
6161 DEFUN (or, args, ,
6162  doc: /* -*- texinfo -*-
6163 @deftypefn {} {@var{z} =} or (@var{x}, @var{y})
6164 @deftypefnx {} {@var{z} =} or (@var{x1}, @var{x2}, @dots{})
6165 Return the logical OR of @var{x} and @var{y}.
6166 
6167