idx-vector.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 1993-2024 The Octave Project Developers
//
// See the file COPYRIGHT.md in the top-level directory of this
// distribution or <https://octave.org/copyright/>.
//
// This file is part of Octave.
//
// Octave is free software: you can redistribute it and/or modify it
// under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// Octave is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with Octave; see the file COPYING.  If not, see
// <https://www.gnu.org/licenses/>.
//
////////////////////////////////////////////////////////////////////////
 
#if defined (HAVE_CONFIG_H)
#  include "config.h"
#endif
 
#include <cinttypes>
#include <cstdlib>
 
#include <ostream>
 
#include "idx-vector.h"
#include "Array.h"
#include "Array-util.h"
#include "Sparse.h"
#include "Range.h"
 
#include "oct-locbuf.h"
#include "lo-error.h"
#include "lo-mappers.h"
 
OCTAVE_BEGIN_NAMESPACE(octave)
 
OCTAVE_NORETURN static void err_invalid_range ()
{
  (*current_liboctave_error_handler) ("invalid range used as index");
}
 
OCTAVE_NORETURN static void
err_index_out_of_range ()
{
  (*current_liboctave_error_handler)
    ("internal error: idx_vector index out of range");
}
 
idx_vector::idx_vector_rep *
idx_vector::nil_rep ()
{
  static idx_vector_rep ivr;
  return &ivr;
}
 
Array<octave_idx_type>
idx_vector::idx_base_rep::as_array ()
{
  (*current_liboctave_error_handler)
    ("internal error: as_array not allowed for this index class");
 
  // Never actually executed, but required to silence compiler warning
  return Array<octave_idx_type> ();
}
 
idx_vector::idx_colon_rep::idx_colon_rep (char c)
  : idx_base_rep ()
{
  if (c != ':')
    (*current_liboctave_error_handler)
      ("internal error: invalid character converted to idx_vector; must be ':'");
}
 
octave_idx_type
idx_vector::idx_colon_rep::checkelem (octave_idx_type i) const
{
  if (i < 0)
    err_index_out_of_range ();
 
  return i;
}
 
idx_vector::idx_base_rep *
idx_vector::idx_colon_rep::sort_idx (Array<octave_idx_type>&)
{
  (*current_liboctave_error_handler)
    ("internal error: idx_colon_rep::sort_idx");
}
 
std::ostream&
idx_vector::idx_colon_rep::print (std::ostream& os) const
{
  return os << ':';
}
 
idx_vector::idx_range_rep::idx_range_rep (octave_idx_type start,
    octave_idx_type limit,
    octave_idx_type step)
  : idx_base_rep (), m_start(start),
    // Round length away from zero to catch incomplete intervals
    m_len (step
           ? std::max ((limit - start + step - (step > 0 ? 1 : -1)) / step,
                       static_cast<octave_idx_type> (0))
           : -1),
    m_step (step)
{
  if (m_len < 0)
    err_invalid_range ();
  if (m_start < 0)
    err_invalid_index (m_start);
  if (m_step < 0 && m_start + (m_len-1)*m_step < 0)
    err_invalid_index (m_start + (m_len-1)*m_step);
}
 
idx_vector::idx_range_rep::idx_range_rep (const range<double>& r)
  : idx_base_rep (), m_start (0), m_len (r.numel ()), m_step (1)
{
  if (m_len < 0)
    err_invalid_range ();
 
  if (m_len > 0)
    {
      if (r.all_elements_are_ints ())
        {
          m_start = static_cast<octave_idx_type> (r.base ()) - 1;
          m_step = static_cast<octave_idx_type> (r.increment ());
          if (m_start < 0)
            err_invalid_index (m_start);
          if (m_step < 0 && m_start + (m_len - 1)*m_step < 0)
            err_invalid_index (m_start + (m_len - 1)*m_step);
        }
      else
        {
          // find first non-integer, then gripe about it
          double b = r.base ();
          double inc = r.increment ();
          err_invalid_index (b != std::trunc (b) ? b : b + inc);
        }
    }
}
 
octave_idx_type
idx_vector::idx_range_rep::checkelem (octave_idx_type i) const
{
  if (i < 0 || i >= m_len)
    err_index_out_of_range ();
 
  return m_start + i*m_step;
}
 
idx_vector::idx_base_rep *
idx_vector::idx_range_rep::sort_uniq_clone (bool)
{
  if (m_step < 0)
    return new idx_range_rep (m_start + (m_len - 1)*m_step, m_len, -m_step, DIRECT);
  else
    {
      m_count++;
      return this;
    }
}
 
idx_vector::idx_base_rep *
idx_vector::idx_range_rep::sort_idx (Array<octave_idx_type>& idx)
{
  if (m_step < 0 && m_len > 0)
    {
      idx.clear (1, m_len);
      for (octave_idx_type i = 0; i < m_len; i++)
        idx.xelem (i) = m_len - 1 - i;
      return new idx_range_rep (m_start + (m_len - 1)*m_step, m_len, -m_step, DIRECT);
    }
  else
    {
      idx.clear (1, m_len);
      for (octave_idx_type i = 0; i < m_len; i++)
        idx.xelem (i) = i;
      m_count++;
      return this;
    }
}
 
std::ostream&
idx_vector::idx_range_rep::print (std::ostream& os) const
{
  os << m_start << ':' << m_step << ':' << m_start + m_len *m_step;
  return os;
}
 
range<double>
idx_vector::idx_range_rep::unconvert () const
{
  return range<double>::make_n_element_range
         (static_cast<double> (m_start+1), static_cast<double> (m_step), m_len);
}
 
Array<octave_idx_type>
idx_vector::idx_range_rep::as_array ()
{
  Array<octave_idx_type> retval (dim_vector (1, m_len));
  for (octave_idx_type i = 0; i < m_len; i++)
    retval.xelem (i) = m_start + i*m_step;
 
  return retval;
}
 
inline octave_idx_type
convert_index (octave_idx_type i, octave_idx_type& ext)
{
  if (i <= 0)
    err_invalid_index (i-1);
 
  if (ext < i)
    ext = i;
 
  return i - 1;
}
 
inline octave_idx_type
convert_index (double x, octave_idx_type& ext)
{
  octave_idx_type i = static_cast<octave_idx_type> (x);
 
  if (static_cast<double> (i) != x)
    err_invalid_index (x-1);
 
  return convert_index (i, ext);
}
 
inline octave_idx_type
convert_index (float x, octave_idx_type& ext)
{
  return convert_index (static_cast<double> (x), ext);
}
 
template <typename T>
inline octave_idx_type
convert_index (octave_int<T> x, octave_idx_type& ext)
{
  octave_idx_type i = octave_int<octave_idx_type> (x).value ();
 
  return convert_index (i, ext);
}
 
template <typename T>
idx_vector::idx_scalar_rep::idx_scalar_rep (T x)
  : idx_base_rep (), m_data (0)
{
  octave_idx_type dummy = 0;
 
  m_data = convert_index (x, dummy);
}
 
idx_vector::idx_scalar_rep::idx_scalar_rep (octave_idx_type i)
  : idx_base_rep (), m_data (i)
{
  if (m_data < 0)
    err_invalid_index (m_data);
}
 
octave_idx_type
idx_vector::idx_scalar_rep::checkelem (octave_idx_type i) const
{
  if (i != 0)
    err_index_out_of_range ();
 
  return m_data;
}
 
idx_vector::idx_base_rep *
idx_vector::idx_scalar_rep::sort_idx (Array<octave_idx_type>& idx)
{
  idx.clear (1, 1);
  idx.fill (0);
  m_count++;
  return this;
}
 
std::ostream&
idx_vector::idx_scalar_rep::print (std::ostream& os) const
{
  return os << m_data;
}
 
double
idx_vector::idx_scalar_rep::unconvert () const
{
  return m_data + 1;
}
 
Array<octave_idx_type>
idx_vector::idx_scalar_rep::as_array ()
{
  return Array<octave_idx_type> (dim_vector (1, 1), m_data);
}
 
template <typename T>
idx_vector::idx_vector_rep::idx_vector_rep (const Array<T>& nda)
  : idx_base_rep (), m_data (nullptr), m_len (nda.numel ()), m_ext (0),
    m_aowner (nullptr), m_orig_dims (nda.dims ())
{
  if (m_len != 0)
    {
      std::unique_ptr<octave_idx_type []> d (new octave_idx_type [m_len]);
 
      for (octave_idx_type i = 0; i < m_len; i++)
        d[i] = convert_index (nda.xelem (i), m_ext);
 
      m_data = d.release ();
    }
}
 
// Note that this makes a shallow copy of the index array.
 
idx_vector::idx_vector_rep::idx_vector_rep (const Array<octave_idx_type>& inda)
  : idx_base_rep (), m_data (inda.data ()), m_len (inda.numel ()), m_ext (0),
    m_aowner (new Array<octave_idx_type> (inda)), m_orig_dims (inda.dims ())
{
  if (m_len != 0)
    {
      octave_idx_type max = -1;
      for (octave_idx_type i = 0; i < m_len; i++)
        {
          octave_idx_type k = inda.xelem (i);
          if (k < 0)
            err_invalid_index (k);
          else if (k > max)
            max = k;
        }
 
      m_ext = max + 1;
    }
}
 
idx_vector::idx_vector_rep::idx_vector_rep (const Array<octave_idx_type>& inda,
    octave_idx_type ext, direct)
  : idx_base_rep (), m_data (inda.data ()), m_len (inda.numel ()),
    m_ext (ext), m_aowner (new Array<octave_idx_type> (inda)),
    m_orig_dims (inda.dims ())
{
  // No checking.
  if (m_ext < 0)
    {
      octave_idx_type max = -1;
      for (octave_idx_type i = 0; i < m_len; i++)
        if (m_data[i] > max)
          max = m_data[i];
 
      m_ext = max + 1;
    }
}
 
idx_vector::idx_vector_rep::idx_vector_rep (bool b)
  : idx_base_rep (), m_data (nullptr), m_len (b ? 1 : 0), m_ext (0),
    m_aowner (nullptr), m_orig_dims (m_len, m_len)
{
  if (m_len != 0)
    {
      octave_idx_type *d = new octave_idx_type [1];
      d[0] = 0;
      m_data = d;
      m_ext = 1;
    }
}
 
idx_vector::idx_vector_rep::idx_vector_rep (const Array<bool>& bnda,
    octave_idx_type nnz)
  : idx_base_rep (), m_data (nullptr), m_len (nnz), m_ext (0),
    m_aowner (nullptr), m_orig_dims ()
{
  if (nnz < 0)
    m_len = bnda.nnz ();
 
  const dim_vector dv = bnda.dims ();
 
  m_orig_dims = dv.make_nd_vector (m_len);
 
  if (m_len != 0)
    {
      octave_idx_type *d = new octave_idx_type [m_len];
 
      octave_idx_type ntot = bnda.numel ();
 
      octave_idx_type k = 0;
      for (octave_idx_type i = 0; i < ntot; i++)
        if (bnda.xelem (i))
          d[k++] = i;
 
      m_data = d;
 
      m_ext = d[k-1] + 1;
    }
}
 
idx_vector::idx_vector_rep::idx_vector_rep (const Sparse<bool>& bnda)
  : idx_base_rep (), m_data (nullptr), m_len (bnda.nnz ()), m_ext (0),
    m_aowner (nullptr), m_orig_dims ()
{
  const dim_vector dv = bnda.dims ();
 
  m_orig_dims = dv.make_nd_vector (m_len);
 
  if (m_len != 0)
    {
      octave_idx_type *d = new octave_idx_type [m_len];
 
      octave_idx_type k = 0;
      octave_idx_type nc = bnda.cols ();
      octave_idx_type nr = bnda.rows ();
 
      for (octave_idx_type j = 0; j < nc; j++)
        for (octave_idx_type i = bnda.cidx (j); i < bnda.cidx (j+1); i++)
          if (bnda.data (i))
            d[k++] = j * nr + bnda.ridx (i);
 
      m_data = d;
 
      m_ext = d[k-1] + 1;
    }
}
 
idx_vector::idx_vector_rep::~idx_vector_rep ()
{
  if (m_aowner)
    delete m_aowner;
  else
    delete [] m_data;
}
 
octave_idx_type
idx_vector::idx_vector_rep::checkelem (octave_idx_type n) const
{
  if (n < 0 || n >= m_len)
    err_invalid_index (n);
 
  return xelem (n);
}
 
idx_vector::idx_base_rep *
idx_vector::idx_vector_rep::sort_uniq_clone (bool uniq)
{
  if (m_len == 0)
    {
      m_count++;
      return this;
    }
 
  // This is wrapped in unique_ptr so that we don't leak on out-of-memory.
  std::unique_ptr<idx_vector_rep> new_rep
  (new idx_vector_rep (nullptr, m_len, m_ext, m_orig_dims, DIRECT));
 
  if (m_ext > m_len*math::log2 (1.0 + m_len))
    {
      // Use standard sort via octave_sort.
      octave_idx_type *new_data = new octave_idx_type [m_len];
      new_rep->m_data = new_data;
 
      std::copy_n (m_data, m_len, new_data);
      octave_sort<octave_idx_type> lsort;
      lsort.set_compare (ASCENDING);
      lsort.sort (new_data, m_len);
 
      if (uniq)
        {
          octave_idx_type new_len = std::unique (new_data, new_data + m_len)
                                    - new_data;
          new_rep->m_len = new_len;
          if (new_rep->m_orig_dims.ndims () == 2 && new_rep->m_orig_dims(0) == 1)
            new_rep->m_orig_dims = dim_vector (1, new_len);
          else
            new_rep->m_orig_dims = dim_vector (new_len, 1);
        }
    }
  else if (uniq)
    {
      // Use two-pass bucket sort (only a mask array needed).
      OCTAVE_LOCAL_BUFFER_INIT (bool, has, m_ext, false);
      for (octave_idx_type i = 0; i < m_len; i++)
        has[m_data[i]] = true;
 
      octave_idx_type new_len = 0;
      for (octave_idx_type i = 0; i < m_ext; i++)
        new_len += has[i];
 
      new_rep->m_len = new_len;
      if (new_rep->m_orig_dims.ndims () == 2 && new_rep->m_orig_dims(0) == 1)
        new_rep->m_orig_dims = dim_vector (1, new_len);
      else
        new_rep->m_orig_dims = dim_vector (new_len, 1);
 
      octave_idx_type *new_data = new octave_idx_type [new_len];
      new_rep->m_data = new_data;
 
      for (octave_idx_type i = 0, j = 0; i < m_ext; i++)
        if (has[i])
          new_data[j++] = i;
    }
  else
    {
      // Use two-pass bucket sort.
      OCTAVE_LOCAL_BUFFER_INIT (octave_idx_type, cnt, m_ext, 0);
      for (octave_idx_type i = 0; i < m_len; i++)
        cnt[m_data[i]]++;
 
      octave_idx_type *new_data = new octave_idx_type [m_len];
      new_rep->m_data = new_data;
 
      for (octave_idx_type i = 0, j = 0; i < m_ext; i++)
        {
          for (octave_idx_type k = 0; k < cnt[i]; k++)
            new_data[j++] = i;
        }
    }
 
  return new_rep.release ();
}
 
idx_vector::idx_base_rep *
idx_vector::idx_vector_rep::sort_idx (Array<octave_idx_type>& idx)
{
  // This is wrapped in unique_ptr so that we don't leak on out-of-memory.
  std::unique_ptr<idx_vector_rep> new_rep
  (new idx_vector_rep (nullptr, m_len, m_ext, m_orig_dims, DIRECT));
 
  if (m_ext > m_len*math::log2 (1.0 + m_len))
    {
      // Use standard sort via octave_sort.
      idx.clear (m_orig_dims);
      octave_idx_type *idx_data = idx.fortran_vec ();
      for (octave_idx_type i = 0; i < m_len; i++)
        idx_data[i] = i;
 
      octave_idx_type *new_data = new octave_idx_type [m_len];
      new_rep->m_data = new_data;
      std::copy_n (m_data, m_len, new_data);
 
      octave_sort<octave_idx_type> lsort;
      lsort.set_compare (ASCENDING);
      lsort.sort (new_data, idx_data, m_len);
    }
  else
    {
      // Use two-pass bucket sort.
      OCTAVE_LOCAL_BUFFER_INIT (octave_idx_type, cnt, m_ext, 0);
 
      for (octave_idx_type i = 0; i < m_len; i++)
        cnt[m_data[i]]++;
 
      idx.clear (m_orig_dims);
      octave_idx_type *idx_data = idx.fortran_vec ();
 
      octave_idx_type *new_data = new octave_idx_type [m_len];
      new_rep->m_data = new_data;
 
      for (octave_idx_type i = 0, k = 0; i < m_ext; i++)
        {
          octave_idx_type j = cnt[i];
          cnt[i] = k;
          k += j;
        }
 
      for (octave_idx_type i = 0; i < m_len; i++)
        {
          octave_idx_type j = m_data[i];
          octave_idx_type k = cnt[j]++;
          new_data[k] = j;
          idx_data[k] = i;
        }
    }
 
  return new_rep.release ();
}
 
std::ostream&
idx_vector::idx_vector_rep::print (std::ostream& os) const
{
  os << '[';
 
  for (octave_idx_type i = 0; i < m_len - 1; i++)
    os << m_data[i] << ',' << ' ';
 
  if (m_len > 0)
    os << m_data[m_len-1];
 
  os << ']';
 
  return os;
}
 
Array<double>
idx_vector::idx_vector_rep::unconvert () const
{
  Array<double> retval (m_orig_dims);
  for (octave_idx_type i = 0; i < m_len; i++)
    retval.xelem (i) = m_data[i] + 1;
  return retval;
}
 
Array<octave_idx_type>
idx_vector::idx_vector_rep::as_array ()
{
  if (m_aowner)
    return *m_aowner;
  else
    {
      Array<octave_idx_type> retval (m_orig_dims);
 
      if (m_data)
        {
          std::memcpy (retval.fortran_vec (), m_data, m_len* sizeof (octave_idx_type));
          // Delete the old copy and share the m_data instead to save memory.
          delete [] m_data;
        }
 
      m_data = retval.fortran_vec ();
      m_aowner = new Array<octave_idx_type> (retval);
 
      return retval;
    }
}
 
idx_vector::idx_mask_rep::idx_mask_rep (bool b)
  : idx_base_rep (), m_data (nullptr), m_len (b ? 1 : 0), m_ext (0),
    m_lsti (-1), m_lste (-1), m_aowner (nullptr), m_orig_dims (m_len, m_len)
{
  if (m_len != 0)
    {
      bool *d = new bool [1];
      d[0] = true;
      m_data = d;
      m_ext = 1;
    }
}
 
idx_vector::idx_mask_rep::idx_mask_rep (const Array<bool>& bnda,
                                        octave_idx_type nnz)
  : idx_base_rep (), m_data (nullptr), m_len (nnz), m_ext (bnda.numel ()),
    m_lsti (-1), m_lste (-1), m_aowner (nullptr), m_orig_dims ()
{
  if (nnz < 0)
    m_len = bnda.nnz ();
 
  // We truncate the extent as much as possible.  For Matlab
  // compatibility, but maybe it's not a bad idea anyway.
  while (m_ext > 0 && ! bnda(m_ext-1))
    m_ext--;
 
  const dim_vector dv = bnda.dims ();
 
  m_orig_dims = dv.make_nd_vector (m_len);
 
  m_aowner = new Array<bool> (bnda);
  m_data = bnda.data ();
}
 
idx_vector::idx_mask_rep::~idx_mask_rep ()
{
  if (m_aowner)
    delete m_aowner;
  else
    delete [] m_data;
}
 
octave_idx_type
idx_vector::idx_mask_rep::xelem (octave_idx_type n) const
{
  if (n == m_lsti + 1)
    {
      m_lsti = n;
      while (! m_data[++m_lste]) ;
    }
  else
    {
      m_lsti = n++;
      m_lste = -1;
      while (n > 0)
        if (m_data[++m_lste]) --n;
    }
  return m_lste;
}
 
octave_idx_type
idx_vector::idx_mask_rep::checkelem (octave_idx_type n) const
{
  if (n < 0 || n >= m_len)
    err_invalid_index (n);
 
  return xelem (n);
}
 
std::ostream&
idx_vector::idx_mask_rep::print (std::ostream& os) const
{
  os << '[';
 
  for (octave_idx_type i = 0; i < m_ext - 1; i++)
    os << m_data[i] << ',' << ' ';
 
  if (m_ext > 0)
    os << m_data[m_ext-1];
 
  os << ']';
 
  return os;
}
 
Array<bool>
idx_vector::idx_mask_rep::unconvert () const
{
  if (m_aowner)
    return *m_aowner;
  else
    {
      Array<bool> retval (dim_vector (m_ext, 1));
      for (octave_idx_type i = 0; i < m_ext; i++)
        retval.xelem (i) = m_data[i];
      return retval;
    }
}
 
Array<octave_idx_type>
idx_vector::idx_mask_rep::as_array ()
{
  if (m_aowner)
    return m_aowner->find ().reshape (m_orig_dims);
  else
    {
      Array<bool> retval (m_orig_dims);
      for (octave_idx_type i = 0, j = 0; i < m_ext; i++)
        if (m_data[i])
          retval.xelem (j++) = i;
 
      return retval;
    }
}
 
idx_vector::idx_base_rep *
idx_vector::idx_mask_rep::sort_idx (Array<octave_idx_type>& idx)
{
  idx.clear (m_len, 1);
  for (octave_idx_type i = 0; i < m_len; i++)
    idx.xelem (i) = i;
 
  m_count++;
  return this;
}
 
const idx_vector idx_vector::colon (new idx_vector::idx_colon_rep ());
 
idx_vector::idx_vector (const Array<bool>& bnda)
  : m_rep (nullptr)
{
  // Convert only if it means saving at least half the memory.
  static const int factor = (2 * sizeof (octave_idx_type));
  octave_idx_type nnz = bnda.nnz ();
  if (nnz <= bnda.numel () / factor)
    m_rep = new idx_vector_rep (bnda, nnz);
  else
    m_rep = new idx_mask_rep (bnda, nnz);
}
 
bool
idx_vector::maybe_reduce (octave_idx_type n, const idx_vector& j,
                          octave_idx_type nj)
{
  bool reduced = false;
 
  // Empty index always reduces.
  if (m_rep->length (n) == 0)
    {
      *this = idx_vector ();
      return true;
    }
 
  // Possibly skip singleton dims.
  if (n == 1 && m_rep->is_colon_equiv (n))
    {
      *this = j;
      return true;
    }
 
  if (nj == 1 && j.is_colon_equiv (nj))
    return true;
 
  switch (j.idx_class ())
    {
    case class_colon:
      switch (m_rep->idx_class ())
        {
        case class_colon:
          // (:,:) reduces to (:)
          reduced = true;
          break;
 
        case class_scalar:
          {
            // (i,:) reduces to a range.
            idx_scalar_rep *r = dynamic_cast<idx_scalar_rep *> (m_rep);
            octave_idx_type k = r->get_data ();
            *this = new idx_range_rep (k, nj, n, DIRECT);
            reduced = true;
          }
          break;
 
        case class_range:
          {
            // (i:k:end,:) reduces to a range if i <= k and k divides n.
            idx_range_rep *r = dynamic_cast<idx_range_rep *> (m_rep);
            octave_idx_type s = r->get_start ();
            octave_idx_type l = r->length (n);
            octave_idx_type t = r->get_step ();
            if (l*t == n)
              {
                *this = new idx_range_rep (s, l * nj, t, DIRECT);
                reduced = true;
              }
          }
          break;
 
        default:
          break;
        }
      break;
 
    case class_range:
      switch (m_rep->idx_class ())
        {
        case class_colon:
          {
            // (:,i:j) reduces to a range (the m_step must be 1)
            idx_range_rep *rj = dynamic_cast<idx_range_rep *> (j.m_rep);
            if (rj->get_step () == 1)
              {
                octave_idx_type sj = rj->get_start ();
                octave_idx_type lj = rj->length (nj);
                *this = new idx_range_rep (sj * n, lj * n, 1, DIRECT);
                reduced = true;
              }
          }
          break;
 
        case class_scalar:
          {
            // (k,i:d:j) reduces to a range.
            idx_scalar_rep *r = dynamic_cast<idx_scalar_rep *> (m_rep);
            idx_range_rep *rj = dynamic_cast<idx_range_rep *> (j.m_rep);
            octave_idx_type k = r->get_data ();
            octave_idx_type sj = rj->get_start ();
            octave_idx_type lj = rj->length (nj);
            octave_idx_type tj = rj->get_step ();
            *this = new idx_range_rep (n * sj + k, lj, n * tj, DIRECT);
            reduced = true;
          }
          break;
 
        case class_range:
          {
            // (i:k:end,p:q) reduces to a range if i <= k and k divides n.
            // (ones (1, m), ones (1, n)) reduces to (ones (1, m*n))
            idx_range_rep *r = dynamic_cast<idx_range_rep *> (m_rep);
            octave_idx_type s = r->get_start ();
            octave_idx_type l = r->length (n);
            octave_idx_type t = r->get_step ();
            idx_range_rep *rj = dynamic_cast<idx_range_rep *> (j.m_rep);
            octave_idx_type sj = rj->get_start ();
            octave_idx_type lj = rj->length (nj);
            octave_idx_type tj = rj->get_step ();
            if ((l*t == n && tj == 1) || (t == 0 && tj == 0))
              {
                *this = new idx_range_rep (s + n * sj, l * lj, t, DIRECT);
                reduced = true;
              }
          }
          break;
 
        default:
          break;
        }
      break;
 
    case class_scalar:
      switch (m_rep->idx_class ())
        {
        case class_scalar:
          {
            // (i,j) reduces to a single index.
            idx_scalar_rep *r = dynamic_cast<idx_scalar_rep *> (m_rep);
            idx_scalar_rep *rj = dynamic_cast<idx_scalar_rep *> (j.m_rep);
            octave_idx_type k = r->get_data () + n * rj->get_data ();
            *this = new idx_scalar_rep (k, DIRECT);
            reduced = true;
          }
          break;
 
        case class_range:
          {
            // (i:d:j,k) reduces to a range.
            idx_range_rep *r = dynamic_cast<idx_range_rep *> (m_rep);
            idx_scalar_rep *rj = dynamic_cast<idx_scalar_rep *> (j.m_rep);
            octave_idx_type s = r->get_start ();
            octave_idx_type l = r->length (nj);
            octave_idx_type t = r->get_step ();
            octave_idx_type k = rj->get_data ();
            *this = new idx_range_rep (n * k + s, l, t, DIRECT);
            reduced = true;
          }
          break;
 
        case class_colon:
          {
            // (:,k) reduces to a range.
            idx_scalar_rep *rj = dynamic_cast<idx_scalar_rep *> (j.m_rep);
            octave_idx_type k = rj->get_data ();
            *this = new idx_range_rep (n * k, n, 1, DIRECT);
            reduced = true;
          }
          break;
 
        default:
          break;
        }
      break;
 
    default:
      break;
    }
 
  return reduced;
}
 
bool
idx_vector::is_cont_range (octave_idx_type n, octave_idx_type& l,
                           octave_idx_type& u) const
{
  bool res = false;
 
  switch (m_rep->idx_class ())
    {
    case class_colon:
      l = 0; u = n;
      res = true;
      break;
 
    case class_range:
      {
        idx_range_rep *r = dynamic_cast<idx_range_rep *> (m_rep);
        if (r->get_step () == 1)
          {
            l = r->get_start ();
            u = l + r->length (n);
            res = true;
          }
      }
      break;
 
    case class_scalar:
      {
        idx_scalar_rep *r = dynamic_cast<idx_scalar_rep *> (m_rep);
        l = r->get_data ();
        u = l + 1;
        res = true;
      }
      break;
 
    case class_mask:
      {
        idx_mask_rep *r = dynamic_cast<idx_mask_rep *> (m_rep);
        octave_idx_type m_ext = r->extent (0);
        octave_idx_type m_len = r->length (0);
        if (m_ext == m_len)
          {
            l = 0;
            u = m_len;
            res = true;
          }
      }
 
    default:
      break;
    }
 
  return res;
}
 
octave_idx_type
idx_vector::increment () const
{
  octave_idx_type retval = 0;
 
  switch (m_rep->idx_class ())
    {
    case class_colon:
      retval = 1;
      break;
 
    case class_range:
      retval = dynamic_cast<idx_range_rep *> (m_rep) -> get_step ();
      break;
 
    case class_vector:
    case class_mask:
      {
        if (length (0) > 1)
          retval = elem (1) - elem (0);
      }
      break;
 
    default:
      break;
    }
 
  return retval;
}
 
const octave_idx_type *
idx_vector::raw ()
{
  if (m_rep->idx_class () != class_vector)
    *this = idx_vector (as_array (), extent (0));
 
  idx_vector_rep *r = dynamic_cast<idx_vector_rep *> (m_rep);
 
  assert (r != nullptr);
 
  return r->get_data ();
}
 
void
idx_vector::copy_data (octave_idx_type *m_data) const
{
  octave_idx_type m_len = m_rep->length (0);
 
  switch (m_rep->idx_class ())
    {
    case class_colon:
      (*current_liboctave_error_handler) ("colon not allowed");
      break;
 
    case class_range:
      {
        idx_range_rep *r = dynamic_cast<idx_range_rep *> (m_rep);
        octave_idx_type m_start = r->get_start ();
        octave_idx_type m_step = r->get_step ();
        octave_idx_type i, j;
        if (m_step == 1)
          for (i = m_start, j = m_start + m_len; i < j; i++) *m_data++ = i;
        else if (m_step == -1)
          for (i = m_start, j = m_start - m_len; i > j; i--) *m_data++ = i;
        else
          for (i = 0, j = m_start; i < m_len; i++, j += m_step) *m_data++ = j;
      }
      break;
 
    case class_scalar:
      {
        idx_scalar_rep *r = dynamic_cast<idx_scalar_rep *> (m_rep);
        *m_data = r->get_data ();
      }
      break;
 
    case class_vector:
      {
        idx_vector_rep *r = dynamic_cast<idx_vector_rep *> (m_rep);
        const octave_idx_type *rdata = r->get_data ();
        std::copy_n (rdata, m_len, m_data);
      }
      break;
 
    case class_mask:
      {
        idx_mask_rep *r = dynamic_cast<idx_mask_rep *> (m_rep);
        const bool *mask = r->get_data ();
        octave_idx_type m_ext = r->extent (0);
        for (octave_idx_type i = 0, j = 0; i < m_ext; i++)
          if (mask[i])
            m_data[j++] = i;
      }
      break;
 
    default:
      assert (false);
      break;
    }
}
 
idx_vector
idx_vector::complement (octave_idx_type n) const
{
  idx_vector retval;
  if (extent (n) > n)
    (*current_liboctave_error_handler)
      ("internal error: out of range complement index requested");
 
  if (idx_class () == class_mask)
    {
      idx_mask_rep *r = dynamic_cast<idx_mask_rep *> (m_rep);
      octave_idx_type nz = r->length (0);
      octave_idx_type m_ext = r->extent (0);
      Array<bool> mask (dim_vector (n, 1));
      const bool *m_data = r->get_data ();
      bool *ndata = mask.fortran_vec ();
      for (octave_idx_type i = 0; i < m_ext; i++)
        ndata[i] = ! m_data[i];
      std::fill_n (ndata + m_ext, n - m_ext, true);
      retval = new idx_mask_rep (mask, n - nz);
    }
  else
    {
      Array<bool> mask (dim_vector (n, 1), true);
      fill (false, length (n), mask.fortran_vec ());
      retval = idx_vector (mask);
    }
 
  return retval;
}
 
bool
idx_vector::is_permutation (octave_idx_type n) const
{
  bool retval = false;
 
  if (is_colon_equiv (n))
    retval = true;
  else if (length(n) == n && extent(n) == n)
    {
      OCTAVE_LOCAL_BUFFER_INIT (bool, left, n, true);
 
      retval = true;
 
      for (octave_idx_type i = 0, m_len = length (); i < m_len; i++)
        {
          octave_idx_type k = xelem (i);
          if (left[k])
            left[k] = false;
          else
            {
              retval = false;
              break;
            }
        }
    }
 
  return retval;
}
 
idx_vector
idx_vector::inverse_permutation (octave_idx_type n) const
{
  assert (n == length (n));
 
  idx_vector retval;
 
  switch (idx_class ())
    {
    case class_range:
      {
        if (increment () == -1)
          retval = sorted ();
        else
          retval = *this;
        break;
      }
    case class_vector:
      {
        idx_vector_rep *r = dynamic_cast<idx_vector_rep *> (m_rep);
        const octave_idx_type *ri = r->get_data ();
        Array<octave_idx_type> idx (orig_dimensions ());
        for (octave_idx_type i = 0; i < n; i++)
          idx.xelem (ri[i]) = i;
        retval = new idx_vector_rep (idx, r->extent (0), DIRECT);
        break;
      }
    default:
      retval = *this;
      break;
    }
 
  return retval;
}
 
idx_vector
idx_vector::unmask () const
{
  if (idx_class () == class_mask)
    {
      idx_mask_rep *r = dynamic_cast<idx_mask_rep *> (m_rep);
      const bool *m_data = r->get_data ();
      octave_idx_type m_ext = r->extent (0);
      octave_idx_type m_len = r->length (0);
      octave_idx_type *idata = new octave_idx_type [m_len];
 
      for (octave_idx_type i = 0, j = 0; i < m_ext; i++)
        if (m_data[i])
          idata[j++] = i;
 
      m_ext = (m_len > 0 ? idata[m_len - 1] + 1 : 0);
 
      return new idx_vector_rep (idata, m_len, m_ext, r->orig_dimensions (),
                                 DIRECT);
    }
  else
    return *this;
}
 
void
idx_vector::unconvert (idx_class_type& iclass,
                       double& scalar, range<double>& range,
                       Array<double>& array, Array<bool>& mask) const
{
  iclass = idx_class ();
  switch (iclass)
    {
    case class_colon:
      break;
 
    case class_range:
      {
        idx_range_rep *r = dynamic_cast<idx_range_rep *> (m_rep);
        range = r->unconvert ();
      }
      break;
 
    case class_scalar:
      {
        idx_scalar_rep *r = dynamic_cast<idx_scalar_rep *> (m_rep);
        scalar = r->unconvert ();
      }
      break;
 
    case class_vector:
      {
        idx_vector_rep *r = dynamic_cast<idx_vector_rep *> (m_rep);
        array = r->unconvert ();
      }
      break;
 
    case class_mask:
      {
        idx_mask_rep *r = dynamic_cast<idx_mask_rep *> (m_rep);
        mask = r->unconvert ();
      }
      break;
 
    default:
      assert (false);
      break;
    }
}
 
Array<octave_idx_type>
idx_vector::as_array () const
{
  return m_rep->as_array ();
}
 
bool
idx_vector::isvector () const
{
  return idx_class () != class_vector || orig_dimensions ().isvector ();
}
 
octave_idx_type
idx_vector::freeze (octave_idx_type z_len, const char *, bool resize_ok)
{
  if (! resize_ok && extent (z_len) > z_len)
    (*current_liboctave_error_handler)
      ("invalid matrix index = %" OCTAVE_IDX_TYPE_FORMAT, extent (z_len));
 
  return length (z_len);
}
 
octave_idx_type
idx_vector::ones_count () const
{
  octave_idx_type n = 0;
 
  if (is_colon ())
    n = 1;
  else
    {
      for (octave_idx_type i = 0; i < length (1); i++)
        if (xelem (i) == 0)
          n++;
    }
 
  return n;
}
 
// Instantiate the octave_int constructors we want.
#define INSTANTIATE_SCALAR_VECTOR_REP_CONST(T)                          \
  template OCTAVE_API idx_vector::idx_scalar_rep::idx_scalar_rep (T);   \
  template OCTAVE_API idx_vector::idx_vector_rep::idx_vector_rep (const Array<T>&);
 
INSTANTIATE_SCALAR_VECTOR_REP_CONST (float)
INSTANTIATE_SCALAR_VECTOR_REP_CONST (double)
INSTANTIATE_SCALAR_VECTOR_REP_CONST (octave_int8)
INSTANTIATE_SCALAR_VECTOR_REP_CONST (octave_int16)
INSTANTIATE_SCALAR_VECTOR_REP_CONST (octave_int32)
INSTANTIATE_SCALAR_VECTOR_REP_CONST (octave_int64)
INSTANTIATE_SCALAR_VECTOR_REP_CONST (octave_uint8)
INSTANTIATE_SCALAR_VECTOR_REP_CONST (octave_uint16)
INSTANTIATE_SCALAR_VECTOR_REP_CONST (octave_uint32)
INSTANTIATE_SCALAR_VECTOR_REP_CONST (octave_uint64)
 
OCTAVE_END_NAMESPACE(octave)
 
/*
 
%!error id=Octave:index-out-of-bounds 1(find ([1,1] != 0))
%!assert ((1:3)(find ([1,0,1] != 0)), [1,3])
 
*/