cube.h

<code cpp>
#ifndef __SW_CUBE_H__
#define __SW_CUBE_H__


#include <iostream>
//#include <iomanip>
//#include <chrono>
#include <ctime>
#include <random>
#include <vector>

// Uncomment to disable assert().
// #define NDEBUG
#include <cassert>


namespace SW {


/*
enum fill_type
{
  fill::zeros,
  fill::ones,
  fill::eye,
  fill::randu,
  fill::randn,
  fill::none
};
*/

namespace ORDER_TYPE {
  enum order_types
  {
    CDR,
    DCR,
    DRC,
    CRD,
    RCD,
    RDC
  };
}


namespace FILL_TYPE {
  enum fill_types
  {
    ZEROS,                             // Set all elements to 0.
    ONES,                              // Set all elements to 1.
    TWOS,                              // Set all elements to 2.
    EYE,                               // Set the elements along the main diagonal to 1 and off-diagonal elements to 0.
    RAND,                              // Set each element to a random value.
    RANDU,                             // Set each element to a random value from a uniform distribution in the [0,1] interval.
    RANDN,                             // Set each element to a random value from a normal/Gaussian distribution with zero mean and unit variance.
    SEQUENTIAL,                        // Set each element to the next sequential value.
    HIGHER,                            // Set each element to a higher value than the previous element.
    LOWER,                             // Set each element to a lower value than the previous element.
    NONE                               // Do not modify the elements.
  };
}


// Default random engine.
int RandRange(int min_value, int max_value)
{
  std::default_random_engine generator;
  std::uniform_int_distribution<int> distribution(min_value, max_value);
  auto random_value = distribution(generator);  // generates number in the range 1..6

  return random_value;
}


// Default random engine.
double RandRangeNormal(double min_value, double max_value)
{
  std::default_random_engine generator;
  std::uniform_real_distribution<double> distribution(min_value, max_value);
  auto random_value = distribution(generator);  // generates number in the range 1..6

  return random_value;
}


// Mersenne Twister 19937 generator.
double RandRangeMersenneTwister(double min_value, double max_value)
{
  std::random_device rd;     // only used once to initialise (seed) engine.
  std::mt19937 rng(rd());    // random-number engine used (Mersenne-Twister in this case).
  std::uniform_real_distribution<double> uni(min_value, max_value); // guaranteed unbiased.

  auto random_value = uni(rng);

  return random_value;
}


// Mersenne Twister 19937 generator.
int RandRange3(int min_value, int max_value)
{
  std::random_device rd;     // only used once to initialise (seed) engine.
  std::mt19937 rng(rd());    // random-number engine used (Mersenne-Twister in this case).
  std::uniform_int_distribution<int> uni(min_value, max_value); // guaranteed unbiased.

  auto random_value = uni(rng);

  return random_value;
}


// C / C++ standard library generator
// The following is the standard C / C++ generator defined in <cstdlib>.Beware: on some machines it has
// an unacceptably short period!!
// Program 1 : http ://www.physics.buffalo.edu/phy410-505/topic2/generators.cpp
/*
double rand(unsigned& seed, bool set)
{
if (set)
srand(seed);
else
seed = std::rand();
return seed / (RAND_MAX + 1.0);
}
*/


/*
// Park - Miller algorithm.
// The following Park - Miller generator is very popular and widely used.It produces sequences of high
// quality and has a long period 231 - 1 = 2147483647.  The following implementation uses Schrage's 
// algorithm, which works on old machines with 16 - bit int's.
// Program 1: http ://www.physics.buffalo.edu/phy410-505/topic2/generators.cpp
double ranf(unsigned& seed, bool set)
{
const int IA = 16807, IC = 2147483647, IQ = 127773, IR = 2836;
int iseed = seed;
if (iseed <= 0)
iseed += IC;
if (!set) {
int h = iseed / IQ;
int l = iseed % IQ;
iseed = IA * l - IR * h;
}
if (iseed <= 0)
iseed += IC;
return (seed = iseed) / double(IC);
}
*/

/*


// Quick and dirty algorithm
// The following implements the Quick and Dirty algorithm from Numerical Recipes[1].
// Program 1: http ://www.physics.buffalo.edu/phy410-505/topic2/generators.cpp
double ranq(unsigned& seed, bool set)
{
static unsigned long idum;
const double TWO_POWER_32 = 4294967296.0;
if (set) {
idum = (unsigned long)seed;
return idum / TWO_POWER_32;
}
idum = 1664525L * idum + 1013904223L;
seed = int(idum);
return idum / TWO_POWER_32;
}
*/

/*
// Marsaglia Xorshift algorithm
// The new Third Edition of Numerical Recipes proposes a high-quality generator based on the Xorshift
// algorithms of Marsaglia[5]. It uses 64-bit unsigned integer arithmetic that is standard on all modern
// machines.  It has an extremely long period ~ 3.13857.  Note the use of C/C++ bitwise xor (^) and right
// shift (>>) operators.
// Program 1: http://www.physics.buffalo.edu/phy410-505/topic2/generators.cpp
#ifdef _MSC_VER // Microsoft C++
typedef unsigned __int64 Ullong; // 64-bit unsigned integer
#else // Macintosh, Linux
typedef unsigned long long Ullong;
#endif
inline Ullong int64(Ullong& u, Ullong& v, Ullong& w) {
u = u * 2862933555777941757LL + 7046029254386353087LL;
v ^= v >> 17; v ^= v << 31; v ^= v >> 8;
w = 4294957665U * (w & 0xffffffff) + (w >> 32);
Ullong x = u ^ (u << 21); x ^= x >> 35; x ^= x << 4;
return (x + v) ^ w;
}
double ranx(unsigned& seed, bool set) {
static Ullong u = 0, v = 4101842887655102017LL, w = 1;
if (set) {
Ullong j = seed;
v = 4101842887655102017LL;
w = 1;
u = j ^ v;
int64(u, v, w);
v = u;
int64(u, v, w);
w = v;
seed = int(int64(u, v, w));
}
else {
6Ullong x = int64(u, v, w);
seed = unsigned(x);
return 5.42101086242752217E-20 * x;
}
}
*/

/*
#define ranf() ((float) rand() / (float) RAND_MAX)
float ranfGauss(int m, float s)
{
  static int pass = 0;
  static float y2;
  float x1, x2, w, y1;

  if (pass)
  {
    y1 = y2;
  }
  else  {
    do {
      x1 = 2.0f * ranf() - 1.0f;
      x2 = 2.0f * ranf() - 1.0f;
      w = x1 * x1 + x2 * x2;
    } while (w >= 1.0f);

    w = (float)sqrt(-2.0 * log(w) / w);
    y1 = x1 * w;
    y2 = x2 * w;
  }
  pass = !pass;

  return ((y1 * s + (float)m));
}
*/

/*
// Gaussian White Noise.
#define PI 3.1415926536f

float R1 = (float) rand() / (float) RAND_MAX;
float R2 = (float) rand() / (float) RAND_MAX;

float X = (float) sqrt( -2.0f * log( R1 )) * cos( 2.0f * PI * R2 );
*/



/*
// set the initial seed to whatever you like
static int RandSeed = 1;

// using rand() (16bit precision)
// takes about 110 seconds for 2 billion calls
float RandFloat1()
{
return ((float)rand()/RAND_MAX) * 2.0f - 1.0f;
}

// direct implementation of rand() (16 bit precision)
// takes about 32 seconds for 2 billion calls
float RandFloat2()
{
return ((float)(((RandSeed = RandSeed * 214013L + 2531011L) >> 16) & 0x7fff)/RAND_MAX) * 2.0f - 1.0f;
}

// fast rand float, using full 32bit precision
// takes about 12 seconds for 2 billion calls
float Fast_RandFloat()
{
RandSeed *= 16807;
return (float)RandSeed * 4.6566129e-010f;
}
*/

/*
// Gaussian random numbers
// This algorithm (adapted from "Natur als fraktale Grafik" by
// Reinhard Scholl) implements a generation method for gaussian
// distributed random numbers with mean=0 and variance=1
// (standard gaussian distribution) mapped to the range of
// -1 to +1 with the maximum at 0.
// For only positive results you might abs() the return value.
// The q variable defines the precision, with q=15 the smallest
// distance between two numbers will be 1/(2^q div 3)=1/10922
// which usually gives good results.

// Note: the random() function used is the standard random
// function from Delphi/Pascal that produces *linear*
// distributed numbers from 0 to parameter-1, the equivalent
// C function is probably rand().

Code :
const q=15;
c1=(1 shl q)-1;
c2=(c1 div 3)+1;
c3=1/c1;

function GRandom:single;
begin
result:=(2*(random(c2)+random(c2)+random(c2))-3*(c2-1))*c3;
end;
*/




// Calculate pseudo-random 32 bit number based on linear congruential method.
unsigned long GenerateRandomNumber(void)
{
  // Change this for different random sequences.
  static unsigned long randSeed = 22222;
  randSeed = (randSeed * 196314165) + 907633515;
  return randSeed;
}



// Mapping of 32768 integers.
int RandN(int min_value, int max_value)
{
  int result = min_value + (rand() % (int)(max_value - min_value + 1));

  return result;
}


unsigned long _Randseed = 1;
int randN(int seed)
{
  _Randseed = _Randseed * 1103515245 + seed;
  return ((unsigned int)(_Randseed >> 16) & RAND_MAX);
}


/*
int RandU(int min_value, int max_value)
{
  return min_value + (int)((double)rand() / (RAND_MAX + 1) * (max_value - min_value + 1));
}
*/
int RandU(int min_value, int max_value)
{
  return ((65539 * rand()) % (2 ^ 31));

  //return ((65539 * seed) & 0x7FFFFFFF);
}


// For 1st result use a random seed.
// For next results use prev result as the seed.
int RandU(int seed)
{
  int result = (25179 * seed + 13849) % 101;
  return result;
}

const int HALF_RAND_MAX = RAND_MAX / 2;
float randu() 
{
  // Generates a random float uniformly between [-1.0, 1.0) using RAND_MAX
  return ((float)rand()) / (HALF_RAND_MAX)-1;
}



int MinMax(int a, int b, int c)
{
  if (b < a)
    return a;
  else if (b > c)
    return c;
  else
    return b;
}


/*
// A simple taylor series, but optimized for integer phase.
// phase is in 0 -> (2^32)-1 range and maps to 0 -> ~2PI
float fastsin(UINT32 phase)
{
  const float frf3 = -1.0f / 6.0f;
  const float frf5 = 1.0f / 120.0f;
  const float frf7 = -1.0f / 5040.0f;
  const float frf9 = 1.0f / 362880.0f;
  const float f0pi5 = 1.570796327f;
  float x, x2, asin;
  UINT32 tmp = 0x3f800000 | (phase >> 7);
  if (phase & 0x40000000)
    tmp ^= 0x007fffff;
  x = (*((float*)&tmp) - 1.0f) * f0pi5;
  x2 = x * x;
  asin = ((((frf9 * x2 + frf7) * x2 + frf5) * x2 + frf3) * x2 + 1.0f) * x;
  return (phase & 0x80000000) ? -asin : asin;
}

// A simple taylor series, but optimized for integer phase.
// fastsin2 is less accurate than fastsin.
// phase is in 0 -> (2^32)-1 range and maps to 0 -> ~2PI
float fastsin2(UINT32 phase)
{
  const float frf3 = -1.0f / 6.0f;
  const float frf5 = 1.0f / 120.0f;
  const float frf7 = -1.0f / 5040.0f;
  const float f0pi5 = 1.570796327f;
  float x, x2, asin;
  UINT32 tmp = 0x3f800000 | (phase >> 7);
  if (phase & 0x40000000)
    tmp ^= 0x007fffff;
  x = (*((float*)&tmp) - 1.0f) * f0pi5;
  x2 = x * x;
  asin = (((frf7 * x2 + frf5) * x2 + frf3) * x2 + 1.0f) * x;
  return (phase & 0x80000000) ? -asin : asin;
}
*/


/*
// C++ gaussian noise generation
#include <cstdlib>
#include <ctime>

// Generate a new random seed from system time - do this once in your constructor.
srand(time(0));

// Setup constants.
const static int q = 15;
const static float c1 = (1 << q) - 1;
const static float c2 = ((int)(c1 / 3)) + 1;
const static float c3 = 1.f / c1;

// random number in range 0 - 1 not including 1.
float random = 0.f;

// the white noise.
float noise = 0.f;

for (int i = 0; i < numSamples; i++)
{
  random = ((float)rand() / (float)(RAND_MAX + 1));
  noise = (2.f * ((random * c2) + (random * c2) + (random * c2)) - 3.f * (c2 - 1.f)) * c3;
}
*/

/*
// intended for 32-bit floating point numbers only and should work a bit faster than the regular ones.

// fastabs() gives you the absolute value of a float
float fastabs(float f)
{int i=((*(int*)&f)&0x7fffffff);return (*(float*)&i);}

// fastneg() gives you the negative number (faster than multiplying with -1)
float fastneg(float f)
{int i=((*(int*)&f)^0x80000000);return (*(float*)&i);}

// fastsgn() gives back +1 for 0 or positive numbers, -1 for negative numbers
int fastsgn(float f)
{return 1+(((*(int*)&f)>>31)<<1);}
*/

/*
// uses IEEE 32-bit floating point representation to
// quickly compute approximations to the log2 of a value.
//
// Both functions return under-estimates of the actual value, although the second flavour is less of an under-estimate than the first.
//
// Try examples of 0.1, 1, 2, 5, 100.
int floorOfLn2( float f ) {
assert( f > 0. );
assert( sizeof(f) == sizeof(int) );
assert( sizeof(f) == 4 );
return (((*(int *)&f)&0x7f800000)>>23)-0x7f;
}

float approxLn2( float f ) {
assert( f > 0. );
assert( sizeof(f) == sizeof(int) );
assert( sizeof(f) == 4 );
int i = (*(int *)&f);
return (((i&0x7f800000)>>23)-0x7f)+(i&0x007fffff)/(float)0x800000;
}
*/

/*
// fastpow(f,n) gives a rather *rough* estimate of a float number f to the power of an integer number n (y=f^n). It is fast but result can be quite a bit off, since we directly mess with the floating point exponent.-> use it only for getting rough estimates of the values and where precision is not that important.
float fastpower(float f,int n)
{
long *lp,l;
lp=(long*)(&f);
l=*lp;l-=0x3F800000l;l<<=(n-1);l+=0x3F800000l;
*lp=l;
return f;
}

// fastroot(f,n) gives the n-th root of f. Same thing concerning precision applies here.
float fastroot(float f,int n)
{
long *lp,l;
lp=(long*)(&f);
l=*lp;l-=0x3F800000l;l>>=(n-1);l+=0x3F800000l;
*lp=l;
return f;
}
*/





template <typename T> class Cube
{
private:
  std::vector<std::vector<std::vector<T> > > data;
  mutable unsigned rows;               // Width.
  mutable unsigned cols;               // Height.
  mutable unsigned depth;              // Depth.

  FILL_TYPE::fill_types fill_type;
  ORDER_TYPE::order_types order_type;

public:
  Cube();
  Cube(unsigned _rows, unsigned _cols, unsigned _depth);
  Cube(unsigned _rows, unsigned _cols, unsigned _depth, FILL_TYPE::fill_types _fill_type);
  Cube(unsigned _rows, unsigned _cols, unsigned _depth, FILL_TYPE::fill_types _fill_type, ORDER_TYPE::order_types _order_type);
  Cube(unsigned _rows, unsigned _cols, unsigned _depth, ORDER_TYPE::order_types _order_type);
  ~Cube();
  
  Cube(const Cube<T>& rhs);

  
  // Arithmetic operators and Relational operators
  const Cube<T>& operator = (const Cube<T> &);  // Assignment operator.

  // Operator overloading, for comparitive operations.
  bool operator==(const Cube<T>& rhs) const;   // Comparison operator.
  bool operator!=(const Cube<T>& rhs) const;   // 
  bool operator<(const Cube<T>& rhs) const;
  bool operator<=(const Cube<T>& rhs) const;
  bool operator>(const Cube<T>& rhs) const;
  bool operator>=(const Cube<T>& rhs) const;

  // Arithmetic operators.
  Cube<T> operator - () const;         // Negate  operator.
  Cube<T> operator ++ ();              // Prefix  increment operator.
  Cube<T> operator ++ (int);           // Postfix increment operator.
  Cube<T> operator -- ();              // Prefix  decrement operator.
  Cube<T> operator -- (int);           // Postfix decrement operator.

  // Mathematical operations.
  Cube<T> operator + (const Cube<T> &);
  Cube<T> operator += (const Cube<T> &);
  Cube<T> operator - (const Cube<T> &);
  Cube<T> operator -= (const Cube<T> &);
  Cube<T> operator * (const Cube<T> &);
  Cube<T> operator *= (const Cube<T> &);
  Cube<T> operator / (const Cube<T> &);
  Cube<T> operator /= (const Cube<T> &);
  Cube<T> operator % (const Cube<T> &);
  Cube<T> operator %= (const Cube<T> &);
  Cube<T> operator ^ (const Cube<T> &);
  Cube<T> operator ^= (const Cube<T> &);

  // Scalar operations.
  Cube<T> operator+(const T& rhs);
  Cube<T> operator+=(const T& rhs);
  Cube<T> operator-(const T& rhs);
  Cube<T> operator-=(const T& rhs);
  Cube<T> operator*(const T& rhs);
  Cube<T> operator*=(const T& rhs);
  Cube<T> operator/(const T& rhs);
  Cube<T> operator/=(const T& rhs);
  Cube<T> operator%(const T& rhs);
  Cube<T> operator%=(const T& rhs);
  Cube<T> operator^(const T& rhs);
  Cube<T> operator^=(const T& rhs);
  
  // Access the individual elements.
  T& operator()(const unsigned& row, const unsigned& col, const unsigned& depth);
  const T& operator()(const unsigned& row, const unsigned& col, const unsigned& depth) const;

  void fill();
  void fill(T value);
  void fill_higher(T seq = 0, T min_inc = 1);
  void fill_lower(T seq = 0, T min_inc = 1);
  void fill_sequentially(T seq = 0, T inc = 1);

  // Access the row, column and depth sizes.
  unsigned getCols() const;
  unsigned getDepth() const;
  FILL_TYPE::fill_types getFillType() const;
  ORDER_TYPE::order_types getOrderType() const;
  unsigned getRows() const;
  unsigned getSize() const;
  unsigned getVolume();
  T& max();
  T& min();
  void print(std::ostream& os = cout);
  void randomize();
  void read(std::istream& is = cin);
  void sequence(T seq, T inc=1);
  T& sum();

  friend Cube<T> operator + (const Cube<T> &, const Cube<T> &);
  friend Cube<T> operator - (const Cube<T> &, const Cube<T> &);
  friend Cube<T> operator * (const Cube<T> &, const Cube<T> &);
  friend Cube<T> operator / (const Cube<T> &, const Cube<T> &);
  friend Cube<T> operator % (const Cube<T> &, const Cube<T> &);
  friend Cube<T> operator ^ (const Cube<T> &, const Cube<T> &);

  friend int operator == (const Cube<T> &, const Cube<T> &);
  friend int operator != (const Cube<T> &, const Cube<T> &);
  friend int operator <  (const Cube<T> &, const Cube<T> &);
  friend int operator <= (const Cube<T> &, const Cube<T> &);
  friend int operator >  (const Cube<T> &, const Cube<T> &);
  friend int operator >= (const Cube<T> &, const Cube<T> &);

  // Other functions
  friend Cube<T> abs(const Cube<T> &);
  friend double div(const Cube<T> &, const Cube<T> &);
  friend unsigned size(const Cube<T> &);
  friend Cube<T> sqrt(const Cube<T> &);
  friend Cube<T> pow(const Cube<T> &, const Cube<T> &);


  // I/O stream functions
  //template<typename T>
  friend std::ostream & operator << (std::ostream &, const Cube<T> &);

  //template<typename T>
  friend std::istream & operator >> (std::istream &, Cube<T> &);
};



// Basic Constructor.
template<typename T>
Cube<T>::Cube()
{
  rows = 0;
  cols = 0;
  depth = 0;
  fill_type = FILL_TYPE::NONE;
  order_type = ORDER_TYPE::CDR;
}


// Parameter Constructor.
template<typename T>
Cube<T>::Cube(unsigned _rows, unsigned _cols, unsigned _depth)
{
  data.resize(_rows);
  for (unsigned i = 0; i<data.size(); i++)
  {
    data[i].resize(_cols);
    for (unsigned j = 0; j < data[i].size(); j++)
    {
      data[i][j].resize(_depth);
    }
  }

  rows = _rows;
  cols = _cols;
  depth = _depth;
  fill_type = FILL_TYPE::NONE;
  order_type = ORDER_TYPE::CDR;
}


// Parameter Constructor.
template<typename T>
Cube<T>::Cube(unsigned _rows, unsigned _cols, unsigned _depth, FILL_TYPE::fill_types _fill_type)
{
  data.resize(_rows);
  for (unsigned i = 0; i<data.size(); i++)
  {
    data[i].resize(_cols);
    for (unsigned j = 0; j < data[i].size(); j++)
    {
      data[i][j].resize(_depth);
    }
  }

  rows = _rows;
  cols = _cols;
  depth = _depth;
  fill_type = _fill_type;
  order_type = ORDER_TYPE::CDR;

  fill();                              // Populate grid depending on FILLS type.
}



// Parameter Constructor.
template<typename T>
Cube<T>::Cube(unsigned _rows, unsigned _cols, unsigned _depth, FILL_TYPE::fill_types _fill_type, ORDER_TYPE::order_types _order_type)
{
  data.resize(_rows);
  for (unsigned i = 0; i<data.size(); i++)
  {
    data[i].resize(_cols);
    for (unsigned j = 0; j < data[i].size(); j++)
    {
      data[i][j].resize(_depth);
    }
  }

  rows = _rows;
  cols = _cols;
  depth = _depth;
  fill_type = _fill_type;
  order_type = _order_type;

  fill();                              // Populate grid depending on FILLS type.
}


// Parameter Constructor.
template<typename T>
Cube<T>::Cube(unsigned _rows, unsigned _cols, unsigned _depth, ORDER_TYPE::order_types _order_type)
{
  data.resize(_rows);
  for (unsigned i = 0; i<data.size(); i++)
  {
    data[i].resize(_cols);
    for (unsigned j = 0; j < data[i].size(); j++)
    {
      data[i][j].resize(_depth);
    }
  }

  rows = _rows;
  cols = _cols;
  depth = _depth;
  fill_type = FILL_TYPE::none;
  order_type = _order_type;
}


// Destructor.
template <typename T>
Cube<T>::~Cube()
{
  // does nothing
}


// Copy Constructor.
template<typename T>
Cube<T>::Cube(const Cube<T>& rhs)
{
  data = rhs.data;
  rows = rhs.getRows();
  cols = rhs.getCols();
  depth = rhs.getDepth();
  fill_type = rhs.getFillType();
  order_type = rhs.getOrderype();
}


// Assignment Operator.
template<typename T>
const Cube<T>& Cube<T>::operator=(const Cube<T>& rhs)
{
  if (&rhs == this)
    return *this;

  unsigned new_cols = rhs.getCols();
  unsigned new_depth = rhs.getDepth();
  unsigned new_rows = rhs.getRows();
  bool new_fill_type = rhs.getFillType();
  bool new_order_type = rhs.getOrderType();

  data.resize(new_rows);
  for (unsigned i = 0; i < data.size(); i++)
  {
    data[i].resize(new_cols);
    for (unsigned i = 0; i < data.size(); i++)
    {
      data[i][j].resize(new_depth);
    }
  }
  /*
  for (unsigned i = 0; i<new_rows; i++)
  {
    for (unsigned j = 0; j<new_cols; j++)
    {
      data[i][j] = rhs(i, j);
    }
  }
  */
  for (unsigned i = 0; i<new_rows; i++)
  {
    for (unsigned j = 0; j<new_cols; j++)
    {
      data[i][j] = rhs(i, j);

      for (unsigned k = 0; k < new_depth; k++)
      {
        data[i][j][k] = rhs(i, j, k);
      }
    }
  }

  rows = new_rows;
  cols = new_cols;
  depth = new_depth;
  fill_type = new_fill_type;
  order_type = new_order_type;

  return *this;
}


// Comparison Operator ==.
template<typename T>
bool Cube<T>::operator==(const Cube<T>& rhs) const
{
  unsigned cols = rhs.getCols();
  unsigned depth = rhs.getDepth();
  unsigned rows = rhs.getRows();
  fill_types fill_type = rhs.getFillType();
  order_types order_type = rhs.getOrderType();

  if (getCols() != cols)
    return false;

  if (getDepth() != depth)
    return false;

  if (getRows() != rows)
    return false;

  if (getFillType() != fill_type)
    return false;

  if (getOrderType() != order_type)
    return false;

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<depth; k++)
      {
        if (this->data[i][j][k] != rhs(i, j, k))
          return false;
      }
    }
  }

  return true;
}


// Comparison Operator !=.
template<typename T>
bool Cube<T>::operator!=(const Cube<T>& rhs) const
{
  unsigned rows = rhs.getRows();
  unsigned cols = rhs.getCols();
  unsigned depth = rhs.getDepth();
  fill_types fill_type = rhs.getFillType();
  order_types order_type = rhs.getOrderType();

  if (getCols() != cols)
    return true;

  if (getDepth() != depth)
    return true;

  if (getRows() != rows)
    return true;

  if (getFillType() != fill_type)
    return true;

  if (getOrderType() != order_type)
    return true;

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<depth; k++)
      {
        if (this->data[i][j][k] == rhs(i, j, k))
          return false;
      }
    }
  }

  return false;
}


// Comparison Operator <.
//
// If every cell is < than rhs.
template<typename T>
bool Cube<T>::operator<(const Cube<T>& rhs) const
{
  unsigned rows = rhs.getRows();
  unsigned cols = rhs.getCols();
  unsigned depth = rhs.getDepth();

  if (getCols() != cols)
    return false;

  if (getDepth() != depth)
    return false;

  if (getRows() != rows)
    return false;

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<depth; k++)
      {
        if (this->data[i][j][k] >= rhs(i, j, k))
          return false;
      }
    }
  }

  return true;
}


// Comparison Operator <=.
//
// If every cell is <= than rhs.
template<typename T>
bool Cube<T>::operator<=(const Cube<T>& rhs) const
{
  unsigned rows = rhs.getRows();
  unsigned cols = rhs.getCols();
  unsigned depth = rhs.getDepth();

  if (getRows() != rows)
    return false;

  if (getCols() != cols)
    return false;

  if (getDepth() != depth)
    return false;

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<depth; k++)
      {
        if (this->data[i][j][k] > rhs(i, j, k))
          return false;
      }
    }
  }

  return true;
}


// Comparison Operator <.
//
// If every cell is > than rhs.
template<typename T>
bool Cube<T>::operator>(const Cube<T>& rhs) const
{
  unsigned rows = rhs.getRows();
  unsigned cols = rhs.getCols();
  unsigned depth = rhs.getDepth();

  if (getRows() != rows)
    return false;

  if (getCols() != cols)
    return false;

  if (getDepth() != depth)
    return false;

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<depth; k++)
      {
        if (this->data[i][j][k] <= rhs(i, j, k))
          return false;
      }
    }
  }

  return true;
}


// Comparison Operator >=.
//
// If every cell is >= than rhs.
template<typename T>
bool Cube<T>::operator>=(const Cube<T>& rhs) const
{
  unsigned rows = rhs.getRows();
  unsigned cols = rhs.getCols();
  unsigned depth = rhs.getDepth();

  if (getRows() != rows)
    return false;

  if (getCols() != cols)
    return false;

  if (getDepth() != depth)
    return false;

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<depth; k++)
      {
        if (this->data[i][j][k] < rhs(i, j, k))
          return false;
      }
    }
  }

  return true;
}


// Unary - operator.
template<typename T>
Cube<T> Cube<T>::operator -() const
{
  Cube result(rows, cols, depth);

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k < depth; k++)
      {
        if (result(i, j, k) != 0)
          result(i, j, k) = -this->data[i][j][k];
      }
    }
  }

  return result;
}


// Prefix increment operator.
template<typename T>
Cube<T> Cube<T>::operator ++ ()
{
  Cube result(rows, cols, depth);

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<depth; k++)
      {
        result(i, j, k) = this->data[i][j][k] + 1;
      }
    }
  }

  return *this = result;
}


// Postfix increment operator.
template<typename T>
Cube<T> Cube<T>::operator ++ (int)
{
  Cube result(rows, cols, depth);

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<depth; k++)
      {
        result(i, j, k) = this->data[i][j][k] + 1;
      }
    }
  }

  return result;
}


// Prefix decrement operator.
template<typename T>
Cube<T> Cube<T>::operator -- ()
{
  Cube result(rows, cols, depth);

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<depth; k++)
      {
        result(i, j, k) = this->data[i][j][k] - 1;
      }
    }
  }

  return *this = result;
}


// Postfix decrement operator.
template<typename T>
Cube<T> Cube<T>::operator -- (int)
{
  Cube result(rows, cols, depth);

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<depth; k++)
      {
        result(i, j, k) = this->data[i][j][k] - 1;
      }
    }
  }

  return result;
}


// Addition of two Cubes.
template<typename T>
Cube<T> Cube<T>::operator+(const Cube<T>& rhs)
{
  float new_rows = (x > rhs.getRows()) ? x : rhs.getRows();
  float new_cols = (y > rhs.getCols()) ? y : rhs.getCols();
  float new_depth = z + rhs.getDepth();
  Cube result(new_rows, new_cols, new_depth);

  for (unsigned i = 0; i<new_rows; i++)
  {
    for (unsigned j = 0; j<new_cols; j++)
    {
      for (unsigned k = 0; k<new_depth; k++)
      {
        //unsigned temp_i = 
        result(i, j, k) = this->data[i][j][k] + rhs(i, j, k);
      }
    }
  }

  return result;
}


// Cumulative addition of this Cube and another.
template<typename T>
Cube<T> Cube<T>::operator+=(const Cube<T>& rhs)
{
  unsigned rows = rhs.getRows();
  unsigned cols = rhs.getCols();
  unsigned depth = rhs.getDepth();

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<depth; k++)
      {
        this->data[i][j][k] += rhs(i, j, k);
      }
    }
  }

  return *this;
}


// Subtraction of this Cube and another.
template<typename T>
Cube<T> Cube<T>::operator-(const Cube<T>& rhs)
{
  unsigned rows = rhs.getRows();
  unsigned cols = rhs.getCols();
  unsigned depth = rhs.getDepth();

  Cube result(rows, cols, depth);

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<depth; k++)
      {
        result(i, j, k) = this->data[i][j][k] - rhs(i, j, k);
      }
    }
  }

  return result;
}


// Cumulative subtraction of this Cube and another.
template<typename T>
Cube<T> Cube<T>::operator-=(const Cube<T>& rhs)
{
  unsigned rows = rhs.getRows();
  unsigned cols = rhs.getCols();
  unsigned depth = rhs.getDepth();

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<depth; k++)
      {
        this->data[i][j][k] -= rhs(i, j, k);
      }
    }
  }

  return *this;
}


// Left multiplication of this Cube and another.
template<typename T>
Cube<T> Cube<T>::operator*(const Cube<T>& rhs)
{
  unsigned rows = rhs.getRows();
  unsigned cols = rhs.getCols();
  unsigned depth = rhs.getDepth();

  Cube result(rows, cols, depth);

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<depth; k++)
      {
        result(i, j, k) += this->data[i][j][k] * rhs(i, j, k);
      }
    }
  }

  return result;
}


// Cumulative left multiplication of this Cube and another.
template<typename T>
Cube<T> Cube<T>::operator*=(const Cube<T>& rhs)
{
  Cube result = (*this) * rhs;
  (*this) = result;
  return *this;
}


// Left division of this Cube and another.
template<typename T>
Cube<T> Cube<T>::operator/(const Cube<T>& rhs)
{
  unsigned rows = rhs.getRows();
  unsigned cols = rhs.getCols();
  unsigned depth = rhs.getDepth();

  Cube result(rows, cols, depth);

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<depth; k++)
      {
        if (rhs(i, j, k) != 0)
          result(i, j, k) += this->data[i][j][k] / rhs(i, j, k);
      }
    }
  }

  return result;
}


// Cumulative left division of this Cube and another.
template<typename T>
Cube<T> Cube<T>::operator/=(const Cube<T>& rhs)
{
  Cube result = (*this) / rhs;
  (*this) = result;
  return *this;
}


// Left mod of this Cube and another.
template<typename T>
Cube<T> Cube<T>::operator%(const Cube<T>& rhs)
{
  unsigned rows = rhs.getRows();
  unsigned cols = rhs.getCols();
  unsigned depth = rhs.getDepth();

  Cube result(rows, cols, depth);

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<rows; k++)
      {
        if (rhs(i, j, k) != 0)
          result(i, j, k) += this->data[i][j][k] % rhs(i, j, k);
      }
    }
  }

  return result;
}


// Cumulative left mod of this Cube and another.
template<typename T>
Cube<T> Cube<T>::operator%=(const Cube<T>& rhs)
{
  Cube result = (*this) % rhs;
  (*this) = result;
  return *this;
}


// Left raise of this Cube and another.
template<typename T>
Cube<T> Cube<T>::operator^(const Cube<T>& rhs)
{
  //	unsigned rows = rhs.get_rows();
  //	unsigned cols = rhs.get_cols();
  Matrix result(rows, cols, depth);

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<depth; k++)
      {
        result(i, j, k) += this->data[i][j][k] ^ rhs(i, j, k);
      }
    }
  }

  return result;
}


// Cumulative left raise of this Cube and another.
template<typename T>
Cube<T> Cube<T>::operator^=(const Cube<T>& rhs)
{
  Cube result = (*this) ^ rhs;
  (*this) = result;
  return *this;
}


// Addition of this Cube and a value.
template<typename T>
Cube<T> Cube<T>::operator+(const T& rhs)
{
  Cube result(rows, cols, depth);

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k < depth; k++)
      {
        result(i, j, k) = this->data[i][j][k] + rhs;
      }
    }
  }

  return result;
}


// Cumulative addition of this Cube and a value.
template<typename T>
Cube<T> Cube<T>::operator+=(const T& rhs)
{
  Cube result = (*this) + rhs;
  (*this) = result;
  return *this;
}


// Subtraction of this Cube and a value.
template<typename T>
Cube<T> Cube<T>::operator-(const T& rhs)
{
  unsigned rows = rhs.getRows();
  unsigned cols = rhs.getCols();
  unsigned depth = rhs.getDepth();

  Cube result(rows, cols, depth);

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k < depth; k++)
      {
        result(i, j, k) = this->data[i][j][k] - rhs;
      }
    }
  }

  return result;
}


// Cumulative subtraction of a valye from this Cube.
template<typename T>
Cube<T> Cube<T>::operator-=(const T& rhs)
{
  Cube result = (*this) - rhs;
  (*this) = result;
  return *this;
}


// Left multiplication of this Cube and a value.
template<typename T>
Cube<T> Cube<T>::operator*(const T& rhs)
{
  unsigned rows = rhs.getRows();
  unsigned cols = rhs.getCols();
  unsigned depth = rhs.getDepth();

  Cube result(rows, cols, depth);

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<depth; k++)
      {
        result(i, j, k) += this->data[i][j][k] * rhs;
      }
    }
  }

  return result;
}


// Cumulative multiplication from this Cube and a value.
template<typename T>
Cube<T> Cube<T>::operator*=(const T& rhs)
{
  Cube result = (*this) * rhs;
  (*this) = result;
  return *this;
}


// Left division of this Cube and a value.
template<typename T>
Cube<T> Cube<T>::operator/(const T& rhs)
{
  unsigned rows = rhs.getRows();
  unsigned cols = rhs.getCols();
  unsigned depth = rhs.getDepth();

  Cube result(rows, cols, depth);

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<depth; k++)
      {
        if (rhs != 0)
          result(i, j, k) += this->data[i][j][k] / rhs;
      }
    }
  }

  return result;
}


// Cumulative division from this Cube and a value.
template<typename T>
Cube<T> Cube<T>::operator/=(const T& rhs)
{
  if (rhs != 0)
  {
    Cube result = (*this) / rhs;
    (*this) = result;
  }

  return *this;
}


// Left mod of this Cube and a value.
template<typename T>
Cube<T> Cube<T>::operator%(const T& rhs)
{
  unsigned rows = rhs.getRows();
  unsigned cols = rhs.getCols();
  unsigned depth = rhs.getDepth();

  Cube result(rows, cols, depth);

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<rows; k++)
      {
        if (rhs != 0)
          result(i, j, k) += this->data[i][j][k] % rhs;
      }
    }
  }

  return result;
}


// Cumulative left mod of this Cube and a value.
template<typename T>
Cube<T> Cube<T>::operator%=(const T& rhs)
{
  if (rhs != 0)
  {
    Cube result = (*this) % rhs;
    (*this) = result;
  }

  return *this;
}


// Left raise of this Cube and a value.
template<typename T>
Cube<T> Cube<T>::operator^(const T& rhs)
{
  //	unsigned rows = rhs.get_rows();
  //	unsigned cols = rhs.get_cols();
  Matrix result(rows, cols, depth);

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<depth; k++)
      {
        result(i, j, k) += this->data[i][j][k] ^ rhs;
      }
    }
  }

  return result;
}


// Cumulative left raise of this Cube and a value..
template<typename T>
Cube<T> Cube<T>::operator^=(const T& rhs)
{
  Cube result = (*this) ^ rhs;
  (*this) = result;
  return *this;
}


// Access the individual elements.
template<typename T>
T& Cube<T>::operator()(const unsigned& row, const unsigned& col, const unsigned& depth)
{
  assert(row >= 0);
  assert(col >= 0);
  assert(depth >= 0);
  assert(row < this->rows);
  assert(col < this->cols);
  assert(depth < this->depth);

  return this->data[row][col][depth];
}


// Access the individual elements (const).
template<typename T>
const T& Cube<T>::operator()(const unsigned& row, const unsigned& col, const unsigned& depth) const
{
  assert(row >= 0);
  assert(col >= 0);
  assert(depth >= 0);
  assert(row < this->rows);
  assert(col < this->cols);
  assert(depth < this->depth);

  return this->data[row][col][depth];
}


// Populates every cell with a specific value which is dependant on the FILL_TYPE.
template<typename T>
void Cube<T>::fill()
{
  std::srand(std::time(0)); // use current time as seed for random generator.

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<depth; k++)
      {
        switch (fill_type)
        {
        case 0: // Zeros.
          data[i][j][k] = 0;
          break;

        case 1: // Ones.
          data[i][j][k] = 1;
          break;

        case 2: // Twos.
          data[i][j][k] = 2;
          break;

        case 3: // Eye.
          break;

        case 4: // Rand.
          //std::srand(std::time(0)); // use current time as seed for random generator
          data[i][j][k] = rand();
          break;

        case 5: // Randu.
          //data[i][j][k] = RandRangeNormal(DBL_MIN, DBL_MAX);
          //data[i][j][k] = RandRange(INT_MIN, INT_MAX);
          //data[i][j][k] = RandRange(DBL_MIN, DBL_MAX);
          //data[i][j][k] = RandU(DBL_MIN, DBL_MAX);
          data[i][j][k] = RandU(INT_MIN, INT_MAX);
          //data[i][j][k] = RandU(0, INT_MAX);
          //data[i][j][k] = RandU(0, INT_MAX);
          //RandRangeMersenneTwister
          break;

        case 6: // Randn.
          //data[i][j][k] = RandRangeMersenneTwister(DBL_MIN, DBL_MAX);
          //data[i][j][k] = RandN(DBL_MIN, DBL_MAX);
          //data[i][j][k] = RandN(INT_MIN, INT_MAX);
          data[i][j][k] = RandN(0, INT_MAX);
          break;

        case 7: // Sequential.
          fill_sequentially();
          return;
          break;

        case 8: // Higher.
          fill_higher();
          return;
          break;

        case 9: // Lower.
          fill_lower();
          return;
          break;

        case 10: // None.
          break;

        Default: // None.
          return;
        }


      }
    }
  }
}


// Populates every cell with the specific value.
template<typename T>
void Cube<T>::fill(T value)
{
  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<depth; k++)
      {
        data[i][j][k] = value;
      }
    }
  }
}


// Populates every cell with the a higher number than the previous cell.
//
// seq is the starting sequence.
// min_inc is the increment to the next number.
template<typename T>
void Cube<T>::fill_higher(T seq = 0, T min_inc = 1)
{
  T tmp;

  switch (order_type)
  {
  case 0: // CDR.
    for (unsigned j = 0; j < cols; j++)
    {
      for (unsigned k = 0; k < depth; k++)
      {
        for (unsigned i = 0; i < rows; i++)
        {
          data[i][j][k] = seq;
          
          tmp = rand();
          while (tmp <= seq + min_inc)
            tmp = rand();
          seq = tmp;
        }
      }
    }
    break;

  case 1: // DCR.
    for (unsigned k = 0; k < depth; k++)
    {
      for (unsigned j = 0; j < cols; j++)
      {
        for (unsigned i = 0; i < rows; i++)
        {
          data[i][j][k] = seq;

          tmp = rand();
          while (tmp <= seq + min_inc)
            tmp = rand();
          seq = tmp;
        }
      }
    }
    break;

  case 2: // DRC.
    for (unsigned k = 0; k < depth; k++)
    {
      for (unsigned i = 0; i < rows; i++)
      {
        for (unsigned j = 0; j < cols; j++)
        {
          data[i][j][k] = seq;

          tmp = rand();
          while (tmp <= seq + min_inc)
            tmp = rand();
          seq = tmp;
        }
      }
    }
    break;

  case 3: // CRD.
    for (unsigned j = 0; j < cols; j++)
    {
      for (unsigned i = 0; i < rows; i++)
      {
        for (unsigned k = 0; k < depth; k++)
        {
          data[i][j][k] = seq;

          tmp = rand();
          while (tmp <= seq + min_inc)
            tmp = rand();
          seq = tmp;
        }
      }
    }
    break;

  case 4: // RCD.
    for (unsigned i = 0; i < rows; i++)
    {
      for (unsigned j = 0; j < cols; j++)
      {
        for (unsigned k = 0; k < depth; k++)
        {
          data[i][j][k] = seq;

          tmp = rand();
          while (tmp <= seq + min_inc)
            tmp = rand();
          seq = tmp;
        }
      }
    }
    break;

  case 5: // RDC.
    for (unsigned i = 0; i < rows; i++)
    {
      for (unsigned k = 0; k < depth; k++)
      {
        for (unsigned j = 0; j < cols; j++)
        {
          data[i][j][k] = seq;

          tmp = rand();
          while (tmp <= seq + min_inc)
            tmp = rand();
          seq = tmp;
        }
      }
    }
    break;

  default: // CDR.
    for (unsigned j = 0; j < cols; j++)
    {
      for (unsigned k = 0; k < depth; k++)
      {
        for (unsigned i = 0; i < rows; i++)
        {
          data[i][j][k] = seq;

          tmp = rand();
          while (tmp <= seq + min_inc)
            tmp = rand();
          seq = tmp;
        }
      }
    }
  }
}


// Populates every cell with the a lower number than the previous cell.
//
// seq is the starting sequence.
// min_inc is the minimum increment to the next number.
template<typename T>
void Cube<T>::fill_lower(T seq = 0, T min_inc = 1)
{
  T tmp;

  switch (order_type)
  {
  case 0: // CDR.
    for (unsigned j = 0; j < cols; j++)
    {
      for (unsigned k = 0; k < depth; k++)
      {
        for (unsigned i = 0; i < rows; i++)
        {
          data[i][j][k] = seq;

          tmp = rand();
          while (tmp >= seq + min_inc)
            tmp = rand();
          seq = tmp;
        }
      }
    }
    break;

  case 1: // DCR.
    for (unsigned k = 0; k < depth; k++)
    {
      for (unsigned j = 0; j < cols; j++)
      {
        for (unsigned i = 0; i < rows; i++)
        {
          data[i][j][k] = seq;

          tmp = rand();
          while (tmp >= seq + min_inc)
            tmp = rand();
          seq = tmp;
        }
      }
    }
    break;

  case 2: // DRC.
    for (unsigned k = 0; k < depth; k++)
    {
      for (unsigned i = 0; i < rows; i++)
      {
        for (unsigned j = 0; j < cols; j++)
        {
          data[i][j][k] = seq;

          tmp = rand();
          while (tmp >= seq + min_inc)
            tmp = rand();
          seq = tmp;
        }
      }
    }
    break;

  case 3: // CRD.
    for (unsigned j = 0; j < cols; j++)
    {
      for (unsigned i = 0; i < rows; i++)
      {
        for (unsigned k = 0; k < depth; k++)
        {
          data[i][j][k] = seq;

          tmp = rand();
          while (tmp >= seq + min_inc)
            tmp = rand();
          seq = tmp;
        }
      }
    }
    break;

  case 4: // RCD.
    for (unsigned i = 0; i < rows; i++)
    {
      for (unsigned j = 0; j < cols; j++)
      {
        for (unsigned k = 0; k < depth; k++)
        {
          data[i][j][k] = seq;

          tmp = rand();
          while (tmp >= seq + min_inc)
            tmp = rand();
          seq = tmp;
        }
      }
    }
    break;

  case 5: // RDC.
    for (unsigned i = 0; i < rows; i++)
    {
      for (unsigned k = 0; k < depth; k++)
      {
        for (unsigned j = 0; j < cols; j++)
        {
          data[i][j][k] = seq;

          tmp = rand();
          while (tmp >= seq + min_inc)
            tmp = rand();
          seq = tmp;
        }
      }
    }
    break;

  default: // CDR.
    for (unsigned j = 0; j < cols; j++)
    {
      for (unsigned k = 0; k < depth; k++)
      {
        for (unsigned i = 0; i < rows; i++)
        {
          data[i][j][k] = seq;

          tmp = rand();
          while (tmp >= seq + min_inc)
            tmp = rand();
          seq = tmp;
        }
      }
    }
  }
}



// Populates every cell with the a sequential number.
//
// seq is the starting sequence.
// inc is the increment to the next number.
template<typename T>
void Cube<T>::fill_sequentially(T seq = 0, T inc = 1)
{
  switch (order_type)
  {
  case 0: // CDR.
    for (unsigned j = 0; j < cols; j++)
    {
      for (unsigned k = 0; k < depth; k++)
      {
        for (unsigned i = 0; i < rows; i++)
        {
          data[i][j][k] = seq;
          seq += inc;
        }
      }
    }
    break;

  case 1: // DCR.
    for (unsigned k = 0; k < depth; k++)
    {
      for (unsigned j = 0; j < cols; j++)
      {
        for (unsigned i = 0; i < rows; i++)
        {
          data[i][j][k] = seq;
          seq += inc;
        }
      }
    }
    break;

  case 2: // DRC.
    for (unsigned k = 0; k < depth; k++)
    {
      for (unsigned i = 0; i < rows; i++)
      {
        for (unsigned j = 0; j < cols; j++)
        {
          data[i][j][k] = seq;
          seq += inc;
        }
      }
    }
    break;

  case 3: // CRD.
    for (unsigned j = 0; j < cols; j++)
    {
      for (unsigned i = 0; i < rows; i++)
      {
        for (unsigned k = 0; k < depth; k++)
        {
          data[i][j][k] = seq;
          seq += inc;
        }
      }
    }
    break;

  case 4: // RCD.
    for (unsigned i = 0; i < rows; i++)
    {
      for (unsigned j = 0; j < cols; j++)
      {
        for (unsigned k = 0; k < depth; k++)
        {
          data[i][j][k] = seq;
          seq += inc;
        }
      }
    }
    break;

  case 5: // RDC.
    for (unsigned i = 0; i < rows; i++)
    {
      for (unsigned k = 0; k < depth; k++)
      {
        for (unsigned j = 0; j < cols; j++)
        {
          data[i][j][k] = seq;
          seq += inc;
        }
      }
    }
    break;

  default: // CDR.
    for (unsigned j = 0; j < cols; j++)
    {
      for (unsigned k = 0; k < depth; k++)
      {
        for (unsigned i = 0; i < rows; i++)
        {
          data[i][j][k] = seq;
          seq += inc;
        }
      }
    }
  }
}


// Get the number of columns.
template<typename T>
unsigned Cube<T>::getCols() const
{
  return this->cols;
}


// Get the number of columns.
template<typename T>
unsigned Cube<T>::getDepth() const
{
  return this->depth;
}


// Gets the fill_type.
//
// Cells are populated using the fill_type.
template<typename T>
FILL_TYPE::fill_types Cube<T>::getFillType() const
{
  return this->fill_type;
}


// Gets the fill_type.
//
// Cells are populated in the order determined by the order_type.
template<typename T>
ORDER_TYPE::order_types Cube<T>::getOrderType() const
{
  return this->order_type;
}


// Get the number of rows.
template<typename T>
unsigned Cube<T>::getRows() const
{
  return this->rows;
}


// Returns the size.
//
// Returns the volume.
template <typename T>
unsigned Cube<T>::getSize() const
{
  return rows * cols * depth;
}


// Returns the volume.
template <typename T>
unsigned Cube<T>::getVolume()
{
  return rows * cols * depth;
}


// Returns the maximum value.
template <typename T>
T& Cube<T>::max()
{
  T result = data[0][0][0];

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<depth; k++)
      {
        if (data[i][j][k] > result)
          result = data[i][j][k];
      }
    }
  }

  return result;
}


// Returns the minimum value.
template <typename T>
T& Cube<T>::min()
{
  T result = data[0][0][0];

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<depth; k++)
      {
        if (data[i][j][k] < result)
          result = data[i][j][k];
      }
    }
  }

  return result;
}


// Prints the cube.
template <typename T>
void Cube<T>::print(std::ostream& os)
{
  unsigned rows = getRows();
  unsigned cols = getCols();
  unsigned depth = getDepth();

  for (unsigned j = 0; j < cols; j++)
  {
    for (unsigned k = 0; k < depth; k++)
    {
      for (unsigned i = 0; i < rows; i++)
      {
        os << data[i][j][k];
        //if ((i*j) < (((rows - 1)*(cols - 1)) - 1))
        //if ((i*j*k) < (((rows - 1)*(cols - 1)*(depth -1)) - 1))
        if ((i < rows - 1) || (j < cols - 1) || (k < depth -1))
          os << ",";
      }

      os << "  ";
    }

    os << std::endl;
  }
}


// Randomizes every cell.
template<typename T>
void Cube<T>::randomize()
{
  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<depth; k++)
      {
        data[i][j][k] = rand();
      }
    }
  }
}


// Reads the input stream into the data element.
template <typename T>
void Cube<T>::read(std::istream& is)
{
  is >> data;
}


// Populates every cell with the a sequential number.
//
// seq is the starting sequence.
// inc is the increment to the next number.
template<typename T>
void Cube<T>::sequence(T seq = 0, T inc = 1)
{
  switch (order_type)
  {
  case 0: // CDR.
    for (unsigned j = 0; j < cols; j++)
    {
      for (unsigned k = 0; k < depth; k++)
      {
        for (unsigned i = 0; i < rows; i++)
        {
          data[i][j][k] = seq;
          seq += inc;
        }
      }
    }
    break;

  case 1: // DCR.
    for (unsigned k = 0; k < depth; k++)
    {
      for (unsigned j = 0; j < cols; j++)
      {
        for (unsigned i = 0; i < rows; i++)
        {
          data[i][j][k] = seq;
          seq += inc;
        }
      }
    }
    break;

  case 2: // DRC.
    for (unsigned k = 0; k < depth; k++)
    {
      for (unsigned i = 0; i < rows; i++)
      {
        for (unsigned j = 0; j < cols; j++)
        {
          data[i][j][k] = seq;
          seq += inc;
        }
      }
    }
    break;

  case 3: // CRD.
    for (unsigned j = 0; j < cols; j++)
    {
      for (unsigned i = 0; i < rows; i++)
      {
        for (unsigned k = 0; k < depth; k++)
        {
          data[i][j][k] = seq;
          seq += inc;
        }
      }
    }
    break;

  case 4: // RCD.
    for (unsigned i = 0; i < rows; i++)
    {
      for (unsigned j = 0; j < cols; j++)
      {
        for (unsigned k = 0; k < depth; k++)
        {
          data[i][j][k] = seq;
          seq += inc;
        }
      }
    }
    break;

  case 5: // RDC.
    for (unsigned i = 0; i < rows; i++)
    {
      for (unsigned k = 0; k < depth; k++)
      {
        for (unsigned j = 0; j < cols; j++)
        {
          data[i][j][k] = seq;
          seq += inc;
        }
      }
    }
    break;

  default: // CDR.
    for (unsigned j = 0; j < cols; j++)
    {
      for (unsigned k = 0; k < depth; k++)
      {
        for (unsigned i = 0; i < rows; i++)
        {
          data[i][j][k] = seq;
          seq += inc;
        }
      }
    }
  }
}


// Returns the sum of all cells.
template <typename T>
T& Cube<T>::sum()
{
  T result = 0;

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<depth; k++)
      {
        if (data[i][j][k] < result)
          result += data[i][j][k];
      }
    }
  }

  return result;
}



//**************************************************************************************************




// Various friendship operators and functions.

template<typename T>
Cube<T> operator + (const Cube<T> &u, const Cube<T> &v)
{
  Cube result(rows, cols, depth);

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<depth; k++)
      {
        result(i, j, k) = u(i, j, k) + v(i, j, k);
      }
    }
  }

  return result;
}


template<typename T>
Cube<T> operator - (const Cube<T> &u, const Cube<T> &v)
{
  Cube result(rows, cols, depth);

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<depth; k++)
      {
        result(i, j, k) = u(i, j, k) - v(i, j, k);
      }
    }
  }

  return result;
}


template<typename T>
Cube<T> operator * (const Cube<T> &u, const Cube<T> &v)
{
  Cube result(rows, cols, depth);

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<depth; k++)
      {
        result(i, j, k) = u(i, j, k) * v(i, j, k);
      }
    }
  }

  return result;
}


//  This algorithm is the long division algorithm.
template<typename T>
Cube<T> operator / (const Cube<T> &u, const Cube<T> &v)
{
  Cube result(rows, cols, depth);

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<depth; k++)
      {
        if (v(i, j, k) != 0)
          result(i, j, k) = u(i, j, k) / v(i, j, k);
      }
    }
  }

  return result;
}


template<typename T>
Cube<T> operator % (const Cube<T> &u, const Cube<T> &v)
{
  Cube result(rows, cols, depth);

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<depth; k++)
      {
        if (v(i, j, k) != 0)
          result(i, j, k) = u(i, j, k) % v(i, j, k);
      }
    }
  }

  return result;
}


template<typename T>
Cube<T> operator ^ (const Cube<T> &u, const Cube<T> &v)
{
  Cube result(rows, cols, depth);

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<depth; k++)
      {
        result(i, j, k) = u(i, j, k) ^ v(i, j, k);
      }
    }
  }

  return result;
}


template<typename T>
int operator == (const Cube<T> &u, const Cube<T> &v)
{
  return (u == v);
}


template<typename T>
int operator != (const Cube<T> &u, const Cube<T> &v)
{
  return !(u == v);
}


template<typename T>
int operator < (const Cube<T> &u, const Cube<T> &v)
{
  return (u < v);
}


template<typename T>
int operator <= (const Cube<T> &u, const Cube<T> &v)
{
  return (u <= v);
}


template<typename T>
int operator >(const Cube<T> &u, const Cube<T> &v)
{
  return (u > v);
}


template<typename T>
int operator >= (const Cube<T> &u, const Cube<T> &v)
{
  return (u >= v);
}


// Calculate the absolute value of a number.
template<typename T>
Cube<T> abs(const Cube<T> &v)
{
  Cube result(rows, cols, depth);

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<depth; k++)
      {
        result(i, j, k) = abs(v(i, j, k));
      }
    }
  }

  return result;
}


// Double division function.
template<typename T>
double div(const Cube<T> &u, const Cube<T> &v)
{
  Cube result(rows, cols, depth);

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<depth; k++)
      {
        if (v(i, j, k) != 0)
          result(i, j, k) = u(i, j, k) / v(i, j, k);
      }
    }
  }

  return result;
}


// Returns the size of the cube.
template<typename T>
unsigned size(const Cube<T> &c)
{
  return c.rows * c.cols * c.depth;
}


// Calculate the integer square root of a number
//    based on the formula (a+b)^2 = a^2 + 2ab + b^2
template<typename T>
Cube<T> sqrt(const Cube<T> &v)
{
  Cube result(rows, cols, depth);

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<depth; k++)
      {
        result(i, j, k) = sqrt(v(i, j, k));
      }
    }
  }

  return result;
}


// Raise a number X to a power of degree
template<typename T>
Cube<T> pow(const Cube<T> &X, const Cube<T> &degree)
{
  Cube result(rows, cols, depth);

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k<depth; k++)
      {
        result(i, j, k) = pow(X(i, j, k), degree);
      }
    }
  }

  return result;
}



template<typename T>
//std::ostream & operator << (std::ostream &s, const Cube<T> &m)
std::ostream & operator << (std::ostream &s, Cube<T> &m)
{
  m.print(s);
  return s;
  /*

  unsigned rows = m.getRows();
  unsigned cols = m.getCols();
  unsigned depth = m.getDepth();

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k < depth; k++)
      {
        s << m(i, j, k);

        if ((i*j) < (((rows - 1)*(cols - 1)) - 1))
          s << ", ";

        //if (())
      }
    }
    s << std::endl;
  }

  return s;
  */
}

/*
template<typename T>
std::ostream & operator << (std::ostream &s, const Cube<T> &m)
{
  unsigned rows = m.getRows();
  unsigned cols = m.getCols();
  unsigned depth = m.getDepth();

  for (unsigned i = 0; i<rows; i++)
  {
    for (unsigned j = 0; j<cols; j++)
    {
      for (unsigned k = 0; k < depth; k++)
      {
        s << m(i, j, k);

        if ((i*j) < (((rows - 1)*(cols - 1)) - 1))
          s << ", ";

        //if (())
      }
    }
    s << std::endl;
  }

  return s;
}
*/
/*
template<typename T>
std::ostream & operator << (std::ostream &s, const Cube<T> &m)
{
  m.print(s);
  return s;
}
*/

/*
template<typename T>
std::istream & operator >> (std::istream &s, Cube<T> &v)
{
  std::string temp(10000, ' ');
  
  s >> temp;
  return s;
}
*/


template<typename T>
std::istream & operator >> (std::istream &s, Cube<T> &m)
{
  m.read(s);
  return s;
}





} // end namespace SW.

#endif
</code>