changelog shortlog tags changeset files revisions annotate raw

src/data.cc

changeset 10289: 4b124317dc38
parent:703038d648f1
author: John W. Eaton <jwe@octave.org>
date: Tue Feb 09 20:58:55 2010 -0500 (68 minutes ago)
permissions: -rw-r--r--
description: base_properties::set_children: account for hidden children
1/*
2
3Copyright (C) 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002,
4 2003, 2004, 2005, 2006, 2007, 2008, 2009 John W. Eaton
5Copyright (C) 2009 Jaroslav Hajek
6Copyright (C) 2009, 2010 VZLU Prague
7
8This file is part of Octave.
9
10Octave is free software; you can redistribute it and/or modify it
11under the terms of the GNU General Public License as published by the
12Free Software Foundation; either version 3 of the License, or (at your
13option) any later version.
14
15Octave is distributed in the hope that it will be useful, but WITHOUT
16ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18for more details.
19
20You should have received a copy of the GNU General Public License
21along with Octave; see the file COPYING. If not, see
22<http://www.gnu.org/licenses/>.
23
24*/
25
26#ifdef HAVE_CONFIG_H
27#include <config.h>
28#endif
29
30#include <sys/types.h>
31
32#ifdef HAVE_SYS_RESOURCE_H
33#include <sys/resource.h>
34#endif
35
36#include <cfloat>
37
38#include <string>
39
40#include "lo-ieee.h"
41#include "lo-math.h"
42#include "oct-time.h"
43#include "str-vec.h"
44#include "quit.h"
45#include "mx-base.h"
46
47#include "Cell.h"
48#include "defun.h"
49#include "error.h"
50#include "gripes.h"
51#include "oct-map.h"
52#include "oct-obj.h"
53#include "ov.h"
54#include "ov-float.h"
55#include "ov-complex.h"
56#include "ov-flt-complex.h"
57#include "ov-cx-mat.h"
58#include "ov-flt-cx-mat.h"
59#include "ov-cx-sparse.h"
60#include "parse.h"
61#include "pt-mat.h"
62#include "utils.h"
63#include "variables.h"
64#include "pager.h"
65#include "xnorm.h"
66
67#if ! defined (HAVE_HYPOTF) && defined (HAVE__HYPOTF)
68#define hypotf _hypotf
69#define HAVE_HYPOTF 1
70#endif
71
72#define ANY_ALL(FCN) \
73 \
74 octave_value retval; \
75 \
76 int nargin = args.length (); \
77 \
78 if (nargin == 1 || nargin == 2) \
79 { \
80 int dim = (nargin == 1 ? -1 : args(1).int_value (true) - 1); \
81 \
82 if (! error_state) \
83 { \
84 if (dim >= -1) \
85 retval = args(0).FCN (dim); \
86 else \
87 error (#FCN ": invalid dimension argument = %d", dim + 1); \
88 } \
89 else \
90 error (#FCN ": expecting dimension argument to be an integer"); \
91 } \
92 else \
93 print_usage (); \
94 \
95 return retval
96
97DEFUN (all, args, ,
98 "-*- texinfo -*-\n\
99@deftypefn {Built-in Function} {} all (@var{x}, @var{dim})\n\
100The function @code{all} behaves like the function @code{any}, except\n\
101that it returns true only if all the elements of a vector, or all the\n\
102elements along dimension @var{dim} of a matrix, are nonzero.\n\
103@end deftypefn")
104{
105 ANY_ALL (all);
106}
107
108/*
109
110%!test
111%! x = ones (3);
112%! x(1,1) = 0;
113%! assert((all (all (rand (3) + 1) == [1, 1, 1]) == 1
114%! && all (all (x) == [0, 1, 1]) == 1
115%! && all (x, 1) == [0, 1, 1]
116%! && all (x, 2) == [0; 1; 1]));
117
118%!test
119%! x = ones (3, 'single');
120%! x(1,1) = 0;
121%! assert((all (all (single (rand (3) + 1)) == [1, 1, 1]) == 1
122%! && all (all (x) == [0, 1, 1]) == 1
123%! && all (x, 1) == [0, 1, 1]
124%! && all (x, 2) == [0; 1; 1]));
125
126%!error <Invalid call to all.*> all ();
127%!error <Invalid call to all.*> all (1, 2, 3);
128
129 */
130
131DEFUN (any, args, ,
132 "-*- texinfo -*-\n\
133@deftypefn {Built-in Function} {} any (@var{x}, @var{dim})\n\
134For a vector argument, return 1 if any element of the vector is\n\
135nonzero.\n\
136\n\
137For a matrix argument, return a row vector of ones and\n\
138zeros with each element indicating whether any of the elements of the\n\
139corresponding column of the matrix are nonzero. For example,\n\
140\n\
141@example\n\
142@group\n\
143any (eye (2, 4))\n\
144 @result{} [ 1, 1, 0, 0 ]\n\
145@end group\n\
146@end example\n\
147\n\
148If the optional argument @var{dim} is supplied, work along dimension\n\
149@var{dim}. For example,\n\
150\n\
151@example\n\
152@group\n\
153any (eye (2, 4), 2)\n\
154 @result{} [ 1; 1 ]\n\
155@end group\n\
156@end example\n\
157@end deftypefn")
158{
159 ANY_ALL (any);
160}
161
162/*
163
164%!test
165%! x = zeros (3);
166%! x(3,3) = 1;
167%! assert((all (any (x) == [0, 0, 1]) == 1
168%! && all (any (ones (3)) == [1, 1, 1]) == 1
169%! && any (x, 1) == [0, 0, 1]
170%! && any (x, 2) == [0; 0; 1]));
171
172%!test
173%! x = zeros (3,'single');
174%! x(3,3) = 1;
175%! assert((all (any (x) == [0, 0, 1]) == 1
176%! && all (any (ones (3, 'single')) == [1, 1, 1]) == 1
177%! && any (x, 1) == [0, 0, 1]
178%! && any (x, 2) == [0; 0; 1]));
179
180%!error <Invalid call to any.*> any ();
181%!error <Invalid call to any.*> any (1, 2, 3);
182
183 */
184
185// These mapping functions may also be useful in other places, eh?
186
187typedef double (*d_dd_fcn) (double, double);
188typedef float (*f_ff_fcn) (float, float);
189
190static NDArray
191map_d_m (d_dd_fcn f, double x, const NDArray& y)
192{
193 NDArray retval (y.dims ());
194 double *r_data = retval.fortran_vec ();
195
196 const double *y_data = y.data ();
197
198 octave_idx_type nel = y.numel ();
199
200 for (octave_idx_type i = 0; i < nel; i++)
201 {
202 octave_quit ();
203 r_data[i] = f (x, y_data[i]);
204 }
205
206 return retval;
207}
208
209static FloatNDArray
210map_f_fm (f_ff_fcn f, float x, const FloatNDArray& y)
211{
212 FloatNDArray retval (y.dims ());
213 float *r_data = retval.fortran_vec ();
214
215 const float *y_data = y.data ();
216
217 octave_idx_type nel = y.numel ();
218
219 for (octave_idx_type i = 0; i < nel; i++)
220 {
221 octave_quit ();
222 r_data[i] = f (x, y_data[i]);
223 }
224
225 return retval;
226}
227
228static NDArray
229map_m_d (d_dd_fcn f, const NDArray& x, double y)
230{
231 NDArray retval (x.dims ());
232 double *r_data = retval.fortran_vec ();
233
234 const double *x_data = x.data ();
235
236 octave_idx_type nel = x.numel ();
237
238 for (octave_idx_type i = 0; i < nel; i++)
239 {
240 octave_quit ();
241 r_data[i] = f (x_data[i], y);
242 }
243
244 return retval;
245}
246
247static FloatNDArray
248map_fm_f (f_ff_fcn f, const FloatNDArray& x, float y)
249{
250 FloatNDArray retval (x.dims ());
251 float *r_data = retval.fortran_vec ();
252
253 const float *x_data = x.data ();
254
255 octave_idx_type nel = x.numel ();
256
257 for (octave_idx_type i = 0; i < nel; i++)
258 {
259 octave_quit ();
260 r_data[i] = f (x_data[i], y);
261 }
262
263 return retval;
264}
265
266static NDArray
267map_m_m (d_dd_fcn f, const NDArray& x, const NDArray& y)
268{
269 assert (x.dims () == y.dims ());
270
271 NDArray retval (x.dims ());
272 double *r_data = retval.fortran_vec ();
273
274 const double *x_data = x.data ();
275 const double *y_data = y.data ();
276
277 octave_idx_type nel = x.numel ();
278
279 for (octave_idx_type i = 0; i < nel; i++)
280 {
281 octave_quit ();
282 r_data[i] = f (x_data[i], y_data[i]);
283 }
284
285 return retval;
286}
287
288static FloatNDArray
289map_fm_fm (f_ff_fcn f, const FloatNDArray& x, const FloatNDArray& y)
290{
291 assert (x.dims () == y.dims ());
292
293 FloatNDArray retval (x.dims ());
294 float *r_data = retval.fortran_vec ();
295
296 const float *x_data = x.data ();
297 const float *y_data = y.data ();
298
299 octave_idx_type nel = x.numel ();
300
301 for (octave_idx_type i = 0; i < nel; i++)
302 {
303 octave_quit ();
304 r_data[i] = f (x_data[i], y_data[i]);
305 }
306
307 return retval;
308}
309
310static SparseMatrix
311map_d_s (d_dd_fcn f, double x, const SparseMatrix& y)
312{
313 octave_idx_type nr = y.rows ();
314 octave_idx_type nc = y.columns ();
315 SparseMatrix retval;
316 double f_zero = f (x, 0.);
317
318 if (f_zero != 0)
319 {
320 retval = SparseMatrix (nr, nc, f_zero);
321
322 for (octave_idx_type j = 0; j < nc; j++)
323 for (octave_idx_type i = y.cidx (j); i < y.cidx (j+1); i++)
324 {
325 octave_quit ();
326 retval.data (y.ridx(i) + j * nr) = f (x, y.data (i));
327 }
328
329 retval.maybe_compress (true);
330 }
331 else
332 {
333 octave_idx_type nz = y.nnz ();
334 retval = SparseMatrix (nr, nc, nz);
335 octave_idx_type ii = 0;
336 retval.cidx (ii) = 0;
337
338 for (octave_idx_type j = 0; j < nc; j++)
339 {
340 for (octave_idx_type i = y.cidx (j); i < y.cidx (j+1); i++)
341 {
342 octave_quit ();
343 double val = f (x, y.data (i));
344
345 if (val != 0.0)
346 {
347 retval.data (ii) = val;
348 retval.ridx (ii++) = y.ridx (i);
349 }
350 }
351 retval.cidx (j + 1) = ii;
352 }
353
354 retval.maybe_compress (false);
355 }
356
357 return retval;
358}
359
360static SparseMatrix
361map_s_d (d_dd_fcn f, const SparseMatrix& x, double y)
362{
363 octave_idx_type nr = x.rows ();
364 octave_idx_type nc = x.columns ();
365 SparseMatrix retval;
366 double f_zero = f (0., y);
367
368 if (f_zero != 0)
369 {
370 retval = SparseMatrix (nr, nc, f_zero);
371
372 for (octave_idx_type j = 0; j < nc; j++)
373 for (octave_idx_type i = x.cidx (j); i < x.cidx (j+1); i++)
374 {
375 octave_quit ();
376 retval.data (x.ridx(i) + j * nr) = f (x.data (i), y);
377 }
378
379 retval.maybe_compress (true);
380 }
381 else
382 {
383 octave_idx_type nz = x.nnz ();
384 retval = SparseMatrix (nr, nc, nz);
385 octave_idx_type ii = 0;
386 retval.cidx (ii) = 0;
387
388 for (octave_idx_type j = 0; j < nc; j++)
389 {
390 for (octave_idx_type i = x.cidx (j); i < x.cidx (j+1); i++)
391 {
392 octave_quit ();
393 double val = f (x.data (i), y);
394
395 if (val != 0.0)
396 {
397 retval.data (ii) = val;
398 retval.ridx (ii++) = x.ridx (i);
399 }
400 }
401 retval.cidx (j + 1) = ii;
402 }
403
404 retval.maybe_compress (false);
405 }
406
407 return retval;
408}
409
410static SparseMatrix
411map_s_s (d_dd_fcn f, const SparseMatrix& x, const SparseMatrix& y)
412{
413 octave_idx_type nr = x.rows ();
414 octave_idx_type nc = x.columns ();
415
416 octave_idx_type y_nr = y.rows ();
417 octave_idx_type y_nc = y.columns ();
418
419 assert (nr == y_nr && nc == y_nc);
420
421 SparseMatrix retval;
422 double f_zero = f (0., 0.);
423
424 if (f_zero != 0)
425 {
426 retval = SparseMatrix (nr, nc, f_zero);
427 octave_idx_type k1 = 0, k2 = 0;
428
429 for (octave_idx_type j = 0; j < nc; j++)
430 {
431 while (k1 < x.cidx(j+1) && k2 < y.cidx(j+1))
432 {
433 octave_quit ();
434 if (k1 >= x.cidx(j+1))
435 {
436 retval.data (y.ridx(k2) + j * nr) = f (0.0, y.data (k2));
437 k2++;
438 }
439 else if (k2 >= y.cidx(j+1))
440 {
441 retval.data (x.ridx(k1) + j * nr) = f (x.data (k1), 0.0);
442 k1++;
443 }
444 else
445 {
446 octave_idx_type rx = x.ridx(k1);
447 octave_idx_type ry = y.ridx(k2);
448
449 if (rx < ry)
450 {
451 retval.data (rx + j * nr) = f (x.data (k1), 0.0);
452 k1++;
453 }
454 else if (rx > ry)
455 {
456 retval.data (ry + j * nr) = f (0.0, y.data (k2));
457 k2++;
458 }
459 else
460 {
461 retval.data (ry + j * nr) = f (x.data (k1), y.data (k2));
462 k1++;
463 k2++;
464 }
465 }
466 }
467 }
468
469 retval.maybe_compress (true);
470 }
471 else
472 {
473 octave_idx_type nz = x.nnz () + y.nnz ();
474 retval = SparseMatrix (nr, nc, nz);
475 octave_idx_type ii = 0;
476 retval.cidx (ii) = 0;
477 octave_idx_type k1 = 0, k2 = 0;
478
479 for (octave_idx_type j = 0; j < nc; j++)
480 {
481 while (k1 < x.cidx(j+1) && k2 < y.cidx(j+1))
482 {
483 octave_quit ();
484 double val;
485 octave_idx_type r;
486 if (k1 >= x.cidx(j+1))
487 {
488 r = y.ridx (k2);
489 val = f (0.0, y.data (k2++));
490 }
491 else if (k2 >= y.cidx(j+1))
492 {
493 r = x.ridx (k1);
494 val = f (x.data (k1++), 0.0);
495 }
496 else
497 {
498 octave_idx_type rx = x.ridx(k1);
499 octave_idx_type ry = y.ridx(k2);
500
501 if (rx < ry)
502 {
503 r = x.ridx (k1);
504 val = f (x.data (k1++), 0.0);
505 }
506 else if (rx > ry)
507 {
508 r = y.ridx (k2);
509 val = f (0.0, y.data (k2++));
510 }
511 else
512 {
513 r = y.ridx (k2);
514 val = f (x.data (k1++), y.data (k2++));
515 }
516 }
517 if (val != 0.0)
518 {
519 retval.data (ii) = val;
520 retval.ridx (ii++) = r;
521 }
522 }
523 retval.cidx (j + 1) = ii;
524 }
525
526 retval.maybe_compress (false);
527 }
528
529 return retval;
530}
531
532DEFUN (atan2, args, ,
533 "-*- texinfo -*-\n\
534@deftypefn {Mapping Function} {} atan2 (@var{y}, @var{x})\n\
535Compute atan (@var{y} / @var{x}) for corresponding elements of @var{y}\n\
536and @var{x}. Signal an error if @var{y} and @var{x} do not match in size\n\
537and orientation.\n\
538@end deftypefn")
539{
540 octave_value retval;
541
542 int nargin = args.length ();
543
544 if (nargin == 2 && args(0).is_defined () && args(1).is_defined ())
545 {
546 if (args(0).is_integer_type () || args(1).is_integer_type ())
547 error ("atan2: not defined for integer types");
548 else
549 {
550 octave_value arg_y = args(0);
551 octave_value arg_x = args(1);
552
553 dim_vector y_dims = arg_y.dims ();
554 dim_vector x_dims = arg_x.dims ();
555
556 bool y_is_scalar = y_dims.all_ones ();
557 bool x_is_scalar = x_dims.all_ones ();
558
559 bool is_float = arg_y.is_single_type () || arg_x.is_single_type ();
560
561 if (y_is_scalar && x_is_scalar)
562 {
563 if (is_float)
564 {
565 float y = arg_y.float_value ();
566
567 if (! error_state)
568 {
569 float x = arg_x.float_value ();
570
571 if (! error_state)
572 retval = atan2f (y, x);
573 }
574 }
575 else
576 {
577 double y = arg_y.double_value ();
578
579 if (! error_state)
580 {
581 double x = arg_x.double_value ();
582
583 if (! error_state)
584 retval = atan2 (y, x);
585 }
586 }
587 }
588 else if (y_is_scalar)
589 {
590 if (is_float)
591 {
592 float y = arg_y.float_value ();
593
594 if (! error_state)
595 {
596 // Even if x is sparse return a full matrix here
597 FloatNDArray x = arg_x.float_array_value ();
598
599 if (! error_state)
600 retval = map_f_fm (atan2f, y, x);
601 }
602 }
603 else
604 {
605 double y = arg_y.double_value ();
606
607 if (! error_state)
608 {
609 // Even if x is sparse return a full matrix here
610 NDArray x = arg_x.array_value ();
611
612 if (! error_state)
613 retval = map_d_m (atan2, y, x);
614 }
615 }
616 }
617 else if (x_is_scalar)
618 {
619 if (arg_y.is_sparse_type ())
620 {
621 SparseMatrix y = arg_y.sparse_matrix_value ();
622
623 if (! error_state)
624 {
625 double x = arg_x.double_value ();
626
627 if (! error_state)
628 retval = map_s_d (atan2, y, x);
629 }
630 }
631 else if (is_float)
632 {
633 FloatNDArray y = arg_y.float_array_value ();
634
635 if (! error_state)
636 {
637 float x = arg_x.float_value ();
638
639 if (! error_state)
640 retval = map_fm_f (atan2f, y, x);
641 }
642 }
643 else
644 {
645 NDArray y = arg_y.array_value ();
646
647 if (! error_state)
648 {
649 double x = arg_x.double_value ();
650
651 if (! error_state)
652 retval = map_m_d (atan2, y, x);
653 }
654 }
655 }
656 else if (y_dims == x_dims)
657 {
658 // Even if y is sparse return a full matrix here
659 if (arg_x.is_sparse_type ())
660 {
661 SparseMatrix y = arg_y.sparse_matrix_value ();
662
663 if (! error_state)
664 {
665 SparseMatrix x = arg_x.sparse_matrix_value ();
666
667 if (! error_state)
668 retval = map_s_s (atan2, y, x);
669 }
670 }
671 else if (is_float)
672 {
673 FloatNDArray y = arg_y.array_value ();
674
675 if (! error_state)
676 {
677 FloatNDArray x = arg_x.array_value ();
678
679 if (! error_state)
680 retval = map_fm_fm (atan2f, y, x);
681 }
682 }
683 else
684 {
685 NDArray y = arg_y.array_value ();
686
687 if (! error_state)
688 {
689 NDArray x = arg_x.array_value ();
690
691 if (! error_state)
692 retval = map_m_m (atan2, y, x);
693 }
694 }
695 }
696 else
697 error ("atan2: nonconformant matrices");
698 }
699 }
700 else
701 print_usage ();
702
703 return retval;
704}
705
706/*
707%!assert (size (atan2 (zeros (0, 2), zeros (0, 2))), [0, 2])
708%!assert (size (atan2 (rand (2, 3, 4), zeros (2, 3, 4))), [2, 3, 4])
709%!assert (size (atan2 (rand (2, 3, 4), 1)), [2, 3, 4])
710%!assert (size (atan2 (1, rand (2, 3, 4))), [2, 3, 4])
711%!assert (size (atan2 (1, 2)), [1, 1])
712
713%!test
714%! rt2 = sqrt (2);
715%! rt3 = sqrt (3);
716%! v = [0, pi/6, pi/4, pi/3, -pi/3, -pi/4, -pi/6, 0];
717%! y = [0, rt3, 1, rt3, -rt3, -1, -rt3, 0];
718%! x = [1, 3, 1, 1, 1, 1, 3, 1];
719%! assert(atan2 (y, x), v, sqrt (eps));
720
721%!test
722%! rt2 = sqrt (2);
723%! rt3 = sqrt (3);
724%! v = single([0, pi/6, pi/4, pi/3, -pi/3, -pi/4, -pi/6, 0]);
725%! y = single([0, rt3, 1, rt3, -rt3, -1, -rt3, 0]);
726%! x = single([1, 3, 1, 1, 1, 1, 3, 1]);
727%! assert(atan2 (y, x), v, sqrt (eps('single')));
728
729%!error <Invalid call to atan2.*> atan2 ();
730%!error <Invalid call to atan2.*> atan2 (1, 2, 3);
731
732*/
733
734
735
736DEFUN (hypot, args, ,
737 "-*- texinfo -*-\n\
738@deftypefn {Built-in Function} {} hypot (@var{x}, @var{y})\n\
739Compute the element-by-element square root of the sum of the squares of\n\
740@var{x} and @var{y}. This is equivalent to\n\
741@code{sqrt (@var{x}.^2 + @var{y}.^2)}, but calculated in a manner that\n\
742avoids overflows for large values of @var{x} or @var{y}.\n\
743@end deftypefn")
744{
745 octave_value retval;
746
747 int nargin = args.length ();
748
749 if (nargin == 2 && args(0).is_defined () && args(1).is_defined ())
750 {
751 if (args(0).is_integer_type () || args(1).is_integer_type ())
752 error ("hypot: not defined for integer types");
753 else
754 {
755 octave_value arg_x = args(0);
756 octave_value arg_y = args(1);
757
758 dim_vector x_dims = arg_x.dims ();
759 dim_vector y_dims = arg_y.dims ();
760
761 bool x_is_scalar = x_dims.all_ones ();
762 bool y_is_scalar = y_dims.all_ones ();
763
764 bool is_float = arg_y.is_single_type () || arg_x.is_single_type ();
765
766 if (y_is_scalar && x_is_scalar)
767 {
768 if (is_float)
769 {
770 float x;
771 if (arg_x.is_complex_type ())
772 x = abs (arg_x.float_complex_value ());
773 else
774 x = arg_x.float_value ();
775
776 if (! error_state)
777 {
778 float y;
779 if (arg_y.is_complex_type ())
780 y = abs (arg_y.float_complex_value ());
781 else
782 y = arg_y.float_value ();
783
784 if (! error_state)
785 retval = hypotf (x, y);
786 }
787 }
788 else
789 {
790 double x;
791 if (arg_x.is_complex_type ())
792 x = abs (arg_x.complex_value ());
793 else
794 x = arg_x.double_value ();
795
796 if (! error_state)
797 {
798 double y;
799 if (arg_y.is_complex_type ())
800 y = abs (arg_y.complex_value ());
801 else
802 y = arg_y.double_value ();
803
804 if (! error_state)
805 retval = hypot (x, y);
806 }
807 }
808 }
809 else if (y_is_scalar)
810 {
811 if (is_float)
812 {
813 FloatNDArray x;
814 if (arg_x.is_complex_type ())
815 x = arg_x.float_complex_array_value ().abs ();
816 else
817 x = arg_x.float_array_value ();
818
819 if (! error_state)
820 {
821 float y;
822 if (arg_y.is_complex_type ())
823 y = abs (arg_y.float_complex_value ());
824 else
825 y = arg_y.float_value ();
826
827 if (! error_state)
828 retval = map_fm_f (hypotf, x, y);
829 }
830 }
831 else
832 {
833 NDArray x;
834 if (arg_x.is_complex_type ())
835 x = arg_x.complex_array_value ().abs ();
836 else
837 x = arg_x.array_value ();
838
839 if (! error_state)
840 {
841 double y;
842 if (arg_y.is_complex_type ())
843 y = abs (arg_y.complex_value ());
844 else
845 y = arg_y.double_value ();
846
847 if (! error_state)
848 retval = map_m_d (hypot, x, y);
849 }
850 }
851 }
852 else if (x_is_scalar)
853 {
854 if (is_float)
855 {
856 float x;
857 if (arg_x.is_complex_type ())
858 x = abs (arg_x.float_complex_value ());
859 else
860 x = arg_x.float_value ();
861
862 if (! error_state)
863 {
864 FloatNDArray y;
865 if (arg_y.is_complex_type ())
866 y = arg_y.float_complex_array_value ().abs ();
867 else
868 y = arg_y.float_array_value ();
869
870 if (! error_state)
871 retval = map_f_fm (hypotf, x, y);
872 }
873 }
874 else
875 {
876 double x;
877 if (arg_x.is_complex_type ())
878 x = abs (arg_x.complex_value ());
879 else
880 x = arg_x.double_value ();
881
882 if (! error_state)
883 {
884 NDArray y;
885 if (arg_y.is_complex_type ())
886 y = arg_y.complex_array_value ().abs ();
887 else
888 y = arg_y.array_value ();
889
890 if (! error_state)
891 retval = map_d_m (hypot, x, y);
892 }
893 }
894 }
895 else if (y_dims == x_dims)
896 {
897 if (arg_x.is_sparse_type () && arg_y.is_sparse_type ())
898 {
899 SparseMatrix x;
900 if (arg_x.is_complex_type ())
901 x = arg_x.sparse_complex_matrix_value ().abs ();
902 else
903 x = arg_x.sparse_matrix_value ();
904
905 if (! error_state)
906 {
907 SparseMatrix y;
908 if (arg_y.is_complex_type ())
909 y = arg_y.sparse_complex_matrix_value ().abs ();
910 else
911 y = arg_y.sparse_matrix_value ();
912
913 if (! error_state)
914 retval = map_s_s (hypot, x, y);
915 }
916 }
917 else if (is_float)
918 {
919 FloatNDArray x;
920 if (arg_x.is_complex_type ())
921 x = arg_x.float_complex_array_value ().abs ();
922 else
923 x = arg_x.float_array_value ();
924
925 if (! error_state)
926 {
927 FloatNDArray y;
928 if (arg_y.is_complex_type ())
929 y = arg_y.float_complex_array_value ().abs ();
930 else
931 y = arg_y.float_array_value ();
932
933 if (! error_state)
934 retval = map_fm_fm (hypotf, x, y);
935 }
936 }
937 else
938 {
939 NDArray x;
940 if (arg_x.is_complex_type ())
941 x = arg_x.complex_array_value ().abs ();
942 else
943 x = arg_x.array_value ();
944
945 if (! error_state)
946 {
947 NDArray y;
948 if (arg_y.is_complex_type ())
949 y = arg_y.complex_array_value ().abs ();
950 else
951 y = arg_y.array_value ();
952
953 if (! error_state)
954 retval = map_m_m (hypot, x, y);
955 }
956 }
957 }
958 else
959 error ("hypot: nonconformant matrices");
960 }
961 }
962 else
963 print_usage ();
964
965 return retval;
966}
967
968/*
969%!assert (size (hypot (zeros (0, 2), zeros (0, 2))), [0, 2])
970%!assert (size (hypot (rand (2, 3, 4), zeros (2, 3, 4))), [2, 3, 4])
971%!assert (size (hypot (rand (2, 3, 4), 1)), [2, 3, 4])
972%!assert (size (hypot (1, rand (2, 3, 4))), [2, 3, 4])
973%!assert (size (hypot (1, 2)), [1, 1])
974%!assert (hypot (1:10, 1:10), sqrt(2) * [1:10], 16*eps)
975%!assert (hypot (single(1:10), single(1:10)), single(sqrt(2) * [1:10]));
976*/
977
978template<typename T, typename ET>
979void
980map_2_xlog2 (const Array<T>& x, Array<T>& f, Array<ET>& e)
981{
982 f = Array<T>(x.dims ());
983 e = Array<ET>(x.dims ());
984 for (octave_idx_type i = 0; i < x.numel (); i++)
985 {
986 int exp;
987 f.xelem (i) = xlog2 (x(i), exp);
988 e.xelem (i) = exp;
989 }
990}
991
992DEFUN (log2, args, nargout,
993 "-*- texinfo -*-\n\
994@deftypefn {Mapping Function} {} log2 (@var{x})\n\
995@deftypefnx {Mapping Function} {[@var{f}, @var{e}] =} log2 (@var{x})\n\
996Compute the base-2 logarithm of each element of @var{x}.\n\
997\n\
998If called with two output arguments, split @var{x} into\n\
999binary mantissa and exponent so that\n\
1000@tex\n\
1001${1 \\over 2} \\le \\left| f \\right| < 1$\n\
1002@end tex\n\
1003@ifnottex\n\
1004@code{1/2 <= abs(f) < 1}\n\
1005@end ifnottex\n\
1006and @var{e} is an integer. If\n\
1007@tex\n\
1008$x = 0$, $f = e = 0$.\n\
1009@end tex\n\
1010@ifnottex\n\
1011@code{x = 0}, @code{f = e = 0}.\n\
1012@end ifnottex\n\
1013@seealso{pow2, log, log10, exp}\n\
1014@end deftypefn")
1015{
1016 octave_value_list retval;
1017
1018 if (args.length () == 1)
1019 {
1020 if (nargout < 2)
1021 retval(0) = args(0).log2 ();
1022 else if (args(0).is_single_type ())
1023 {
1024 if (args(0).is_real_type ())
1025 {
1026 FloatNDArray f;
1027 FloatNDArray x = args(0).float_array_value ();
1028 // FIXME -- should E be an int value?
1029 FloatMatrix e;
1030 map_2_xlog2 (x, f, e);
1031 retval (1) = e;
1032 retval (0) = f;
1033 }
1034 else if (args(0).is_complex_type ())
1035 {
1036 FloatComplexNDArray f;
1037 FloatComplexNDArray x = args(0).float_complex_array_value ();
1038 // FIXME -- should E be an int value?
1039 FloatNDArray e;
1040 map_2_xlog2 (x, f, e);
1041 retval (1) = e;
1042 retval (0) = f;
1043 }
1044 }
1045 else if (args(0).is_real_type ())
1046 {
1047 NDArray f;
1048 NDArray x = args(0).array_value ();
1049 // FIXME -- should E be an int value?
1050 Matrix e;
1051 map_2_xlog2 (x, f, e);
1052 retval (1) = e;
1053 retval (0) = f;
1054 }
1055 else if (args(0).is_complex_type ())
1056 {
1057 ComplexNDArray f;
1058 ComplexNDArray x = args(0).complex_array_value ();
1059 // FIXME -- should E be an int value?
1060 NDArray e;
1061 map_2_xlog2 (x, f, e);
1062 retval (1) = e;
1063 retval (0) = f;
1064 }
1065 else
1066 gripe_wrong_type_arg ("log2", args(0));
1067 }
1068 else
1069 print_usage ();
1070
1071 return retval;
1072}
1073
1074/*
1075%!assert(log2 ([1/4, 1/2, 1, 2, 4]), [-2, -1, 0, 1, 2]);
1076%!assert(log2(Inf), Inf);
1077%!assert(isnan(log2(NaN)));
1078%!assert(log2(4*i), 2 + log2(1*i));
1079%!assert(log2(complex(0,Inf)), Inf + log2(i));
1080
1081%!test
1082%! [f, e] = log2 ([0,-1; 2,-4; Inf,-Inf]);
1083%! assert (f, [0,-0.5; 0.5,-0.5; Inf,-Inf]);
1084%! assert (e(1:2,:), [0,1;2,3])
1085
1086%!test
1087%! [f, e] = log2 (complex (zeros (3, 2), [0,-1; 2,-4; Inf,-Inf]));
1088%! assert (f, complex (zeros (3, 2), [0,-0.5; 0.5,-0.5; Inf,-Inf]));
1089%! assert (e(1:2,:), [0,1; 2,3]);
1090*/
1091
1092DEFUN (fmod, args, ,
1093 "-*- texinfo -*-\n\
1094@deftypefn {Mapping Function} {} fmod (@var{x}, @var{y})\n\
1095Compute the floating point remainder of dividing @var{x} by @var{y}\n\
1096using the C library function @code{fmod}. The result has the same\n\
1097sign as @var{x}. If @var{y} is zero, the result is implementation-dependent.\n\
1098@seealso{mod, rem}\n\
1099@end deftypefn")
1100{
1101 octave_value retval;
1102
1103 int nargin = args.length ();
1104
1105 if (nargin == 2 && args(0).is_defined () && args(1).is_defined ())
1106 {
1107 octave_value arg_x = args(0);
1108 octave_value arg_y = args(1);
1109
1110 dim_vector y_dims = arg_y.dims ();
1111 dim_vector x_dims = arg_x.dims ();
1112
1113 bool y_is_scalar = y_dims.all_ones ();
1114 bool x_is_scalar = x_dims.all_ones ();
1115
1116 bool is_float = arg_y.is_single_type () || arg_x.is_single_type ();
1117
1118 if (y_is_scalar && x_is_scalar)
1119 {
1120 if (is_float)
1121 {
1122 float y = arg_y.float_value ();
1123
1124 if (! error_state)
1125 {
1126 float x = arg_x.float_value ();
1127
1128 if (! error_state)
1129 retval = fmod (x, y);
1130 }
1131 }
1132 else
1133 {
1134 double y = arg_y.double_value ();
1135
1136 if (! error_state)
1137 {
1138 double x = arg_x.double_value ();
1139
1140 if (! error_state)
1141 retval = fmod (x, y);
1142 }
1143 }
1144 }
1145 else if (y_is_scalar)
1146 {
1147 if (is_float)
1148 {
1149 float y = arg_y.float_value ();
1150
1151 if (! error_state)
1152 {
1153 FloatNDArray x = arg_x.float_array_value ();
1154
1155 if (! error_state)
1156 retval = map_fm_f (fmodf, x, y);
1157 }
1158 }
1159 else
1160 {
1161 double y = arg_y.double_value ();
1162
1163 if (! error_state)
1164 {
1165 if (arg_x.is_sparse_type ())
1166 {
1167 SparseMatrix x = arg_x.sparse_matrix_value ();
1168
1169 if (! error_state)
1170 retval = map_s_d (fmod, x, y);
1171 }
1172 else
1173 {
1174 NDArray x = arg_x.array_value ();
1175
1176 if (! error_state)
1177 retval = map_m_d (fmod, x, y);
1178 }
1179 }
1180 }
1181 }
1182 else if (x_is_scalar)
1183 {
1184 if (arg_y.is_sparse_type ())
1185 {
1186 SparseMatrix y = arg_y.sparse_matrix_value ();
1187
1188 if (! error_state)
1189 {
1190 double x = arg_x.double_value ();
1191
1192 if (! error_state)
1193 retval = map_d_s (fmod, x, y);
1194 }
1195 }
1196 else if (is_float)
1197 {
1198 FloatNDArray y = arg_y.float_array_value ();
1199
1200 if (! error_state)
1201 {
1202 float x = arg_x.float_value ();
1203
1204 if (! error_state)
1205 retval = map_f_fm (fmodf, x, y);
1206 }
1207 }
1208 else
1209 {
1210 NDArray y = arg_y.array_value ();
1211
1212 if (! error_state)
1213 {
1214 double x = arg_x.double_value ();
1215
1216 if (! error_state)
1217 retval = map_d_m (fmod, x, y);
1218 }
1219 }
1220 }
1221 else if (y_dims == x_dims)
1222 {
1223 if (arg_y.is_sparse_type () || arg_x.is_sparse_type ())
1224 {
1225 SparseMatrix y = arg_y.sparse_matrix_value ();
1226
1227 if (! error_state)
1228 {
1229 SparseMatrix x = arg_x.sparse_matrix_value ();
1230
1231 if (! error_state)
1232 retval = map_s_s (fmod, x, y);
1233 }
1234 }
1235 else if (is_float)
1236 {
1237 FloatNDArray y = arg_y.float_array_value ();
1238
1239 if (! error_state)
1240 {
1241 FloatNDArray x = arg_x.float_array_value ();
1242
1243 if (! error_state)
1244 retval = map_fm_fm (fmodf, x, y);
1245 }
1246 }
1247 else
1248 {
1249 NDArray y = arg_y.array_value ();
1250
1251 if (! error_state)
1252 {
1253 NDArray x = arg_x.array_value ();
1254
1255 if (! error_state)
1256 retval = map_m_m (fmod, x, y);
1257 }
1258 }
1259 }
1260 else
1261 error ("fmod: nonconformant matrices");
1262 }
1263 else
1264 print_usage ();
1265
1266 return retval;
1267}
1268
1269/*
1270%!assert (size (fmod (zeros (0, 2), zeros (0, 2))), [0, 2])
1271%!assert (size (fmod (rand (2, 3, 4), zeros (2, 3, 4))), [2, 3, 4])
1272%!assert (size (fmod (rand (2, 3, 4), 1)), [2, 3, 4])
1273%!assert (size (fmod (1, rand (2, 3, 4))), [2, 3, 4])
1274%!assert (size (fmod (1, 2)), [1, 1])
1275*/
1276
1277// FIXME Need to convert the reduction functions of this file for single precision
1278
1279#define NATIVE_REDUCTION_1(FCN, TYPE, DIM) \
1280 (arg.is_ ## TYPE ## _type ()) \
1281 { \
1282 TYPE ## NDArray tmp = arg. TYPE ##_array_value (); \
1283 \
1284 if (! error_state) \
1285 { \
1286 octave_ ## TYPE::clear_conv_flags (); \
1287 retval = tmp.FCN (DIM); \
1288 if (octave_ ## TYPE::get_trunc_flag ()) \
1289 { \
1290 gripe_native_integer_math_truncated (#FCN, \
1291 octave_ ## TYPE::type_name ()); \
1292 octave_ ## TYPE::clear_conv_flags (); \
1293 } \
1294 } \
1295 }
1296
1297#define NATIVE_REDUCTION(FCN, BOOL_FCN) \
1298 \
1299 octave_value retval; \
1300 \
1301 int nargin = args.length (); \
1302 \
1303 bool isnative = false; \
1304 bool isdouble = false; \
1305 \
1306 if (nargin > 1 && args(nargin - 1).is_string ()) \
1307 { \
1308 std::string str = args(nargin - 1).string_value (); \
1309 \
1310 if (! error_state) \
1311 { \
1312 if (str == "native") \
1313 isnative = true; \
1314 else if (str == "double") \
1315 isdouble = true; \
1316 else \
1317 error ("sum: unrecognized string argument"); \
1318 nargin --; \
1319 } \
1320 } \
1321 \
1322 if (nargin == 1 || nargin == 2) \
1323 { \
1324 octave_value arg = args(0); \
1325 \
1326 int dim = (nargin == 1 ? -1 : args(1).int_value (true) - 1); \
1327 \
1328 if (! error_state) \
1329 { \
1330 if (dim >= -1) \
1331 { \
1332 if (arg.is_sparse_type ()) \
1333 { \
1334 if (arg.is_real_type ()) \
1335 { \
1336 SparseMatrix tmp = arg.sparse_matrix_value (); \
1337 \
1338 if (! error_state) \
1339 retval = tmp.FCN (dim); \
1340 } \
1341 else \
1342 { \
1343 SparseComplexMatrix tmp = arg.sparse_complex_matrix_value (); \
1344 \
1345 if (! error_state) \
1346 retval = tmp.FCN (dim); \
1347 } \
1348 } \
1349 else \
1350 { \
1351 if (isnative) \
1352 { \
1353 if NATIVE_REDUCTION_1 (FCN, uint8, dim) \
1354 else if NATIVE_REDUCTION_1 (FCN, uint16, dim) \
1355 else if NATIVE_REDUCTION_1 (FCN, uint32, dim) \
1356 else if NATIVE_REDUCTION_1 (FCN, uint64, dim) \
1357 else if NATIVE_REDUCTION_1 (FCN, int8, dim) \
1358 else if NATIVE_REDUCTION_1 (FCN, int16, dim) \
1359 else if NATIVE_REDUCTION_1 (FCN, int32, dim) \
1360 else if NATIVE_REDUCTION_1 (FCN, int64, dim) \
1361 else if (arg.is_bool_type ()) \
1362 { \
1363 boolNDArray tmp = arg.bool_array_value (); \
1364 if (! error_state) \
1365 retval = boolNDArray (tmp.BOOL_FCN (dim)); \
1366 } \
1367 else if (arg.is_char_matrix ()) \
1368 { \
1369 error (#FCN, ": invalid char type"); \
1370 } \
1371 else if (!isdouble && arg.is_single_type ()) \
1372 { \
1373 if (arg.is_complex_type ()) \
1374 { \
1375 FloatComplexNDArray tmp = \
1376 arg.float_complex_array_value (); \
1377 \
1378 if (! error_state) \
1379 retval = tmp.FCN (dim); \
1380 } \
1381 else if (arg.is_real_type ()) \
1382 { \
1383 FloatNDArray tmp = arg.float_array_value (); \
1384 \
1385 if (! error_state) \
1386 retval = tmp.FCN (dim); \
1387 } \
1388 } \
1389 else if (arg.is_complex_type ()) \
1390 { \
1391 ComplexNDArray tmp = arg.complex_array_value (); \
1392 \
1393 if (! error_state) \
1394 retval = tmp.FCN (dim); \
1395 } \
1396 else if (arg.is_real_type ()) \
1397 { \
1398 NDArray tmp = arg.array_value (); \
1399 \
1400 if (! error_state) \
1401 retval = tmp.FCN (dim); \
1402 } \
1403 else \
1404 { \
1405 gripe_wrong_type_arg (#FCN, arg); \
1406 return retval; \
1407 } \
1408 } \
1409 else if (arg.is_bool_type ()) \
1410 { \
1411 boolNDArray tmp = arg.bool_array_value (); \
1412 if (! error_state) \
1413 retval = tmp.FCN (dim); \
1414 } \
1415 else if (!isdouble && arg.is_single_type ()) \
1416 { \
1417 if (arg.is_real_type ()) \
1418 { \
1419 FloatNDArray tmp = arg.float_array_value (); \
1420 \
1421 if (! error_state) \
1422 retval = tmp.FCN (dim); \
1423 } \
1424 else if (arg.is_complex_type ()) \
1425 { \
1426 FloatComplexNDArray tmp = \
1427 arg.float_complex_array_value (); \
1428 \
1429 if (! error_state) \
1430 retval = tmp.FCN (dim); \
1431 } \
1432 } \
1433 else if (arg.is_real_type ()) \
1434 { \
1435 NDArray tmp = arg.array_value (); \
1436 \
1437 if (! error_state) \
1438 retval = tmp.FCN (dim); \
1439 } \
1440 else if (arg.is_complex_type ()) \
1441 { \
1442 ComplexNDArray tmp = arg.complex_array_value (); \
1443 \
1444 if (! error_state) \
1445 retval = tmp.FCN (dim); \
1446 } \
1447 else \
1448 { \
1449 gripe_wrong_type_arg (#FCN, arg); \
1450 return retval; \
1451 } \
1452 } \
1453 } \
1454 else \
1455 error (#FCN ": invalid dimension argument = %d", dim + 1); \
1456 } \
1457 \
1458 } \
1459 else \
1460 print_usage (); \
1461 \
1462 return retval
1463
1464#define DATA_REDUCTION(FCN) \
1465 \
1466 octave_value retval; \
1467 \
1468 int nargin = args.length (); \
1469 \
1470 if (nargin == 1 || nargin == 2) \
1471 { \
1472 octave_value arg = args(0); \
1473 \
1474 int dim = (nargin == 1 ? -1 : args(1).int_value (true) - 1); \
1475 \
1476 if (! error_state) \
1477 { \
1478 if (dim >= -1) \
1479 { \
1480 if (arg.is_real_type ()) \
1481 { \
1482 if (arg.is_sparse_type ()) \
1483 { \
1484 SparseMatrix tmp = arg.sparse_matrix_value (); \
1485 \
1486 if (! error_state) \
1487 retval = tmp.FCN (dim); \
1488 } \
1489 else if (arg.is_single_type ()) \
1490 { \
1491 FloatNDArray tmp = arg.float_array_value (); \
1492 \
1493 if (! error_state) \
1494 retval = tmp.FCN (dim); \
1495 } \
1496 else \
1497 { \
1498 NDArray tmp = arg.array_value (); \
1499 \
1500 if (! error_state) \
1501 retval = tmp.FCN (dim); \
1502 } \
1503 } \
1504 else if (arg.is_complex_type ()) \
1505 { \
1506 if (arg.is_sparse_type ()) \
1507 { \
1508 SparseComplexMatrix tmp = arg.sparse_complex_matrix_value (); \
1509 \
1510 if (! error_state) \
1511 retval = tmp.FCN (dim); \
1512 } \
1513 else if (arg.is_single_type ()) \
1514 { \
1515 FloatComplexNDArray tmp = arg.float_complex_array_value (); \
1516 \
1517 if (! error_state) \
1518 retval = tmp.FCN (dim); \
1519 } \
1520 else \
1521 { \
1522 ComplexNDArray tmp = arg.complex_array_value (); \
1523 \
1524 if (! error_state) \
1525 retval = tmp.FCN (dim); \
1526 } \
1527 } \
1528 else \
1529 { \
1530 gripe_wrong_type_arg (#FCN, arg); \
1531 return retval; \
1532 } \
1533 } \
1534 else \
1535 error (#FCN ": invalid dimension argument = %d", dim + 1); \
1536 } \
1537 } \
1538 else \
1539 print_usage (); \
1540 \
1541 return retval
1542
1543DEFUN (cumprod, args, ,
1544 "-*- texinfo -*-\n\
1545@deftypefn {Built-in Function} {} cumprod (@var{x})\n\
1546@deftypefnx {Built-in Function} {} cumprod (@var{x}, @var{dim})\n\
1547Cumulative product of elements along dimension @var{dim}. If\n\
1548@var{dim} is omitted, it defaults to the first non-singleton dimension.\n\
1549\n\
1550@seealso{prod, cumsum}\n\
1551@end deftypefn")
1552{
1553 DATA_REDUCTION (cumprod);
1554}
1555
1556/*
1557
1558%!assert (cumprod ([1, 2, 3]), [1, 2, 6]);
1559%!assert (cumprod ([-1; -2; -3]), [-1; 2; -6]);
1560%!assert (cumprod ([i, 2+i, -3+2i, 4]), [i, -1+2i, -1-8i, -4-32i]);
1561%!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]);
1562
1563%!assert (cumprod (single([1, 2, 3])), single([1, 2, 6]));
1564%!assert (cumprod (single([-1; -2; -3])), single([-1; 2; -6]));
1565%!assert (cumprod (single([i, 2+i, -3+2i, 4])), single([i, -1+2i, -1-8i, -4-32i]));
1566%!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]));
1567
1568%!error <Invalid call to cumprod.*> cumprod ();
1569
1570%!assert (cumprod ([2, 3; 4, 5], 1), [2, 3; 8, 15]);
1571%!assert (cumprod ([2, 3; 4, 5], 2), [2, 6; 4, 20]);
1572
1573%!assert (cumprod (single([2, 3; 4, 5]), 1), single([2, 3; 8, 15]));
1574%!assert (cumprod (single([2, 3; 4, 5]), 2), single([2, 6; 4, 20]));
1575
1576 */
1577
1578DEFUN (cumsum, args, ,
1579 "-*- texinfo -*-\n\
1580@deftypefn {Built-in Function} {} cumsum (@var{x})\n\
1581@deftypefnx {Built-in Function} {} cumsum (@var{x}, @var{dim})\n\
1582@deftypefnx {Built-in Function} {} cumsum (@dots{}, 'native')\n\
1583Cumulative sum of elements along dimension @var{dim}. If @var{dim}\n\
1584is omitted, it defaults to the first non-singleton dimension.\n\
1585\n\
1586The \"native\" argument implies the summation is performed in native type.\n\
1587 See @code{sum} for a complete description and example of the use of\n\
1588\"native\".\n\
1589@seealso{sum, cumprod}\n\
1590@end deftypefn")
1591{
1592 octave_value retval;
1593
1594 int nargin = args.length ();
1595
1596 bool isnative = false;
1597 bool isdouble = false;
1598
1599 if (nargin > 1 && args(nargin - 1).is_string ())
1600 {
1601 std::string str = args(nargin - 1).string_value ();
1602
1603 if (! error_state)
1604 {
1605 if (str == "native")
1606 isnative = true;
1607 else if (str == "double")
1608 isdouble = true;
1609 else
1610 error ("sum: unrecognized string argument");
1611 nargin --;
1612 }
1613 }
1614
1615 if (error_state)
1616 return retval;
1617
1618 if (nargin == 1 || nargin == 2)
1619 {
1620 octave_value arg = args(0);
1621
1622 int dim = -1;
1623 if (nargin == 2)
1624 {
1625 dim = args(1).int_value () - 1;
1626 if (dim < 0)
1627 error ("cumsum: invalid dimension argument = %d", dim + 1);
1628 }
1629
1630 if (! error_state)
1631 {
1632 switch (arg.builtin_type ())
1633 {
1634 case btyp_double:
1635 if (arg.is_sparse_type ())
1636 retval = arg.sparse_matrix_value ().cumsum (dim);
1637 else
1638 retval = arg.array_value ().cumsum (dim);
1639 break;
1640 case btyp_complex:
1641 if (arg.is_sparse_type ())
1642 retval = arg.sparse_complex_matrix_value ().cumsum (dim);
1643 else
1644 retval = arg.complex_array_value ().cumsum (dim);
1645 break;
1646 case btyp_float:
1647 if (isdouble)
1648 retval = arg.array_value ().cumsum (dim);
1649 else
1650 retval = arg.float_array_value ().cumsum (dim);
1651 break;
1652 case btyp_float_complex:
1653 if (isdouble)
1654 retval = arg.complex_array_value ().cumsum (dim);
1655 else
1656 retval = arg.float_complex_array_value ().cumsum (dim);
1657 break;
1658
1659#define MAKE_INT_BRANCH(X) \
1660 case btyp_ ## X: \
1661 if (isnative) \
1662 retval = arg.X ## _array_value ().cumsum (dim); \
1663 else \
1664 retval = arg.array_value ().cumsum (dim); \
1665 break
1666 MAKE_INT_BRANCH (int8);
1667 MAKE_INT_BRANCH (int16);
1668 MAKE_INT_BRANCH (int32);
1669 MAKE_INT_BRANCH (int64);
1670 MAKE_INT_BRANCH (uint8);
1671 MAKE_INT_BRANCH (uint16);
1672 MAKE_INT_BRANCH (uint32);
1673 MAKE_INT_BRANCH (uint64);
1674#undef MAKE_INT_BRANCH
1675
1676 case btyp_bool:
1677 if (arg.is_sparse_type ())
1678 {
1679 SparseMatrix cs = arg.sparse_matrix_value ().cumsum (dim);
1680 if (isnative)
1681 retval = cs != 0.0;
1682 else
1683 retval = cs;
1684 }
1685 else
1686 {
1687 NDArray cs = arg.bool_array_value ().cumsum (dim);
1688 if (isnative)
1689 retval = cs != 0.0;
1690 else
1691 retval = cs;
1692 }
1693 break;
1694
1695 default:
1696 gripe_wrong_type_arg ("cumsum", arg);
1697 }
1698 }
1699 }
1700 else
1701 print_usage ();
1702
1703 return retval;
1704}
1705
1706/*
1707
1708%!assert (cumsum ([1, 2, 3]), [1, 3, 6]);
1709%!assert (cumsum ([-1; -2; -3]), [-1; -3; -6]);
1710%!assert (cumsum ([i, 2+i, -3+2i, 4]), [i, 2+2i, -1+4i, 3+4i]);
1711%!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]);
1712
1713%!assert (cumsum (single([1, 2, 3])), single([1, 3, 6]));
1714%!assert (cumsum (single([-1; -2; -3])), single([-1; -3; -6]));
1715%!assert (cumsum (single([i, 2+i, -3+2i, 4])), single([i, 2+2i, -1+4i, 3+4i]));
1716%!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]));
1717
1718%!error <Invalid call to cumsum.*> cumsum ();
1719
1720%!assert (cumsum ([1, 2; 3, 4], 1), [1, 2; 4, 6]);
1721%!assert (cumsum ([1, 2; 3, 4], 2), [1, 3; 3, 7]);
1722
1723%!assert (cumsum (single([1, 2; 3, 4]), 1), single([1, 2; 4, 6]));
1724%!assert (cumsum (single([1, 2; 3, 4]), 2), single([1, 3; 3, 7]));
1725
1726 */
1727
1728DEFUN (diag, args, ,
1729 "-*- texinfo -*-\n\
1730@deftypefn {Built-in Function} {} diag (@var{v}, @var{k})\n\
1731Return a diagonal matrix with vector @var{v} on diagonal @var{k}. The\n\
1732second argument is optional. If it is positive, the vector is placed on\n\
1733the @var{k}-th super-diagonal. If it is negative, it is placed on the\n\
1734@var{-k}-th sub-diagonal. The default value of @var{k} is 0, and the\n\
1735vector is placed on the main diagonal. For example,\n\
1736\n\
1737@example\n\
1738@group\n\
1739diag ([1, 2, 3], 1)\n\
1740 @result{} 0 1 0 0\n\
1741 0 0 2 0\n\
1742 0 0 0 3\n\
1743 0 0 0 0\n\
1744@end group\n\
1745@end example\n\
1746\n\
1747@noindent\n\
1748Given a matrix argument, instead of a vector, @code{diag} extracts the\n\
1749@var{k}-th diagonal of the matrix.\n\
1750@end deftypefn")
1751{
1752 octave_value retval;
1753
1754 int nargin = args.length ();
1755
1756 if (nargin == 1 && args(0).is_defined ())
1757 retval = args(0).diag();
1758 else if (nargin == 2 && args(0).is_defined () && args(1).is_defined ())
1759 {
1760 octave_idx_type k = args(1).int_value ();
1761
1762 if (error_state)
1763 error ("diag: invalid second argument");
1764 else
1765 retval = args(0).diag(k);
1766 }
1767 else if (nargin == 3)
1768 {
1769 octave_value arg0 = args(0);
1770 if (arg0.ndims () == 2 && (args(0).rows () == 1 || args(0).columns () == 1))
1771 {
1772 octave_idx_type m = args(1).int_value (), n = args(2).int_value ();
1773 if (! error_state)
1774 retval = arg0.diag ().resize (dim_vector (m, n));
1775 else
1776 error ("diag: invalid dimensions");
1777 }
1778 else
1779 error ("diag: first argument must be a vector");
1780 }
1781 else
1782 print_usage ();
1783
1784 return retval;
1785}
1786
1787/*
1788
1789%!assert(full (diag ([1; 2; 3])), [1, 0, 0; 0, 2, 0; 0, 0, 3]);
1790%!assert(diag ([1; 2; 3], 1), [0, 1, 0, 0; 0, 0, 2, 0; 0, 0, 0, 3; 0, 0, 0, 0]);
1791%!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]);
1792%!assert(diag ([1; 2; 3],-1), [0, 0, 0, 0; 1, 0, 0, 0; 0, 2, 0, 0; 0, 0, 3, 0]);
1793%!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]);
1794
1795%!assert(diag ([1, 0, 0; 0, 2, 0; 0, 0, 3]), [1; 2; 3]);
1796%!assert(diag ([0, 1, 0, 0; 0, 0, 2, 0; 0, 0, 0, 3; 0, 0, 0, 0], 1), [1; 2; 3]);
1797%!assert(diag ([0, 0, 0, 0; 1, 0, 0, 0; 0, 2, 0, 0; 0, 0, 3, 0], -1), [1; 2; 3]);
1798
1799%!assert(full (diag (single([1; 2; 3]))), single([1, 0, 0; 0, 2, 0; 0, 0, 3]));
1800%!assert(diag (single([1; 2; 3]), 1), single([0, 1, 0, 0; 0, 0, 2, 0; 0, 0, 0, 3; 0, 0, 0, 0]));
1801%!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]));
1802%!assert(diag (single([1; 2; 3]),-1), single([0, 0, 0, 0; 1, 0, 0, 0; 0, 2, 0, 0; 0, 0, 3, 0]));
1803%!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]));
1804
1805%!assert(diag (single([1, 0, 0; 0, 2, 0; 0, 0, 3])), single([1; 2; 3]));
1806%!assert(diag (single([0, 1, 0, 0; 0, 0, 2, 0; 0, 0, 0, 3; 0, 0, 0, 0]), 1), single([1; 2; 3]));
1807%!assert(diag (single([0, 0, 0, 0; 1, 0, 0, 0; 0, 2, 0, 0; 0, 0, 3, 0]), -1), single([1; 2; 3]));
1808
1809%!assert(diag (int8([1; 2; 3])), int8([1, 0, 0; 0, 2, 0; 0, 0, 3]));
1810%!assert(diag (int8([1; 2; 3]), 1), int8([0, 1, 0, 0; 0, 0, 2, 0; 0, 0, 0, 3; 0, 0, 0, 0]));
1811%!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]));
1812%!assert(diag (int8([1; 2; 3]),-1), int8([0, 0, 0, 0; 1, 0, 0, 0; 0, 2, 0, 0; 0, 0, 3, 0]));
1813%!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]));
1814
1815%!assert(diag (int8([1, 0, 0; 0, 2, 0; 0, 0, 3])), int8([1; 2; 3]));
1816%!assert(diag (int8([0, 1, 0, 0; 0, 0, 2, 0; 0, 0, 0, 3; 0, 0, 0, 0]), 1), int8([1; 2; 3]));
1817%!assert(diag (int8([0, 0, 0, 0; 1, 0, 0, 0; 0, 2, 0, 0; 0, 0, 3, 0]), -1), int8([1; 2; 3]));
1818
1819%!error <Invalid call to diag.*> diag ();
1820
1821 */
1822
1823DEFUN (prod, args, ,
1824 "-*- texinfo -*-\n\
1825@deftypefn {Built-in Function} {} prod (@var{x})\n\
1826@deftypefnx {Built-in Function} {} prod (@var{x}, @var{dim})\n\
1827Product of elements along dimension @var{dim}. If @var{dim} is\n\
1828omitted, it defaults to the first non-singleton dimension.\n\
1829@seealso{cumprod, sum}\n\
1830@end deftypefn")
1831{
1832 DATA_REDUCTION (prod);
1833}
1834
1835/*
1836
1837%!assert (prod ([1, 2, 3]), 6);
1838%!assert (prod ([-1; -2; -3]), -6);
1839%!assert (prod ([i, 2+i, -3+2i, 4]), -4 - 32i);
1840%!assert (prod ([1, 2, 3; i, 2i, 3i; 1+i, 2+2i, 3+3i]), [-1+i, -8+8i, -27+27i]);
1841
1842%!assert (prod (single([1, 2, 3])), single(6));
1843%!assert (prod (single([-1; -2; -3])), single(-6));
1844%!assert (prod (single([i, 2+i, -3+2i, 4])), single(-4 - 32i));
1845%!assert (prod (single([1, 2, 3; i, 2i, 3i; 1+i, 2+2i, 3+3i])), single([-1+i, -8+8i, -27+27i]));
1846
1847%!error <Invalid call to prod.*> prod ();
1848
1849%!assert (prod ([1, 2; 3, 4], 1), [3, 8]);
1850%!assert (prod ([1, 2; 3, 4], 2), [2; 12]);
1851%!assert (prod (zeros (1, 0)), 1);
1852%!assert (prod (zeros (1, 0), 1), zeros (1, 0));
1853%!assert (prod (zeros (1, 0), 2), 1);
1854%!assert (prod (zeros (0, 1)), 1);
1855%!assert (prod (zeros (0, 1), 1), 1);
1856%!assert (prod (zeros (0, 1), 2), zeros (0, 1));
1857%!assert (prod (zeros (2, 0)), zeros (1, 0));
1858%!assert (prod (zeros (2, 0), 1), zeros (1, 0));
1859%!assert (prod (zeros (2, 0), 2), [1; 1]);
1860%!assert (prod (zeros (0, 2)), [1, 1]);
1861%!assert (prod (zeros (0, 2), 1), [1, 1]);
1862%!assert (prod (zeros (0, 2), 2), zeros(0, 1));
1863
1864%!assert (prod (single([1, 2; 3, 4]), 1), single([3, 8]));
1865%!assert (prod (single([1, 2; 3, 4]), 2), single([2; 12]));
1866%!assert (prod (zeros (1, 0, 'single')), single(1));
1867%!assert (prod (zeros (1, 0, 'single'), 1), zeros (1, 0, 'single'));
1868%!assert (prod (zeros (1, 0, 'single'), 2), single(1));
1869%!assert (prod (zeros (0, 1, 'single')), single(1));
1870%!assert (prod (zeros (0, 1, 'single'), 1), single(1));
1871%!assert (prod (zeros (0, 1, 'single'), 2), zeros (0, 1, 'single'));
1872%!assert (prod (zeros (2, 0, 'single')), zeros (1, 0, 'single'));
1873%!assert (prod (zeros (2, 0, 'single'), 1), zeros (1, 0, 'single'));
1874%!assert (prod (zeros (2, 0, 'single'), 2), single([1; 1]));
1875%!assert (prod (zeros (0, 2, 'single')), single([1, 1]));
1876%!assert (prod (zeros (0, 2, 'single'), 1), single([1, 1]));
1877%!assert (prod (zeros (0, 2, 'single'), 2), zeros(0, 1, 'single'));
1878
1879 */
1880
1881#define SINGLE_TYPE_CONCAT(TYPE, EXTRACTOR) \
1882 do \
1883 { \
1884 int dv_len = dv.length (); \
1885 Array<octave_idx_type> ra_idx (dv_len > 1 ? dv_len : 2, 0); \
1886 \
1887 for (int j = 1; j < n_args; j++) \
1888 { \
1889 octave_quit (); \
1890 \
1891 TYPE ra = args(j).EXTRACTOR (); \
1892 \
1893 if (! error_state) \
1894 { \
1895 result.insert (ra, ra_idx); \
1896 \
1897 if (error_state) \
1898 return retval; \
1899 \
1900 dim_vector dv_tmp = args (j).dims (); \
1901 \
1902 if (dim >= dv_len) \
1903 { \
1904 if (j > 1) \
1905 error ("%s: indexing error", fname.c_str ()); \
1906 break; \
1907 } \
1908 else \
1909 ra_idx (dim) += (dim < dv_tmp.length () ? dv_tmp (dim) : 1); \
1910 } \
1911 } \
1912 } \
1913 while (0)
1914
1915#define DO_SINGLE_TYPE_CONCAT(TYPE, EXTRACTOR) \
1916 do \
1917 { \
1918 TYPE result (dv); \
1919 \
1920 SINGLE_TYPE_CONCAT(TYPE, EXTRACTOR); \
1921 \
1922 retval = result; \
1923 } \
1924 while (0)
1925
1926static octave_value
1927do_cat (const octave_value_list& args, std::string fname)
1928{
1929 octave_value retval;
1930
1931 int n_args = args.length ();
1932
1933 if (n_args == 1)
1934 retval = Matrix ();
1935 else if (n_args == 2)
1936 retval = args(1);
1937 else if (n_args > 2)
1938 {
1939 octave_idx_type dim = args(0).int_value () - 1;
1940
1941 if (error_state)
1942 {
1943 error ("cat: expecting first argument to be a integer");
1944 return retval;
1945 }
1946
1947 if (dim >= 0)
1948 {
1949
1950 dim_vector dv = args(1).dims ();
1951 std::string result_type = args(1).class_name ();
1952
1953 bool all_sq_strings_p = args(1).is_sq_string ();
1954 bool all_dq_strings_p = args(1).is_dq_string ();
1955 bool all_real_p = args(1).is_real_type ();
1956 bool any_sparse_p = args(1).is_sparse_type();
1957
1958 for (int i = 2; i < args.length (); i++)
1959 {
1960 // add_dims constructs a dimension vector which holds the
1961 // dimensions of the final array after concatenation.
1962
1963 if (! dv.concat (args(i).dims (), dim))
1964 {
1965 // Dimensions do not match.
1966 error ("cat: dimension mismatch");
1967 return retval;
1968 }
1969
1970 result_type =
1971 get_concat_class (result_type, args(i).class_name ());
1972
1973 if (all_sq_strings_p && ! args(i).is_sq_string ())
1974 all_sq_strings_p = false;
1975 if (all_dq_strings_p && ! args(i).is_dq_string ())
1976 all_dq_strings_p = false;
1977 if (all_real_p && ! args(i).is_real_type ())
1978 all_real_p = false;
1979 if (!any_sparse_p && args(i).is_sparse_type ())
1980 any_sparse_p = true;
1981 }
1982
1983 if (result_type == "double")
1984 {
1985 if (any_sparse_p)
1986 {
1987 if (all_real_p)
1988 DO_SINGLE_TYPE_CONCAT (SparseMatrix, sparse_matrix_value);
1989 else
1990 DO_SINGLE_TYPE_CONCAT (SparseComplexMatrix, sparse_complex_matrix_value);
1991 }
1992 else
1993 {
1994 if (all_real_p)
1995 DO_SINGLE_TYPE_CONCAT (NDArray, array_value);
1996 else
1997 DO_SINGLE_TYPE_CONCAT (ComplexNDArray, complex_array_value);
1998 }
1999 }
2000 else if (result_type == "single")
2001 {
2002 if (all_real_p)
2003 DO_SINGLE_TYPE_CONCAT (FloatNDArray, float_array_value);
2004 else
2005 DO_SINGLE_TYPE_CONCAT (FloatComplexNDArray,
2006 float_complex_array_value);
2007 }
2008 else if (result_type == "char")
2009 {
2010 char type = all_dq_strings_p ? '"' : '\'';
2011
2012 maybe_warn_string_concat (all_dq_strings_p, all_sq_strings_p);
2013
2014 charNDArray result (dv, Vstring_fill_char);
2015
2016 SINGLE_TYPE_CONCAT (charNDArray, char_array_value);
2017
2018 retval = octave_value (result, type);
2019 }
2020 else if (result_type == "logical")
2021 {
2022 if (any_sparse_p)
2023 DO_SINGLE_TYPE_CONCAT (SparseBoolMatrix, sparse_bool_matrix_value);
2024 else
2025 DO_SINGLE_TYPE_CONCAT (boolNDArray, bool_array_value);
2026 }
2027 else if (result_type == "int8")
2028 DO_SINGLE_TYPE_CONCAT (int8NDArray, int8_array_value);
2029 else if (result_type == "int16")
2030 DO_SINGLE_TYPE_CONCAT (int16NDArray, int16_array_value);
2031 else if (result_type == "int32")
2032 DO_SINGLE_TYPE_CONCAT (int32NDArray, int32_array_value);
2033 else if (result_type == "int64")
2034 DO_SINGLE_TYPE_CONCAT (int64NDArray, int64_array_value);
2035 else if (result_type == "uint8")
2036 DO_SINGLE_TYPE_CONCAT (uint8NDArray, uint8_array_value);
2037 else if (result_type == "uint16")
2038 DO_SINGLE_TYPE_CONCAT (uint16NDArray, uint16_array_value);
2039 else if (result_type == "uint32")
2040 DO_SINGLE_TYPE_CONCAT (uint32NDArray, uint32_array_value);
2041 else if (result_type == "uint64")
2042 DO_SINGLE_TYPE_CONCAT (uint64NDArray, uint64_array_value);
2043 else
2044 {
2045 // The lines below might seem crazy, since we take a copy
2046 // of the first argument, resize it to be empty and then resize
2047 // it to be full. This is done since it means that there is no
2048 // recopying of data, as would happen if we used a single resize.
2049 // It should be noted that resize operation is also significantly
2050 // slower than the do_cat_op function, so it makes sense to have
2051 // an empty matrix and copy all data.
2052 //
2053 // We might also start with a empty octave_value using
2054 // tmp = octave_value_typeinfo::lookup_type
2055 // (args(1).type_name());
2056 // and then directly resize. However, for some types there might
2057 // be some additional setup needed, and so this should be avoided.
2058
2059 octave_value tmp = args (1);
2060 tmp = tmp.resize (dim_vector (0,0)).resize (dv);
2061
2062 if (error_state)
2063 return retval;
2064
2065 int dv_len = dv.length ();
2066 Array<octave_idx_type> ra_idx (dv_len, 0);
2067
2068 for (int j = 1; j < n_args; j++)
2069 {
2070 // Can't fast return here to skip empty matrices as something
2071 // like cat(1,[],single([])) must return an empty matrix of
2072 // the right type.
2073 tmp = do_cat_op (tmp, args (j), ra_idx);
2074
2075 if (error_state)
2076 return retval;
2077
2078 dim_vector dv_tmp = args (j).dims ();
2079
2080 if (dim >= dv_len)
2081 {
2082 if (j > 1)
2083 error ("%s: indexing error", fname.c_str ());
2084 break;
2085 }
2086 else
2087 ra_idx (dim) += (dim < dv_tmp.length () ?
2088 dv_tmp (dim) : 1);
2089 }
2090 retval = tmp;
2091 }
2092
2093 if (! error_state)
2094 {
2095 // Reshape, chopping trailing singleton dimensions
2096 dv.chop_trailing_singletons ();
2097 retval = retval.reshape (dv);
2098 }
2099 }
2100 else
2101 error ("%s: invalid dimension argument", fname.c_str ());
2102 }
2103 else
2104 print_usage ();
2105
2106 return retval;
2107}
2108
2109DEFUN (horzcat, args, ,
2110 "-*- texinfo -*-\n\
2111@deftypefn {Built-in Function} {} horzcat (@var{array1}, @var{array2}, @dots{}, @var{arrayN})\n\
2112Return the horizontal concatenation of N-d array objects, @var{array1},\n\
2113@var{array2}, @dots{}, @var{arrayN} along dimension 2.\n\
2114@seealso{cat, vertcat}\n\
2115@end deftypefn")
2116{
2117 octave_value_list args_tmp = args;
2118
2119 int dim = 2;
2120
2121 octave_value d (dim);
2122
2123 args_tmp.prepend (d);
2124
2125 return do_cat (args_tmp, "horzcat");
2126}
2127
2128DEFUN (vertcat, args, ,
2129 "-*- texinfo -*-\n\
2130@deftypefn {Built-in Function} {} vertcat (@var{array1}, @var{array2}, @dots{}, @var{arrayN})\n\
2131Return the vertical concatenation of N-d array objects, @var{array1},\n\
2132@var{array2}, @dots{}, @var{arrayN} along dimension 1.\n\
2133@seealso{cat, horzcat}\n\
2134@end deftypefn")
2135{
2136 octave_value_list args_tmp = args;
2137
2138 int dim = 1;
2139
2140 octave_value d (dim);
2141
2142 args_tmp.prepend (d);
2143
2144 return do_cat (args_tmp, "vertcat");
2145}
2146
2147DEFUN (cat, args, ,
2148 "-*- texinfo -*-\n\
2149@deftypefn {Built-in Function} {} cat (@var{dim}, @var{array1}, @var{array2}, @dots{}, @var{arrayN})\n\
2150Return the concatenation of N-d array objects, @var{array1},\n\
2151@var{array2}, @dots{}, @var{arrayN} along dimension @var{dim}.\n\
2152\n\
2153@example\n\
2154@group\n\
2155A = ones (2, 2);\n\
2156B = zeros (2, 2);\n\
2157cat (2, A, B)\n\
2158@result{} ans =\n\
2159\n\
2160 1 1 0 0\n\
2161 1 1 0 0\n\
2162@end group\n\
2163@end example\n\
2164\n\
2165Alternatively, we can concatenate @var{A} and @var{B} along the\n\
2166second dimension the following way:\n\
2167\n\
2168@example\n\
2169@group\n\
2170[A, B].\n\
2171@end group\n\
2172@end example\n\
2173\n\
2174@var{dim} can be larger than the dimensions of the N-d array objects\n\
2175and the result will thus have @var{dim} dimensions as the\n\
2176following example shows:\n\
2177@example\n\
2178@group\n\
2179cat (4, ones(2, 2), zeros (2, 2))\n\
2180@result{} ans =\n\
2181\n\
2182 ans(:,:,1,1) =\n\
2183\n\
2184 1 1\n\
2185 1 1\n\
2186\n\
2187 ans(:,:,1,2) =\n\
2188 0 0\n\
2189 0 0\n\
2190@end group\n\
2191@end example\n\
2192@seealso{horzcat, vertcat}\n\
2193@end deftypefn")
2194{
2195 return do_cat (args, "cat");
2196}
2197
2198/*
2199
2200%!function ret = testcat (t1, t2, tr, cmplx)
2201%! assert (cat (1, cast ([], t1), cast([], t2)), cast ([], tr));
2202%!
2203%! assert (cat (1, cast (1, t1), cast (2, t2)), cast ([1; 2], tr));
2204%! assert (cat (1, cast (1, t1), cast ([2; 3], t2)), cast ([1; 2; 3], tr));
2205%! assert (cat (1, cast ([1; 2], t1), cast (3, t2)), cast ([1; 2; 3], tr));
2206%! assert (cat (1, cast ([1; 2], t1), cast ([3; 4], t2)), cast ([1; 2; 3; 4], tr));
2207%! assert (cat (2, cast (1, t1), cast (2, t2)), cast ([1, 2], tr));
2208%! assert (cat (2, cast (1, t1), cast ([2, 3], t2)), cast ([1, 2, 3], tr));
2209%! assert (cat (2, cast ([1, 2], t1), cast (3, t2)), cast ([1, 2, 3], tr));
2210%! assert (cat (2, cast ([1, 2], t1), cast ([3, 4], t2)), cast ([1, 2, 3, 4], tr));
2211%!
2212%! assert ([cast(1, t1); cast(2, t2)], cast ([1; 2], tr));
2213%! assert ([cast(1, t1); cast([2; 3], t2)], cast ([1; 2; 3], tr));
2214%! assert ([cast([1; 2], t1); cast(3, t2)], cast ([1; 2; 3], tr));
2215%! assert ([cast([1; 2], t1); cast([3; 4], t2)], cast ([1; 2; 3; 4], tr));
2216%! assert ([cast(1, t1), cast(2, t2)], cast ([1, 2], tr));
2217%! assert ([cast(1, t1), cast([2, 3], t2)], cast ([1, 2, 3], tr));
2218%! assert ([cast([1, 2], t1), cast(3, t2)], cast ([1, 2, 3], tr));
2219%! assert ([cast([1, 2], t1), cast([3, 4], t2)], cast ([1, 2, 3, 4], tr));
2220%!
2221%! if (nargin == 3 || cmplx)
2222%! assert (cat (1, cast (1i, t1), cast (2, t2)), cast ([1i; 2], tr));
2223%! assert (cat (1, cast (1i, t1), cast ([2; 3], t2)), cast ([1i; 2; 3], tr));
2224%! assert (cat (1, cast ([1i; 2], t1), cast (3, t2)), cast ([1i; 2; 3], tr));
2225%! assert (cat (1, cast ([1i; 2], t1), cast ([3; 4], t2)), cast ([1i; 2; 3; 4], tr));
2226%! assert (cat (2, cast (1i, t1), cast (2, t2)), cast ([1i, 2], tr));
2227%! assert (cat (2, cast (1i, t1), cast ([2, 3], t2)), cast ([1i, 2, 3], tr));
2228%! assert (cat (2, cast ([1i, 2], t1), cast (3, t2)), cast ([1i, 2, 3], tr));
2229%! assert (cat (2, cast ([1i, 2], t1), cast ([3, 4], t2)), cast ([1i, 2, 3, 4], tr));
2230%!
2231%! assert ([cast(1i, t1); cast(2, t2)], cast ([1i; 2], tr));
2232%! assert ([cast(1i, t1); cast([2; 3], t2)], cast ([1i; 2; 3], tr));
2233%! assert ([cast([1i; 2], t1); cast(3, t2)], cast ([1i; 2; 3], tr));
2234%! assert ([cast([1i; 2], t1); cast([3; 4], t2)], cast ([1i; 2; 3; 4], tr));
2235%! assert ([cast(1i, t1), cast(2, t2)], cast ([1i, 2], tr));
2236%! assert ([cast(1i, t1), cast([2, 3], t2)], cast ([1i, 2, 3], tr));
2237%! assert ([cast([1i, 2], t1), cast(3, t2)], cast ([1i, 2, 3], tr));
2238%! assert ([cast([1i, 2], t1), cast([3, 4], t2)], cast ([1i, 2, 3, 4], tr));
2239%!
2240%! assert (cat (1, cast (1, t1), cast (2i, t2)), cast ([1; 2i], tr));
2241%! assert (cat (1, cast (1, t1), cast ([2i; 3], t2)), cast ([1; 2i; 3], tr));
2242%! assert (cat (1, cast ([1; 2], t1), cast (3i, t2)), cast ([1; 2; 3i], tr));
2243%! assert (cat (1, cast ([1; 2], t1), cast ([3i; 4], t2)), cast ([1; 2; 3i; 4], tr));
2244%! assert (cat (2, cast (1, t1), cast (2i, t2)), cast ([1, 2i], tr));
2245%! assert (cat (2, cast (1, t1), cast ([2i, 3], t2)), cast ([1, 2i, 3], tr));
2246%! assert (cat (2, cast ([1, 2], t1), cast (3i, t2)), cast ([1, 2, 3i], tr));
2247%! assert (cat (2, cast ([1, 2], t1), cast ([3i, 4], t2)), cast ([1, 2, 3i, 4], tr));
2248%!
2249%! assert ([cast(1, t1); cast(2i, t2)], cast ([1; 2i], tr));
2250%! assert ([cast(1, t1); cast([2i; 3], t2)], cast ([1; 2i; 3], tr));
2251%! assert ([cast([1; 2], t1); cast(3i, t2)], cast ([1; 2; 3i], tr));
2252%! assert ([cast([1; 2], t1); cast([3i; 4], t2)], cast ([1; 2; 3i; 4], tr));
2253%! assert ([cast(1, t1), cast(2i, t2)], cast ([1, 2i], tr));
2254%! assert ([cast(1, t1), cast([2i, 3], t2)], cast ([1, 2i, 3], tr));
2255%! assert ([cast([1, 2], t1), cast(3i, t2)], cast ([1, 2, 3i], tr));
2256%! assert ([cast([1, 2], t1), cast([3i, 4], t2)], cast ([1, 2, 3i, 4], tr));
2257%!
2258%! assert (cat (1, cast (1i, t1), cast (2i, t2)), cast ([1i; 2i], tr));
2259%! assert (cat (1, cast (1i, t1), cast ([2i; 3], t2)), cast ([1i; 2i; 3], tr));
2260%! assert (cat (1, cast ([1i; 2], t1), cast (3i, t2)), cast ([1i; 2; 3i], tr));
2261%! assert (cat (1, cast ([1i; 2], t1), cast ([3i; 4], t2)), cast ([1i; 2; 3i; 4], tr));
2262%! assert (cat (2, cast (1i, t1), cast (2i, t2)), cast ([1i, 2i], tr));
2263%! assert (cat (2, cast (1i, t1), cast ([2i, 3], t2)), cast ([1i, 2i, 3], tr));
2264%! assert (cat (2, cast ([1i, 2], t1), cast (3i, t2)), cast ([1i, 2, 3i], tr));
2265%! assert (cat (2, cast ([1i, 2], t1), cast ([3i, 4], t2)), cast ([1i, 2, 3i, 4], tr));
2266%!
2267%! assert ([cast(1i, t1); cast(2i, t2)], cast ([1i; 2i], tr));
2268%! assert ([cast(1i, t1); cast([2i; 3], t2)], cast ([1i; 2i; 3], tr));
2269%! assert ([cast([1i; 2], t1); cast(3i, t2)], cast ([1i; 2; 3i], tr));
2270%! assert ([cast([1i; 2], t1); cast([3i; 4], t2)], cast ([1i; 2; 3i; 4], tr));
2271%! assert ([cast(1i, t1), cast(2i, t2)], cast ([1i, 2i], tr));
2272%! assert ([cast(1i, t1), cast([2i, 3], t2)], cast ([1i, 2i, 3], tr));
2273%! assert ([cast([1i, 2], t1), cast(3i, t2)], cast ([1i, 2, 3i], tr));
2274%! assert ([cast([1i, 2], t1), cast([3i, 4], t2)], cast ([1i, 2, 3i, 4], tr));
2275%! endif
2276%! ret = true;
2277
2278%!assert (testcat('double', 'double', 'double'));
2279%!assert (testcat('single', 'double', 'single'));
2280%!assert (testcat('double', 'single', 'single'));
2281%!assert (testcat('single', 'single', 'single'));
2282
2283%!assert (testcat('double', 'int8', 'int8', false));
2284%!assert (testcat('int8', 'double', 'int8', false));
2285%!assert (testcat('single', 'int8', 'int8', false));
2286%!assert (testcat('int8', 'single', 'int8', false));
2287%!assert (testcat('int8', 'int8', 'int8', false));
2288%!assert (testcat('double', 'int16', 'int16', false));
2289%!assert (testcat('int16', 'double', 'int16', false));
2290%!assert (testcat('single', 'int16', 'int16', false));
2291%!assert (testcat('int16', 'single', 'int16', false));
2292%!assert (testcat('int16', 'int16', 'int16', false));
2293%!assert (testcat('double', 'int32', 'int32', false));
2294%!assert (testcat('int32', 'double', 'int32', false));
2295%!assert (testcat('single', 'int32', 'int32', false));
2296%!assert (testcat('int32', 'single', 'int32', false));
2297%!assert (testcat('int32', 'int32', 'int32', false));
2298%!assert (testcat('double', 'int64', 'int64', false));
2299%!assert (testcat('int64', 'double', 'int64', false));
2300%!assert (testcat('single', 'int64', 'int64', false));
2301%!assert (testcat('int64', 'single', 'int64', false));
2302%!assert (testcat('int64', 'int64', 'int64', false));
2303
2304%!assert (testcat('double', 'uint8', 'uint8', false));
2305%!assert (testcat('uint8', 'double', 'uint8', false));
2306%!assert (testcat('single', 'uint8', 'uint8', false));
2307%!assert (testcat('uint8', 'single', 'uint8', false));
2308%!assert (testcat('uint8', 'uint8', 'uint8', false));
2309%!assert (testcat('double', 'uint16', 'uint16', false));
2310%!assert (testcat('uint16', 'double', 'uint16', false));
2311%!assert (testcat('single', 'uint16', 'uint16', false));
2312%!assert (testcat('uint16', 'single', 'uint16', false));
2313%!assert (testcat('uint16', 'uint16', 'uint16', false));
2314%!assert (testcat('double', 'uint32', 'uint32', false));
2315%!assert (testcat('uint32', 'double', 'uint32', false));
2316%!assert (testcat('single', 'uint32', 'uint32', false));
2317%!assert (testcat('uint32', 'single', 'uint32', false));
2318%!assert (testcat('uint32', 'uint32', 'uint32', false));
2319%!assert (testcat('double', 'uint64', 'uint64', false));
2320%!assert (testcat('uint64', 'double', 'uint64', false));
2321%!assert (testcat('single', 'uint64', 'uint64', false));
2322%!assert (testcat('uint64', 'single', 'uint64', false));
2323%!assert (testcat('uint64', 'uint64', 'uint64', false));
2324
2325*/
2326
2327static octave_value
2328do_permute (const octave_value_list& args, bool inv)
2329{
2330 octave_value retval;
2331
2332 if (args.length () == 2 && args(1).length () >= args(1).ndims ())
2333 {
2334 Array<int> vec = args(1).int_vector_value ();
2335
2336 // FIXME -- maybe we should create an idx_vector object
2337 // here and pass that to permute?
2338
2339 int n = vec.length ();
2340
2341 for (int i = 0; i < n; i++)
2342 vec(i)--;
2343
2344 octave_value ret = args(0).permute (vec, inv);
2345
2346 if (! error_state)
2347 retval = ret;
2348 }
2349 else
2350 print_usage ();
2351
2352 return retval;
2353}
2354
2355DEFUN (permute, args, ,
2356 "-*- texinfo -*-\n\
2357@deftypefn {Built-in Function} {} permute (@var{a}, @var{perm})\n\
2358Return the generalized transpose for an N-d array object @var{a}.\n\
2359The permutation vector @var{perm} must contain the elements\n\
2360@code{1:ndims(a)} (in any order, but each element must appear just once).\n\
2361@seealso{ipermute}\n\
2362@end deftypefn")
2363{
2364 return do_permute (args, false);
2365}
2366
2367DEFUN (ipermute, args, ,
2368 "-*- texinfo -*-\n\
2369@deftypefn {Built-in Function} {} ipermute (@var{a}, @var{iperm})\n\
2370The inverse of the @code{permute} function. The expression\n\
2371\n\
2372@example\n\
2373ipermute (permute (a, perm), perm)\n\
2374@end example\n\
2375returns the original array @var{a}.\n\
2376@seealso{permute}\n\
2377@end deftypefn")
2378{
2379 return do_permute (args, true);
2380}
2381
2382DEFUN (length, args, ,
2383 "-*- texinfo -*-\n\
2384@deftypefn {Built-in Function} {} length (@var{a})\n\
2385Return the `length' of the object @var{a}. For matrix objects, the\n\
2386length is the number of rows or columns, whichever is greater (this\n\
2387odd definition is used for compatibility with @sc{matlab}).\n\
2388@end deftypefn")
2389{
2390 octave_value retval;
2391
2392 if (args.length () == 1)
2393 retval = args(0).length ();
2394 else
2395 print_usage ();
2396
2397 return retval;
2398}
2399
2400DEFUN (ndims, args, ,
2401 "-*- texinfo -*-\n\
2402@deftypefn {Built-in Function} {} ndims (@var{a})\n\
2403Returns the number of dimensions of array @var{a}.\n\
2404For any array, the result will always be larger than or equal to 2.\n\
2405Trailing singleton dimensions are not counted.\n\
2406@end deftypefn")
2407{
2408 octave_value retval;
2409
2410 if (args.length () == 1)
2411 retval = args(0).ndims ();
2412 else
2413 print_usage ();
2414
2415 return retval;
2416}
2417
2418DEFUN (numel, args, ,
2419 "-*- texinfo -*-\n\
2420@deftypefn {Built-in Function} {} numel (@var{a})\n\
2421@deftypefnx {Built-in Function} {} numel (@var{a}, @var{idx1}, @var{idx2}, @dots{})\n\
2422Returns the number of elements in the object @var{a}.\n\
2423Optionally, if indices @var{idx1}, @var{idx2}, @dots{} are supplied,\n\
2424return the number of elements that would result from the indexing\n\
2425@example\n\
2426 @var{a}(@var{idx1}, @var{idx2}, @dots{})\n\
2427@end example\n\
2428This method is also called when an object appears as lvalue with cs-list\n\
2429indexing, i.e., @code{object@{@dots{}@}} or @code{object(@dots{}).field}.\n\
2430@seealso{size}\n\
2431@end deftypefn")
2432{
2433 octave_value retval;
2434 octave_idx_type nargin = args.length ();
2435
2436 if (nargin == 1)
2437 retval = args(0).numel ();
2438 else if (nargin > 1)
2439 {
2440 // Don't use numel (const octave_value_list&) here as that corresponds to
2441 // an overloaded call, not to builtin!
2442 retval = dims_to_numel (args(0).dims (), args.slice (1, nargin-1));
2443 }
2444 else
2445 print_usage ();
2446
2447 return retval;
2448}
2449
2450DEFUN (size, args, nargout,
2451 "-*- texinfo -*-\n\
2452@deftypefn {Built-in Function} {} size (@var{a}, @var{n})\n\
2453Return the number rows and columns of @var{a}.\n\
2454\n\
2455With one input argument and one output argument, the result is returned\n\
2456in a row vector. If there are multiple output arguments, the number of\n\
2457rows is assigned to the first, and the number of columns to the second,\n\
2458etc. For example,\n\
2459\n\
2460@example\n\
2461@group\n\
2462size ([1, 2; 3, 4; 5, 6])\n\
2463 @result{} [ 3, 2 ]\n\
2464\n\
2465[nr, nc] = size ([1, 2; 3, 4; 5, 6])\n\
2466 @result{} nr = 3\n\
2467 @result{} nc = 2\n\
2468@end group\n\
2469@end example\n\
2470\n\
2471If given a second argument, @code{size} will return the size of the\n\
2472corresponding dimension. For example\n\
2473\n\
2474@example\n\
2475@group\n\
2476size ([1, 2; 3, 4; 5, 6], 2)\n\
2477 @result{} 2\n\
2478@end group\n\
2479@end example\n\
2480\n\
2481@noindent\n\
2482returns the number of columns in the given matrix.\n\
2483@seealso{numel}\n\
2484@end deftypefn")
2485{
2486 octave_value_list retval;
2487
2488 int nargin = args.length ();
2489
2490 if (nargin == 1)
2491 {
2492 const dim_vector dimensions = args(0).dims ();
2493
2494 if (nargout > 1)
2495 {
2496 const dim_vector rdims = dimensions.redim (nargout);
2497 retval.resize (nargout);
2498 for (int i = 0; i < nargout; i++)
2499 retval(i) = rdims(i);
2500 }
2501 else
2502 {
2503 int ndims = dimensions.length ();
2504
2505 NoAlias<Matrix> m (1, ndims);
2506
2507 for (int i = 0; i < ndims; i++)
2508 m(i) = dimensions(i);
2509
2510 retval(0) = m;
2511 }
2512 }
2513 else if (nargin == 2 && nargout < 2)
2514 {
2515 octave_idx_type nd = args(1).int_value (true);
2516
2517 if (error_state)
2518 error ("size: expecting scalar as second argument");
2519 else
2520 {
2521 const dim_vector dv = args(0).dims ();
2522
2523 if (nd > 0)
2524 {
2525 if (nd <= dv.length ())
2526 retval(0) = dv(nd-1);
2527 else
2528 retval(0) = 1;
2529 }
2530 else
2531 error ("size: requested dimension (= %d) out of range", nd);
2532 }
2533 }
2534 else
2535 print_usage ();
2536
2537 return retval;
2538}
2539
2540DEFUN (size_equal, args, ,
2541 "-*- texinfo -*-\n\
2542@deftypefn {Built-in Function} {} size_equal (@var{a}, @var{b}, @dots{})\n\
2543Return true if the dimensions of all arguments agree.\n\
2544Trailing singleton dimensions are ignored.\n\
2545Called with a single or no argument, size_equal returns true.\n\
2546@seealso{size, numel}\n\
2547@end deftypefn")
2548{
2549 octave_value retval;
2550
2551 int nargin = args.length ();
2552
2553 retval = true;
2554
2555 if (nargin >= 1)
2556 {
2557 dim_vector a_dims = args(0).dims ();
2558
2559 for (int i = 1; i < nargin; ++i)
2560 {
2561 dim_vector b_dims = args(i).dims ();
2562
2563 if (a_dims != b_dims)
2564 {
2565 retval = false;
2566 break;
2567 }
2568 }
2569 }
2570
2571 return retval;
2572}
2573
2574DEFUN (nnz, args, ,
2575 "-*- texinfo -*-\n\
2576@deftypefn {Built-in Function} {@var{scalar} =} nnz (@var{a})\n\
2577Returns the number of non zero elements in @var{a}.\n\
2578@seealso{sparse}\n\
2579@end deftypefn")
2580{
2581 octave_value retval;
2582
2583 if (args.length () == 1)
2584 retval = args(0).nnz ();
2585 else
2586 print_usage ();
2587
2588 return retval;
2589}
2590
2591DEFUN (nzmax, args, ,
2592 "-*- texinfo -*-\n\
2593@deftypefn {Built-in Function} {@var{scalar} =} nzmax (@var{SM})\n\
2594Return the amount of storage allocated to the sparse matrix @var{SM}.\n\
2595Note that Octave tends to crop unused memory at the first opportunity\n\
2596for sparse objects. There are some cases of user created sparse objects\n\
2597where the value returned by @dfn{nzmax} will not be the same as @dfn{nnz},\n\
2598but in general they will give the same result.\n\
2599@seealso{sparse, spalloc}\n\
2600@end deftypefn")
2601{
2602 octave_value retval;
2603
2604 if (args.length() == 1)
2605 retval = args(0).nzmax ();
2606 else
2607 print_usage ();
2608
2609 return retval;
2610}
2611
2612DEFUN (rows, args, ,
2613 "-*- texinfo -*-\n\
2614@deftypefn {Built-in Function} {} rows (@var{a})\n\
2615Return the number of rows of @var{a}.\n\
2616@seealso{size, numel, columns, length, isscalar, isvector, ismatrix}\n\
2617@end deftypefn")
2618{
2619 octave_value retval;
2620
2621 if (args.length () == 1)
2622 retval = args(0).rows ();
2623 else
2624 print_usage ();
2625
2626 return retval;
2627}
2628
2629DEFUN (columns, args, ,
2630 "-*- texinfo -*-\n\
2631@deftypefn {Built-in Function} {} columns (@var{a})\n\
2632Return the number of columns of @var{a}.\n\
2633@seealso{size, numel, rows, length, isscalar, isvector, ismatrix}\n\
2634@end deftypefn")
2635{
2636 octave_value retval;
2637
2638 if (args.length () == 1)
2639 retval = args(0).columns ();
2640 else
2641 print_usage ();
2642
2643 return retval;
2644}
2645
2646DEFUN (sum, args, ,
2647 "-*- texinfo -*-\n\
2648@deftypefn {Built-in Function} {} sum (@var{x})\n\
2649@deftypefnx {Built-in Function} {} sum (@var{x}, @var{dim})\n\
2650@deftypefnx {Built-in Function} {} sum (@dots{}, 'native')\n\
2651@deftypefnx {Built-in Function} {} sum (@dots{}, 'double')\n\
2652@deftypefnx {Built-in Function} {} sum (@dots{}, 'extra')\n\
2653Sum of elements along dimension @var{dim}. If @var{dim} is\n\
2654omitted, it defaults to the first non-singleton dimension.\n\
2655\n\
2656If the optional argument 'native' is given, then the sum is performed\n\
2657in the same type as the original argument, rather than in the default\n\
2658double type. For example\n\
2659\n\
2660@example\n\
2661@group\n\
2662sum ([true, true])\n\
2663 @result{} 2\n\
2664sum ([true, true], 'native')\n\
2665 @result{} true\n\
2666@end group\n\
2667@end example\n\
2668 \n\
2669On the contrary, if 'double' is given, the sum is performed in double precision\n\
2670even for single precision inputs.\n\
2671\n\
2672For double precision inputs, 'extra' indicates that a more accurate algorithm\n\
2673than straightforward summation is to be used. For single precision inputs, 'extra' is\n\
2674the same as 'double'. Otherwise, 'extra' has no effect.\n\
2675@seealso{cumsum, sumsq, prod}\n\
2676@end deftypefn")
2677{
2678 octave_value retval;
2679
2680 int nargin = args.length ();
2681
2682 bool isnative = false;
2683 bool isdouble = false;
2684 bool isextra = false;
2685
2686 if (nargin > 1 && args(nargin - 1).is_string ())
2687 {
2688 std::string str = args(nargin - 1).string_value ();
2689
2690 if (! error_state)
2691 {
2692 if (str == "native")
2693 isnative = true;
2694 else if (str == "double")
2695 isdouble = true;
2696 else if (str == "extra")
2697 isextra = true;
2698 else
2699 error ("sum: unrecognized string argument");
2700 nargin --;
2701 }
2702 }
2703
2704 if (error_state)
2705 return retval;
2706
2707 if (nargin == 1 || nargin == 2)
2708 {
2709 octave_value arg = args(0);
2710
2711 int dim = -1;
2712 if (nargin == 2)
2713 {
2714 dim = args(1).int_value () - 1;
2715 if (dim < 0)
2716 error ("sum: invalid dimension argument = %d", dim + 1);
2717 }
2718
2719 if (! error_state)
2720 {
2721 switch (arg.builtin_type ())
2722 {
2723 case btyp_double:
2724 if (arg.is_sparse_type ())
2725 {
2726 if (isextra)
2727 warning ("sum: 'extra' not yet implemented for sparse matrices");
2728 retval = arg.sparse_matrix_value ().sum (dim);
2729 }
2730 else if (isextra)
2731 retval = arg.array_value ().xsum (dim);
2732 else
2733 retval = arg.array_value ().sum (dim);
2734 break;
2735 case btyp_complex:
2736 if (arg.is_sparse_type ())
2737 {
2738 if (isextra)
2739 warning ("sum: 'extra' not yet implemented for sparse matrices");
2740 retval = arg.sparse_complex_matrix_value ().sum (dim);
2741 }
2742 else if (isextra)
2743 retval = arg.complex_array_value ().xsum (dim);
2744 else
2745 retval = arg.complex_array_value ().sum (dim);
2746 break;
2747 case btyp_float:
2748 if (isdouble || isextra)
2749 retval = arg.float_array_value ().dsum (dim);
2750 else
2751 retval = arg.float_array_value ().sum (dim);
2752 break;
2753 case btyp_float_complex:
2754 if (isdouble || isextra)
2755 retval = arg.float_complex_array_value ().dsum (dim);
2756 else
2757 retval = arg.float_complex_array_value ().sum (dim);
2758 break;
2759
2760#define MAKE_INT_BRANCH(X) \
2761 case btyp_ ## X: \
2762 if (isnative) \
2763 retval = arg.X ## _array_value ().sum (dim); \
2764 else \
2765 retval = arg.X ## _array_value ().dsum (dim); \
2766 break
2767 MAKE_INT_BRANCH (int8);
2768 MAKE_INT_BRANCH (int16);
2769 MAKE_INT_BRANCH (int32);
2770 MAKE_INT_BRANCH (int64);
2771 MAKE_INT_BRANCH (uint8);
2772 MAKE_INT_BRANCH (uint16);
2773 MAKE_INT_BRANCH (uint32);
2774 MAKE_INT_BRANCH (uint64);
2775#undef MAKE_INT_BRANCH
2776
2777 case btyp_bool:
2778 if (arg.is_sparse_type ())
2779 {
2780 if (isnative)
2781 retval = arg.sparse_bool_matrix_value ().any (dim);
2782 else
2783 retval = arg.sparse_matrix_value ().sum (dim);
2784 }
2785 else if (isnative)
2786 retval = arg.bool_array_value ().any (dim);
2787 else
2788 retval = arg.bool_array_value ().sum (dim);
2789 break;
2790
2791 default:
2792 gripe_wrong_type_arg ("sum", arg);
2793 }
2794 }
2795 }
2796 else
2797 print_usage ();
2798
2799 return retval;
2800}
2801
2802/*
2803
2804%!assert (sum([true,true]), 2)
2805%!assert (sum([true,true],'native'), true)
2806%!assert (sum(int8([127,10,-20])), 117);
2807%!assert (sum(int8([127,10,-20]),'native'), int8(107));
2808
2809%!assert(sum ([1, 2, 3]), 6)
2810%!assert(sum ([-1; -2; -3]), -6);
2811%!assert(sum ([i, 2+i, -3+2i, 4]), 3+4i);
2812%!assert(sum ([1, 2, 3; i, 2i, 3i; 1+i, 2+2i, 3+3i]), [2+2i, 4+4i, 6+6i]);
2813
2814%!assert(sum (single([1, 2, 3])), single(6))
2815%!assert(sum (single([-1; -2; -3])), single(-6));
2816%!assert(sum (single([i, 2+i, -3+2i, 4])), single(3+4i));
2817%!assert(sum (single([1, 2, 3; i, 2i, 3i; 1+i, 2+2i, 3+3i])), single([2+2i, 4+4i, 6+6i]));
2818
2819%!error <Invalid call to sum.*> sum ();
2820
2821%!assert (sum ([1, 2; 3, 4], 1), [4, 6]);
2822%!assert (sum ([1, 2; 3, 4], 2), [3; 7]);
2823%!assert (sum (zeros (1, 0)), 0);
2824%!assert (sum (zeros (1, 0), 1), zeros(1, 0));
2825%!assert (sum (zeros (1, 0), 2), 0);
2826%!assert (sum (zeros (0, 1)), 0);
2827%!assert (sum (zeros (0, 1), 1), 0);
2828%!assert (sum (zeros (0, 1), 2), zeros(0, 1));
2829%!assert (sum (zeros (2, 0)), zeros(1, 0));
2830%!assert (sum (zeros (2, 0), 1), zeros(1, 0));
2831%!assert (sum (zeros (2, 0), 2), [0; 0]);
2832%!assert (sum (zeros (0, 2)), [0, 0]);
2833%!assert (sum (zeros (0, 2), 1), [0, 0]);
2834%!assert (sum (zeros (0, 2), 2), zeros(0, 1));
2835%!assert (sum (zeros (2, 2, 0, 3)), zeros(1, 2, 0, 3));
2836%!assert (sum (zeros (2, 2, 0, 3), 2), zeros(2, 1, 0, 3));
2837%!assert (sum (zeros (2, 2, 0, 3), 3), zeros(2, 2, 1, 3));
2838%!assert (sum (zeros (2, 2, 0, 3), 4), zeros(2, 2, 0));
2839%!assert (sum (zeros (2, 2, 0, 3), 7), zeros(2, 2, 0, 3));
2840
2841%!assert (sum (single([1, 2; 3, 4]), 1), single([4, 6]));
2842%!assert (sum (single([1, 2; 3, 4]), 2), single([3; 7]));
2843%!assert (sum (zeros (1, 0, 'single')), single(0));
2844%!assert (sum (zeros (1, 0, 'single'), 1), zeros(1, 0, 'single'));
2845%!assert (sum (zeros (1, 0, 'single'), 2), single(0));
2846%!assert (sum (zeros (0, 1, 'single')), single(0));
2847%!assert (sum (zeros (0, 1, 'single'), 1), single(0));
2848%!assert (sum (zeros (0, 1, 'single'), 2), zeros(0, 1, 'single'));
2849%!assert (sum (zeros (2, 0, 'single')), zeros(1, 0, 'single'));
2850%!assert (sum (zeros (2, 0, 'single'), 1), zeros(1, 0, 'single'));
2851%!assert (sum (zeros (2, 0, 'single'), 2), single([0; 0]));
2852%!assert (sum (zeros (0, 2, 'single')), single([0, 0]));
2853%!assert (sum (zeros (0, 2, 'single'), 1), single([0, 0]));
2854%!assert (sum (zeros (0, 2, 'single'), 2), zeros(0, 1, 'single'));
2855%!assert (sum (zeros (2, 2, 0, 3, 'single')), zeros(1, 2, 0, 3, 'single'));
2856%!assert (sum (zeros (2, 2, 0, 3, 'single'), 2), zeros(2, 1, 0, 3, 'single'));
2857%!assert (sum (zeros (2, 2, 0, 3, 'single'), 3), zeros(2, 2, 1, 3, 'single'));
2858%!assert (sum (zeros (2, 2, 0, 3, 'single'), 4), zeros(2, 2, 0, 'single'));
2859%!assert (sum (zeros (2, 2, 0, 3, 'single'), 7), zeros(2, 2, 0, 3, 'single'));
2860
2861*/
2862
2863DEFUN (sumsq, args, ,
2864 "-*- texinfo -*-\n\
2865@deftypefn {Built-in Function} {} sumsq (@var{x})\n\
2866@deftypefnx {Built-in Function} {} sumsq (@var{x}, @var{dim})\n\
2867Sum of squares of elements along dimension @var{dim}. If @var{dim}\n\
2868is omitted, it defaults to the first non-singleton dimension.\n\
2869\n\
2870This function is conceptually equivalent to computing\n\
2871@example\n\
2872sum (x .* conj (x), dim)\n\
2873@end example\n\
2874but it uses less memory and avoids calling @code{conj} if @var{x} is real.\n\
2875@seealso{sum}\n\
2876@end deftypefn")
2877{
2878 DATA_REDUCTION (sumsq);
2879}
2880
2881/*
2882
2883%!assert(sumsq ([1, 2, 3]), 14)
2884%!assert(sumsq ([-1; -2; 4i]), 21);
2885%!assert(sumsq ([1, 2, 3; 2, 3, 4; 4i, 6i, 2]), [21, 49, 29]);
2886
2887%!assert(sumsq (single([1, 2, 3])), single(14))
2888%!assert(sumsq (single([-1; -2; 4i])), single(21));
2889%!assert(sumsq (single([1, 2, 3; 2, 3, 4; 4i, 6i, 2])), single([21, 49, 29]));
2890
2891%!error <Invalid call to sumsq.*> sumsq ();
2892
2893%!assert (sumsq ([1, 2; 3, 4], 1), [10, 20]);
2894%!assert (sumsq ([1, 2; 3, 4], 2), [5; 25]);
2895
2896%!assert (sumsq (single([1, 2; 3, 4]), 1), single([10, 20]));
2897%!assert (sumsq (single([1, 2; 3, 4]), 2), single([5; 25]));
2898
2899 */
2900
2901DEFUN (islogical, args, ,
2902 "-*- texinfo -*-\n\
2903@deftypefn {Built-in Function} {} islogical (@var{x})\n\
2904Return true if @var{x} is a logical object.\n\
2905@end deftypefn")
2906{
2907 octave_value retval;
2908
2909 if (args.length () == 1)
2910 retval = args(0).is_bool_type ();
2911 else
2912 print_usage ();
2913
2914 return retval;
2915}
2916
2917DEFALIAS (isbool, islogical);
2918
2919/*
2920
2921%!assert (islogical(true), true)
2922%!assert (islogical(false), true)
2923%!assert (islogical([true, false]), true)
2924%!assert (islogical(1), false)
2925%!assert (islogical(1i), false)
2926%!assert (islogical([1,1]), false)
2927%!assert (islogical(single(1)), false)
2928%!assert (islogical(single(1i)), false)
2929%!assert (islogical(single([1,1])), false)
2930
2931 */
2932
2933DEFUN (isinteger, args, ,
2934 "-*- texinfo -*-\n\
2935@deftypefn {Built-in Function} {} isinteger (@var{x})\n\
2936Return true if @var{x} is an integer object (int8, uint8, int16, etc.).\n\
2937Note that @code{isinteger (14)} is false because numeric constants in\n\
2938Octave are double precision floating point values.\n\
2939@seealso{isreal, isnumeric, class, isa}\n\
2940@end deftypefn")
2941{
2942 octave_value retval;
2943
2944 if (args.length () == 1)
2945 retval = args(0).is_integer_type ();
2946 else
2947 print_usage ();
2948
2949 return retval;
2950}
2951
2952DEFUN (iscomplex, args, ,
2953 "-*- texinfo -*-\n\
2954@deftypefn {Built-in Function} {} iscomplex (@var{x})\n\
2955Return true if @var{x} is a complex-valued numeric object.\n\
2956@end deftypefn")
2957{
2958 octave_value retval;
2959
2960 if (args.length () == 1)
2961 retval = args(0).is_complex_type ();
2962 else
2963 print_usage ();
2964
2965 return retval;
2966}
2967
2968DEFUN (isfloat, args, ,
2969 "-*- texinfo -*-\n\
2970@deftypefn {Built-in Function} {} isfloat (@var{x})\n\
2971Return true if @var{x} is a floating-point numeric object.\n\
2972@end deftypefn")
2973{
2974 octave_value retval;
2975
2976 if (args.length () == 1)
2977 retval = args(0).is_float_type ();
2978 else
2979 print_usage ();
2980
2981 return retval;
2982}
2983
2984// FIXME -- perhaps this should be implemented with an
2985// octave_value member function?
2986
2987DEFUN (complex, args, ,
2988 "-*- texinfo -*-\n\
2989@deftypefn {Built-in Function} {} complex (@var{x})\n\
2990@deftypefnx {Built-in Function} {} complex (@var{re}, @var{im})\n\
2991Return a complex result from real arguments. With 1 real argument @var{x},\n\
2992return the complex result @code{@var{x} + 0i}. With 2 real arguments,\n\
2993return the complex result @code{@var{re} + @var{im}}. @code{complex} can\n\
2994often be more convenient than expressions such as @code{a + i*b}.\n\
2995For example:\n\
2996\n\
2997@example\n\
2998@group\n\
2999complex ([1, 2], [3, 4])\n\
3000@result{}\n\
3001 1 + 3i 2 + 4i\n\
3002@end group\n\
3003@end example\n\
3004@seealso{real, imag, iscomplex}\n\
3005@end deftypefn")
3006{
3007 octave_value retval;
3008
3009 int nargin = args.length ();
3010
3011 if (nargin == 1)
3012 {
3013 octave_value arg = args(0);
3014
3015 if (arg.is_complex_type ())
3016 retval = arg;
3017 else
3018 {
3019 if (arg.is_sparse_type ())
3020 {
3021 SparseComplexMatrix val = arg.sparse_complex_matrix_value ();
3022
3023 if (! error_state)
3024 retval = octave_value (new octave_sparse_complex_matrix (val));
3025 }
3026 else if (arg.is_single_type ())
3027 {
3028 if (arg.numel () == 1)
3029 {
3030 FloatComplex val = arg.float_complex_value ();
3031
3032 if (! error_state)
3033 retval = octave_value (new octave_float_complex (val));
3034 }
3035 else
3036 {
3037 FloatComplexNDArray val = arg.float_complex_array_value ();
3038
3039 if (! error_state)
3040 retval = octave_value (new octave_float_complex_matrix (val));
3041 }
3042 }
3043 else
3044 {
3045 if (arg.numel () == 1)
3046 {
3047 Complex val = arg.complex_value ();
3048
3049 if (! error_state)
3050 retval = octave_value (new octave_complex (val));
3051 }
3052 else
3053 {
3054 ComplexNDArray val = arg.complex_array_value ();
3055
3056 if (! error_state)
3057 retval = octave_value (new octave_complex_matrix (val));
3058 }
3059 }
3060
3061 if (error_state)
3062 error ("complex: invalid conversion");
3063 }
3064 }
3065 else if (nargin == 2)
3066 {
3067 octave_value re = args(0);
3068 octave_value im = args(1);
3069
3070 if (re.is_sparse_type () && im.is_sparse_type ())
3071 {
3072 const SparseMatrix re_val = re.sparse_matrix_value ();
3073 const SparseMatrix im_val = im.sparse_matrix_value ();
3074
3075 if (!error_state)
3076 {
3077 if (re.numel () == 1)
3078 {
3079 SparseComplexMatrix result;
3080 if (re_val.nnz () == 0)
3081 result = Complex(0, 1) * SparseComplexMatrix (im_val);
3082 else
3083 {
3084 result = SparseComplexMatrix (im_val.dims (), re_val (0));
3085 octave_idx_type nr = im_val.rows ();
3086 octave_idx_type nc = im_val.cols ();
3087
3088 for (octave_idx_type j = 0; j < nc; j++)
3089 {
3090 octave_idx_type off = j * nr;
3091 for (octave_idx_type i = im_val.cidx(j);
3092 i < im_val.cidx(j + 1); i++)
3093 result.data (im_val.ridx(i) + off) =
3094 result.data (im_val.ridx(i) + off) +
3095 Complex (0, im_val.data (i));
3096 }
3097 }
3098 retval = octave_value (new octave_sparse_complex_matrix (result));
3099 }
3100 else if (im.numel () == 1)
3101 {
3102 SparseComplexMatrix result;
3103 if (im_val.nnz () == 0)
3104 result = SparseComplexMatrix (re_val);
3105 else
3106 {
3107 result = SparseComplexMatrix (re_val.rows(), re_val.cols(), Complex(0, im_val (0)));
3108 octave_idx_type nr = re_val.rows ();
3109 octave_idx_type nc = re_val.cols ();
3110
3111 for (octave_idx_type j = 0; j < nc; j++)
3112 {
3113 octave_idx_type off = j * nr;
3114 for (octave_idx_type i = re_val.cidx(j);
3115 i < re_val.cidx(j + 1); i++)
3116 result.data (re_val.ridx(i) + off) =
3117 result.data (re_val.ridx(i) + off) +
3118 re_val.data (i);
3119 }
3120 }
3121 retval = octave_value (new octave_sparse_complex_matrix (result));
3122 }
3123 else
3124 {
3125 if (re_val.dims () == im_val.dims ())
3126 {
3127 SparseComplexMatrix result = SparseComplexMatrix(re_val)
3128 + Complex(0, 1) * SparseComplexMatrix (im_val);
3129 retval = octave_value (new octave_sparse_complex_matrix (result));
3130 }
3131 else
3132 error ("complex: dimension mismatch");
3133 }
3134 }
3135 }
3136 else if (re.is_single_type () || im.is_single_type ())
3137 {
3138 if (re.numel () == 1)
3139 {
3140 float re_val = re.float_value ();
3141
3142 if (im.numel () == 1)
3143 {
3144 float im_val = im.double_value ();
3145
3146 if (! error_state)
3147 retval = octave_value (new octave_float_complex (FloatComplex (re_val, im_val)));
3148 }
3149 else
3150 {
3151 const FloatNDArray im_val = im.float_array_value ();
3152
3153 if (! error_state)
3154 {
3155 FloatComplexNDArray result (im_val.dims (), FloatComplex ());
3156
3157 for (octave_idx_type i = 0; i < im_val.numel (); i++)
3158 result.xelem (i) = FloatComplex (re_val, im_val(i));
3159
3160 retval = octave_value (new octave_float_complex_matrix (result));
3161 }
3162 }
3163 }
3164 else
3165 {
3166 const FloatNDArray re_val = re.float_array_value ();
3167
3168 if (im.numel () == 1)
3169 {
3170 float im_val = im.float_value ();
3171
3172 if (! error_state)
3173 {
3174 FloatComplexNDArray result (re_val.dims (), FloatComplex ());
3175
3176 for (octave_idx_type i = 0; i < re_val.numel (); i++)
3177 result.xelem (i) = FloatComplex (re_val(i), im_val);
3178
3179 retval = octave_value (new octave_float_complex_matrix (result));
3180 }
3181 }
3182 else
3183 {
3184 const FloatNDArray im_val = im.float_array_value ();
3185
3186 if (! error_state)
3187 {
3188 if (re_val.dims () == im_val.dims ())
3189 {
3190 FloatComplexNDArray result (re_val.dims (), FloatComplex ());
3191
3192 for (octave_idx_type i = 0; i < re_val.numel (); i++)
3193 result.xelem (i) = FloatComplex (re_val(i), im_val(i));
3194
3195 retval = octave_value (new octave_float_complex_matrix (result));
3196 }
3197 else
3198 error ("complex: dimension mismatch");
3199 }
3200 }
3201 }
3202 }
3203 else if (re.numel () == 1)
3204 {
3205 double re_val = re.double_value ();
3206
3207 if (im.numel () == 1)
3208 {
3209 double im_val = im.double_value ();
3210
3211 if (! error_state)
3212 retval = octave_value (new octave_complex (Complex (re_val, im_val)));
3213 }
3214 else
3215 {
3216 const NDArray im_val = im.array_value ();
3217
3218 if (! error_state)
3219 {
3220 ComplexNDArray result (im_val.dims (), Complex ());
3221
3222 for (octave_idx_type i = 0; i < im_val.numel (); i++)
3223 result.xelem (i) = Complex (re_val, im_val(i));
3224
3225 retval = octave_value (new octave_complex_matrix (result));
3226 }
3227 }
3228 }
3229 else
3230 {
3231 const NDArray re_val = re.array_value ();
3232
3233 if (im.numel () == 1)
3234 {
3235 double im_val = im.double_value ();
3236
3237 if (! error_state)
3238 {
3239 ComplexNDArray result (re_val.dims (), Complex ());
3240
3241 for (octave_idx_type i = 0; i < re_val.numel (); i++)
3242 result.xelem (i) = Complex (re_val(i), im_val);
3243
3244 retval = octave_value (new octave_complex_matrix (result));
3245 }
3246 }
3247 else
3248 {
3249 const NDArray im_val = im.array_value ();
3250
3251 if (! error_state)
3252 {
3253 if (re_val.dims () == im_val.dims ())
3254 {
3255 ComplexNDArray result (re_val.dims (), Complex ());
3256
3257 for (octave_idx_type i = 0; i < re_val.numel (); i++)
3258 result.xelem (i) = Complex (re_val(i), im_val(i));
3259
3260 retval = octave_value (new octave_complex_matrix (result));
3261 }
3262 else
3263 error ("complex: dimension mismatch");
3264 }
3265 }
3266 }
3267
3268 if (error_state)
3269 error ("complex: invalid conversion");
3270 }
3271 else
3272 print_usage ();
3273
3274 return retval;
3275}
3276
3277DEFUN (isreal, args, ,
3278 "-*- texinfo -*-\n\
3279@deftypefn {Built-in Function} {} isreal (@var{x})\n\
3280Return true if @var{x} is a real-valued numeric object.\n\
3281@end deftypefn")
3282{
3283 octave_value retval;
3284
3285 if (args.length () == 1)
3286 retval = args(0).is_real_type ();
3287 else
3288 print_usage ();
3289
3290 return retval;
3291}
3292
3293DEFUN (isempty, args, ,
3294 "-*- texinfo -*-\n\
3295@deftypefn {Built-in Function} {} isempty (@var{a})\n\
3296Return 1 if @var{a} is an empty matrix (either the number of rows, or\n\
3297the number of columns, or both are zero). Otherwise, return 0.\n\
3298@end deftypefn")
3299{
3300 octave_value retval = false;
3301
3302 if (args.length () == 1)
3303 retval = args(0).is_empty ();
3304 else
3305 print_usage ();
3306
3307 return retval;
3308}
3309
3310DEFUN (isnumeric, args, ,
3311 "-*- texinfo -*-\n\
3312@deftypefn {Built-in Function} {} isnumeric (@var{x})\n\
3313Return nonzero if @var{x} is a numeric object.\n\
3314@end deftypefn")
3315{
3316 octave_value retval;
3317
3318 if (args.length () == 1)
3319 retval = args(0).is_numeric_type ();
3320 else
3321 print_usage ();
3322
3323 return retval;
3324}
3325
3326DEFUN (islist, args, ,
3327 "-*- texinfo -*-\n\
3328@deftypefn {Built-in Function} {} islist (@var{x})\n\
3329Return nonzero if @var{x} is a list.\n\
3330@end deftypefn")
3331{
3332 octave_value retval;
3333
3334 if (args.length () == 1)
3335 retval = args(0).is_list ();
3336 else
3337 print_usage ();
3338
3339 return retval;
3340}
3341
3342DEFUN (ismatrix, args, ,
3343 "-*- texinfo -*-\n\
3344@deftypefn {Built-in Function} {} ismatrix (@var{a})\n\
3345Return 1 if @var{a} is a matrix. Otherwise, return 0.\n\
3346@end deftypefn")
3347{
3348 octave_value retval = false;
3349
3350 if (args.length () == 1)
3351 {
3352 octave_value arg = args(0);
3353
3354 retval = arg.is_matrix_type () || arg.is_scalar_type () || arg.is_range ();
3355 }
3356 else
3357 print_usage ();
3358
3359 return retval;
3360}
3361
3362/*
3363
3364%!assert(ismatrix ([]));
3365%!assert(ismatrix (1));
3366%!assert(ismatrix ([1, 2, 3]));
3367%!assert(ismatrix ([1, 2; 3, 4]));
3368%!assert(ismatrix (zeros (3, 2, 4)));
3369
3370%!assert(ismatrix (single([])));
3371%!assert(ismatrix (single(1)));
3372%!assert(ismatrix (single([1, 2, 3])));
3373%!assert(ismatrix (single([1, 2; 3, 4])));
3374
3375%!assert(ismatrix ("t"));
3376%!assert(ismatrix ("test"));
3377%!assert(ismatrix (["test"; "ing"]));
3378
3379%!test
3380%! s.a = 1;
3381%! assert(ismatrix (s), false);
3382
3383%!error <Invalid call to ismatrix.*> ismatrix ();
3384%!error <Invalid call to ismatrix.*> ismatrix ([1, 2; 3, 4], 2);
3385
3386 */
3387
3388static octave_value
3389fill_matrix (const octave_value_list& args, int val, const char *fcn)
3390{
3391 octave_value retval;
3392
3393 int nargin = args.length ();
3394
3395 oct_data_conv::data_type dt = oct_data_conv::dt_double;
3396
3397 dim_vector dims (1, 1);
3398
3399 if (nargin > 0 && args(nargin-1).is_string ())
3400 {
3401 std::string nm = args(nargin-1).string_value ();
3402 nargin--;
3403
3404 dt = oct_data_conv::string_to_data_type (nm);
3405
3406 if (error_state)
3407 return retval;
3408 }
3409
3410 switch (nargin)
3411 {
3412 case 0:
3413 break;
3414
3415 case 1:
3416 get_dimensions (args(0), fcn, dims);
3417 break;
3418
3419 default:
3420 {
3421 dims.resize (nargin);
3422
3423 for (int i = 0; i < nargin; i++)
3424 {
3425 dims(i) = args(i).is_empty () ? 0 : args(i).idx_type_value ();
3426
3427 if (error_state)
3428 {
3429 error ("%s: expecting scalar integer arguments", fcn);
3430 break;
3431 }
3432 }
3433 }
3434 break;
3435 }
3436
3437 if (! error_state)
3438 {
3439 dims.chop_trailing_singletons ();
3440
3441 check_dimensions (dims, fcn);
3442
3443 // FIXME -- perhaps this should be made extensible by
3444 // using the class name to lookup a function to call to create
3445 // the new value.
3446
3447 // Note that automatic narrowing will handle conversion from
3448 // NDArray to scalar.
3449
3450 if (! error_state)
3451 {
3452 switch (dt)
3453 {
3454 case oct_data_conv::dt_int8:
3455 retval = int8NDArray (dims, val);
3456 break;
3457
3458 case oct_data_conv::dt_uint8:
3459 retval = uint8NDArray (dims, val);
3460 break;
3461
3462 case oct_data_conv::dt_int16:
3463 retval = int16NDArray (dims, val);
3464 break;
3465
3466 case oct_data_conv::dt_uint16:
3467 retval = uint16NDArray (dims, val);
3468 break;
3469
3470 case oct_data_conv::dt_int32:
3471 retval = int32NDArray (dims, val);
3472 break;
3473
3474 case oct_data_conv::dt_uint32:
3475 retval = uint32NDArray (dims, val);
3476 break;
3477
3478 case oct_data_conv::dt_int64:
3479 retval = int64NDArray (dims, val);
3480 break;
3481
3482 case oct_data_conv::dt_uint64:
3483 retval = uint64NDArray (dims, val);
3484 break;
3485
3486 case oct_data_conv::dt_single:
3487 retval = FloatNDArray (dims, val);
3488 break;
3489
3490 case oct_data_conv::dt_double:
3491 {
3492 if (val == 1 && dims.length () == 2 && dims (0) == 1)
3493 retval = Range (1.0, 0.0, dims (1)); // packed form
3494 else
3495 retval = NDArray (dims, val);
3496 }
3497 break;
3498
3499 case oct_data_conv::dt_logical:
3500 retval = boolNDArray (dims, val);
3501 break;
3502
3503 default:
3504 error ("%s: invalid class name", fcn);
3505 break;
3506 }
3507 }
3508 }
3509
3510 return retval;
3511}
3512
3513static octave_value
3514fill_matrix (const octave_value_list& args, double val, float fval,
3515 const char *fcn)
3516{
3517 octave_value retval;
3518
3519 int nargin = args.length ();
3520
3521 oct_data_conv::data_type dt = oct_data_conv::dt_double;
3522
3523 dim_vector dims (1, 1);
3524
3525 if (nargin > 0 && args(nargin-1).is_string ())
3526 {
3527 std::string nm = args(nargin-1).string_value ();
3528 nargin--;
3529
3530 dt = oct_data_conv::string_to_data_type (nm);
3531
3532 if (error_state)
3533 return retval;
3534 }
3535
3536 switch (nargin)
3537 {
3538 case 0:
3539 break;
3540
3541 case 1:
3542 get_dimensions (args(0), fcn, dims);
3543 break;
3544
3545 default:
3546 {
3547 dims.resize (nargin);
3548
3549 for (int i = 0; i < nargin; i++)
3550 {
3551 dims(i) = args(i).is_empty () ? 0 : args(i).idx_type_value ();
3552
3553 if (error_state)
3554 {
3555 error ("%s: expecting scalar integer arguments", fcn);
3556 break;
3557 }
3558 }
3559 }
3560 break;
3561 }
3562
3563 if (! error_state)
3564 {
3565 dims.chop_trailing_singletons ();
3566
3567 check_dimensions (dims, fcn);
3568
3569 // Note that automatic narrowing will handle conversion from
3570 // NDArray to scalar.
3571
3572 if (! error_state)
3573 {
3574 switch (dt)
3575 {
3576 case oct_data_conv::dt_single:
3577 retval = FloatNDArray (dims, fval);
3578 break;
3579
3580 case oct_data_conv::dt_double:
3581 retval = NDArray (dims, val);
3582 break;
3583
3584 default:
3585 error ("%s: invalid class name", fcn);
3586 break;
3587 }
3588 }
3589 }
3590
3591 return retval;
3592}
3593
3594static octave_value
3595fill_matrix (const octave_value_list& args, double val, const char *fcn)
3596{
3597 octave_value retval;
3598
3599 int nargin = args.length ();
3600
3601 oct_data_conv::data_type dt = oct_data_conv::dt_double;
3602
3603 dim_vector dims (1, 1);
3604
3605 if (nargin > 0 && args(nargin-1).is_string ())
3606 {
3607 std::string nm = args(nargin-1).string_value ();
3608 nargin--;
3609
3610 dt = oct_data_conv::string_to_data_type (nm);
3611
3612 if (error_state)
3613 return retval;
3614 }
3615
3616 switch (nargin)
3617 {
3618 case 0:
3619 break;
3620
3621 case 1:
3622 get_dimensions (args(0), fcn, dims);
3623 break;
3624
3625 default:
3626 {
3627 dims.resize (nargin);
3628
3629 for (int i = 0; i < nargin; i++)
3630 {
3631 dims(i) = args(i).is_empty () ? 0 : args(i).idx_type_value ();
3632
3633 if (error_state)
3634 {
3635 error ("%s: expecting scalar integer arguments", fcn);
3636 break;
3637 }
3638 }
3639 }
3640 break;
3641 }
3642
3643 if (! error_state)
3644 {
3645 dims.chop_trailing_singletons ();
3646
3647 check_dimensions (dims, fcn);
3648
3649 // Note that automatic narrowing will handle conversion from
3650 // NDArray to scalar.
3651
3652 if (! error_state)
3653 {
3654 switch (dt)
3655 {
3656 case oct_data_conv::dt_single:
3657 retval = FloatNDArray (dims, static_cast <float> (val));
3658 break;
3659
3660 case oct_data_conv::dt_double:
3661 retval = NDArray (dims, val);
3662 break;
3663
3664 default:
3665 error ("%s: invalid class name", fcn);
3666 break;
3667 }
3668 }
3669 }
3670
3671 return retval;
3672}
3673
3674static octave_value
3675fill_matrix (const octave_value_list& args, const Complex& val,
3676 const char *fcn)
3677{
3678 octave_value retval;
3679
3680 int nargin = args.length ();
3681
3682 oct_data_conv::data_type dt = oct_data_conv::dt_double;
3683
3684 dim_vector dims (1, 1);
3685
3686 if (nargin > 0 && args(nargin-1).is_string ())
3687 {
3688 std::string nm = args(nargin-1).string_value ();
3689 nargin--;
3690
3691 dt = oct_data_conv::string_to_data_type (nm);
3692
3693 if (error_state)
3694 return retval;
3695 }
3696
3697 switch (nargin)
3698 {
3699 case 0:
3700 break;
3701
3702 case 1:
3703 get_dimensions (args(0), fcn, dims);
3704 break;
3705
3706 default:
3707 {
3708 dims.resize (nargin);
3709
3710 for (int i = 0; i < nargin; i++)
3711 {
3712 dims(i) = args(i).is_empty () ? 0 : args(i).idx_type_value ();
3713
3714 if (error_state)
3715 {
3716 error ("%s: expecting scalar integer arguments", fcn);
3717 break;
3718 }
3719 }
3720 }
3721 break;
3722 }
3723
3724 if (! error_state)
3725 {
3726 dims.chop_trailing_singletons ();
3727
3728 check_dimensions (dims, fcn);
3729
3730 // Note that automatic narrowing will handle conversion from
3731 // NDArray to scalar.
3732
3733 if (! error_state)
3734 {
3735 switch (dt)
3736 {
3737 case oct_data_conv::dt_single:
3738 retval = FloatComplexNDArray (dims, static_cast<FloatComplex> (val));
3739 break;
3740
3741 case oct_data_conv::dt_double:
3742 retval = ComplexNDArray (dims, val);
3743 break;
3744
3745 default:
3746 error ("%s: invalid class name", fcn);
3747 break;
3748 }
3749 }
3750 }
3751
3752 return retval;
3753}
3754
3755static octave_value
3756fill_matrix (const octave_value_list& args, bool val, const char *fcn)
3757{
3758 octave_value retval;
3759
3760 int nargin = args.length ();
3761
3762 dim_vector dims (1, 1);
3763
3764 switch (nargin)
3765 {
3766 case 0:
3767 break;
3768
3769 case 1:
3770 get_dimensions (args(0), fcn, dims);
3771 break;
3772
3773 default:
3774 {
3775 dims.resize (nargin);
3776
3777 for (int i = 0; i < nargin; i++)
3778 {
3779 dims(i) = args(i).is_empty () ? 0 : args(i).idx_type_value ();
3780
3781 if (error_state)
3782 {
3783 error ("%s: expecting scalar integer arguments", fcn);
3784 break;
3785 }
3786 }
3787 }
3788 break;
3789 }
3790
3791 if (! error_state)
3792 {
3793 dims.chop_trailing_singletons ();
3794
3795 check_dimensions (dims, fcn);
3796
3797 // Note that automatic narrowing will handle conversion from
3798 // NDArray to scalar.
3799
3800 if (! error_state)
3801 retval = boolNDArray (dims, val);
3802 }
3803
3804 return retval;
3805}
3806
3807DEFUN (ones, args, ,
3808 "-*- texinfo -*-\n\
3809@deftypefn {Built-in Function} {} ones (@var{x})\n\
3810@deftypefnx {Built-in Function} {} ones (@var{n}, @var{m})\n\
3811@deftypefnx {Built-in Function} {} ones (@var{n}, @var{m}, @var{k}, @dots{})\n\
3812@deftypefnx {Built-in Function} {} ones (@dots{}, @var{class})\n\
3813Return a matrix or N-dimensional array whose elements are all 1.\n\
3814If invoked with a single scalar integer argument, return a square\n\
3815matrix of the specified size. If invoked with two or more scalar\n\
3816integer arguments, or a vector of integer values, return an array with\n\
3817given dimensions.\n\
3818\n\
3819If you need to create a matrix whose values are all the same, you should\n\
3820use an expression like\n\
3821\n\
3822@example\n\
3823val_matrix = val * ones (n, m)\n\
3824@end example\n\
3825\n\
3826The optional argument @var{class}, allows @code{ones} to return an array of\n\
3827the specified type, for example\n\
3828\n\
3829@example\n\
3830val = ones (n,m, \"uint8\")\n\
3831@end example\n\
3832@end deftypefn")
3833{
3834 return fill_matrix (args, 1, "ones");
3835}
3836
3837/*
3838
3839%!assert(ones (3), [1, 1, 1; 1, 1, 1; 1, 1, 1]);
3840%!assert(ones (2, 3), [1, 1, 1; 1, 1, 1]);
3841%!assert(ones (3, 2), [1, 1; 1, 1; 1, 1]);
3842%!assert(size (ones (3, 4, 5)), [3, 4, 5]);
3843
3844%!assert(ones (3,'single'), single([1, 1, 1; 1, 1, 1; 1, 1, 1]));
3845%!assert(ones (2, 3,'single'), single([1, 1, 1; 1, 1, 1]));
3846%!assert(ones (3, 2,'single'), single([1, 1; 1, 1; 1, 1]));
3847%!assert(size (ones (3, 4, 5, 'single')), [3, 4, 5]);
3848
3849%!assert(ones (3,'int8'), int8([1, 1, 1; 1, 1, 1; 1, 1, 1]));
3850%!assert(ones (2, 3,'int8'), int8([1, 1, 1; 1, 1, 1]));
3851%!assert(ones (3, 2,'int8'), int8([1, 1; 1, 1; 1, 1]));
3852%!assert(size (ones (3, 4, 5, 'int8')), [3, 4, 5]);
3853
3854 */
3855
3856DEFUN (zeros, args, ,
3857 "-*- texinfo -*-\n\
3858@deftypefn {Built-in Function} {} zeros (@var{x})\n\
3859@deftypefnx {Built-in Function} {} zeros (@var{n}, @var{m})\n\
3860@deftypefnx {Built-in Function} {} zeros (@var{n}, @var{m}, @var{k}, @dots{})\n\
3861@deftypefnx {Built-in Function} {} zeros (@dots{}, @var{class})\n\
3862Return a matrix or N-dimensional array whose elements are all 0.\n\
3863The arguments are handled the same as the arguments for @code{ones}.\n\
3864\n\
3865The optional argument @var{class}, allows @code{zeros} to return an array of\n\
3866the specified type, for example\n\
3867\n\
3868@example\n\
3869val = zeros (n,m, \"uint8\")\n\
3870@end example\n\
3871@end deftypefn")
3872{
3873 return fill_matrix (args, 0, "zeros");
3874}
3875
3876/*
3877
3878%!assert(zeros (3), [0, 0, 0; 0, 0, 0; 0, 0, 0]);
3879%!assert(zeros (2, 3), [0, 0, 0; 0, 0, 0]);
3880%!assert(zeros (3, 2), [0, 0; 0, 0; 0, 0]);
3881%!assert(size (zeros (3, 4, 5)), [3, 4, 5]);
3882
3883%!assert(zeros (3,'single'), single([0, 0, 0; 0, 0, 0; 0, 0, 0]));
3884%!assert(zeros (2, 3,'single'), single([0, 0, 0; 0, 0, 0]));
3885%!assert(zeros (3, 2,'single'), single([0, 0; 0, 0; 0, 0]));
3886%!assert(size (zeros (3, 4, 5, 'single')), [3, 4, 5]);
3887
3888%!assert(zeros (3,'int8'), int8([0, 0, 0; 0, 0, 0; 0, 0, 0]));
3889%!assert(zeros (2, 3,'int8'), int8([0, 0, 0; 0, 0, 0]));
3890%!assert(zeros (3, 2,'int8'), int8([0, 0; 0, 0; 0, 0]));
3891%!assert(size (zeros (3, 4, 5, 'int8')), [3, 4, 5]);
3892
3893 */
3894
3895DEFUN (Inf, args, ,
3896 "-*- texinfo -*-\n\
3897@deftypefn {Built-in Function} {} Inf\n\
3898@deftypefnx {Built-in Function} {} Inf (@var{n})\n\
3899@deftypefnx {Built-in Function} {} Inf (@var{n}, @var{m})\n\
3900@deftypefnx {Built-in Function} {} Inf (@var{n}, @var{m}, @var{k}, @dots{})\n\
3901@deftypefnx {Built-in Function} {} Inf (@dots{}, @var{class})\n\
3902Return a scalar, matrix or N-dimensional array whose elements are all equal\n\
3903to the IEEE representation for positive infinity.\n\
3904\n\
3905Infinity is produced when results are too large to be represented using the\n\
3906the IEEE floating point format for numbers. Two common examples which\n\
3907produce infinity are division by zero and overflow.\n\
3908@example\n\
3909@group\n\
3910[1/0 e^800]\n\
3911@result{}\n\
3912Inf Inf\n\
3913@end group\n\
3914@end example\n\
3915\n\
3916When called with no arguments, return a scalar with the value @samp{Inf}.\n\
3917When called with a single argument, return a square matrix with the dimension\n\
3918specified. When called with more than one scalar argument the first two\n\
3919arguments are taken as the number of rows and columns and any further\n\
3920arguments specify additional matrix dimensions.\n\
3921The optional argument @var{class} specifies the return type and may be\n\
3922either \"double\" or \"single\".\n\
3923@seealso{isinf}\n\
3924@end deftypefn")
3925{
3926 return fill_matrix (args, lo_ieee_inf_value (),
3927 lo_ieee_float_inf_value (), "Inf");
3928}
3929
3930DEFALIAS (inf, Inf);
3931
3932/*
3933
3934%!assert(inf (3), [Inf, Inf, Inf; Inf, Inf, Inf; Inf, Inf, Inf]);
3935%!assert(inf (2, 3), [Inf, Inf, Inf; Inf, Inf, Inf]);
3936%!assert(inf (3, 2), [Inf, Inf; Inf, Inf; Inf, Inf]);
3937%!assert(size (inf (3, 4, 5)), [3, 4, 5]);
3938
3939%!assert(inf (3,'single'), single([Inf, Inf, Inf; Inf, Inf, Inf; Inf, Inf, Inf]));
3940%!assert(inf (2, 3,'single'), single([Inf, Inf, Inf; Inf, Inf, Inf]));
3941%!assert(inf (3, 2,'single'), single([Inf, Inf; Inf, Inf; Inf, Inf]));
3942%!assert(size (inf (3, 4, 5, 'single')), [3, 4, 5]);
3943
3944%!error(inf (3,'int8'));
3945%!error(inf (2, 3,'int8'));
3946%!error(inf (3, 2,'int8'));
3947%!error(inf (3, 4, 5, 'int8'));
3948
3949 */
3950
3951DEFUN (NaN, args, ,
3952 "-*- texinfo -*-\n\
3953@deftypefn {Built-in Function} {} NaN\n\
3954@deftypefnx {Built-in Function} {} NaN (@var{n})\n\
3955@deftypefnx {Built-in Function} {} NaN (@var{n}, @var{m})\n\
3956@deftypefnx {Built-in Function} {} NaN (@var{n}, @var{m}, @var{k}, @dots{})\n\
3957@deftypefnx {Built-in Function} {} NaN (@dots{}, @var{class})\n\
3958Return a scalar, matrix, or N-dimensional array whose elements are all equal\n\
3959to the IEEE symbol NaN (Not a Number).\n\
3960NaN is the result of operations which do not produce a well defined numerical\n\
3961result. Common operations which produce a NaN are arithmetic with infinity\n\
3962@tex\n\
3963($\\infty - \\infty$), zero divided by zero ($0/0$),\n\
3964@end tex\n\
3965@ifnottex\n\
3966(Inf - Inf), zero divided by zero (0/0),\n\
3967@end ifnottex\n\
3968and any operation involving another NaN value (5 + NaN).\n\
3969\n\
3970Note that NaN always compares not equal to NaN (NaN != NaN). This behavior\n\
3971is specified by the IEEE standard for floating point arithmetic. To\n\
3972find NaN values, use the @code{isnan} function.\n\
3973\n\
3974When called with no arguments, return a scalar with the value @samp{NaN}.\n\
3975When called with a single argument, return a square matrix with the dimension\n\
3976specified. When called with more than one scalar argument the first two\n\
3977arguments are taken as the number of rows and columns and any further\n\
3978arguments specify additional matrix dimensions.\n\
3979The optional argument @var{class} specifies the return type and may be\n\
3980either \"double\" or \"single\".\n\
3981@seealso{isnan}\n\
3982@end deftypefn")
3983{
3984 return fill_matrix (args, lo_ieee_nan_value (),
3985 lo_ieee_float_nan_value (), "NaN");
3986}
3987
3988DEFALIAS (nan, NaN);
3989
3990/*
3991%!assert(NaN (3), [NaN, NaN, NaN; NaN, NaN, NaN; NaN, NaN, NaN]);
3992%!assert(NaN (2, 3), [NaN, NaN, NaN; NaN, NaN, NaN]);
3993%!assert(NaN (3, 2), [NaN, NaN; NaN, NaN; NaN, NaN]);
3994%!assert(size (NaN (3, 4, 5)), [3, 4, 5]);
3995
3996%!assert(NaN (3,'single'), single([NaN, NaN, NaN; NaN, NaN, NaN; NaN, NaN, NaN]));
3997%!assert(NaN (2, 3,'single'), single([NaN, NaN, NaN; NaN, NaN, NaN]));
3998%!assert(NaN (3, 2,'single'), single([NaN, NaN; NaN, NaN; NaN, NaN]));
3999%!assert(size (NaN (3, 4, 5, 'single')), [3, 4, 5]);
4000
4001%!error(NaN (3,'int8'));
4002%!error(NaN (2, 3,'int8'));
4003%!error(NaN (3, 2,'int8'));
4004%!error(NaN (3, 4, 5, 'int8'));
4005
4006 */
4007
4008DEFUN (e, args, ,
4009 "-*- texinfo -*-\n\
4010@deftypefn {Built-in Function} {} e\n\
4011@deftypefnx {Built-in Function} {} e (@var{n})\n\
4012@deftypefnx {Built-in Function} {} e (@var{n}, @var{m})\n\
4013@deftypefnx {Built-in Function} {} e (@var{n}, @var{m}, @var{k}, @dots{})\n\
4014@deftypefnx {Built-in Function} {} e (@dots{}, @var{class})\n\
4015Return a scalar, matrix, or N-dimensional array whose elements are all equal\n\
4016to the base of natural logarithms. The constant\n\
4017@tex\n\
4018$e$ satisfies the equation $\\log (e) = 1$.\n\
4019@end tex\n\
4020@ifnottex\n\
4021@samp{e} satisfies the equation @code{log} (e) = 1.\n\
4022@end ifnottex\n\
4023\n\
4024When called with no arguments, return a scalar with the value @math{e}. When\n\
4025called with a single argument, return a square matrix with the dimension\n\
4026specified. When called with more than one scalar argument the first two\n\
4027arguments are taken as the number of rows and columns and any further\n\
4028arguments specify additional matrix dimensions.\n\
4029The optional argument @var{class} specifies the return type and may be\n\
4030either \"double\" or \"single\".\n\
4031@end deftypefn")
4032{
4033#if defined (M_E)
4034 double e_val = M_E;
4035#else
4036 double e_val = exp (1.0);
4037#endif
4038
4039 return fill_matrix (args, e_val, "e");
4040}
4041
4042DEFUN (eps, args, ,
4043 "-*- texinfo -*-\n\
4044@deftypefn {Built-in Function} {} eps\n\
4045@deftypefnx {Built-in Function} {} eps (@var{x})\n\
4046@deftypefnx {Built-in Function} {} eps (@var{n}, @var{m})\n\
4047@deftypefnx {Built-in Function} {} eps (@var{n}, @var{m}, @var{k}, @dots{})\n\
4048@deftypefnx {Built-in Function} {} eps (@dots{}, @var{class})\n\
4049Return a scalar, matrix or N-dimensional array whose elements are all eps,\n\
4050the machine precision. More precisely, @code{eps} is the relative spacing\n\
4051between any two adjacent numbers in the machine's floating point system.\n\
4052This number is obviously system dependent. On machines that support IEEE\n\
4053floating point arithmetic, @code{eps} is approximately\n\
4054@tex\n\
4055$2.2204\\times10^{-16}$ for double precision and $1.1921\\times10^{-7}$\n\
4056@end tex\n\
4057@ifnottex\n\
40582.2204e-16 for double precision and 1.1921e-07\n\
4059@end ifnottex\n\
4060for single precision.\n\
4061\n\
4062When called with no arguments, return a scalar with the value\n\
4063@code{eps(1.0)}.\n\
4064Given a single argument @var{x}, return the distance between @var{x} and\n\
4065the next largest value.\n\
4066When called with more than one argument the first two arguments are taken as\n\
4067the number of rows and columns and any further\n\
4068arguments specify additional matrix dimensions.\n\
4069The optional argument @var{class} specifies the return type and may be\n\
4070either \"double\" or \"single\".\n\
4071@end deftypefn")
4072{
4073 int nargin = args.length ();
4074 octave_value retval;
4075
4076 if (nargin == 1 && ! args(0).is_string ())
4077 {
4078 if (args(0).is_single_type ())
4079 {
4080 float val = args(0).float_value ();
4081
4082 if (! error_state)
4083 {
4084 val = ::fabsf(val);
4085 if (xisnan (val) || xisinf (val))
4086 retval = fill_matrix (octave_value ("single"),
4087 lo_ieee_nan_value (),
4088 lo_ieee_float_nan_value (), "eps");
4089 else if (val < FLT_MIN)
4090 retval = fill_matrix (octave_value ("single"), 0e0,
4091 powf (2.0, -149e0), "eps");
4092 else
4093 {
4094 int expon;
4095 frexpf (val, &expon);
4096 val = std::pow (static_cast <float> (2.0),
4097 static_cast <float> (expon - 24));
4098 retval = fill_matrix (octave_value ("single"), DBL_EPSILON,
4099 val, "eps");
4100 }
4101 }
4102 }
4103 else
4104 {
4105 double val = args(0).double_value ();
4106
4107 if (! error_state)
4108 {
4109 val = ::fabs(val);
4110 if (xisnan (val) || xisinf (val))
4111 retval = fill_matrix (octave_value_list (),
4112 lo_ieee_nan_value (),
4113 lo_ieee_float_nan_value (), "eps");
4114 else if (val < DBL_MIN)
4115 retval = fill_matrix (octave_value_list (),
4116 pow (2.0, -1074e0), 0e0, "eps");
4117 else
4118 {
4119 int expon;
4120 frexp (val, &expon);
4121 val = std::pow (static_cast <double> (2.0),
4122 static_cast <double> (expon - 53));
4123 retval = fill_matrix (octave_value_list (), val,
4124 FLT_EPSILON, "eps");
4125 }
4126 }
4127 }
4128 }
4129 else
4130 retval = fill_matrix (args, DBL_EPSILON, FLT_EPSILON, "eps");
4131
4132 return retval;
4133}
4134
4135/*
4136
4137%!assert(eps(1/2),2^(-53))
4138%!assert(eps(1),2^(-52))
4139%!assert(eps(2),2^(-51))
4140%!assert(eps(realmax),2^971)
4141%!assert(eps(0),2^(-1074))
4142%!assert(eps(realmin/2),2^(-1074))
4143%!assert(eps(realmin/16),2^(-1074))
4144%!assert(eps(Inf),NaN)
4145%!assert(eps(NaN),NaN)
4146%!assert(eps(single(1/2)),single(2^(-24)))
4147%!assert(eps(single(1)),single(2^(-23)))
4148%!assert(eps(single(2)),single(2^(-22)))
4149%!assert(eps(realmax('single')),single(2^104))
4150%!assert(eps(single(0)),single(2^(-149)))
4151%!assert(eps(realmin('single')/2),single(2^(-149)))
4152%!assert(eps(realmin('single')/16),single(2^(-149)))
4153%!assert(eps(single(Inf)),single(NaN))
4154%!assert(eps(single(NaN)),single(NaN))
4155
4156*/
4157
4158
4159DEFUN (pi, args, ,
4160 "-*- texinfo -*-\n\
4161@deftypefn {Built-in Function} {} pi\n\
4162@deftypefnx {Built-in Function} {} pi (@var{n})\n\
4163@deftypefnx {Built-in Function} {} pi (@var{n}, @var{m})\n\
4164@deftypefnx {Built-in Function} {} pi (@var{n}, @var{m}, @var{k}, @dots{})\n\
4165@deftypefnx {Built-in Function} {} pi (@dots{}, @var{class})\n\
4166Return a scalar, matrix, or N-dimensional array whose elements are all equal\n\
4167to the ratio of the circumference of a circle to its\n\
4168@tex\n\
4169diameter($\\pi$).\n\
4170@end tex\n\
4171@ifnottex\n\
4172diameter.\n\
4173@end ifnottex\n\
4174Internally, @code{pi} is computed as @samp{4.0 * atan (1.0)}.\n\
4175\n\
4176When called with no arguments, return a scalar with the value of\n\
4177@tex\n\
4178$\\pi$.\n\
4179@end tex\n\
4180@ifnottex\n\
4181pi.\n\
4182@end ifnottex\n\
4183When called with a single argument, return a square matrix with the dimension\n\
4184specified. When called with more than one scalar argument the first two\n\
4185arguments are taken as the number of rows and columns and any further\n\
4186arguments specify additional matrix dimensions.\n\
4187The optional argument @var{class} specifies the return type and may be\n\
4188either \"double\" or \"single\".\n\
4189@end deftypefn")
4190{
4191#if defined (M_PI)
4192 double pi_val = M_PI;
4193#else
4194 double pi_val = 4.0 * atan (1.0);
4195#endif
4196
4197 return fill_matrix (args, pi_val, "pi");
4198}
4199
4200DEFUN (realmax, args, ,
4201 "-*- texinfo -*-\n\
4202@deftypefn {Built-in Function} {} realmax\n\
4203@deftypefnx {Built-in Function} {} realmax (@var{n})\n\
4204@deftypefnx {Built-in Function} {} realmax (@var{n}, @var{m})\n\
4205@deftypefnx {Built-in Function} {} realmax (@var{n}, @var{m}, @var{k}, @dots{})\n\
4206@deftypefnx {Built-in Function} {} realmax (@dots{}, @var{class})\n\
4207Return a scalar, matrix or N-dimensional array whose elements are all equal\n\
4208to the largest floating point number that is representable. The actual\n\
4209value is system dependent. On machines that support IEEE\n\
4210floating point arithmetic, @code{realmax} is approximately\n\
4211@tex\n\
4212$1.7977\\times10^{308}$ for double precision and $3.4028\\times10^{38}$\n\
4213@end tex\n\
4214@ifnottex\n\
42151.7977e+308 for double precision and 3.4028e+38\n\
4216@end ifnottex\n\
4217for single precision.\n\
4218\n\
4219When called with no arguments, return a scalar with the value\n\
4220@code{realmax(\"double\")}.\n\
4221When called with a single argument, return a square matrix with the dimension\n\
4222specified. When called with more than one scalar argument the first two\n\
4223arguments are taken as the number of rows and columns and any further\n\
4224arguments specify additional matrix dimensions.\n\
4225The optional argument @var{class} specifies the return type and may be\n\
4226either \"double\" or \"single\".\n\
4227@seealso{realmin, intmax, bitmax}\n\
4228@end deftypefn")
4229{
4230 return fill_matrix (args, DBL_MAX, FLT_MAX, "realmax");
4231}
4232
4233DEFUN (realmin, args, ,
4234 "-*- texinfo -*-\n\
4235@deftypefn {Built-in Function} {} realmin\n\
4236@deftypefnx {Built-in Function} {} realmin (@var{n})\n\
4237@deftypefnx {Built-in Function} {} realmin (@var{n}, @var{m})\n\
4238@deftypefnx {Built-in Function} {} realmin (@var{n}, @var{m}, @var{k}, @dots{})\n\
4239@deftypefnx {Built-in Function} {} realmin (@dots{}, @var{class})\n\
4240Return a scalar, matrix or N-dimensional array whose elements are all equal\n\
4241to the smallest normalized floating point number that is representable.\n\
4242The actual value is system dependent. On machines that support\n\
4243IEEE floating point arithmetic, @code{realmin} is approximately\n\
4244@tex\n\
4245$2.2251\\times10^{-308}$ for double precision and $1.1755\\times10^{-38}$\n\
4246@end tex\n\
4247@ifnottex\n\
42482.2251e-308 for double precision and 1.1755e-38\n\
4249@end ifnottex\n\
4250for single precision.\n\
4251\n\
4252When called with no arguments, return a scalar with the value\n\
4253@code{realmin(\"double\")}.\n\
4254When called with a single argument, return a square matrix with the dimension\n\
4255specified. When called with more than one scalar argument the first two\n\
4256arguments are taken as the number of rows and columns and any further\n\
4257arguments specify additional matrix dimensions.\n\
4258The optional argument @var{class} specifies the return type and may be\n\
4259either \"double\" or \"single\".\n\
4260@seealso{realmax, intmin}\n\
4261@end deftypefn")
4262{
4263 return fill_matrix (args, DBL_MIN, FLT_MIN, "realmin");
4264}
4265
4266DEFUN (I, args, ,
4267 "-*- texinfo -*-\n\
4268@deftypefn {Built-in Function} {} I\n\
4269@deftypefnx {Built-in Function} {} I (@var{n})\n\
4270@deftypefnx {Built-in Function} {} I (@var{n}, @var{m})\n\
4271@deftypefnx {Built-in Function} {} I (@var{n}, @var{m}, @var{k}, @dots{})\n\
4272@deftypefnx {Built-in Function} {} I (@dots{}, @var{class})\n\
4273Return a scalar, matrix, or N-dimensional array whose elements are all equal\n\
4274to the pure imaginary unit, defined as\n\
4275@tex\n\
4276$\\sqrt{-1}$.\n\
4277@end tex\n\
4278@ifnottex\n\
4279@code{sqrt (-1)}.\n\
4280@end ifnottex\n\
4281 I, and its equivalents i, J, and j, are functions so any of the names may\n\
4282be reused for other purposes (such as i for a counter variable).\n\
4283\n\
4284When called with no arguments, return a scalar with the value @math{i}. When\n\
4285called with a single argument, return a square matrix with the dimension\n\
4286specified. When called with more than one scalar argument the first two\n\
4287arguments are taken as the number of rows and columns and any further\n\
4288arguments specify additional matrix dimensions.\n\
4289The optional argument @var{class} specifies the return type and may be\n\
4290either \"double\" or \"single\".\n\
4291@end deftypefn")
4292{
4293 return fill_matrix (args, Complex (0.0, 1.0), "I");
4294}
4295
4296DEFALIAS (i, I);
4297DEFALIAS (J, I);
4298DEFALIAS (j, I);
4299
4300DEFUN (NA, args, ,
4301 "-*- texinfo -*-\n\
4302@deftypefn {Built-in Function} {} NA\n\
4303@deftypefnx {Built-in Function} {} NA (@var{n})\n\
4304@deftypefnx {Built-in Function} {} NA (@var{n}, @var{m})\n\
4305@deftypefnx {Built-in Function} {} NA (@var{n}, @var{m}, @var{k}, @dots{})\n\
4306@deftypefnx {Built-in Function} {} NA (@dots{}, @var{class})\n\
4307Return a scalar, matrix, or N-dimensional array whose elements are all equal\n\
4308to the special constant used to designate missing values.\n\
4309\n\
4310Note that NA always compares not equal to NA (NA != NA).\n\
4311To find NA values, use the @code{isna} function.\n\
4312\n\
4313When called with no arguments, return a scalar with the value @samp{NA}.\n\
4314When called with a single argument, return a square matrix with the dimension\n\
4315specified. When called with more than one scalar argument the first two\n\
4316arguments are taken as the number of rows and columns and any further\n\
4317arguments specify additional matrix dimensions.\n\
4318The optional argument @var{class} specifies the return type and may be\n\
4319either \"double\" or \"single\".\n\
4320@seealso{isna}\n\
4321@end deftypefn")
4322{
4323 return fill_matrix (args, lo_ieee_na_value (),
4324 lo_ieee_float_na_value (), "NA");
4325}
4326
4327/*
4328
4329%!assert(single(NA('double')),NA('single'))
4330%!assert(double(NA('single')),NA('double'))
4331
4332 */
4333
4334DEFUN (false, args, ,
4335 "-*- texinfo -*-\n\
4336@deftypefn {Built-in Function} {} false (@var{x})\n\
4337@deftypefnx {Built-in Function} {} false (@var{n}, @var{m})\n\
4338@deftypefnx {Built-in Function} {} false (@var{n}, @var{m}, @var{k}, @dots{})\n\
4339Return a matrix or N-dimensional array whose elements are all logical 0.\n\
4340The arguments are handled the same as the arguments for @code{ones}.\n\
4341@end deftypefn")
4342{
4343 return fill_matrix (args, false, "false");
4344}
4345
4346DEFUN (true, args, ,
4347 "-*- texinfo -*-\n\
4348@deftypefn {Built-in Function} {} true (@var{x})\n\
4349@deftypefnx {Built-in Function} {} true (@var{n}, @var{m})\n\
4350@deftypefnx {Built-in Function} {} true (@var{n}, @var{m}, @var{k}, @dots{})\n\
4351Return a matrix or N-dimensional array whose elements are all logical 1.\n\
4352The arguments are handled the same as the arguments for @code{ones}.\n\
4353@end deftypefn")
4354{
4355 return fill_matrix (args, true, "true");
4356}
4357
4358template <class MT>
4359octave_value
4360identity_matrix (int nr, int nc)
4361{
4362 octave_value retval;
4363
4364 typename MT::element_type one (1);
4365
4366 if (nr == 1 && nc == 1)
4367 retval = one;
4368 else
4369 {
4370 dim_vector dims (nr, nc);
4371
4372 typename MT::element_type zero (0);
4373
4374 MT m (dims, zero);
4375
4376 if (nr > 0 && nc > 0)
4377 {
4378 int n = std::min (nr, nc);
4379
4380 for (int i = 0; i < n; i++)
4381 m(i,i) = one;
4382 }
4383
4384 retval = m;
4385 }
4386
4387 return retval;
4388}
4389
4390#define INSTANTIATE_EYE(T) \
4391 template octave_value identity_matrix<T> (int, int)
4392
4393INSTANTIATE_EYE (int8NDArray);
4394INSTANTIATE_EYE (uint8NDArray);
4395INSTANTIATE_EYE (int16NDArray);
4396INSTANTIATE_EYE (uint16NDArray);
4397INSTANTIATE_EYE (int32NDArray);
4398INSTANTIATE_EYE (uint32NDArray);
4399INSTANTIATE_EYE (int64NDArray);
4400INSTANTIATE_EYE (uint64NDArray);
4401INSTANTIATE_EYE (FloatNDArray);
4402INSTANTIATE_EYE (NDArray);
4403INSTANTIATE_EYE (boolNDArray);
4404
4405static octave_value
4406identity_matrix (int nr, int nc, oct_data_conv::data_type dt)
4407{
4408 octave_value retval;
4409
4410 // FIXME -- perhaps this should be made extensible by using
4411 // the class name to lookup a function to call to create the new
4412 // value.
4413
4414 if (! error_state)
4415 {
4416 switch (dt)
4417 {
4418 case oct_data_conv::dt_int8:
4419 retval = identity_matrix<int8NDArray> (nr, nc);
4420 break;
4421
4422 case oct_data_conv::dt_uint8:
4423 retval = identity_matrix<uint8NDArray> (nr, nc);
4424 break;
4425
4426 case oct_data_conv::dt_int16:
4427 retval = identity_matrix<int16NDArray> (nr, nc);
4428 break;
4429
4430 case oct_data_conv::dt_uint16:
4431 retval = identity_matrix<uint16NDArray> (nr, nc);
4432 break;
4433
4434 case oct_data_conv::dt_int32:
4435 retval = identity_matrix<int32NDArray> (nr, nc);
4436 break;
4437
4438 case oct_data_conv::dt_uint32:
4439 retval = identity_matrix<uint32NDArray> (nr, nc);
4440 break;
4441
4442 case oct_data_conv::dt_int64:
4443 retval = identity_matrix<int64NDArray> (nr, nc);
4444 break;
4445
4446 case oct_data_conv::dt_uint64:
4447 retval = identity_matrix<uint64NDArray> (nr, nc);
4448 break;
4449
4450 case oct_data_conv::dt_single:
4451 retval = FloatDiagMatrix (nr, nc, 1.0f);
4452 break;
4453
4454 case oct_data_conv::dt_double:
4455 retval = DiagMatrix (nr, nc, 1.0);
4456 break;
4457
4458 case oct_data_conv::dt_logical:
4459 retval = identity_matrix<boolNDArray> (nr, nc);
4460 break;
4461
4462 default:
4463 error ("eye: invalid class name");
4464 break;
4465 }
4466 }
4467
4468 return retval;
4469}
4470
4471#undef INT_EYE_MATRIX
4472
4473DEFUN (eye, args, ,
4474 "-*- texinfo -*-\n\
4475@deftypefn {Built-in Function} {} eye (@var{x})\n\
4476@deftypefnx {Built-in Function} {} eye (@var{n}, @var{m})\n\
4477@deftypefnx {Built-in Function} {} eye (@dots{}, @var{class})\n\
4478Return an identity matrix. If invoked with a single scalar argument,\n\
4479@code{eye} returns a square matrix with the dimension specified. If you\n\
4480supply two scalar arguments, @code{eye} takes them to be the number of\n\
4481rows and columns. If given a vector with two elements, @code{eye} uses\n\
4482the values of the elements as the number of rows and columns,\n\
4483respectively. For example,\n\
4484\n\
4485@example\n\
4486@group\n\
4487eye (3)\n\
4488 @result{} 1 0 0\n\
4489 0 1 0\n\
4490 0 0 1\n\
4491@end group\n\
4492@end example\n\
4493\n\
4494The following expressions all produce the same result:\n\
4495\n\
4496@example\n\
4497@group\n\
4498eye (2)\n\
4499@equiv{}\n\
4500eye (2, 2)\n\
4501@equiv{}\n\
4502eye (size ([1, 2; 3, 4])\n\
4503@end group\n\
4504@end example\n\
4505\n\
4506The optional argument @var{class}, allows @code{eye} to return an array of\n\
4507the specified type, like\n\
4508\n\
4509@example\n\
4510val = zeros (n,m, \"uint8\")\n\
4511@end example\n\
4512\n\
4513Calling @code{eye} with no arguments is equivalent to calling it\n\
4514with an argument of 1. This odd definition is for compatibility\n\
4515with @sc{matlab}.\n\
4516@end deftypefn")
4517{
4518 octave_value retval;
4519
4520 int nargin = args.length ();
4521
4522 oct_data_conv::data_type dt = oct_data_conv::dt_double;
4523
4524 // Check for type information.
4525
4526 if (nargin > 0 && args(nargin-1).is_string ())
4527 {
4528 std::string nm = args(nargin-1).string_value ();
4529 nargin--;
4530
4531 dt = oct_data_conv::string_to_data_type (nm);
4532
4533 if (error_state)
4534 return retval;
4535 }
4536
4537 switch (nargin)
4538 {
4539 case 0:
4540 retval = identity_matrix (1, 1, dt);
4541 break;
4542
4543 case 1:
4544 {
4545 octave_idx_type nr, nc;
4546 get_dimensions (args(0), "eye", nr, nc);
4547
4548 if (! error_state)
4549 retval = identity_matrix (nr, nc, dt);
4550 }
4551 break;
4552
4553 case 2:
4554 {
4555 octave_idx_type nr, nc;
4556 get_dimensions (args(0), args(1), "eye", nr, nc);
4557
4558 if (! error_state)
4559 retval = identity_matrix (nr, nc, dt);
4560 }
4561 break;
4562
4563 default:
4564 print_usage ();
4565 break;
4566 }
4567
4568 return retval;
4569}
4570
4571
4572/*
4573
4574%!assert (full (eye(3)), [1, 0, 0; 0, 1, 0; 0, 0, 1]);
4575%!assert (full (eye(2, 3)), [1, 0, 0; 0, 1, 0]);
4576
4577%!assert (full (eye(3,'single')), single([1, 0, 0; 0, 1, 0; 0, 0, 1]));
4578%!assert (full (eye(2, 3,'single')), single([1, 0, 0; 0, 1, 0]));
4579
4580%!assert (eye(3,'int8'), int8([1, 0, 0; 0, 1, 0; 0, 0, 1]));
4581%!assert (eye(2, 3,'int8'), int8([1, 0, 0; 0, 1, 0]));
4582
4583%!error <Invalid call to eye.*> eye (1, 2, 3);
4584
4585 */
4586
4587template <class MT>
4588static octave_value
4589do_linspace (const octave_value& base, const octave_value& limit,
4590 octave_idx_type n)
4591{
4592 typedef typename MT::column_vector_type CVT;
4593 typedef typename MT::element_type T;
4594
4595 octave_value retval;
4596
4597 if (base.is_scalar_type ())
4598 {
4599 T bs = octave_value_extract<T> (base);
4600 if (limit.is_scalar_type ())
4601 {
4602 T ls = octave_value_extract<T> (limit);
4603 retval = linspace (bs, ls, n);
4604 }
4605 else
4606 {
4607 CVT lv = octave_value_extract<CVT> (limit);
4608 CVT bv (lv.length (), bs);
4609 retval = linspace (bv, lv, n);
4610 }
4611 }
4612 else
4613 {
4614 CVT bv = octave_value_extract<CVT> (base);
4615 if (limit.is_scalar_type ())
4616 {
4617 T ls = octave_value_extract<T> (limit);
4618 CVT lv (bv.length (), ls);
4619 retval = linspace (bv, lv, n);
4620 }
4621 else
4622 {
4623 CVT lv = octave_value_extract<CVT> (limit);
4624 retval = linspace (bv, lv, n);
4625 }
4626 }
4627
4628 return retval;
4629}
4630
4631DEFUN (linspace, args, ,
4632 "-*- texinfo -*-\n\
4633@deftypefn {Built-in Function} {} linspace (@var{base}, @var{limit}, @var{n})\n\
4634Return a row vector with @var{n} linearly spaced elements between\n\
4635@var{base} and @var{limit}. If the number of elements is greater than one,\n\
4636then the @var{base} and @var{limit} are always included in\n\
4637the range. If @var{base} is greater than @var{limit}, the elements are\n\
4638stored in decreasing order. If the number of points is not specified, a\n\
4639value of 100 is used.\n\
4640\n\
4641The @code{linspace} function always returns a row vector if both\n\
4642@var{base} and @var{limit} are scalars. If one of them or both are column\n\
4643vectors, @code{linspace} returns a matrix.\n\
4644\n\
4645For compatibility with @sc{matlab}, return the second argument if\n\
4646fewer than two values are requested.\n\
4647@end deftypefn")
4648{
4649 octave_value retval;
4650
4651 int nargin = args.length ();
4652
4653 octave_idx_type npoints = 100;
4654
4655 if (nargin != 2 && nargin != 3)
4656 {
4657 print_usage ();
4658 return retval;
4659 }
4660
4661 if (nargin == 3)
4662 npoints = args(2).idx_type_value ();
4663
4664 if (! error_state)
4665 {
4666 octave_value arg_1 = args(0);
4667 octave_value arg_2 = args(1);
4668
4669 if (arg_1.is_single_type () || arg_2.is_single_type ())
4670 {
4671 if (arg_1.is_complex_type () || arg_2.is_complex_type ())
4672 retval = do_linspace<FloatComplexMatrix> (arg_1, arg_2, npoints);
4673 else
4674 retval = do_linspace<FloatMatrix> (arg_1, arg_2, npoints);
4675
4676 }
4677 else
4678 {
4679 if (arg_1.is_complex_type () || arg_2.is_complex_type ())
4680 retval = do_linspace<ComplexMatrix> (arg_1, arg_2, npoints);
4681 else
4682 retval = do_linspace<Matrix> (arg_1, arg_2, npoints);
4683 }
4684 }
4685 else
4686 error ("linspace: expecting third argument to be an integer");
4687
4688 return retval;
4689}
4690
4691
4692/*
4693
4694%!test
4695%! x1 = linspace (1, 2);
4696%! x2 = linspace (1, 2, 10);
4697%! x3 = linspace (1, -2, 10);
4698%! assert((size (x1) == [1, 100] && x1(1) == 1 && x1(100) == 2
4699%! && size (x2) == [1, 10] && x2(1) == 1 && x2(10) == 2
4700%! && size (x3) == [1, 10] && x3(1) == 1 && x3(10) == -2));
4701
4702
4703% assert(linspace ([1, 2; 3, 4], 5, 6), linspace (1, 5, 6));
4704
4705%!error <Invalid call to linspace.*> linspace ();
4706%!error <Invalid call to linspace.*> linspace (1, 2, 3, 4);
4707
4708%!test
4709%! fail("linspace ([1, 2; 3, 4], 5, 6)","warning");
4710
4711*/
4712
4713// FIXME -- should accept dimensions as separate args for N-d
4714// arrays as well as 1-d and 2-d arrays.
4715
4716DEFUN (resize, args, ,
4717 "-*- texinfo -*-\n\
4718@deftypefn {Built-in Function} {} resize (@var{x}, @var{m})\n\
4719@deftypefnx {Built-in Function} {} resize (@var{x}, @var{m}, @var{n})\n\
4720@deftypefnx {Built-in Function} {} resize (@var{x}, @var{m}, @var{n}, @dots{})\n\
4721Resize @var{x} cutting off elements as necessary.\n\
4722\n\
4723In the result, element with certain indices is equal to the corresponding\n\
4724element of @var{x} if the indices are within the bounds of @var{x};\n\
4725otherwise, the element is set to zero.\n\
4726\n\
4727In other words, the statement\n\
4728\n\
4729@example\n\
4730 y = resize (x, dv);\n\
4731@end example\n\
4732\n\
4733@noindent\n\
4734is equivalent to the following code:\n\
4735\n\
4736@example\n\
4737@group\n\
4738 y = zeros (dv, class (x));\n\
4739 sz = min (dv, size (x));\n\
4740 for i = 1:length (sz), idx@{i@} = 1:sz(i); endfor\n\
4741 y(idx@{:@}) = x(idx@{:@});\n\
4742@end group\n\
4743@end example\n\
4744\n\
4745@noindent\n\
4746but is performed more efficiently.\n\
4747\n\
4748If only @var{m} is supplied and it is a scalar, the dimension of the\n\
4749result is @var{m}-by-@var{m}. If @var{m} is a vector, then the\n\
4750dimensions of the result are given by the elements of @var{m}.\n\
4751If both @var{m} and @var{n} are scalars, then the dimensions of\n\
4752the result are @var{m}-by-@var{n}.\n\
4753\n\
4754An object can be resized to more dimensions than it has;\n\
4755in such case the missing dimensions are assumed to be 1.\n\
4756Resizing an object to fewer dimensions is not possible.\n\
4757@seealso{reshape, postpad}\n\
4758@end deftypefn")
4759{
4760 octave_value retval;
4761 int nargin = args.length ();
4762
4763 if (nargin == 2)
4764 {
4765 Array<double> vec = args(1).vector_value ();
4766 int ndim = vec.length ();
4767 if (ndim == 1)
4768 {
4769 octave_idx_type m = static_cast<octave_idx_type> (vec(0));
4770 retval = args(0);
4771 retval = retval.resize (dim_vector (m, m), true);
4772 }
4773 else
4774 {
4775 dim_vector dv;
4776 dv.resize (ndim);
4777 for (int i = 0; i < ndim; i++)
4778 dv(i) = static_cast<octave_idx_type> (vec(i));
4779 retval = args(0);
4780 retval = retval.resize (dv, true);
4781 }
4782 }
4783 else if (nargin > 2)
4784 {
4785 dim_vector dv;
4786 dv.resize (nargin - 1);
4787 for (octave_idx_type i = 1; i < nargin; i++)
4788 dv(i-1) = static_cast<octave_idx_type> (args(i).scalar_value ());
4789 if (!error_state)
4790 {
4791 retval = args(0);
4792 retval = retval.resize (dv, true);
4793 }
4794
4795 }
4796 else
4797 print_usage ();
4798 return retval;
4799}
4800
4801// FIXME -- should use octave_idx_type for dimensions.
4802
4803DEFUN (reshape, args, ,
4804 "-*- texinfo -*-\n\
4805@deftypefn {Built-in Function} {} reshape (@var{a}, @var{m}, @var{n}, @dots{})\n\
4806@deftypefnx {Built-in Function} {} reshape (@var{a}, @var{size})\n\
4807Return a matrix with the given dimensions whose elements are taken\n\
4808from the matrix @var{a}. The elements of the matrix are accessed in\n\
4809column-major order (like Fortran arrays are stored).\n\
4810\n\
4811For example,\n\
4812\n\
4813@example\n\
4814@group\n\
4815reshape ([1, 2, 3, 4], 2, 2)\n\
4816 @result{} 1 3\n\
4817 2 4\n\
4818@end group\n\
4819@end example\n\
4820\n\
4821@noindent\n\
4822Note that the total number of elements in the original\n\
4823matrix must match the total number of elements in the new matrix.\n\
4824\n\
4825A single dimension of the return matrix can be unknown and is flagged\n\
4826by an empty argument.\n\
4827@end deftypefn")
4828{
4829 octave_value retval;
4830
4831 int nargin = args.length ();
4832
4833 Array<int> new_size;
4834
4835 if (nargin == 2)
4836 new_size = args(1).int_vector_value ();
4837 else if (nargin > 2)
4838 {
4839 new_size.resize (nargin-1);