AbsArr.h

Go to the documentation of this file.

#ifndef ABSARR_H
#define ABSARR_H
/*
Simple arrays with variable length.
DynLinArr - "Dynamic Linear Array" - it is what it is.
The name is choosen to distinguish this template class from every known
other libraries.
The information is kept in one block of memory and accessible with
indexes. The array is characterized by the only parameter - its physical
size. If you need a logicl size different from physical, keep the
logical size in a separated place.
Increasing of the array size leads to allocation of
another block of memory and copying all elements to it by
operator= (which can be correctly defined in the class definitions).
This is difference with STL, in which vector have logical size
which can be smaller than the volume of allocated memory.
Currently the array boundaries are checked at each access, this is
also difference with STL.
These checks can be switched out by undefining macro ALR_CHECK_BOUND.
The copying of array leads to copying all elements, descruction - to
desctruction of elements.
The array DynLinArr can keep other array DynLinArr as elements
and constitute multy-dimensional array  of not "parallelogram" shape,
that is for each first index, the dimension corresponding to
the second index can be different from that of the other first indexes.
This is very powerfull feature, but it is not always desirable.
If you need box-like array or "parallelogram" shape, use DynArr -
"Dynamic Array".
 
In DynArr the number of dimensions is arbitrary, but size of this
array along each dimension is fixed and independent on
the indexes of other dimensions. So this is rectangular array.
The access is performed by indexes,
but not through operators []. Since it is not trivial to provide
such access for this array, I preferred the simplest solutions
and arranged this access through little functions "ac":
  T& DynArr::ac(long i)  // for 1-dimensional array
  T& DynArr::ac(long i1, long i2)  // for 2-dimensional array
  T& DynArr::ac(const DynLinArr<long>& ind)
    // for arbitrary number of dimensions
    // but the number of them in array should be equal to size of ind.
    // ind is array of indexes. Its first element if the first index,
    // second is second, etc.
class DynArr is constructed with the help of DynLinArr.
DynArr can keep as elements another DynArr's or DynLinArr's.
DynLinArr can also keep DynArr as elements.
 
There are many various auxiliary utilites assosiated with these arrays.
Some of them look little bit antiquated, since they was created in
stone age and were not revised after innovation of electricity.:)
This happened because the author needed to port his programs on
very various computers, not all of them had very advanced C++-related
software. Today the situation with C++ is being changed rapidly
and some of these pearls may be arranged by more modern way.
But the main components of this file are actual and very convenient.
 
Copyright (c) 1999-2005 I. B. Smirnov
 
The file can be used, copied, modified, and distributed
according to the terms of GNU Lesser General Public License version 2.1
as published by the Free Software Foundation,
and provided that the above copyright notice, this permission notice,
and notices about any modifications of the original text
appear in all copies and in supporting documentation.
The file is provided "as is" without express or implied warranty.
*/
#include <iostream>
#include <iomanip>
#include "wcpplib/stream/prstream.h"
#include "wcpplib/util/FunNameStack.h"
#include "wcpplib/util/emul_new_stand.h"
#include "wcpplib/util/String.h"
 
// Here there is a good place to switch on the bound check in all programs
#ifndef ALR_CHECK_BOUND
#define ALR_CHECK_BOUND
#endif
#ifndef ALR_CHECK_EACH_BOUND
#define ALR_CHECK_EACH_BOUND
#endif
 
#ifdef USE_REPLACE_ALLOC
#include "wcpplib/safetl/ReplaceAlloc.h"
#endif
 
//#define DEBUG_DYNLINARR  // make some print
//#define DEBUG_DYNARR  // make some print and in addition
// functions from DynArr make some formally unnecessary checks, which are
// hiwever useful for debug.
 
// Here there is a good place to switch off the bound check in all programs
//#ifdef ALR_CHECK_BOUND
//#undef ALR_CHECK_BOUND
//#endif
//#ifdef ALR_CHECK_BOUND
//error in program
//#endif
//#define FUNCTINH
 
//#define DONT_USE_ABSPTR
#ifndef DONT_USE_ABSPTR  // in oder to supply some programs without
                         // smart pointers
#include "wcpplib/safetl/AbsPtr.h"
 
// Note: at current setup, if DONT_USE_ABSPTR is not activated,
// DynLinArr and DynArr are both addressable from ActivePtr,
// since they are derived from RegPassivePtr, which is itself addressable,
// and since they have virtual copy functions and constructors themselves.
 
#else
// Some necessary repetitions from AbsPtr.h:
enum Pilfer {
  steal
};
 
// See AnsPtr.h for details regarding the following:
#define PILF_CONST const
#define PILF_MUTABLE mutable
//#define PILF_CONST
//#define PILF_MUTABLE
 
#endif
#include "wcpplib/math/minmax.h"
 
extern long max_qel_DynLinArr;  // increase it if need
// Helps to detect access to not inited DynLinArr,
// what may happen at initializetion of class members
// and use of not inited member for initialization of the other.
 
template <class T> class DynArr;
 
//enum ArgInterp { interp_as_arr };  // these "enums" are replaced to
// class types because it was noticed that unexpected conversions
// from and to int may spoil everything, also see comment just below.
class ArgInterp_Arr {
};
//enum PutArgInterp { interp_as_adr };
class ArgInterp_SingleAdr  // for put_qel()
    {
};
//class PutArgInterp_Arr
//{ };
//enum ArgInterpVal { interp_as_val };  // No, it can be converted
// to int and leads to call of incorrect contructors.
// Trying to use empty classes.
class ArgInterp_Val {
};
 
#ifndef DONT_USE_ABSPTR
template <class T>
class DynLinArr : public RegPassivePtr
#else
                  template <class T>
                  class DynLinArr
#endif
                  {
 public:
  // Constructors
  DynLinArr(void) : qel(0), el(NULL) { ; }
  explicit DynLinArr(long fqel) : qel(fqel), el(NULL) {
    if (qel > max_qel_DynLinArr) {
      mcerr << "ERROR in DynLinArr(long fqel):\n";
      mcerr << "qel > max_qel_DynLinArr:\n";
      Iprint2n(mcout, qel, max_qel_DynLinArr);
      mcerr << "Type of T is (in internal notations) " << typeid(T).name()
            << '\n';
      spexit(mcerr);
    }
    if (qel < 0) {
      mcerr << "ERROR in DynLinArr(long fqel):\n";
      mcerr << "qel < 0:\n";
      Iprintn(mcout, qel);
      mcerr << "Type of T is (in internal notations) " << typeid(T).name()
            << '\n';
      spexit(mcerr);
    }
#ifdef USE_REPLACE_ALLOC
    if (fqel > 0) {
      el = (T*)malloc(sizeof(T) * fqel);
      long n;
      for (n = 0; n < fqel; n++)
        new (&(el[n])) T;
    }
#else
    el = (fqel > 0) ? (new T[fqel]) : (T*)NULL;
#endif
  }
 
  DynLinArr(long fqel, const T& val) : qel(fqel), el(NULL) {
    if (qel > max_qel_DynLinArr) {
      mcerr << "ERROR in DynLinArr(long fqel, const T& val):\n";
      mcerr << "qel > max_qel_DynLinArr:\n";
      Iprint2n(mcout, qel, max_qel_DynLinArr);
      mcerr << "Type of T is (in internal notations) " << typeid(T).name()
            << '\n';
      spexit(mcerr);
    }
    if (qel < 0) {
      mcerr << "ERROR in DynLinArr(long fqel, const T& val):\n";
      mcerr << "qel < 0:\n";
      Iprintn(mcout, qel);
      mcerr << "Type of T is (in internal notations) " << typeid(T).name()
            << '\n';
      spexit(mcerr);
    }
#ifdef USE_REPLACE_ALLOC
    if (fqel > 0) {
      el = (T*)malloc(sizeof(T) * fqel);
      long n;
      for (n = 0; n < fqel; n++)
        new (&(el[n])) T;
    }
#else
    el = (fqel > 0) ? (new T[fqel]) : (T*)NULL;
#endif
    assignAll(val);
  }
  DynLinArr(long fqel, const T* ar, ArgInterp_Arr /*t*/) : qel(fqel), el(NULL) {
    if (qel > max_qel_DynLinArr) {
      mcerr << "ERROR in DynLinArr(long fqel, const T* ar, ArgInterp_Arr):\n";
      mcerr << "qel > max_qel_DynLinArr:\n";
      Iprint2n(mcout, qel, max_qel_DynLinArr);
      mcerr << "Type of T is (in internal notations) " << typeid(T).name()
            << '\n';
      spexit(mcerr);
    }
    if (qel < 0) {
      mcerr << "ERROR in DynLinArr(long fqel, const T* ar, ArgInterp_Arr):\n";
      mcerr << "qel < 0:\n";
      Iprintn(mcout, qel);
      mcerr << "Type of T is (in internal notations) " << typeid(T).name()
            << '\n';
      spexit(mcerr);
    }
#ifdef USE_REPLACE_ALLOC
    if (fqel > 0) {
      el = (T*)malloc(sizeof(T) * fqel);
      long n;
      for (n = 0; n < fqel; n++)
        new (&(el[n])) T;
    }
#else
    el = (fqel > 0) ? (new T[fqel]) : (T*)NULL;
#endif
    long n;
    for (n = 0; n < qel; n++)
      el[n] = ar[n];
  }
  // const T* ar is array here (of course).
  // ArgInterp_Arr serves to distinguish this
  // from previous constructor when second argument
  // does not coincide with T but can be converted to it.
  // Typical example is when T is double and second argument is 0.
  // It is assumed that  0 should be converted to 0.0 and previous
  // constructor should be called. But in practice (Red Hat Linux 6.0)
  // the compiler says
  // call of overloaded `DynLinArr(int &, int)' is ambiguous.
  // I don't know whether this is error of particular compiler or
  // general problem. But the third dummy argument is anyway convenient
  // to distringuish these cases.
 
  DynLinArr<T>& operator=(const DynLinArr<T>& f);
  template <class D> DynLinArr<T>& operator=(const DynLinArr<D>& f);
 
  void pass(long fqel, T* fel) {
    // Do not call directly! Is to be used only
    // from assignment operator above
    clear();
    qel = fqel;
    el = fel;
  }
 
  inline DynLinArr(const DynLinArr<T>& f);
  DynLinArr(PILF_CONST DynLinArr<T>& f, Pilfer) : qel(f.qel), el(f.el) {
#ifdef DEBUG_DYNLINARR
    mcout << "DynLinArr( DynLinArr<T>& f, Pilfer) is working\n";
#endif
    f.qel = 0;
    f.el = 0;
  }
  DynLinArr(const DynArr<T>& f);  // works only if f has one dimension
                                  // otherwise calls spexit.
 
  DynLinArr(const DynArr<T>& f, int n_of_dim,
            // 0 - first dim) 1 - second dim)
            long roc_number);
  // takes only mentioned raw or column.
 
  DynLinArr& assignAll(const T& f) {
    check();
    long n;
    for (n = 0; n < qel; n++)
      el[n] = f;
    return *this;
  }
  template <class X>
  DynLinArr<T>& assignAll1(const X& f)  // assumes that
      // element is object
      // which also accepts assignAll, which is called for it.
      {
    check();
    long n;
    for (n = 0; n < qel; n++)
      el[n].assignAll(f);
    return *this;
  }
  //DynLinArr& operator=(const T& f) { int n; for( n=0; n<qel; n++) el[n]=f; }
  inline T& operator[](long n) {
#ifdef ALR_CHECK_BOUND
    if (n >= 0 && n < qel) {
      return el[n];
    } else {
      mcerr << "ERROR in const T& DynLinArr::operator[](long n) const: "
            << "n is out of bounds, n=" << n << " qel=" << qel << '\n';
      mcerr << "Type of T is (in internal notations) " << typeid(T).name()
            << '\n';
      spexit(mcerr);
#ifdef INS_CRETURN
      return el[0];  // to quiet Microsoft Visial Studio compiler against
                     // "not all control paths return a value"
#endif
    }
#else
    return el[n];
#endif
    /*
#ifdef ALR_CHECK_BOUND
      if(n<0 || n>=qel)
      {
    mcerr<<"ERROR in T& DynLinArr::operator[](long n): "
         <<"n is out of bounds, n="<<n
         <<" qel="<<qel<<'\n';
    mcerr<<"Type of T is (in internal notations) "<<typeid(T).name()<<'\n';
    spexit(mcerr);
      }
#endif
      return el[n];
      */
  }
  inline const T& operator[](long n) const {
#ifdef ALR_CHECK_BOUND
    if (n >= 0 && n < qel) {
      return el[n];
    } else {
      mcerr << "ERROR in const T& DynLinArr::operator[](long n) const: "
            << "n is out of bounds, n=" << n << " qel=" << qel << '\n';
      mcerr << "Type of T is (in internal notations) " << typeid(T).name()
            << '\n';
      spexit(mcerr);
#ifdef INS_CRETURN   // Insert calming return
      return el[0];  // to quiet Microsoft Visial Studio compiler against
                     // "not all control paths return a value"
#endif
    }
#else
    return el[n];
#endif
    /* // it is more simple to read but might be slower to execute
#ifdef ALR_CHECK_BOUND
    if(n<0 || n>=qel)
    {
      mcerr<<"ERROR in const T& DynLinArr::operator[](long n) const: "
       <<"n is out of bounds, n="<<n
       <<" qel="<<qel<<'\n';
      mcerr<<"Type of T is (in internal notations) "<<typeid(T).name()<<'\n';
      spexit(mcerr);
    }
#endif
    return el[n];
    */
  }
  inline T& acu(long n)  // unchecked access
      {
    return el[n];
  }
  inline const T& acu(long n) const  // unchecked access
      {
    return el[n];
  }
  inline T& last_el(void) {
#ifdef ALR_CHECK_BOUND
    if (qel > 0) {
      return el[qel - 1];
    } else {
      mcerr << "ERROR in const T& DynLinArr::last_el(void) const: qel <=0:"
            << " qel" << qel << '\n';
      mcerr << "Type of T is (in internal notations) " << typeid(T).name()
            << '\n';
      spexit(mcerr);
#ifdef INS_CRETURN   // Insert calming return
      return el[0];  // to quiet Microsoft Visial Studio compiler against
                     // "not all control paths return a value"
#endif
    }
#else
    return el[qel];
#endif
  }
 
  inline const T& last_el(void) const {
#ifdef ALR_CHECK_BOUND
    if (qel > 0) {
      return el[qel - 1];
    } else {
      mcerr << "ERROR in const T& DynLinArr::last_el(void) const: qel <=0:"
            << " qel" << qel << '\n';
      mcerr << "Type of T is (in internal notations) " << typeid(T).name()
            << '\n';
      spexit(mcerr);
#ifdef INS_CRETURN   // Insert calming return
      return el[0];  // to quiet Microsoft Visial Studio compiler against
                     // "not all control paths return a value"
#endif
    }
#else
    return el[qel];
#endif
  }
 
  long get_qel(void) const { return qel; }
  void put_qel(long fqel);
  void put_qel(long fqel, const T* val, ArgInterp_SingleAdr t);
  // creates array with size fqel
  // If old array existed, then
  //   If it was less than fqel, it all is copied to new array
  //       and the other elements either assigned *val or
  //       remains not inited.
  //   else its fqel part is copyed to new array.
  // ArgInterp_SingleAdr serves to distinguish this
  // from next function when second argument
  // does not coincide with T but can be converted to it
  // (for example NULL to int and it is not clear which function
  // should be called.
  // Attention!, val is an element, this is assimetry with contructor,
  // in which ar is an array.
 
  void put_qel(long fqel, const T& val) {
    put_qel(fqel, &val, ArgInterp_SingleAdr());
  }
  void increment(const T* val = NULL) {
    check();
    long q = qel + 1;
    put_qel(q, *val);
  }
  void increment(const T& val) {
    check();
    long q = qel + 1;
    put_qel(q, val);
  }
  void clear(void) {
    put_qel(0);
  }  // Not only clears the content,
     // but makes zero dimension.
  void pilfer(PILF_CONST DynLinArr<T>& f) {
#ifdef DEBUG_DYNLINARR
    mcout << "DynLinArr::pilfer is called\n";
#endif
    if (this != &f) {
      if (qel != 0) {
        if (f.qel != 0) {
          mcerr << "ERROR in DynLinArr::pilfer(...):\n";
          mcerr << "Both the destination and source arrays are not empty\n";
          // For explanations why it is dengerous, see similar function
          // of ActivePtr.
          spexit(mcerr);
        } else {
#ifdef USE_REPLACE_ALLOC
          {
            long n;
            for (n = 0; n < qel; n++)
              el[n].~T();
            free(el);
          }
#else
          delete[] el;
#endif
          qel = 0;
        }
      }
      el = f.el;
      qel = f.qel;
      f.el = NULL;
      f.qel = 0;
    }
  }
 
  void check(void) const;
 
  /*
  void print(std::ostream& file, long qpr) const
    {
      Ifile<<"DynLinArr<T>: qel="<<get_qel()<<" qpr="<<qpr<<'\n';
      long n;
      indn.n+=2;
      for( n=0; n<qpr; n++)
      {
    Ifile<<"n="<<n<<" el[n]="<<this->DynLinArr<T>::operator[](n)<<'\n';
      }
      indn.n-=2;
    }
  */
  DynArr<T> top(void);  // transpose the vector, rotate it
  // from vertical colunm to horisontal line for the purpose
  // of linear algebra calculations.
 
  //friend void DLA_sort<T>(DynLinArr<T>& f);
  // Apply any function of one argument (of type T) to each element.
  template <class P> friend void apply1(DynLinArr<P>& ar, void (*fun)(P& f));
  //template<class P> friend void apply1m(DynLinArr<P>& ar,
  //                    void (*P::fun)(void));
  // Apply any function of two arguments
  // (first of which is of type T and the second is of type of address
  // of another function (possibly apply1) to each element.
  template <class P, class X>
  friend void apply2(DynLinArr<P>& ar, void (*fun1)(P& f, void (*fun21)(X& f)),
                     void (*fun2)(X& f));
 
  //template<class P, class X> friend void apply2m(DynLinArr<P>& ar,
  //            void (*fun1)(P& f, void (*fun21)(X& f) ),
  //            void (*X::fun2)(void) );
  //void apply(void (*fun)(T& f))
  //{ long n; for(n=0; n<qel; n++) (*fun)(el[n]); }
 
  // Attention: the both sorts below at large N are much less efficient
  // then sort for standard vectors.
  // At N of the order of 50000 thay take several seconds,
  // whereas sort for vectors takes several tenths of second.
 
  void sort(long q_to_sort = 0);  // sorts first q_to_sort to increasing order.
  // If q_to_sort = 0, it sorts all get_qel() elements.
 
  void sort(DynLinArr<long>& sort_ind, long q_to_sort = 0) const;
  // to increasing order
  // Does not change array, but creates index array which gives access
  // in sorted order. For example:
  // DynLinArr< long > sort_ind;
  // collection.sort(sort_ind);
  // for(n=0; n<q; n++)
  // {
  //   collection_ordered[n] = collection[sort_ind[n]];
  // }
  // Thus content of each n-th element of sort_ind gives
  // position of n-th (by size) element of array.
 
  // finds and sorts minimal q_to_sort elements from all.
  void sort_select_increasing(DynLinArr<long>& sort_ind,
                              long q_to_sort = 0) const;
  void sort_select_decreasing(DynLinArr<long>& sort_ind,
                              long q_to_sort = 0) const;
 
//void sort(DynLinArr< long >& sort_ind, long q_to_sort = 0) const;
 
/*
  The following does not work.
  It will attempt to instantiate it even for non-class types like long etc.
  virtual void print(std::ostream& file, int l) const;
*/
#ifndef DONT_USE_ABSPTR
  macro_copy_header(DynLinArr);
#endif
 
  virtual ~DynLinArr() {
    check();
    if (el != NULL)
#ifdef USE_REPLACE_ALLOC
        {
      long n;
      for (n = 0; n < qel; n++) {
        (&(el[n]))->~T();
      }
      free(el);
    }
#else
    delete[] el;
#endif
  }
 
 private:
  PILF_MUTABLE long qel;  // number of elements, mutable only for pilfer
  PILF_MUTABLE T* el;     // array of qel elements, mutable only for pilfer
//(regarding mutable and pilfer see ActivePtr for more comments).
#ifdef USE_REPLACE_ALLOC
  macro_alloc
#endif
};
#ifndef DONT_USE_ABSPTR
template <class T>
macro_copy_body(DynLinArr<T>)
#endif
    template <class T>
void apply1(DynLinArr<T>& ar, void (*fun)(T& f)) {
  long n;
  for (n = 0; n < ar.qel; n++)
    (*fun)(ar.el[n]);
}
 
//template<class T>
//void apply1m(DynLinArr<T>& ar, void (*T::fun)(void))
//{ long n; for(n=0; n<ar.qel; n++) ar.el[n].(*fun)(); }
 
template <class T, class X>
void apply2(DynLinArr<T>& ar, void (*fun1)(T& f, void (*fun21)(X& f)),
            void (*fun2)(X& f)) {
  long n;
  for (n = 0; n < ar.qel; n++)
    (*fun1)(ar.el[n], fun2);
}
 
//template<class T, class X>
//void apply2m(DynLinArr<T>& ar,
//       void (*fun1)(T& f, void (*fun21)(X& f) ),
//       void (*X::fun2)(void) )
//{ long n; for(n=0; n<ar.qel; n++) (*fun1)(ar.el[n], fun2); }
 
template <class T> void DynLinArr<T>::check(void) const {
  if (qel < 0) {
    mcerr << "ERROR in template<class T> void DynLinArr<T>::check(void):\n";
    mcerr << "qel < 0, qel=" << qel << '\n';
    mcerr << "Type of T is (in internal notations) " << typeid(T).name()
          << '\n';
    spexit(mcerr);
  }
  if (qel == 0 && el != NULL) {
    mcerr << "ERROR in template<class T> void DynLinArr<T>::check(void):\n";
    //mcerr<<"qel == 0 && el != NULL: long(el)="<<long(el)<<'\n';
    mcerr << "qel == 0 && el != NULL: el=" << el << '\n';
    mcerr << "Type of T is (in internal notations) " << typeid(T).name()
          << '\n';
    spexit(mcerr);
  }
  if (qel > 0) {
    if (el == NULL) {
      mcerr << "ERROR in template<class T> void DynLinArr<T>::check(void):\n";
      mcerr << "qel > 0 && el == NULL: qel=" << qel << '\n';
      mcerr << "Type of T is (in internal notations) " << typeid(T).name()
            << '\n';
      spexit(mcerr);
    }
    if (qel > max_qel_DynLinArr) {
      mcerr << "ERROR in template<class T> void DynLinArr<T>::check(void):\n";
      mcerr << "qel > max_qel_DynLinArr: \n";
      Iprint2n(mcout, qel, max_qel_DynLinArr);
      mcerr << "Type of T is (in internal notations) " << typeid(T).name()
            << '\n';
      spexit(mcerr);
    }
  }
}
 
template <class T>
DynLinArr<T>& DynLinArr<T>::operator=(const DynLinArr<T>& f) {
#ifdef DEBUG_DYNLINARR
  mcout << "DynLinArr<T>& DynLinArr<T>::operator=(const DynLinArr<T>& f) is "
           "working\n";
#endif
  if (this != &f) {
    //mcout<<"DynLinArr<T>& operator=(const DynLinArr<T>& f): long(el)="
    //<<long(el)<<'\n';
    check();
    f.check();
    // First of all we allocate memory and copy to there
    // all elements from f, because f can be part of this array,
    // for example, one of its elements.
    long q = f.get_qel();
#ifdef USE_REPLACE_ALLOC
    T* temp_el;
    if (q > 0) {
      temp_el = (T*)malloc(sizeof(T) * q);
      //long n; for(n=0; n<q; n++) new(&(temp_el[n])) T;
      long n;
      for (n = 0; n < q; n++)
        new (&(temp_el[n])) T(f.el[n]);  // here
      // we call copy constructor instead of creating default as for above
      // and then to copy to it.
    } else
      temp_el = NULL;
#else
    //mcout<<"operator=:  sizeof(T)="<<std::setw(10)<<sizeof(T)<<" q="<<q<<'\n';
    //T* temp_el = (q > 0) ? (new T[q]) : (T*)NULL;
    T* temp_el = (T*)NULL;
    if (q > 0) {
      //mcout<<"operator=:  sizeof(T)="<<std::setw(10)<<sizeof(T)<<"
      //q="<<q<<'\n';
      temp_el = new T[q];
      long n;
      for (n = 0; n < q; n++)
        temp_el[n] = f.el[n];
    }
#endif
    //long n; for( n=0; n<q; n++) temp_el[n]=f.el[n];
    /*
#ifdef USE_REPLACE_ALLOC
    for(n=0; n<q; n++) el[n].~T();
    free(el);
#else
    delete[] el;
#endif
    qel = q;
    el = temp_el;
    */
    pass(q, temp_el);
    //el = (qel > 0) ? (new T[qel]) : (T*)NULL;
    //long n; for( n=0; n<qel; n++) el[n]=f.el[n];
  }
  return *this;
}
 
template <class T>
template <class D>
DynLinArr<T>& DynLinArr<T>::operator=(const DynLinArr<D>& f) {
#ifdef DEBUG_DYNLINARR
  mcout << "DynLinArr<T>& DynLinArr<T>::operator=(const DynLinArr<D>& f) is "
           "working\n";
#endif
  //if(this != &f)
  //{
  //mcout<<"DynLinArr<T>& operator=(const DynLinArr<T>& f): long(el)="
  //<<long(el)<<'\n';
  check();
  f.check();
  // First of all we allocate memory and copy to there
  // all elements from f, because f can be part of this array,
  // for example, one of its elements.
  long q = f.get_qel();
#ifdef USE_REPLACE_ALLOC
  T* temp_el;
  if (q > 0) {
    temp_el = (T*)malloc(sizeof(T) * q);
    //long n; for(n=0; n<q; n++) T* t = new(&(temp_el[n])) T;
    // What is this (above)? Perhaps I did not debug this.
 
    long n;
    for (n = 0; n < q; n++)
      new (&(temp_el[n])) T(f.el[n]);  // here
    // we call copy constructor instead of creating default as for above
    // and then to copy to it.
  } else
    temp_el = NULL;
#else
  //T* temp_el = (q > 0) ? (new T[q]) : (T*)NULL;
  T* temp_el = (T*)NULL;
  if (q > 0) {
    //mcout<<"operator=:  sizeof(T)="<<std::setw(10)<<sizeof(T)<<" q="<<q<<'\n';
    temp_el = new T[q];
    //long n; for( n=0; n<q; n++) temp_el[n]=f.el[n];
    long n;
    for (n = 0; n < q; n++)
      temp_el[n] = f[n];
  }
 
#endif
  //long n; for( n=0; n<q; n++) temp_el[n] = f.acu(n);
  pass(q, temp_el);
  //el = (qel > 0) ? (new T[qel]) : (T*)NULL;
  //long n; for( n=0; n<qel; n++) el[n]=f.el[n];
  //}
  return *this;
}
 
template <class T>
inline DynLinArr<T>::DynLinArr(const DynLinArr<T>& f)
    :
#ifndef DONT_USE_ABSPTR
      RegPassivePtr(),
#endif
      qel(0),
      el(NULL) {
#ifdef DEBUG_DYNLINARR
  mcout << "DynLinArr(const DynLinArr<T>& f) is working\n";
#endif
  *this = f;
}
 
template <class T> void DynLinArr<T>::put_qel(long fqel) {
  //
  // creates array with size fqel
  // If old array existed, then
  //   If it was less than fqel, it all is copied to new array
  //       and the other elements either assigned *val or
  //       remains not inited.
  //   else its fqel part is copyed to new array.
  //mcout<<"put_qel: *this="<<(*this);
  if (fqel < 0) {
    mcerr << "ERROR in template<class T> void DynLinArr<T>::put_qel(long "
             "fqel):\n";
    mcerr << "fqel < 0, fqel=" << fqel << '\n';
    mcerr << "Type of T is (in internal notations) " << typeid(T).name()
          << '\n';
    spexit(mcerr);
  }
  check();
  if (el == NULL) {
    qel = fqel;
    if (qel > 0) {
#ifdef USE_REPLACE_ALLOC
      el = (T*)malloc(sizeof(T) * fqel);
      for (long n = 0; n < fqel; ++n)
        new (&(el[n])) T;
#else
      el = new T[fqel];
#endif
    }
  } else {
    if (qel != fqel) {
      if (fqel <= 0) {
        qel = 0;
#ifdef USE_REPLACE_ALLOC
        for (long n = 0; n < qel; ++n)
          el[n].~T();
        free(el);
#else
        delete[] el;
#endif
        el = NULL;
      } else {
        T* elh;
#ifdef USE_REPLACE_ALLOC
        elh = (T*)malloc(sizeof(T) * fqel);
        for (n = 0; n < fqel; ++n)
          new (&(elh[n])) T;
#else
        elh = new T[fqel];  // long q = find_min(qel,fqel);
#endif
        for (long n = 0; n < fqel; ++n) {
          if (n < qel) elh[n] = el[n];
        }
#ifdef USE_REPLACE_ALLOC
        for (long n = 0; n < qel; ++n)
          el[n].~T();
        free(el);
#else
        delete[] el;
#endif
        el = elh;
        qel = fqel;
      }
    }
  }
}
 
template <class T>
void DynLinArr<T>::put_qel(long fqel, const T* val, ArgInterp_SingleAdr t) {
  // By default val == NULL
  // creates array with size fqel
  // If old array existed, then
  //   If it was less than fqel, it all is copied to new array
  //       and the other elements either assigned *val or
  //       remains not inited.
  //   else its fqel part is copyed to new array.
  //mcout<<"put_qel: *this="<<(*this);
  if (fqel < 0) {
    mcerr << "ERROR in template<class T> void DynLinArr<T>::put_qel(long fqel, "
             "const T* val, ArgInterp_SingleAdr):\n";
    mcerr << "fqel < 0, fqel=" << fqel << '\n';
    mcerr << "Type of T is (in internal notations) " << typeid(T).name()
          << '\n';
    spexit(mcerr);
    // Avoid compiler warning because of unused variable t (HS).
    mcerr << sizeof(t) << "\n";
  }
  check();
  if (el == NULL) {
    qel = fqel;
    if (qel > 0) {
#ifdef USE_REPLACE_ALLOC
      el = (T*)malloc(sizeof(T) * fqel);
      for (long n = 0; n < fqel; ++n)
        new (&(el[n])) T;
#else
      el = new T[fqel];
#endif
    }
    if (val != NULL) for (long n = 0; n < qel; ++n)
        el[n] = *val;
  } else {
    if (qel != fqel) {
      if (fqel <= 0) {
        qel = 0;
#ifdef USE_REPLACE_ALLOC
        for (long n = 0; n < qel; ++n)
          el[n].~T();
        free(el);
#else
        delete[] el;
#endif
        el = NULL;
      } else {
        T* elh;
#ifdef USE_REPLACE_ALLOC
        elh = (T*)malloc(sizeof(T) * fqel);
        for (long n = 0; n < fqel; ++n)
          new (&(elh[n])) T;
#else
        elh = new T[fqel];  // long q = find_min(qel,fqel);
#endif
        for (long n = 0; n < fqel; ++n) {
          if (n < qel)
            elh[n] = el[n];
          else if (val != NULL)
            elh[n] = *val;
        }
#ifdef USE_REPLACE_ALLOC
        for (long n = 0; n < qel; ++n)
          el[n].~T();
        free(el);
#else
        delete[] el;
#endif
        el = elh;
        qel = fqel;
      }
    }
  }
}
 
// Simple sorting routine optimized for DynLiArr. )ptimized means
//not doing checks at each indexing, but it is not the best for the large
//number of N.
// Note that it can be made better if at choosing n_possible_next
// to arrange going backward. But this is much more complicated so currently
// I am not going to do this.
 
template <class T> void DynLinArr<T>::sort(long q_to_sort) {
  mfunnamep("void DynLinArr<T>::sort(long q_to_sort = 0)");
 
  check_econd12(q_to_sort, >, qel, mcerr);
  if (q_to_sort <= 0) q_to_sort = qel;
  if (q_to_sort <= 1) return;
 
  long n_possible_next = 1;
  long q_comp = 0;
  long n, m;
  for (n = 0; n < q_to_sort - 1; n++) {
    //Iprint2n(mcout, n, n_possible_next);
    // first it finds the minimum along the rest and replaces if it is less
    //long nmin = n+1;
    long nmin = n_possible_next;
    T el_min = el[nmin];
    int s_change_possible_next = 0;
 
    //for(m=n+2; m<q_to_sort; m++)
    for (m = n_possible_next + 1; m < q_to_sort; m++) {
      q_comp++;
      //if(el[nmin] > el[m])
      if (el_min > el[m]) {
        n_possible_next = nmin;
        s_change_possible_next = 1;
        nmin = m;
        el_min = el[nmin];
      }
    }
    //Iprint4n(mcout, n_possible_next, s_change_possible_next, nmin, el_min);
    if (s_change_possible_next == 0 || n_possible_next < n + 2) {
      n_possible_next = n + 2;
    }
    //if(el[n] > el[nmin])
    //{
    //  T t = el[nmin];
    //  el[nmin] = el[n];
    //  el[n] = t;
    //}
    //Iprintn(mcout, n_possible_next);
    //Iprint2n(mcout, el[n], el_min);
    if (el[n] > el_min) {
      if (s_change_possible_next == 1) {
        //if(n_possible_next < q_to_sort)
        if (n_possible_next < q_to_sort && el[n] < el[n_possible_next]) {
          n_possible_next = nmin;
        }
      }
      //mcout<<"replacing el[n] and el_min\n";
      T t = el_min;
      el[nmin] = el[n];
      el[n] = t;
    }
    //Iprintn(mcout, (*this));
 
  }
  //Iprintn(mcout, q_comp);
}
 
// New variant, should be faster, the old is below.
 
template <class T>
void DynLinArr<T>::sort(DynLinArr<long>& sort_ind, long q_to_sort) const {
  mfunnamep("void DynLinArr<T>::sort(DynLinArr< long >& sort_ind, long "
            "q_to_sort = 0) const");
 
  check_econd12(q_to_sort, >, qel, mcerr);
  if (q_to_sort <= 0) q_to_sort = qel;
  //if(q_to_sort <= 1) return;
 
  sort_ind.clear();
  sort_ind.pilfer(DynLinArr<long>(q_to_sort));
  long n, m;
  for (n = 0; n < q_to_sort; n++) {
    sort_ind.acu(n) = n;
  }
  if (q_to_sort <= 1) return;
 
  long n_possible_next = 1;
 
  for (n = 0; n < q_to_sort - 1; n++) {
    // first it finds the minimum along the rest and replaces if it is less
    long nmin = n_possible_next;
    long ind_nmin = sort_ind.acu(nmin);
    int s_change_possible_next = 0;
 
    //for(m=n+2; m<q_to_sort; m++)
    for (m = n_possible_next + 1; m < q_to_sort; m++) {
      if (el[ind_nmin] > el[sort_ind.acu(m)])
          //if(el[sort_ind.acu(nmin)] > el[sort_ind.acu(m)])
          {
        n_possible_next = nmin;
        s_change_possible_next = 1;
        nmin = m;
        ind_nmin = sort_ind.acu(nmin);
      }
    }
    if (s_change_possible_next == 0 || n_possible_next < n + 2) {
      n_possible_next = n + 2;
    }
    if (el[sort_ind.acu(n)] > el[ind_nmin])
        //if(el[sort_ind.acu(n)] > el[sort_ind.acu(nmin)])
        {
      if (s_change_possible_next == 1) {
        if (n_possible_next < q_to_sort &&
            el[sort_ind.acu(n)] < el[sort_ind.acu(n_possible_next)]) {
          n_possible_next = nmin;
        }
      }
      //long t = sort_ind.acu(nmin);
      sort_ind.acu(nmin) = sort_ind.acu(n);
      sort_ind.acu(n) = ind_nmin;
    }
  }
}
 
template <class T>
void DynLinArr<T>::sort_select_increasing(DynLinArr<long>& sort_ind,
                                          long q_to_sort) const {
  mfunnamep("void DynLinArr<T>::sort_select_increasing(DynLinArr< long >& "
            "sort_ind, long q_to_sort = 0) const");
 
  check_econd12(q_to_sort, >, qel, mcerr);
  long s_last_noninc = 0;
  if (q_to_sort <= 0) {
    q_to_sort = qel;
    s_last_noninc = 1;
  } else if (q_to_sort == qel) {
    s_last_noninc = 1;
  }
 
  //if(qel <= 1) return;
  sort_ind.clear();
  sort_ind.pilfer(DynLinArr<long>(qel));
  long n, m;
  for (n = 0; n < qel; n++) {
    sort_ind[n] = n;
  }
  if (q_to_sort <= 1) return;
 
  long n_possible_next = 1;
 
  for (n = 0; n < q_to_sort - s_last_noninc; n++) {
    // first it finds the minimum along the rest and replaces if it is less
    long nmin = n_possible_next;
    long ind_nmin = sort_ind[nmin];
    int s_change_possible_next = 0;
 
    //for(m=n+2; m<q_to_sort; m++)
    for (m = n_possible_next + 1; m < qel; m++) {
      if (el[ind_nmin] > el[sort_ind[m]])
          //if(el[sort_ind.acu(nmin)] > el[sort_ind.acu(m)])
          {
        n_possible_next = nmin;
        s_change_possible_next = 1;
        nmin = m;
        ind_nmin = sort_ind[nmin];
      }
    }
    if (s_change_possible_next == 0 || n_possible_next < n + 2) {
      n_possible_next = n + 2;
    }
    if (el[sort_ind[n]] > el[ind_nmin])
        //if(el[sort_ind.acu(n)] > el[sort_ind.acu(nmin)])
        {
      if (s_change_possible_next == 1) {
        if (n_possible_next < q_to_sort &&
            el[sort_ind[n]] < el[sort_ind[n_possible_next]]) {
          n_possible_next = nmin;
        }
      }
      //long t = sort_ind.acu(nmin);
      sort_ind[nmin] = sort_ind[n];
      sort_ind[n] = ind_nmin;
    }
  }
  sort_ind.put_qel(qel);
}
 
template <class T>
void DynLinArr<T>::sort_select_decreasing(DynLinArr<long>& sort_ind,
                                          long q_to_sort) const {
  mfunnamep("void DynLinArr<T>::sort_select_decreasing(DynLinArr< long >& "
            "sort_ind, long q_to_sort = 0) const");
 
  check_econd12(q_to_sort, >, qel, mcerr);
  long s_last_noninc = 0;
  if (q_to_sort <= 0) {
    q_to_sort = qel;
    s_last_noninc = 1;
  } else if (q_to_sort == qel) {
    s_last_noninc = 1;
  }
  //Iprintn(mcout, q_to_sort);
 
  //if(qel <= 1) return;
  sort_ind.clear();
  sort_ind.pilfer(DynLinArr<long>(qel));
  long n, m;
  for (n = 0; n < qel; n++) {
    sort_ind[n] = n;
  }
  if (q_to_sort <= 1) return;
 
  long n_possible_next = 1;
 
  for (n = 0; n < q_to_sort - s_last_noninc; n++) {
    //Iprintn(mcout, n);
    // first it finds the minimum along the rest and replaces if it is less
    long nmax = n_possible_next;
    //Iprintn(mcout, nmax);
    long ind_nmax = sort_ind[nmax];
    int s_change_possible_next = 0;
 
    //for(m=n+2; m<q_to_sort; m++)
    for (m = n_possible_next + 1; m < qel; m++) {
      //Iprint3n(mcout, ind_nmax, m, sort_ind[m]);
      if (el[ind_nmax] < el[sort_ind[m]]) {
        n_possible_next = nmax;
        s_change_possible_next = 1;
        nmax = m;
        //Iprintn(mcout, nmax);
        ind_nmax = sort_ind[nmax];
      }
    }
    if (s_change_possible_next == 0 || n_possible_next < n + 2) {
      n_possible_next = n + 2;
    }
    //Iprint4n(mcout, n, sort_ind[n], ind_nmax, el[ind_nmax]);
    if (el[sort_ind[n]] < el[ind_nmax]) {
      if (s_change_possible_next == 1) {
        if (n_possible_next < q_to_sort &&
            el[sort_ind[n]] > el[sort_ind[n_possible_next]]) {
          n_possible_next = nmax;
        }
      }
      //long t = sort_ind.acu(nmin);
      sort_ind[nmax] = sort_ind[n];
      sort_ind[n] = ind_nmax;
    }
  }
  sort_ind.put_qel(qel);
}
 
/*
template<class T>
void DynLinArr<T>::sort(DynLinArr< long >& sort_ind, long q_to_sort) const
{
  mfunnamep("void DynLinArr<T>::sort(DynLinArr< long >& sort_ind, long q_to_sort
= 0) const");
 
  check_econd12(q_to_sort , > , qel , mcerr);
  if(q_to_sort <= 0) q_to_sort = qel;
  if(q_to_sort <= 1) return;
 
  //if(qel <= 1) return;
  sort_ind.clear();
  sort_ind.pilfer(DynLinArr< long >(q_to_sort));
  long n,m;
  for(n=0; n<q_to_sort; n++)
  {
    sort_ind.acu(n) = n;
  }
 
  for(n=0; n<q_to_sort-1; n++)
  {
    // first it finds the minimum along the rest and replaces if it is less
    long nmin = n+1;
    long ind_nmin = sort_ind.acu(nmin);
    for(m=n+2; m<q_to_sort; m++)
    {
      if(el[ind_nmin] > el[sort_ind.acu(m)])
      //if(el[sort_ind.acu(nmin)] > el[sort_ind.acu(m)])
      {
    nmin = m;
    ind_nmin = sort_ind.acu(nmin);
      }
    }
    if(el[sort_ind.acu(n)] > el[ind_nmin])
      //if(el[sort_ind.acu(n)] > el[sort_ind.acu(nmin)])
    {
      //long t = sort_ind.acu(nmin);
      sort_ind.acu(nmin) = sort_ind.acu(n);
      sort_ind.acu(n) = ind_nmin;
    }
  }
}
*/
 
/*
template<T> void DLA_sort(DynLinArr<T>& f)
{
  mfunnamep("template<T> void DLA_sort(DynLinArr<T>& f)");
 
  long q = f.get_qel();
  if(q <= 1) return;
  long n,m;
  T* a = &(f[0]);
  for(n=0; n<q-1; n++)
  {
    for(m=n+1; m<q; m++)
    {
      if(a[n] > a[m])
      {
    T t = a[m];
    a[m] = a[n];
    a[n] = t;
      }
    }
  }
}
*/
/*
//Somewhy this doen not work prolerly if T is template class itself
//AbsCont_ts.c:167: type unification failed for function template
//`template <class T> long int append_TreeNode(TreeNode<...> &,
//TreeNode<...> &, long int &, T * = 0, long int = 0)'
*/
/*
This is not compiled in Solaris:
Error: Function templates may not have default parameters.
*/
#ifndef BAN_DEFAULT_PAR_FUN_TEMPL
template <class T>
long append(const T& t,         // value to assign
            DynLinArr<T>& dla,  // array to whose elements this value is to be
            // assigned
            long& qael,  //input: index of element to which it will be assigned
            // if qael == get_qel(), the qel will be increased
            // to either new_qel or by 3 times.
            // But if qael > get_qel(), it is considered as error.
            //output: the next index
            T* tempt = NULL,  // data filled to extra elements
            long new_qel = 0  // recommended new size
            ) {
  if (dla.get_qel() < qael) {
    mcerr << "ERROR in long append(class DynLinArr& dla, ...): dla.get_qel() < "
             "qael\n"
          << "dla.get_qel()=" << dla.get_qel() << " qael=" << qael << '\n';
    mcerr << "Type of T is (in internal notations) " << typeid(T).name()
          << '\n';
    spexit(mcerr);
  }
  if (dla.get_qel() == qael) {
    if (new_qel <= qael) new_qel = find_max(3 * qael, long(3));
    dla.put_qel(new_qel, tempt, ArgInterp_SingleAdr());
  }
  dla[qael++] = t;
  return qael;
}
 
#else
 
//long append(const T& t, DynLinArr<T>& dla, long& qael,
//              T* tempt=NULL, long new_qel=0  )
#define append_to_DynLinArr(t, dla, qael, tempt, new_qel)                     \
  {                                                                           \
    if (dla.get_qel() < qael) {                                               \
      mcerr << "append(class DynLinArr& dla, ...): dla.get_qel() < qael, "    \
            << "dla.get_qel()=" << dla.get_qel() << " qael=" << qael << '\n'; \
      spexit(mcerr);                                                          \
    }                                                                         \
    int nn = new_qel;                                                         \
    if (dla.get_qel() == qael) {                                              \
      if (new_qel <= qael) nn = find_max(3 * qael, long(3));                  \
      dla.put_qel(nn, tempt);                                                 \
    }                                                                         \
    dla[qael++] = t;                                                          \
  }
#define append_to_DynLinArr_short(t, dla, qael)                               \
  {                                                                           \
    if (dla.get_qel() < qael) {                                               \
      mcerr << "append(class DynLinArr& dla, ...): dla.get_qel() < qael, "    \
            << "dla.get_qel()=" << dla.get_qel() << " qael=" << qael << '\n'; \
      spexit(mcerr);                                                          \
    }                                                                         \
    int nn = 0;                                                               \
    if (dla.get_qel() == qael) {                                              \
      nn = find_max(3 * qael, long(3));                                       \
      dla.put_qel(nn, NULL);                                                  \
    }                                                                         \
    dla[qael++] = t;                                                          \
  }
#endif
 
/*
template<class T, int >
DynLinArr<T>& convert_Stat_to_Dyn(const StatLinArr<T,fqel>& fs,
                  DynLinArr<T>& fd)
{
  fd.put_qel(fqel);
  long n; for( n=0; n<fqel; n++) el[n] = f[n];
  return fd;
}
*/
 
template <class T>
    std::ostream& operator<<(std::ostream& file, const DynLinArr<T>& f) {
  //mfunnamep("template<class T> std::ostream& operator<<(std::ostream& file,
  //const DynLinArr<T>& f)");
  //mcout<<"operator<<(std::ostream& file, const DynLinArr<T>& f) is started\n";
  Ifile << "DynLinArr<T>: qel=" << f.get_qel() << '\n';
  f.check();
  long n;
  indn.n += 2;
  for (n = 0; n < f.get_qel(); n++) {
    //Ifile<<"n="<<n<<" el[n]="<<noindent<<f[n]<<yesindent<<'\n';
    if (s_short_output == 0) {
      Ifile << "n=" << n << " el[n]=";
    }
    std::ostringstream ost;
    ost << indn << noindent << f[n] << yesindent;
    put_one_n(ost);
    file << ost.str();
  }
  //file<<yesindent;
  indn.n -= 2;
  return file;
}
 
#include "wcpplib/stream/definp.h"
 
template <class T>
    std::istream& operator>>(std::istream& file, DynLinArr<T>& f) {
  mfunnamep(
      "template<class T> istream& operator>>(istream& file, DynLinArr<T>& f)");
  //mcout<<"operator<<(std::ostream& file, const DynLinArr<T>& f) is started\n";
  definp_endpar dep(&file, 0, 1, 0, s_short_output);
  long qel = 0;
  dep.s_short = 0;
  DEFINPAP(qel);
  dep.s_short = s_short_output;
  check_econd11(qel, < 0, mcerr);
  f.clear();
  f.put_qel(qel);
  long fn;
  for (fn = 0; fn < qel; fn++) {
    if (s_short_output == 0) {
      long n = 0;
      DEFINPAP(n);
      check_econd12(fn, !=, n, mcerr);
    }
    //set_position("el[n]=", *dep.istrm, dep.s_rewind, dep.s_req_sep);
    //file >> f[n];
    definp_any_par(f[fn], "el[n]=", dep);
  }
  return file;
}
 
// Commented out unused function (hschindl)
/*
template<class T>
void output_DynLinArr(std::ostream& file, const DynLinArr<T>& f, int l, long q)
{
  //mfunnamep("template<class T> void output_DynLinArr(std::ostream& file, const
DynLinArr<T>& f, int l, long q)");
  Ifile<<"DynLinArr<T>: qel="<<f.get_qel()<<" q to print is "<<q<<'\n';
  f.check();
  if(q>f.get_qel())
  {
    mcerr<<"output_DynLinArr(...): q>f.get_qel(), q="<<q
     <<" f.get_qel()="<<f.get_qel()<<'\n';
    mcerr<<"Type of T is (in internal notations) "<<typeid(T).name()<<'\n';
    spexit(mcerr);
  }
  long n;
  indn.n+=2;
  for( n=0; n<q; n++)
  {
    //Ifile<<"n="<<n<<" el[n]="<<noindent<<f[n]<<yesindent<<'\n';
    Ifile<<"n="<<n<<" el[n]="<<noindent;
    std::ostringstream ost;
    ost<<f[n]<<yesindent;
    put_one_n(ost);
    file<<ost.str();
  }
  //file<<yesindent;
  indn.n-=2;
}
*/
 
template <class T>
void print_DynLinArr(std::ostream& file, const DynLinArr<T>& f, int l) {
  //mfunnamep("template<class T> void print_DynLinArr(std::ostream& file, const
  //DynLinArr<T>& f, int l)");
  Ifile << "DynLinArr<T>: qel=" << f.get_qel() << '\n';
  f.check();
  long n;
  indn.n += 2;
  for (n = 0; n < f.get_qel(); n++) {
    //Ifile<<"n="<<n<<" el[n]="<<noindent; f[n].print(file, l);
    //file<<yesindent;
    Ifile << "n=" << n << " el[n]=" << noindent;
    std::ostringstream ost;
    f[n].print(ost, l);
    ost << yesindent;
    put_one_n(ost);
    file << ost.str();
  }
  indn.n -= 2;
}
/*
Does not work.
It will attempt to instantiate it even for non-class types like long etc.
template<class T>
void DynLinArr<T>::print(std::ostream& file, int l) const
{
  mfunnamep("template<class T> void DynLinArr::print(std::ostream& file, int
l)");
  Ifile<<"DynLinArr<T>: qel="<<get_qel()<<'\n';
  check();
  long n;
  indn.n+=2;
  for( n=0; n<get_qel(); n++)
  {
    Ifile<<"n="<<n<<" el[n]="<<noindent; this->operator[](n).print(file, l);
    file<<yesindent;
  }
  indn.n-=2;
}
*/
template <class T>
void print_DynLinArr(std::ostream& file, const DynLinArr<T>& f, int l, long q) {
  //mfunnamep("template<class T> void print_DynLinArr(std::ostream& file, const
  //DynLinArr<T>& f, int l, long q)");
  Ifile << "DynLinArr<T>: qel=" << f.get_qel() << " q to print is " << q
        << '\n';
  f.check();
  if (q > f.get_qel()) {
    mcerr << "print_DynLinArr(...): q>f.get_qel(), q=" << q
          << " f.get_qel()=" << f.get_qel() << '\n';
    mcerr << "Type of T is (in internal notations) " << typeid(T).name()
          << '\n';
    spexit(mcerr);
  }
  long n;
  indn.n += 2;
  for (n = 0; n < q; n++) {
    //Ifile<<"n="<<n<<" el[n]="<<noindent; f[n].print(file, l);
    //file<<yesindent;
    Ifile << "n=" << n << " el[n]=" << noindent;
    std::ostringstream ost;
    f[n].print(ost, l);
    ost << yesindent;
    put_one_n(ost);
    file << ost.str();
  }
  indn.n -= 2;
}
 
template <class T>
void print_adr_DynLinArr(std::ostream& file, const DynLinArr<T>& f, int l,
                         long q) {
  //mfunnamep("template<class T> void print_adr_DynLinArr(std::ostream& file,
  //const DynLinArr<T>& f, int l, long q)");
  Ifile << "DynLinArr<T>: qel=" << f.get_qel() << " q to print is " << q
        << '\n';
  f.check();
  if (q > f.get_qel()) {
    mcerr << "print_adr_DynLinArr(...): q>f.get_qel(), q=" << q
          << " f.get_qel()=" << f.get_qel() << '\n';
    mcerr << "Type of T is (in internal notations) " << typeid(T).name()
          << '\n';
    spexit(mcerr);
  }
  long n;
  indn.n += 2;
  for (n = 0; n < q; n++) {
    //Ifile<<"n="<<n<<" el[n]="<<noindent; f[n]->print(file, l);
    //file<<yesindent;
    Ifile << "n=" << n << " el[n]=" << noindent;
    std::ostringstream ost;
    f[n]->print(ost, l);
    ost << yesindent;
    put_one_n(ost);
    file << ost.str();
  }
  indn.n -= 2;
}
 
const int pq_arrelem_in_line = 5;
 
void print_DynLinArr_int(std::ostream& file, const DynLinArr<int>& f);
void print_DynLinArr_long(std::ostream& file, const DynLinArr<long>& f);
void print_DynLinArr_float(std::ostream& file, const DynLinArr<float>& f);
void print_DynLinArr_double(std::ostream& file, const DynLinArr<double>& f);
// See AbsArrD for similar function with DoubleAc
 
void print_DynLinArr_double2(std::ostream& file, const DynLinArr<double>& f1,
                             const DynLinArr<double>& f2);
// Print two arrays in two colums side-by-side.
// Good for arrays of equal dimensions.
 
void print_DynLinArr_int_double(std::ostream& file, const DynLinArr<int>& iar,
                                const DynLinArr<double>& dar);
 
void print_DynLinArr_int_double3(std::ostream& file, const DynLinArr<int>& iar,
                                 const DynLinArr<double>& dar1,
                                 const DynLinArr<double>& dar2,
                                 const DynLinArr<double>& dar3);
// Print 4 arrays in two colums side-by-side.
// Good for arrays of equal dimensions.
 
#define Iprintdla_int(file, name)           \
  file << indn << #name << "=" << noindent; \
  print_DynLinArr_int(file, name);
#define Iprintdla_long(file, name)          \
  file << indn << #name << "=" << noindent; \
  print_DynLinArr_long(file, name);
#define Iprintdla_float(file, name)         \
  file << indn << #name << "=" << noindent; \
  print_DynLinArr_float(file, name);
#define Iprintdla_double(file, name)        \
  file << indn << #name << "=" << noindent; \
  print_DynLinArr_double(file, name);
// See AbsArrD for similar function with DoubleAc
 
template <class T, class X> void copy_DynLinArr(const T& s, X& d) {
  mfunnamep("template<class T, class X> void copy_DynLinArr(const T& s, X& d)");
  s.check();
  d.check();
  long q = s.get_qel();
  d.put_qel(q);
  long n;
  for (n = 0; n < q; n++) {
    d[n] = s[n];
  }
}
 
// Covert to another compatible type
template <class T, class X> void convert_DynLinArr(const T& s, X& d) {
  long q = s.get_qel();
  d.put_qel(q);
  long n;
  for (n = 0; n < q; n++) {
    d[n] = X(s[n]);
  }
}
 
template <class T>
void put_qel_1(DynLinArr<T>& f, long fq)  // change dimensions of arrays
    // which are elements of the main one
    {
  long q = f.get_qel();
  long n;
  for (n = 0; n < q; n++)
    f[n].put_qel(fq);
}
 
template <class T, class T1>
void assignAll_1(DynLinArr<T>& f, const T1& ft)  // assign ft to all elements
    // of arrays which are elements of the main one
    {
  long q = f.get_qel();
  long n;
  for (n = 0; n < q; n++)
    f[n].assignAll(ft);
}
 
template <class T> class IterDynLinArr {
 public:
  IterDynLinArr(DynLinArr<T>* fdar) : dar(fdar), ncur(-1) { ; }
  T* more(void) {
    if (ncur < dar->get_qel() - 1)
      return &((*dar)[++ncur]);
    else
      ncur = dar->get_qel();
    return NULL;
  }
  T* current(void) {
    if (ncur >= 0 || ncur < dar->get_qel())
      return &((*dar)[ncur]);
    else
      return NULL;
  }
  T* less(void)  // switch current to previous
      {
    if (ncur >= 1)
      return &((*dar)[--ncur]);
    else
      ncur = -1;
    return NULL;
  }
  long get_ncur(void) { return ncur; }
 
 private:
  long ncur;
  DynLinArr<T>* dar;
};
 
int gconfirm_ind(const DynLinArr<long>& qel, const DynLinArr<long>& ind);
int gconfirm_ind_ext(const DynLinArr<long>& qel, const DynLinArr<long>& ind);
 
// The following function checks whether the arrays or their first qfirst
// elements are equal.
// In the first case this means that their dimensions and all elements
// should be equal to each other.
// In the second case if either of the dimensions is less than qfirst,
// the result is negative (0), so the arrays are considered different.
// If the dimensions are equal or more than qfirst,
// the function compares that number of first elements.
// If to put here just equal, the program may not be compiled since
// instead of this function the program will insert
//template <class InputIterator1, class InputIterator2>
//inline bool equal(InputIterator1 first1, InputIterator1 last1,
//        InputIterator2 first2) {
// Therefore I substituted it to ifequal.
template <class T>
int ifequal(const DynLinArr<T>& fd1, const DynLinArr<T>& fd2,
            long qfirst = -1) {
  long n;
  if (qfirst == -1) {
    if ((qfirst = fd1.get_qel()) != fd2.get_qel()) return 0;
    //qfirst=fd1.get_qel();
  } else {
    if (qfirst > fd1.get_qel() || qfirst > fd2.get_qel()) return 0;
  }
  //Iprintn(mcout, qfirst);
  for (n = 0; n < qfirst; n++) {
    //Iprint3n(mcout, n, fd1[n], fd2[n]);
    if (!(fd1[n] == fd2[n])) return 0;
  }
  return 1;
}
 
// The same as above, but compares DynLinArr with ordinary array.
template <class T>
int ifequal(const DynLinArr<T>& fd1, const T* fd2, long qfirst = -1) {
  long n;
  if (qfirst == -1) {
    qfirst = fd1.get_qel();
  } else {
    if (qfirst > fd1.get_qel()) return 0;
  }
  for (n = 0; n < qfirst; n++)
    if (!(fd1[n] == fd2[n])) return 0;
  return 1;
}
 
template <class T> int ifequal(T* fd1, T* fd2, long qfirst) {
  long n;
  for (n = 0; n < qfirst; n++)
    if (!(fd1[n] == fd2[n])) return 0;
  return 1;
}
 
template <class T>
DynLinArr<T> merge(const DynLinArr<T>& fd1, long qfd1, const DynLinArr<T>& fd2,
                   long qfd2) {
  long n;
  if (qfd1 < 0) {
    qfd1 = fd1.get_qel();
  }
  if (qfd2 < 0) {
    qfd2 = fd2.get_qel();
  }
  DynLinArr<T> ret(qfd1 + qfd2);
  if (qfd1 + qfd2 == 0) return ret;
  for (n = 0; n < qfd1; n++) {
    ret[n] = fd1[n];
  }
  for (n = 0; n < qfd2; n++) {
    ret[qfd1 + n] = fd2[n];
  }
  return ret;
}
 
// The following might be in String.h, but String.h is included into this
// file, and the whole system cannot work.
String DynLinArr_char_we_to_String(DynLinArr<char>& ar);
 
//template<class T> class IndexingProvider;
 
#ifndef DONT_USE_ABSPTR
template <class T>
class DynArr : public RegPassivePtr
#else
               template <class T>
               class DynArr
#endif
               {
 public:
  // Constructors
  DynArr(void) { ; }
  // For one-dimensional array:
  explicit DynArr(long fqel, T* val = NULL)
      : qel(DynLinArr<long>(1)),
        cum_qel(DynLinArr<long>(1)),
        el(DynLinArr<T>(fqel)) {
#ifdef DEBUG_DYNARR
    mcout << "explicit DynArr(long fqel, T* val=NULL) is called\n";
#endif
    qel[0] = fqel;
    cum_qel[0] = 1;
    if (val != NULL) assignAll(*val);
  }
  DynArr(long fqel, T val, ArgInterp_Val /*t*/)
      : qel(DynLinArr<long>(1)),
        cum_qel(DynLinArr<long>(1)),
        el(DynLinArr<T>(fqel)) {
#ifdef DEBUG_DYNARR
    mcout << "explicit DynArr(long fqel, T* val=NULL) is called\n";
#endif
    qel[0] = fqel;
    cum_qel[0] = 1;
    assignAll(val);
  }
  DynArr(long fqel, const T* ar, ArgInterp_Arr /*t*/)
      : qel(DynLinArr<long>(1)), cum_qel(DynLinArr<long>(1)), el(fqel) {
    qel[0] = fqel;
    cum_qel[0] = 1;
    long n;
    for (n = 0; n < fqel; n++)
      el.acu(n) = ar[n];
  }
  // Another variant for one-dimensional array
  // Attention: if the T is long, this might be mixed with array of dimensions.
  // To avoid this the latter should be accompanied by address, see below.
  DynArr(const DynLinArr<T>& f)
      : qel(DynLinArr<long>(1)), cum_qel(DynLinArr<long>(1)), el(f.get_qel()) {
    qel[0] = f.get_qel();
    cum_qel[0] = 1;
    for (long n = 0; n < qel[0]; n++)
      ac(n) = f[n];
  }
 
  // For two-dimensional array:
  DynArr(long fqel1, long fqel2, T* val = NULL)
      : qel(DynLinArr<long>(2)),
        cum_qel(DynLinArr<long>(2)),
        el(fqel1 * fqel2) {
#ifdef DEBUG_DYNARR
    mcout << "DynArr(long fqel1, long fqel2, T* val=NULL) is called\n";
#endif
    qel[0] = fqel1;
    qel[1] = fqel2;
    cum_qel[0] = fqel2;
    cum_qel[1] = 1;
    if (val != NULL) assignAll(*val);
  }
 
  DynArr(long fqel1, long fqel2, T val, ArgInterp_Val /*t*/)
      : qel(DynLinArr<long>(2)),
        cum_qel(DynLinArr<long>(2)),
        el(fqel1 * fqel2) {
#ifdef DEBUG_DYNARR
    mcout
        << "DynArr(long fqel1, long fqel2, T val, ArgInterp_Val t) is called\n";
#endif
    qel[0] = fqel1;
    qel[1] = fqel2;
    cum_qel[0] = fqel2;
    cum_qel[1] = 1;
    assignAll(val);
  }
 
  // For three-dimensional array:
  DynArr(long fqel1, long fqel2, long fqel3, T* val = NULL)
      : qel(DynLinArr<long>(3)),
        cum_qel(DynLinArr<long>(3)),
        el(fqel1 * fqel2 * fqel3) {
    qel[0] = fqel1;
    qel[1] = fqel2;
    qel[2] = fqel3;
    cum_qel[0] = fqel2 * fqel3;
    cum_qel[1] = fqel3;
    cum_qel[2] = 1;
    if (val != NULL) assignAll(*val);
  }
 
  // For four-dimensional array:
  DynArr(long fqel1, long fqel2, long fqel3, long fqel4, T* val = NULL)
      : qel(DynLinArr<long>(4)),
        cum_qel(DynLinArr<long>(4)),
        el(fqel1 * fqel2 * fqel3 * fqel4) {
    qel[0] = fqel1;
    qel[1] = fqel2;
    qel[2] = fqel3;
    qel[3] = fqel4;
    cum_qel[0] = fqel2 * fqel3 * fqel4;
    cum_qel[1] = fqel3 * fqel4;
    cum_qel[2] = fqel4;
    cum_qel[3] = 1;
    if (val != NULL) assignAll(*val);
  }
 
  // Default value is removed in order to avoid confusions with copy
  // constructor from DynLinArr at T = long.
  // If initialization of values is not necessary, just put NULL as argument.
  // It creates array with structure determined by fqel.
 
  DynArr(const DynLinArr<long>& fqel, T val, ArgInterp_Val /*t*/) : qel(fqel) {
    long qdim = qel.get_qel();
    if (qdim <= 0) return;
    cum_qel.put_qel(qdim);
    long ndim;
    long size = qel[0];
    for (ndim = 1; ndim < qdim; ndim++)
      size *= qel[ndim];
    el.put_qel(size);
    cum_qel[qdim - 1] = 1;
    for (ndim = qdim - 2; ndim >= 0; ndim--)
      cum_qel[ndim] = qel[ndim + 1] * cum_qel[ndim + 1];
    assignAll(val);
  }
 
  explicit DynArr(const DynLinArr<long>& fqel, T* val) : qel(fqel) {
    // fqel: array of dimensions
    // val: address of value to fill, may be NULL
    long qdim = qel.get_qel();
    if (qdim <= 0) return;
    cum_qel.put_qel(qdim);
    long ndim;
    long size = qel[0];
    for (ndim = 1; ndim < qdim; ndim++)
      size *= qel[ndim];
    el.put_qel(size);
    cum_qel[qdim - 1] = 1;
    for (ndim = qdim - 2; ndim >= 0; ndim--)
      cum_qel[ndim] = qel[ndim + 1] * cum_qel[ndim + 1];
    if (val != NULL) assignAll(*val);
  }
 
  DynArr(const DynArr<T>& f)
#ifndef DONT_USE_ABSPTR
      : RegPassivePtr()
#endif
        {
#ifdef DEBUG_DYNARR
    mcout << "DynArr(const DynArr<T>& f) is working\n";
#endif
    *this = f;
  }
  DynArr(PILF_CONST DynArr<T>& f, Pilfer)
      : qel(f.qel, steal), cum_qel(f.cum_qel, steal), el(f.el, steal) {
#ifdef DEBUG_DYNARR
    mcout << "DynArr( DynArr<T>& f, Pilfer) is working\n";
#endif
  }
 
  void pilfer(PILF_CONST DynArr<T>& f) {
#ifdef DEBUG_DYNARR
    mcout << "DynArr::pilfer is called\n";
#endif
    if (this != &f) {
      if (qel.get_qel() != 0) {
        if (f.qel.get_qel() != 0) {
          mcerr << "ERROR in DynArr::pilfer(...):\n";
          mcerr << "Both the destination and source arrays are not empty\n";
          // For explanations why it is dangerous, see similar function
          // of ActivePtr.
          spexit(mcerr);
        } else {
          qel.clear();
          cum_qel.clear();
          el.clear();
        }
      }
      qel.pilfer(f.qel);
      cum_qel.pilfer(f.cum_qel);
      el.pilfer(f.el);
      f.qel.clear();
      f.cum_qel.clear();
      f.el.clear();
    }
  }
 
  DynArr<T>& operator=(const DynArr<T>& f);
  template <class D> DynArr<T>& operator=(const DynArr<D>& f);
 
  void pass(long q, DynLinArr<long> fqel, DynLinArr<long> fcum_qel, T* fel)
      // Do not call directly! Is to be used only
      // from assignment operator above
      {
    clear();
    qel = fqel;
    cum_qel = fcum_qel;
    el.pass(q, fel);
  }
 
  // Auxiliary class that provides indexing through ordinary way
  // It is perhaps quite slow.
  template <class D> class IndexingProvider {
   public:
    DynArr<D>& arr;
    mutable long q_deref_ind;
    mutable long current_pos;
    IndexingProvider(DynArr<D>& farr, long fq_deref_ind, long fcurrent_pos)
        : arr(farr), q_deref_ind(fq_deref_ind), current_pos(fcurrent_pos) {}
    operator D&() const {
      if (q_deref_ind != arr.qel.get_qel()) {
        mcerr << "ERROR in IndexingProvider::operator D& (): q_deref_ind != "
                 "qel.get_qel()\n";
        Iprint2n(mcerr, q_deref_ind, arr.qel.get_qel());
        mcerr << "Type of T is (in internal notations) " << typeid(D).name()
              << '\n';
        spexit(mcerr);
      }
#ifdef ALR_CHECK_BOUND
      return arr.el[current_pos];
#else
      return arr.el.acu(current_pos);
#endif
    }
    D& operator=(const D& f) {
      if (q_deref_ind != arr.get_qel().get_qel()) {
        mcerr << "ERROR in T& IndexingProvider::operator=(T& f): q_deref_ind "
                 "!= arr.get_qel().get_qel()\n";
        Iprint2n(mcerr, q_deref_ind, arr.get_qel().get_qel());
        mcerr << "Type of T is (in internal notations) " << typeid(D).name()
              << '\n';
        spexit(mcerr);
      }
#ifdef ALR_CHECK_BOUND
      arr.ac_lin(current_pos) = f;
#else
      arr.acu_lin(current_pos) = f;
#endif
      //return arr.el[current_pos];
    }
 
    inline IndexingProvider<D>& operator[](long n) {
      if (q_deref_ind >= arr.get_qel().get_qel()) {
        mcerr << "ERROR in DynArr::IndexingProvider& DynArr::operator[](long "
                 "n): q_deref_ind >= arr.get_qel().get_qel()\n";
        Iprint2n(mcerr, q_deref_ind, arr.get_qel().get_qel());
        mcerr << "Type of T is (in internal notations) " << typeid(D).name()
              << '\n';
        spexit(mcerr);
      }
#ifdef ALR_CHECK_EACH_BOUND
      if (n >= 0 && n < arr.qel.acu(q_deref_ind)) {
#endif
        current_pos += n * arr.get_cum_qel().acu(q_deref_ind);
        q_deref_ind++;
        return *this;
#ifdef ALR_CHECK_EACH_BOUND
      } else {
        mcerr << "Error in IndexingProvider<D>& "
                 "IndexingProvider::operator[](long n): n < 0 || n >= "
                 "qel.acu(q_deref_ind)\n";
        Iprint2n(mcout, n, arr.qel.acu(q_deref_ind));
        mcerr << "Type of T is (in internal notations) " << typeid(T).name()
              << '\n';
        spexit(mcerr);
      }
#endif
    }
 
    inline const IndexingProvider<D>& operator[](long n) const {
      if (q_deref_ind >= arr.get_qel().get_qel()) {
        mcerr << "ERROR in DynArr::IndexingProvider& DynArr::operator[](long "
                 "n): q_deref_ind >= arr.get_qel().get_qel()\n";
        Iprint2n(mcerr, q_deref_ind, arr.get_qel().get_qel());
        mcerr << "Type of T is (in internal notations) " << typeid(D).name()
              << '\n';
        spexit(mcerr);
      }
#ifdef ALR_CHECK_EACH_BOUND
      if (n >= 0 && n < arr.qel.acu(q_deref_ind)) {
#endif
        current_pos += n * arr.get_cum_qel().acu(q_deref_ind);
        q_deref_ind++;
        return *this;
#ifdef ALR_CHECK_EACH_BOUND
      } else {
        mcerr << "Error in IndexingProvider<D>& "
                 "IndexingProvider::operator[](long n): n < 0 || n >= "
                 "qel.acu(q_deref_ind)\n";
        Iprint2n(mcout, n, arr.qel.acu(q_deref_ind));
        mcerr << "Type of T is (in internal notations) " << typeid(T).name()
              << '\n';
        spexit(mcerr);
      }
#endif
    }
 
  };
 
  // This operator can be used to provide ordinary indexing [][][]... .
  // It is perhaps quite slow compared with ac() functions.
  // It is provided only for the use in contexts in which [][][]...
  // is needed.
  inline IndexingProvider<T> operator[](long n) {
    if (qel.get_qel() < 1) {
      mcerr << "ERROR in DynArr::IndexingProvider DynArr::operator[](long n): "
               "qel.get_qel()< 1, qel.get_qel()=" << qel.get_qel() << '\n';
      mcerr << "Type of T is (in internal notations) " << typeid(T).name()
            << '\n';
      spexit(mcerr);
    }
#ifdef ALR_CHECK_EACH_BOUND
    if (n >= 0 && n < qel.acu(0)) {
#endif
      return IndexingProvider<T>(*this, 1, n * cum_qel.acu(0));
#ifdef ALR_CHECK_EACH_BOUND
    } else {
      mcerr << "Error in IndexingProvider<T> DynArr::operator[](long n): n < 0 "
               "|| n >= qel.acu(0)\n";
      Iprint2n(mcout, n, qel.acu(0));
      mcerr << "Type of T is (in internal notations) " << typeid(T).name()
            << '\n';
      spexit(mcerr);
    }
#endif
  }
 
  inline const IndexingProvider<T> operator[](long n) const {
    if (qel.get_qel() < 1) {
      mcerr << "ERROR in DynArr::IndexingProvider DynArr::operator[](long n): "
               "qel.get_qel()< 1, qel.get_qel()=" << qel.get_qel() << '\n';
      mcerr << "Type of T is (in internal notations) " << typeid(T).name()
            << '\n';
      spexit(mcerr);
    }
#ifdef ALR_CHECK_EACH_BOUND
    if (n >= 0 && n < qel.acu(0)) {
#endif
      //DynArr<T>* temp = static_cast<DynArr<T>* >(this);
      DynArr<T>* temp = const_cast<DynArr<T>*>(this);
      //    static_cast<const IndexingProvider<T> >(*this);
      return IndexingProvider<T>(*temp, 1, n * cum_qel.acu(0));
#ifdef ALR_CHECK_EACH_BOUND
    } else {
      mcerr << "Error in IndexingProvider<T> DynArr::operator[](long n): n < 0 "
               "|| n >= qel.acu(0)\n";
      Iprint2n(mcout, n, qel.acu(0));
      mcerr << "Type of T is (in internal notations) " << typeid(T).name()
            << '\n';
      spexit(mcerr);
    }
#endif
  }
 
  T& ac(long i)  // for 1-dimensional array
      {
    if (qel.get_qel() == 1)
      return el[i];
    else {
      mcerr << "ERROR in DynArr::ac(long i): qel.get_qel()!= 1, qel.get_qel()="
            << qel.get_qel() << '\n';
      mcerr << "Type of T is (in internal notations) " << typeid(T).name()
            << '\n';
      spexit(mcerr);
#ifdef INS_CRETURN   // Insert calming return
      return el[0];  // to quiet Microsoft Visual Studio compiler against
                     // "not all control paths return a value"
#endif
    }
    //return el[i];
  }
  const T& ac(long i) const  // for 1-dimensional array
      {
    if (qel.get_qel() == 1)
      return el[i];
    else {
      mcerr << "ERROR in DynArr::ac(long i): qel.get_qel()!= 1, qel.get_qel()="
            << qel.get_qel() << '\n';
      mcerr << "Type of T is (in internal notations) " << typeid(T).name()
            << '\n';
      spexit(mcerr);
#ifdef INS_CRETURN   // Insert calming return
      return el[0];  // to quiet Microsoft Visual Studio compiler against
                     // "not all control paths return a value"
#endif
    }
    /*
      if(qel.get_qel() != 1)
      { mcerr<<"ERROR in DynArr::ac(long i): qel.get_qel()!= 1, qel.get_qel()="
         <<qel.get_qel()<<'\n';
      mcerr<<"Type of T is (in internal notations) "<<typeid(T).name()<<'\n';
      spexit(mcerr); }
      return el[i];
      */
  }
  inline T& acu(long i1)  // for 1-dimensional array, completely unchecked
      {
    return el.acu(i1);
  }
  inline const T& acu(
      long i1) const  // for 1-dimensional array, completely unchecked
      {
    return el.acu(i1);
  }
 
  T& ac(const DynLinArr<long>& ind)  // for arbitrary number of dimensions
      // but the number of them in array should be equal to size of ind.
      // ind is array of indexes. Its first element if the first index,
      // second is second, etc.
      {
    long q;
    if ((q = qel.get_qel()) != ind.get_qel()) {
      mcerr << "ERROR in DynArr::ac(const DynLinArr<long>& ind): "
            << "qel.get_qel()!= ind.get_qel()\n"
            << "qel.get_qel()=" << qel.get_qel()
            << " ind.get_qel()=" << ind.get_qel() << '\n';
      mcerr << "Type of T is (in internal notations) " << typeid(T).name()
            << '\n';
      spexit(mcerr);
    }
#ifdef ALR_CHECK_EACH_BOUND
#ifdef DEBUG_DYNARR
    if (q == 1)  // faster for this case
      return el[ind.acu(0)];
    else
      return el[calc_lin_ind(ind)];
#else
    if (q == 1)
      return el[ind.acu(0)];
    else
      return el.acu(calc_lin_ind(ind));
#endif
#else
    if (q == 1)
      return el[ind.acu(0)];
    else
      return el[calc_lin_ind(ind)];
#endif
  }
  const T& ac(const DynLinArr<long>& ind) const  // the same as above
      {
    long q;
    if ((q = qel.get_qel()) != ind.get_qel()) {
      mcerr << "ERROR in DynArr::ac(const DynLinArr<long>& ind): "
            << "qel.get_qel()!= ind.get_qel()\n"
            << "qel.get_qel()=" << qel.get_qel()
            << " ind.get_qel()=" << ind.get_qel() << '\n';
      mcerr << "Type of T is (in internal notations) " << typeid(T).name()
            << '\n';
      spexit(mcerr);
    }
/*
#ifdef ALR_CHECK_EACH_BOUND
#ifdef DEBUG_DYNARR
      return el[calc_lin_ind(ind)];
#else
      return el.acu(calc_lin_ind(ind));
#endif
#else
      return el[calc_lin_ind(ind)];
#endif
      //return el[calc_lin_ind(ind)];
      */
 
#ifdef ALR_CHECK_EACH_BOUND
#ifdef DEBUG_DYNARR
    if (q == 1)  // faster for this case
      return el[ind.acu(0)];
    else
      return el[calc_lin_ind(ind)];
#else
    if (q == 1)
      return el[ind.acu(0)];
    else
      return el.acu(calc_lin_ind(ind));
#endif
#else
    if (q == 1)
      return el[ind.acu(0)];
    else
      return el[calc_lin_ind(ind)];
#endif
  }
 
  T& acp(const DynLinArr<long>& ind)  // the same as above, but
      // the size of ind can be more than the number of indexes
      // (the rest is unused)
      {
    long q;
    if ((q = qel.get_qel()) > ind.get_qel()) {
      mcerr << "ERROR in DynArr::acp(const DynLinArr<long>& ind): "
            << "qel.get_qel()!= ind.get_qel()\n"
            << "qel.get_qel()=" << qel.get_qel()
            << " ind.get_qel()=" << ind.get_qel() << '\n';
      mcerr << "Type of T is (in internal notations) " << typeid(T).name()
            << '\n';
      spexit(mcerr);
    }
 
#ifdef ALR_CHECK_EACH_BOUND
#ifdef DEBUG_DYNARR
    if (q == 1)  // faster for this case
      return el[ind.acu(0)];
    else
      return el[calc_lin_ind(ind)];
#else
    if (q == 1)
      return el[ind.acu(0)];
    else
      return el.acu(calc_lin_ind(ind));
#endif
#else
    if (q == 1)
      return el[ind.acu(0)];
    else
      return el[calc_lin_ind(ind)];
#endif
  }
 
  const T& acp(const DynLinArr<long>& ind) const {
    long q;
    if ((q = qel.get_qel()) > ind.get_qel()) {
      mcerr << "ERROR in DynArr::acp(const DynLinArr<long>& ind): "
            << "qel.get_qel()!= ind.get_qel()\n"
            << "qel.get_qel()=" << qel.get_qel()
            << " ind.get_qel()=" << ind.get_qel() << '\n';
      mcerr << "Type of T is (in internal notations) " << typeid(T).name()
            << '\n';
      spexit(mcerr);
    }
#ifdef ALR_CHECK_EACH_BOUND
#ifdef DEBUG_DYNARR
    if (q == 1)  // faster for this case
      return el[ind.acu(0)];
    else
      return el[calc_lin_ind(ind)];
#else
    if (q == 1)
      return el[ind.acu(0)];
    else
      return el.acu(calc_lin_ind(ind));
#endif
#else
    if (q == 1)
      return el[ind.acu(0)];
    else
      return el[calc_lin_ind(ind)];
#endif
  }
 
  T& acu(const DynLinArr<long>& ind)  // unchecked
      {
    if (qel.get_qel() == 1)
      return el.acu(ind.acu(0));
    else
      return el.acu(calc_lin_ind(ind));
  }
  const T& acu(const DynLinArr<long>& ind) const  // unchecked
      {
    if (qel.get_qel() == 1)
      return el.acu(ind.acu(0));
    else
      return el.acu(calc_lin_ind(ind));
  }
 
  T& ac(long i1, long i2)  // for 2-dimensional array
      {
    if (qel.get_qel() == 2) {
#ifdef ALR_CHECK_EACH_BOUND
      if (i1 >= 0 && i1 < qel.acu(0)) {
        if (i2 >= 0 && i2 < qel.acu(1)) {
#ifdef DEBUG_DYNARR
          return el[i1 * cum_qel.acu(0) + i2];
#else
          return el.acu(i1 * cum_qel.acu(0) + i2);
#endif
        } else {
          mcerr << "Error in DynArr::ac(long i1, long i2): i2 < 0 || i2 >= "
                   "qel.acu(1)\n";
          Iprint2n(mcout, i2, qel[1]);
        }
      } else {
        mcerr << "Error in DynArr::ac(long i1, long i2): i1 < 0 || i1 >= "
                 "qel.acu(0)\n";
        Iprint2n(mcout, i1, qel[0]);
      }
    } else {
      mcerr << "ERROR in DynArr::ac(long i1, long i2): qel.get_qel()!= 2,"
            << " qel.get_qel()=" << qel.get_qel() << '\n';
    }
    mcerr << "Type of T is (in internal notations) " << typeid(T).name()
          << '\n';
    spexit(mcerr);
#ifdef INS_CRETURN  // Insert calming return
    return el[0];   // to quiet Microsoft Visual Studio compiler against
                    // "not all control paths return a value"
#endif
#else  // for ifdef ALR_CHECK_EACH_BOUND
    return el[i1 * cum_qel.acu(0) + i2];
 
#endif
  }
 
  /*
  {
    if(qel.get_qel() != 2)
    { mcerr<<"ERROR in DynArr::ac(long i1, long i2): qel.get_qel()!= 2,"
       <<" qel.get_qel()=" <<qel.get_qel()<<'\n';
    mcerr<<"Type of T is (in internal notations) "<<typeid(T).name()<<'\n';
    spexit(mcerr); }
#ifdef ALR_CHECK_EACH_BOUND
    if(i1 < 0 || i1 >= qel.acu(0))
    {
      mcerr<<"Error in DynArr::ac(long i1, long i2): i1 < 0 || i1 >=
qel.acu(0)\n";
      Iprint2n(mcout, i1, qel[0]);
      spexit(mcerr);
    }
    if(i2 < 0 || i2 >= qel.acu(1))
    {
      mcerr<<"Error in DynArr::ac(long i1, long i2): i2 < 0 || i2 >=
qel.acu(1)\n";
      Iprint2n(mcout, i2, qel[1]);
      spexit(mcerr);
    }
#ifdef DEBUG_DYNARR
    return el[i1*cum_qel.acu(0) + i2];
#else
    return el.acu(i1*cum_qel.acu(0) + i2);
#endif
#else
    return el[i1*cum_qel.acu(0) + i2];
#endif
  }
  */
  const T& ac(long i1, long i2) const  // for 2-dimensional array
      {
    if (qel.get_qel() == 2) {
#ifdef ALR_CHECK_EACH_BOUND
      if (i1 >= 0 && i1 < qel.acu(0)) {
        if (i2 >= 0 && i2 < qel.acu(1)) {
#ifdef DEBUG_DYNARR
          return el[i1 * cum_qel.acu(0) + i2];
#else
          return el.acu(i1 * cum_qel.acu(0) + i2);
#endif
        } else {
          mcerr << "Error in DynArr::ac(long i1, long i2): i2 < 0 || i2 >= "
                   "qel.acu(1)\n";
          Iprint2n(mcout, i2, qel[1]);
        }
      } else {
        mcerr << "Error in DynArr::ac(long i1, long i2): i1 < 0 || i1 >= "
                 "qel.acu(0)\n";
        Iprint2n(mcout, i1, qel[0]);
      }
    } else {
      mcerr << "ERROR in DynArr::ac(long i1, long i2): qel.get_qel()!= 2,"
            << " qel.get_qel()=" << qel.get_qel() << '\n';
    }
    mcerr << "Type of T is (in internal notations) " << typeid(T).name()
          << '\n';
    spexit(mcerr);
#ifdef INS_CRETURN  // Insert calming return
    return el[0];   // to quiet Microsoft Visual Studio compiler against
                    // "not all control paths return a value"
#endif
#else  // for ifdef ALR_CHECK_EACH_BOUND
    return el[i1 * cum_qel.acu(0) + i2];
 
#endif
  }
  /*
  {
    if(qel.get_qel() != 2)
    { mcerr<<"ERROR in DynArr::ac(long i1, long i2): qel.get_qel()!= 2,"
       <<" qel.get_qel()=" <<qel.get_qel()<<'\n';
    mcerr<<"Type of T is (in internal notations) "<<typeid(T).name()<<'\n';
    spexit(mcerr); }
#ifdef ALR_CHECK_EACH_BOUND
    if(i1 < 0 || i1 >= qel.acu(0))
    {
      mcerr<<"Error in DynArr::ac(long i1, long i2): i1 < 0 || i1 >=
qel.acu(0)\n";
      Iprint2n(mcout, i1, qel[0]);
      spexit(mcerr);
    }
    if(i2 < 0 || i2 >= qel.acu(1))
    {
      mcerr<<"Error in DynArr::ac(long i1, long i2): i2 < 0 || i2 >=
qel.acu(1)\n";
      Iprint2n(mcout, i2, qel[1]);
      spexit(mcerr);
    }
#ifdef DEBUG_DYNARR
    return el.acu(i1*cum_qel.acu(0) + i2);
#else
    return el[i1*cum_qel.acu(0) + i2];
#endif
#else
    return el[i1*cum_qel.acu(0) + i2];
#endif
    //return el[i1*cum_qel[0] + i2];
  }
    */
  inline T& acu(long i1,
                long i2)  // for 2-dimensional array, completely unchecked
      {
    return el.acu(i1 * cum_qel.acu(0) + i2);
  }
  inline const T& acu(long i1, long i2) const  // for 2-dimensional array
      {
    return el.acu(i1 * cum_qel.acu(0) + i2);
  }
 
  T& ac(long i1, long i2, long i3)  // for 3-dimensional array
      {
    if (qel.get_qel() != 3) {
      mcerr << "ERROR in DynArr::ac(long i1, long i2, long i3): "
               "qel.get_qel()!= 3,"
            << " qel.get_qel()=" << qel.get_qel() << '\n';
      mcerr << "Type of T is (in internal notations) " << typeid(T).name()
            << '\n';
      spexit(mcerr);
    }
#ifdef ALR_CHECK_EACH_BOUND
    if (i1 < 0 || i1 >= qel.acu(0)) {
      mcerr << "Error in DynArr::ac(long i1, long i2, long i3): i1 < 0 || i1 "
               ">= qel.acu(0)\n";
      Iprint2n(mcout, i1, qel[0]);
      spexit(mcerr);
    }
    if (i2 < 0 || i2 >= qel.acu(1)) {
      mcerr << "Error in DynArr::ac(long i1, long i2, long i3): i2 < 0 || i2 "
               ">= qel.acu(1)\n";
      Iprint2n(mcout, i2, qel[1]);
      spexit(mcerr);
    }
    if (i3 < 0 || i3 >= qel.acu(2)) {
      mcerr << "Error in DynArr::ac(long i1, long i2, long i3): i3 < 0 || i3 "
               ">= qel.acu(2)\n";
      Iprint2n(mcout, i3, qel[2]);
      spexit(mcerr);
    }
#ifdef DEBUG_DYNARR
    return el.acu(i1 * cum_qel.acu(0) + i2 * cum_qel.acu(1) + i3);
#else
    return el[i1 * cum_qel.acu(0) + i2 * cum_qel[1] + i3];
#endif
#else
    return el[i1 * cum_qel.acu(0) + i2 * cum_qel[1] + i3];
#endif
    //return el[i1*cum_qel[0] + i2*cum_qel[1] + i3];
  }
 
  const T& ac(long i1, long i2, long i3) const  // for 3-dimensional array
      {
    if (qel.get_qel() != 3) {
      mcerr << "ERROR in DynArr::ac(long i1, long i2, long i3): "
               "qel.get_qel()!= 3,"
            << " qel.get_qel()=" << qel.get_qel() << '\n';
      mcerr << "Type of T is (in internal notations) " << typeid(T).name()
            << '\n';
      spexit(mcerr);
    }
#ifdef ALR_CHECK_EACH_BOUND
    if (i1 < 0 || i1 >= qel.acu(0)) {
      mcerr << "Error in DynArr::ac(long i1, long i2, long i3): i1 < 0 || i1 "
               ">= qel.acu(0)\n";
      Iprint2n(mcout, i1, qel[0]);
      spexit(mcerr);
    }
    if (i2 < 0 || i2 >= qel.acu(1)) {
      mcerr << "Error in DynArr::ac(long i1, long i2, long i3): i2 < 0 || i2 "
               ">= qel.acu(1)\n";
      Iprint2n(mcout, i2, qel[1]);
      spexit(mcerr);
    }
    if (i3 < 0 || i3 >= qel.acu(2)) {
      mcerr << "Error in DynArr::ac(long i1, long i2, long i3): i3 < 0 || i3 "
               ">= qel.acu(2)\n";
      Iprint2n(mcout, i3, qel[2]);
      spexit(mcerr);
    }
#ifdef DEBUG_DYNARR
    return el.acu(i1 * cum_qel.acu(0) + i2 * cum_qel.acu(1) + i3);
#else
    return el[i1 * cum_qel.acu(0) + i2 * cum_qel[1] + i3];
#endif
#else
    return el[i1 * cum_qel.acu(0) + i2 * cum_qel[1] + i3];
#endif
    //return el[i1*cum_qel[0] + i2*cum_qel[1] + i3];
  }
 
  long get_qel_lin(void) const { return el.get_qel(); }
 
  // access to array as linear array
  inline T& ac_lin(long n) {
    long qelln = el.get_qel();
#ifdef ALR_CHECK_BOUND
    if (n >= 0 && n < qelln) {
      return el[n];
    } else {
      mcerr << "ERROR in T& DynArr::ac_lin(long n): "
            << "n is out of bounds, n=" << n << " qelln=" << qelln << '\n';
      mcerr << "Type of T is (in internal notations) " << typeid(T).name()
            << '\n';
      spexit(mcerr);
#ifdef INS_CRETURN
      return el[0];  // to quiet Microsoft Visial Studio compiler against
                     // "not all control paths return a value"
#endif
    }
#else
    return el[n];
#endif
  }
  inline const T& ac_lin(long n) const {
    long qelln = el.get_qel();
#ifdef ALR_CHECK_BOUND
    if (n >= 0 && n < qelln) {
      return el[n];
    } else {
      mcerr << "ERROR in T& DynArr::ac_lin(long n): "
            << "n is out of bounds, n=" << n << " qelln=" << qelln << '\n';
      mcerr << "Type of T is (in internal notations) " << typeid(T).name()
            << '\n';
      spexit(mcerr);
#ifdef INS_CRETURN
      return el[0];  // to quiet Microsoft Visial Studio compiler against
                     // "not all control paths return a value"
#endif
    }
#else
    return el[n];
#endif
  }
  // access to array as linear array always without check
  inline T& acu_lin(long n) { return el[n]; }
  inline const T& acu_lin(long n) const { return el[n]; }
 
  void assignAll(const T& val);
 
  long get_qdim(void) const { return qel.get_qel(); }
  const DynLinArr<long>& get_qel(void) const { return qel; }
  const DynLinArr<T>& get_el(void) const { return el; }
 
  // The following is mainly for debug (diagnostic print):
  const DynLinArr<long>& get_cum_qel(void) const { return cum_qel; }
 
  //void put_qel(const DynLinArr<long>& fqel, T* val=NULL);
  void put_qel(T* val = NULL);
  // 25.10.2006: Today I do not understand these following  comments.
  // They looks like simple copy from these in DynLinArr.
  // creates array with size fqel
  // If old array existed, then
  //   If it was less than fqel, it all is copied to new array
  //       and the other elements are either remains not inited
  //       or assignned *val.
  //   else its fqel part is copyed to new array.
 
  void clear(void) {
    qel.clear();
    cum_qel.clear();
    el.clear();
  }
 
  int confirm_ind(const DynLinArr<long>& ind) { return gconfirm_ind(qel, ind); }
  int confirm_ind_ext(const DynLinArr<long>& ind) {
    return gconfirm_ind_ext(qel, ind);
  }
  DynArr<T> top(void);
 
  // Apply any function of one argument (of type T) to each element.
  template <class P> friend void apply1(DynArr<P>& ar, void (*fun)(P& f));
 
  // Apply any function of two arguments
  // (first of which is of type T and the second is of type of address
  // of another function (possibly apply1) to each element.
  template <class P, class X>
  friend void apply2(DynArr<P>& ar, void (*fun1)(P& f, void (*fun21)(X& f)),
                     void (*fun2)(X& f));
 
  void check(void) const {
    qel.check();
    cum_qel.check();
    el.check();
  }
  int get_s_non_emplty(void) const {
    long q = qel.get_qel();
    if (q == 0) return 0;
    long n;
    for (n = 0; n < q; n++) {
      if (qel[n] <= 0) return 0;
    }
    return 1;
  }
#ifndef DONT_USE_ABSPTR
  macro_copy_total(DynArr);
#endif
  virtual ~DynArr(void) {}
 
 private:
  PILF_MUTABLE DynLinArr<long> qel;
  // Linear array with number of elements by each dimension
  PILF_MUTABLE DynLinArr<long> cum_qel;
  // "cumulative qel": each element contains product
  // of next elements of qel.
  // The last element is always 1.
  // Used for fast search of proper element in el.
 
  PILF_MUTABLE DynLinArr<T> el;
  // Contains all elements.
  // The last index varies faster.
 
  long calc_lin_ind(const DynLinArr<long>& ind) const {
    long i = 0;
    long n;
    long qdim1 = qel.get_qel() - 1;
    /*
      // This check is not necessary by two reasons
      // 1. The checks of this condition are made in calling functions.
      // 2. The correct condition is not != but >.
      if(qdim1 != ind.get_qel() - 1)
      {
    mcerr<<"ERROR in long DynArr::calc_lin_ind(const DynLinArr<long>& ind)
const\n";
    mcerr<<"qdim1 != ind.get_qel() - 1\n";
    Iprint2n(mcerr, qdim1, (ind.get_qel() - 1));
    spexit(mcerr);
      }
      */
    for (n = 0; n < qdim1; n++) {
#ifdef ALR_CHECK_EACH_BOUND
#ifdef DEBUG_DYNARR
      if (ind[n] < 0 || ind[n] >= qel[n]) {
        mcerr << "ERROR in long DynArr::calc_lin_ind(const DynLinArr<long>& "
                 "ind) const\n";
        mcerr << "ind[n] < 0 || ind[n] >= qel[n]\n";
        Iprint3n(mcout, n, ind[n], qel[n]);
        Iprintn(mcout, qel);
        spexit(mcerr);
      }
#else
      if (ind.acu(n) < 0 || ind.acu(n) >= qel.acu(n)) {
        mcerr << "ERROR in long DynArr::calc_lin_ind(const DynLinArr<long>& "
                 "ind) const\n";
        mcerr << "ind.acu(n) < 0 || ind.acu(n) >= qel.acu(n)\n";
        Iprint3n(mcout, n, ind[n], qel[n]);
        Iprintn(mcout, qel);
        spexit(mcerr);
      }
#endif
#endif
      i += ind[n] * cum_qel[n];
    }
#ifdef ALR_CHECK_EACH_BOUND
#ifdef DEBUG_DYNARR
    if (ind[qdim1] < 0 || ind[qdim1] >= qel[qdim1]) {
      mcerr << "ERROR in long DynArr::calc_lin_ind(const DynLinArr<long>& ind) "
               "const\n";
      mcerr << "ind[qdim1] < 0 || ind[qdim1] >= qel[qdim1]\n";
      Iprint3n(mcout, n, ind[qdim1], qel[qdim1]);
      Iprintn(mcout, qel);
      spexit(mcerr);
    }
#else
    if (ind.acu(qdim1) < 0 || ind.acu(qdim1) >= qel.acu(qdim1)) {
      mcerr << "ERROR in long DynArr::calc_lin_ind(const DynLinArr<long>& ind) "
               "const\n";
      mcerr << "ind.acu(qdim1) < 0 || ind.acu(qdim1) >= qel.acu(qdim1)\n";
      Iprint3n(mcout, n, ind[qdim1], qel[qdim1]);
      Iprintn(mcout, qel);
      spexit(mcerr);
    }
#endif
#endif
    i += ind[qdim1];  // the last index, for which multiplication by cum_qel
                      // is not necessary.
    return i;
  }
#ifdef USE_REPLACE_ALLOC
  macro_alloc
#endif
 
};
 
template <class T> DynArr<T>& DynArr<T>::operator=(const DynArr<T>& f) {
#ifdef DEBUG_DYNARR
  mcout << "DynArr<T>& DynArr<T>::operator=(const DynArr<T>& f)\n";
#endif
  if (this != &f) {
    //mcout<<"DynLinArr<T>& operator=(const DynLinArr<T>& f): long(el)="
    //<<long(el)<<'\n';
    check();
    f.check();
    qel = f.qel;
    cum_qel = f.cum_qel;
    el = f.el;
  }
  return *this;
}
 
template <class T>
template <class D>
DynArr<T>& DynArr<T>::operator=(const DynArr<D>& f) {
#ifdef DEBUG_DYNLINARR
  mcout << "DynArr<T>& DynArr<T>::operator=(const DynArr<D>& f)\n";
#endif
  //if(this != &f)
  //{
  //mcout<<"DynLinArr<T>& operator=(const DynLinArr<T>& f): long(el)="
  //<<long(el)<<'\n';
  check();
  f.check();
  DynLinArr<long> qel = f.get_qel();
  DynLinArr<long> cum_qel = f.get_cum_qel();
  // for example, one of its elements.
  long q = f.get_qel_lin();
#ifdef USE_REPLACE_ALLOC
  T* temp_el;
  if (q > 0) {
    temp_el = (T*)malloc(sizeof(T) * q);
    long n;
    for (n = 0; n < q; n++)
      new (&(temp_el[n])) T;
  } else
    temp_el = NULL;
#else
  T* temp_el = (q > 0) ? (new T[q]) : (T*)NULL;
#endif
  long n;
  for (n = 0; n < q; n++)
    temp_el[n] = f.acu_lin(n);
  pass(q, qel, cum_qel, temp_el);
  //}
  return *this;
}
 
template <class T> void apply1(DynArr<T>& ar, void (*fun)(T& f)) {
  long q = ar.el.get_qel();
  long n;
  for (n = 0; n < q; n++)
    (*fun)(ar.el[n]);
}
template <class T, class X>
void apply2(DynArr<T>& ar, void (*fun1)(T& f, void (*fun21)(X& f)),
            void (*fun2)(X& f)) {
  long q = ar.el.get_qel();
  long n;
  for (n = 0; n < q; n++)
    (*fun1)(ar.el[n], fun2);
}
 
int find_next_comb(const DynLinArr<long>& qel, DynLinArr<long>& f);
int find_prev_comb(const DynLinArr<long>& qel, DynLinArr<long>& f);
int find_next_comb_not_less(const DynLinArr<long>& qel, DynLinArr<long>& f);
 
template <class T> class IterDynArr {
 public:
  IterDynArr(const DynArr<T>* fdar)
      : ncur(fdar->get_qdim()), dar((DynArr<T>*)fdar) {
    long n;
    long qdim1 = ncur.get_qel() - 1;
    if (qdim1 >= 0) {
      for (n = 0; n < qdim1; n++)
        ncur[n] = 0;
      ncur[qdim1] = -1;
    }
  }
  T* more(void)  // Just next element. Why not "next", do not remember.
      {
    if (find_next_comb(dar->get_qel(), ncur))
      return &(dar->ac(ncur));
    else
      return NULL;
  }
 
  T* current(void)  // Element currently pointed by ncut.
      {
    long n;
    long qdim = ncur.get_qel();
    if (qdim <= 0) return NULL;
    for (n = 0; n < qdim; n++) {
      if (ncur[n] < 0 || ncur[n] >= dar->get_qel()[n]) return NULL;
    }
    return &(dar->ac(ncur));
  }
 
  T* less(void)  // Switch current to previous.
      // Again do not remember why not" prev"
      {
    if (find_prev_comb(dar->get_qel(), ncur))
      return &(dar->ac(ncur));
    else
      return NULL;
  }
 
  const DynLinArr<long>& get_ncur(void) const { return ncur; }
 
 private:
  DynLinArr<long> ncur;
  DynArr<T>* dar;
};
 
extern DynLinArr<long> qel_communicat;
 
template <class T>
void DynArr<T>::put_qel(T* val)
    //  by default           val=NULL
    {
  check();
  if (qel.get_qel() == 0) {
    *this = DynArr<T>(qel_communicat, val);
    return;
  }
 
  DynArr<T> datemp(qel_communicat, val);  // all init to val
  IterDynArr<T> iter(&datemp);
  T* at;
  while ((at = iter.more()) != NULL) {
    if (confirm_ind_ext(iter.get_ncur()))
      *at = acp(iter.get_ncur());  // change to old values where they were
  }
  *this = datemp;
}
 
template <class T>
    int operator==(const DynLinArr<T>& f1, const DynLinArr<T>& f2) {
  if (f1.get_qel() != f2.get_qel()) return 0;
  long q = f1.get_qel();
  long n;
  for (n = 0; n < q; n++) {
    if (!(f1.acu(n) == f2.acu(n))) return 0;
  }
  return 1;
}
 
template <class T, class P>
int apeq_mant(const DynLinArr<T>& f1, const DynLinArr<T>& f2, P prec) {
  if (f1.get_qel() != f2.get_qel()) return 0;
  long q = f1.get_qel();
  long n;
  for (n = 0; n < q; n++) {
    if (!apeq_mant(f1.acu(n), f2.acu(n), prec)) return 0;
  }
  return 1;
}
 
template <class T>
    int operator!=(const DynLinArr<T>& f1, const DynLinArr<T>& f2) {
  if (f1 == f2)
    return 0;
  else
    return 1;
}
 
template <class T> int operator==(const DynArr<T>& f1, const DynArr<T>& f2) {
  if (f1.get_qel() != f2.get_qel()) return 0;
  if (f1.get_el() != f2.get_el()) return 0;
  return 1;
}
 
template <class T, class P>
int apeq_mant(const DynArr<T>& f1, const DynArr<T>& f2, P prec) {
  if (f1.get_qel() != f2.get_qel()) return 0;
  if (!apeq_mant(f1.get_el(), f2.get_el(), prec)) return 0;
  return 1;
}
 
template <class T> int operator!=(const DynArr<T>& f1, const DynArr<T>& f2) {
  if (f1.get_qel() == f2.get_qel()) return 0;
  if (f1.get_el() == f2.get_el()) return 0;
  return 1;
}
 
template <class T> void DynArr<T>::assignAll(const T& val) {
  check();
  // try faster and simpler way (30.10.2006):
  el.assignAll(val);
  /*
  IterDynArr<T> iter(this);
  T* at;
  while( (at=iter.more()) != NULL )
  {
    *at=val;
  }
  */
}
 
template <class T, class X> void copy_DynArr(const DynArr<T>& s, DynArr<X>& d) {
  mfunnamep("template<class T, class X> void copy_DynArr(const DynArr<T>& s, "
            "DynArr<X>& d)");
  s.check();
  d.check();
  d = DynArr<X>(s.get_qel(), NULL);
  IterDynArr<T> iter(&s);
  T* at;
  while ((at = iter.more()) != NULL) {
    const DynLinArr<long>& ncur = iter.get_ncur();
    d.ac(ncur) = *at;
  }
}
 
// Convert types of elements from T to X.
template <class T, class X>
void convert_DynArr(const DynArr<T>& s, DynArr<X>& d) {
  mfunnamep("template<class T, class X> void convert_DynArr(const DynArr<T>& "
            "s, DynArr<X>& d)");
  s.check();
  d.check();
  d = DynArr<X>(s.get_qel(), NULL);
  IterDynArr<T> iter(&s);
  T* at;
  while ((at = iter.more()) != NULL) {
    const DynLinArr<long>& ncur = iter.get_ncur();
    d.ac(ncur) = X(*at);
  }
}
 
template <class T> DynArr<T> DynLinArr<T>::top(void) {
  check();
  DynArr<T> r(1, qel);
  long n;
  for (n = 0; n < qel; n++) {
    r.ac(0, n) = el[n];
  }
  return r;
}
template <class T> DynArr<T> DynArr<T>::top(void) {
  mfunnamep("template<class T> DynArr<T> DynArr<T>::top(void)");
  check();
  long qdim = get_qdim();
  check_econd11(qdim, != 2, mcerr);
  long n1, n2;
  DynArr<T> r(qel[1], qel[0]);
  for (n1 = 0; n1 < qel[0]; n1++) {
    for (n2 = 0; n2 < qel[1]; n2++) {
      r.ac(n2, n1) = ac(n1, n2);
    }
  }
  return r;
}
 
template <class T> DynLinArr<T>::DynLinArr(const DynArr<T>& f) {
  //mcout<<"template<class T> DynLinArr<T>::DynLinArr(const DynArr<T>& f):\n";
  f.check();
  long qdim = f.get_qdim();
  if (qdim != 1) {
    mcerr << "ERROR in DynLinArr<T>::DynLinArr(const DynArr<T>& f):\n"
          << "f.get_qdim() != 1, f.get_qdim()=" << f.get_qdim() << '\n';
    mcerr << "Type of T is (in internal notations) " << typeid(T).name()
          << '\n';
    spexit(mcerr);
  }
  const DynLinArr<long>& qelem = f.get_qel();
  qel = qelem[0];
  //mcout<<"qel="<<qel<<'\n';
  if (qel > 0) {
#ifdef USE_REPLACE_ALLOC
    el = (T*)malloc(sizeof(T) * qel);
    long n;
    for (n = 0; n < qel; n++)
      new (&(el[n])) T;
#else
    el = new T[qel];
#endif
    long n;
    for (n = 0; n < qel; n++)
      el[n] = f.ac(n);
  } else
    el = NULL;
}
 
template <class T>
DynLinArr<T>::DynLinArr(const DynArr<T>& f, int n_of_dim,
                        // 0 - first dim) 1 - second dim)
                        long roc_number)
    // takes only mentioned raw or column.
    {
  f.check();
  long qdim = f.get_qdim();
  if (qdim != 2) {
    mcerr << "ERROR in DynLinArr<T>::DynLinArr(const DynArr<T>& f, int "
             "n_of_dim,long roc_number):\n"
          << "f.get_qdim() != 2, f.get_qdim()=" << f.get_qdim() << '\n';
    mcerr << "Type of T is (in internal notations) " << typeid(T).name()
          << '\n';
    spexit(mcerr);
  }
  const DynLinArr<long>& qelem = f.get_qel();
  if (n_of_dim == 0) {
    qel = qelem[1];
  } else {
    qel = qelem[0];
  }
  //mcout<<"qel="<<qel<<'\n';
  if (qel > 0) {
#ifdef USE_REPLACE_ALLOC
    el = (T*)malloc(sizeof(T) * qel);
    long n;
    for (n = 0; n < qel; n++)
      new (&(el[n])) T;
#else
    el = new T[qel];
#endif
    long n;
    if (n_of_dim == 0) {
      for (n = 0; n < qel; n++)
        el[n] = f.ac(roc_number, n);
    } else {
      for (n = 0; n < qel; n++)
        el[n] = f.ac(n, roc_number);
    }
  }
}
 
template <class T>
    std::ostream& operator<<(std::ostream& file, const DynArr<T>& f) {
  //mfunnamep("template<class T> std::ostream& operator<<(std::ostream& file,
  //const DynArr<T>& f)");
  f.check();
  Ifile << "DynArr<T>: qdim=" << f.get_qdim() << '\n';
  indn.n += 2;
  if (s_short_output > 0) {
    Ifile << noindent << f.get_qel() << yesindent;
  } else {
    Ifile << "qel=" << noindent << f.get_qel() << yesindent;
    Ifile << "cum_qel=" << noindent << f.get_cum_qel() << yesindent;
  }
  if (f.get_s_non_emplty() == 1) {
    if (s_short_output == 0) {
      Ifile << "Content element by element:\n";
      Ifile << "(The first number is sequencial number, then there are "
               "indexes, the last is the element)\n";
      //DynArr<T>& ff(f);
    }
    long nseq = 0;
    IterDynArr<T> iter_f(&((DynArr<T>&)f));
    T* at;
    while ((at = iter_f.more()) != NULL) {
      std::ostringstream ost;
      if (s_short_output == 0) {
        //Ifile<<"ncur="<<noindent<<iter_f.get_ncur()<<yesindent;
        Ifile << "nseq=" << std::setw(5) << nseq << " ncur=";
        long n;
        for (n = 0; n < iter_f.get_ncur().get_qel(); n++) {
          file << ' ' << std::setw(5) << iter_f.get_ncur()[n];
        }
        ost << indn << " element=" << noindent << (*at) << yesindent;
      } else {
        ost << indn << noindent << (*at) << yesindent;
      }
      put_one_n(ost);
      file << ost.str();
      nseq++;
    }
    file << yesindent;
  } else {
    if (s_short_output == 0) {
      Ifile << "Content is empty.\n";
    }
  }
  indn.n -= 2;
  return file;
}
/*
template<class T>
void DybArr<T>::short_output(std::ostream& file)
{
  mfunnamep("template<class T> void DybArr<T>::short_output(std::ostream&
file))");
  f.check();
  Ifile<<"DynArr<T>: qdim="<<f.get_qdim()<<'\n';
  indn.n+=2;
  qel.short_output(file);
  if(f.get_s_non_emplty() == 1)
  {
    long nseq=0;
    IterDynArr<T> iter_f( &((DynArr<T>&) f));
    T* at;
    while( (at=iter_f.more()) != NULL )
    {
      std::ostringstream ost;
      ost<<indn<<noindent<<(*at)<<yesindent;
      put_one_n(ost);
      file<<ost.str();
      nseq++;
    }
    file<<yesindent;
  }
  indn.n-=2;
  return file;
}
*/
template <class T> std::istream& operator>>(std::istream& file, DynArr<T>& f) {
  mfunnamep(
      "template<class T> istream& operator>>(istream& file, DynArr<T>& f)");
  definp_endpar dep(&file, 0, 1, 0, s_short_output);
  long qdim = 0;
  dep.s_short = 0;
  DEFINPAP(qdim);
  dep.s_short = s_short_output;
  check_econd11(qdim, < 0, mcerr);
  if (s_short_output == 0) {
    set_position("qel=DynLinArr<T>:", *dep.istrm, dep.s_rewind, dep.s_req_sep);
  } else {
    set_position("DynLinArr<T>:", *dep.istrm, dep.s_rewind, dep.s_req_sep);
  }
  DynLinArr<long> qel_loc;
  //mcout<<"now will read qel_loc\n";
  file >> qel_loc;
  //mcout<<"qel_loc is read\n";
  if (s_short_output == 0) {
    set_position("cum_qel=DynLinArr<T>:", *dep.istrm, dep.s_rewind,
                 dep.s_req_sep);
    DynLinArr<long> cum_qel_loc;
    file >> cum_qel_loc;  // this is in general unnecessary
  }
  if (qel_loc.get_qel() > 0) {
    f.pilfer(DynArr<T>(qel_loc, NULL));
    long nseq;
    long n;
    long qseq = qel_loc[0];
    for (n = 1; n < qel_loc.get_qel(); n++) {
      qseq *= qel_loc[n];
    }
    for (n = 0; n < qseq; n++) {
      if (s_short_output == 0) {
        DEFINPAP(nseq);
        check_econd12(nseq, !=, n, mcerr);
        DynLinArr<long> ncur(qel_loc.get_qel());
        set_position("ncur=", *dep.istrm, dep.s_rewind, dep.s_req_sep);
        long m;
        for (m = 0; m < qel_loc.get_qel(); m++) {
          file >> ncur[m];
        }
        //T element;
        //DEFINPAP(element);
        //f.ac(ncur) = element;
        set_position("element=", *dep.istrm, dep.s_rewind, dep.s_req_sep);
      }
      file >> f.ac_lin(n);
    }
  } else {
    if (s_short_output == 0) {
      // just pass to end
      set_position("Content is empty.", *dep.istrm, dep.s_rewind,
                   dep.s_req_sep);
    }
    f.pilfer(DynArr<T>());
 
  }
  return file;
}
 
/*
template<class T>
void DynArr<T>::short_read(istream& file)
{
  mfunnamep("template<class T> void DynArr<T>::short_read(istream& file)");
  definp_endpar dep(&file, 0, 1, 0);
  long qdim=0;
  DEFINPAP(qdim);
  check_econd11(qdim, < 0 , mcerr);
  DynLinArr<long> qel_loc;
  qel_loc.short_read(file);
  // generate cum
  if(qel_loc.get_qel() > 0 )
  {
    f.pilfer(DynArr<T>(qel_loc, NULL));
    long nseq;
    long n;
    long qseq=qel[0];
    for(n=1; n<el.get_qel(); n++)
    {
      qseq*=el[n];
    }
    for(n=0; n<qseq; n++)
    {
      file>>f.ac(ncur);
    }
  }
  else
  {
    // just pass to end
    set_position("Content is empty.",
         *dep.istrm, dep.s_rewind, dep.s_req_sep);
  }
  return file;
}
*/
 
template <class T>
void print_DynArr(std::ostream& file, const DynArr<T>& f, int l) {
  //mfunnamep("template<class T> oid print_DynArr(std::ostream& file, const
  //DynArr<T>& f, int l)");
  f.check();
  //Ifile<<"DynArr<T>: qdim="<<f.get_qdim()
  //     <<" qel="<<noindent<<f.get_qel()<<yesindent<<'\n';
  Ifile << "DynArr<T>: qdim=" << f.get_qdim() << '\n';
  indn.n += 2;
  Ifile << "qel=" << noindent << f.get_qel() << yesindent;
  Ifile << "cum_qel=" << noindent << f.get_cum_qel() << yesindent;
  Ifile << "Content element by element:\n";
  Ifile << "(The first number is sequencial number, then there are indexes, "
           "the last is the element)\n";
  //DynArr<T>& ff(f);
  long nseq = 0;
  IterDynArr<T> iter_f(&((DynArr<T>&)f));
  T* at;
  while ((at = iter_f.more()) != NULL) {
    //Ifile<<"ncur="<<noindent<<iter_f.get_ncur()<<yesindent;
    Ifile << "nseq=" << std::setw(5) << nseq << " ncur=";
    long n;
    for (n = 0; n < iter_f.get_ncur().get_qel(); n++) {
      file << ' ' << std::setw(5) << iter_f.get_ncur()[n];
    }
    //file<<'\n';
    //Ifile<<"element="<<noindent; at->print(file, l);
    //file<<yesindent<<'\n';
    std::ostringstream ost;
    ost << indn << " element=" << noindent;
    at->print(ost, l);
    ost << yesindent;
    put_one_n(ost);
    file << ost.str();
  }
  file << yesindent;
  indn.n -= 2;
 
}
 
// New experimental approach.
// give the width of field, which is put in setw().
// Whether the array will be printed in single lines
// or by columns, is determined by whether 80 symbols are enough
void print_DynArr_int_w(std::ostream& file, const DynArr<int>& f, int w);
 
void print_DynArr_float(std::ostream& file, const DynArr<float>& f);
void print_DynArr_double(std::ostream& file, const DynArr<double>& f);
// ^Identical functions
// See AbsArrD for similar function with DoubleAc
 
#define Iprintda_double(file, name)         \
  file << indn << #name << "=" << noindent; \
  print_DynArr_double(file, name);
// See AbsArrD for similar function with DoubleAc
 
#endif