#include "Matrix.h" #include //#include #define _CHECK_INDEX //////////////////////// // Array Constructors // //////////////////////// tcT /*inline*/ Array::Array(int l, int m, int n, T* d/**=0**/) : L(l), M(m), N(n), X(d) { } tcT /*inline*/ Array::Array(int m, int n, T* d/**=0**/) : L(n), M(m), N(n), X(d) { } tcT /*inline*/ Array::Array(int n, T* d/**=0**/) : L(n), M(1), N(n), X(d) { } tcT /*inline*/ Array::Array(const Array& a) : L(a.L), M(a.M), N(a.N), X(a.X) { } tcT /*inline*/ Array::Array() : L(0), M(0), N(0), X(0) { } tcT /*inline*/ /**virtual**/ Array::~Array() { } ///////////////////// // Array selectors // ///////////////////// tcT /*inline*/ int Array::l() const { return L; } tcT /*inline*/ int Array::m() const { return M; } tcT /*inline*/ int Array::n() const { return N; } tcT /*inline*/ T* Array::x() const { return X; } // subArray (without space reallocation) tcT /*inline*/ Array Array::s(int i, int m, int j, int n) const { #ifdef _CHECK_INDEX if ( i<0 ) throw whereOutOfRange( "Array Array::s(int,int,int,int) const","i",i,0,M-1); if ( m<0 || i+m>M ) throw whereOutOfRange( "Array Array::s(int,int,int,int) const","i+m",i+m,0,M-1); if ( j<0 ) throw whereOutOfRange( "Array Array::s(int,int,int,int) const","j",j,0,N-1); if ( n<0 || j+n>N ) throw whereOutOfRange( "Array Array::s(int,int,int,int) const","j+n",j+n,0,N-1); if ( !X ) throw nullArray("Array Array::s(int,int,int,int) const"); #endif return Array(L, m, n, &X[i*L+j]); } tcT /*inline*/ Array Array::s(int i/**=0**/,int m/**=1**/,int j/**=0**/) const { #ifdef _CHECK_INDEX if ( i<0 ) throw whereOutOfRange( "Array Array::s(int=0,int=1,int=0) const","i",i,0,M-1); if ( m<0 || i+m>M ) throw whereOutOfRange( "Array Array::s(int=0,int=1,int=0) const","i+m",i+m,0,M-1); if ( j<0 || j>=N ) throw whereOutOfRange( "Array Array::s(int=0,int=1,int=0) const","j",j,0,N-1); if ( !X ) throw nullArray("Array Array::s(int,int,int,int) const"); #endif return Array(L, m, N-j, &X[i*L+j]); } // subMatrix (with space reallocation) tcT /*inline*/ Matrix Array::_(int i, int m, int j, int n) const { return Matrix(s(i, m, j, n)); } tcT /*inline*/ Matrix Array::_(int i/**=0**/,int m/**=1**/,int j/**=0**/)const { return Matrix(s(i, m, j, N-j)); } // subscripting tcT /*inline*/ T* Array::operator [] (int i) const { #ifdef _CHECK_INDEX if ( i<0 || i>=M ) throw whereOutOfRange( "T* Array::operator[](int) const","i",i,0,M-1); if ( !X ) throw nullArray("T* Array::operator[](int) const"); #endif return &X[i*L]; } tcT /*inline*/ T& Array::operator () (int i) const { #ifdef _CHECK_INDEX if ( i<0 || i>=N ) throw whereOutOfRange( "T& Array::operator()(int) const","i",i,0,N-1); if ( !X ) throw nullArray("T& Array::operator()(int) const"); #endif return X[i]; } tcT /*inline*/ T& Array::operator () (int i, int j) const { #ifdef _CHECK_INDEX if ( i<0 || i>=M ) throw whereOutOfRange( "T& Array::operator()(int,int) const","i",i,0,M-1); if ( j<0 || j>=N ) throw whereOutOfRange( "T& Array::operator()(int,int) const","j",j,0,N-1); if ( !X ) throw nullArray("T& Array::operator()(int,int) const"); #endif return X[i*L+j]; } ////////////////////// // Arrays modifiers // ////////////////////// tcT /*inline*/ Array& Array::resize(const Array& a) { L = a.L; M = a.M; N = a.N; X = a.X; return *this; } tcT /*inline*/ Array& Array::resize(int m, int n, T* d) { L = n; M = m; N = n; X = d; return *this; } tcT /*inline*/ Array& Array::resize(int n, T* d) { L = n; M = 1; N = n; X = d; return *this; } tcT /*inline*/ int& Array::l(int i) { L = i; return L; } tcT /*inline*/ int& Array::m(int j) { M = j; return M; } tcT /*inline*/ int& Array::n(int k) { N = k; return N; } tcT /*inline*/ T*& Array::x(T* u) { X = u; return X; } ////////////////////////// // Assignment of arrays // ////////////////////////// // simple assignment tcT Array& Array::operator = (const T& b) { if ( !N || !M ) return *this; if ( !X ) throw nullArray("Array& Array::operator=(const T&)"); T* s = X; T* t = X + N; do { while (s < t) *s++ = b; s += L - N; } while ((t += L) <= &X[M*L]); return *this; } tcT Array& Array::operator = (const Array& a) { if ( !N || !M ) { if ( N != a.N || M != a.M ) throw whereDimMismatch("Array& Array::operator=(const Array&)",M,N,a.M,a.N); return *this; } if ( !X || !a.X ) throw nullArray("Array& Array::operator=(const Array&)"); if (N == a.N) { if (M == a.M) { // N == a.N && M == a.M T* t = a.X; T* u = X; T* v = X + N; do { while (u < v) *u++ = *t++; t += a.L - a.N; u += L - N; } while ((v += L) <= &X[M*L]); return *this; }; if (a.M == 1) { // N == a.N && M != a.M == 1 T* t; T* u = X; T* v = X + N; do { t = a.X; while (u < v) *u++ = *t++; u += L - N; } while ((v += L) <= &X[M*L]); return *this; } } if (a.N == 1) { if (M == a.M) { // N != a.N == 1 && M == a.M T* t = a.X; T* u = X; T* v = X + N; do { while (u < v) *u++ = *t; t++; u += L - N; } while ((v += L) <= &X[M*L]); return *this; } if (a.M == 1) { // N != a.N == 1 && M != a.M == 1 T* t = a.X; T* u = X; T* v = X + N; do { while (u < v) *u++ = *t; u += L - N; } while ((v += L) <= &X[M*L]); return *this; } } throw whereDimMismatch("Array& Array::operator=(const Array&)",M,N,a.M,a.N); } // add and assign tcT Array& Array::operator += (const T& b) { if ( !N || !M ) return *this; if ( !X ) throw nullArray("Array& Array::operator+=(const T&)"); T* s = X; T* t = X + N; do { while (s < t) *s++ += b; s += L - N; } while ((t += L) <= &X[M*L]); return *this; } tcT Array& Array::operator += (const Array& a) { if ( !N || !M ) { if ( N != a.N || M != a.M ) throw whereDimMismatch( "Array& Array::operator+=(const Array&)",M,N,a.M,a.N); return *this; } if ( !X || !a.X ) throw nullArray("Array& Array::operator+=(const Array&)"); if (N == a.N) { if (M == a.M) { // N == a.N T* t = a.X; T* u = X; T* v = X + N; do { while (u < v) *u++ += *t++; t += a.L - a.N; u += L - N; } while ((v += L) <= &X[M*L]); return *this; }; if (a.M == 1) { // N == a.N && M != a.M T* t; T* u = X; T* v = X + N; do { t = a.X; while (u < v) *u++ += *t++; u += L - N; } while ((v += L) <= &X[M*L]); return *this; }; }; if (a.N == 1) { // N != a.N if (M == a.M) { // N != a.N == 1 T* t = a.X; T* u = X; T* v = X + N; do { while (u < v) *u++ += *t; t++; u += L - N; } while ((v += L) <= &X[M*L]); return *this; }; if (a.M == 1) { // N != a.N == 1 && M != a.M T* t = a.X; T* u = X; T* v = X + N; do { while (u < v) *u++ += *t; u += L - N; } while ((v += L) <= &X[M*L]); return *this; }; }; throw whereDimMismatch( "Array& Array::operator+=(const Array&)",M,N,a.M,a.N); } // substract and assign tcT Array& Array::operator -= (const T& b) { if ( !N || !M ) return *this; if ( !X ) throw nullArray("Array& Array::operator-=(const T&)"); T* s = X; T* t = X + N; do { while (s < t) *s++ -= b; s += L - N; } while ((t += L) <= &X[M*L]); return *this; } tcT Array& Array::operator -= (const Array& a) { if ( !N || !M ) { if ( N != a.N || M != a.M ) throw whereDimMismatch( "Array& Array::operator-=(const Array&)",M,N,a.M,a.N); return *this; } if ( !X || !a.X ) throw nullArray("Array& Array::operator-=(const Array&)"); if (N == a.N) { if (M == a.M) { // N == a.N T* t = a.X; T* u = X; T* v = X + N; do { while (u < v) *u++ -= *t++; t += a.L - a.N; u += L - N; } while ((v += L) <= &X[M*L]); return *this; }; if (a.M == 1) { // N == a.N && M != a.M T* t; T* u = X; T* v = X + N; do { t = a.X; while (u < v) *u++ -= *t++; u += L - N; } while ((v += L) <= &X[M*L]); return *this; }; }; if (a.N == 1) { // N != a.N if (M == a.M) { // N != a.N == 1 T* t = a.X; T* u = X; T* v = X + N; do { while (u < v) *u++ -= *t; t++; u += L - N; } while ((v += L) <= &X[M*L]); return *this; }; if (a.M == 1) { // N != a.N == 1 && M != a.M T* t = a.X; T* u = X; T* v = X + N; do { while (u < v) *u++ -= *t; u += L - N; } while ((v += L) <= &X[M*L]); return *this; }; }; throw whereDimMismatch( "Array& Array::operator-=(const Array&)",M,N,a.M,a.N); } // multiply and assign tcT Array& Array::operator *= (const T& b) { if ( !N || !M ) return *this; if ( !X ) throw nullArray("Array& Array::operator*=(const T&)"); T* s = X; T* t = X + N; do { while (s < t) *s++ *= b; s += L - N; } while ((t += L) <= &X[M*L]); return *this; } tcT Array& Array::operator *= (const Array& a) { if ( !N || !M ) { if ( N != a.N || M != a.M ) throw whereDimMismatch( "Array& Array::operator*=(const Array&)",M,N,a.M,a.N); return *this; } if ( !X || !a.X ) throw nullArray("Array& Array::operator*=(const Array&)"); if (N == a.N) { if (M == a.M) { // N == a.N T* t = a.X; T* u = X; T* v = X + N; do { while (u < v) *u++ *= *t++; t += a.L - a.N; u += L - N; } while ((v += L) <= &X[M*L]); return *this; }; if (a.M == 1) { // N == a.N && M != a.M T* t; T* u = X; T* v = X + N; do { t = a.X; while (u < v) *u++ *= *t++; u += L - N; } while ((v += L) <= &X[M*L]); return *this; }; }; if (a.N == 1) { // N != a.N if (M == a.M) { // N != a.N == 1 T* t = a.X; T* u = X; T* v = X + N; do { while (u < v) *u++ *= *t; t++; u += L - N; } while ((v += L) <= &X[M*L]); return *this; }; if (a.M == 1) { // N != a.N == 1 && M != a.M T* t = a.X; T* u = X; T* v = X + N; do { while (u < v) *u++ *= *t; u += L - N; } while ((v += L) <= &X[M*L]); return *this; }; }; throw whereDimMismatch( "Array& Array::operator*=(const Array&)",M,N,a.M,a.N); } // divide and assign tcT Array& Array::operator /= (const T& b) { if ( !N || !M ) return *this; if ( !X ) throw nullArray("Array& Array::operator/=(const T&)"); T* s = X; T* t = X + N; do { while (s < t) *s++ /= b; s += L - N; } while ((t += L) <= &X[M*L]); return *this; } tcT Array& Array::operator /= (const Array& a) { if ( !N || !M ) { if ( N != a.N || M != a.M ) throw whereDimMismatch( "Array& Array::operator/=(const Array&)",M,N,a.M,a.N); return *this; } if ( !X || !a.X ) throw nullArray("Array& Array::operator/=(const Array&)"); if (N == a.N) { if (M == a.M) { // N == a.N T* t = a.X; T* u = X; T* v = X + N; do { while (u < v) *u++ /= *t++; t += a.L - a.N; u += L - N; } while ((v += L) <= &X[M*L]); return *this; }; if (a.M == 1) { // N == a.N && M != a.M T* t; T* u = X; T* v = X + N; do { t = a.X; while (u < v) *u++ /= *t++; u += L - N; } while ((v += L) <= &X[M*L]); return *this; }; }; if (a.N == 1) { // N != a.N if (M == a.M) { // N != a.N == 1 T* t = a.X; T* u = X; T* v = X + N; do { while (u < v) *u++ /= *t; t++; u += L - N; } while ((v += L) <= &X[M*L]); return *this; }; if (a.M == 1) { // N != a.N == 1 && M != a.M T* t = a.X; T* u = X; T* v = X + N; do { while (u < v) *u++ /= *t; u += L - N; } while ((v += L) <= &X[M*L]); return *this; }; }; throw whereDimMismatch( "Array& Array::operator/=(const Array&)",M,N,a.M,a.N); } // shift right and assign tcT Array& Array::operator >>= (int n) { if ( !N || !M ) return *this; if ( !X ) throw nullArray("Array& Array::operator>>=(int)"); T* t = &X[M*L]; T* u = &X[M*L]; T* v = &X[M*L]; while (X <= (u -= L)) { t -= L - N + n; v -= L - N; while (t > u) *--v = *--t; while (v > u) *--v = *t; }; return *this; } // shift left and assign tcT Array& Array::operator <<= (int n) { if ( !N || !M ) return *this; if ( !X ) throw nullArray("Array& Array::operator<<=(int)"); T* t = X + n; T* u = X; T* v = X + N; do { while (t < v) *u++ = *t++; while (u < v) *u++ = (T)0; t += L - N + n; u += L - N; } while ((v += L) <= &X[M*L]); return *this; } ////////////////////////// // Comparison of arrays // ////////////////////////// // less than tcT boolean Array::operator < (const T& b) const { if ( !X ) throw nullArray("boolean Array::operator<(const T&)"); T* u = X; T* v = X + N; do { while (u < v && *u < b) u++; if (u < v) return false; u += L - N; } while ((v += L) <= &X[M*L]); return true; } tcT /*inline*/ boolean operator < (const T& a, const Array& b) { return b > a; } tcT boolean Array::operator < (const Array& b) const { if ( !X || !b.X ) throw nullArray("boolean Array::operator<(const Array&)"); if (N == b.N) { if (M == b.M) { // N == b.N T* t = X; T* u = b.X; T* v = b.X + b.N; do { while (u < v && *t < *u) { t++; u++; }; if (u < v) return false; t += L - N; u += b.L - b.N; } while ((v += b.L) <= &b.X[b.M*b.L]); return true; }; if (M == 1) { // N == b.N && M != b.M T* u = b.X; T* v = b.X + b.N; do { T* t = X; while (u < v && *t < *u) { t++; u++; }; if (u < v) return false; u += b.L - b.N; } while ((v += b.L) <= &b.X[b.M*b.L]); return true; }; if (b.M == 1) { // N == b.N && M != b.M && M != 1 T* u = X; T* v = X + N; do { T* t = b.X; while (u < v && *u < *t) { t++; u++; }; if (u < v) return false; t += b.L - b.N; u += L - N; } while ((v += L) <= &X[M*L]); return true; }; }; if (N == 1) { // b.N != N if (M == b.M) { // b.N != N == 1 T* t = X; T* u = b.X; T* v = b.X + b.N; do { while (u < v && *t < *u) u++; if (u < v) return false; t++; u += b.L - b.N; } while ((v += b.L) <= &b.X[b.M*b.L]); return true; }; if (M == 1) { // b.N != N == 1 && M != b.M T* u = b.X; T* v = b.X + b.N; do { while (u < v && *X < *u) u++; if (u < v) return false; u += b.L - b.N; } while ((v += b.L) <= &b.X[b.M*b.L]); return true; }; }; if (b.N == 1) { // N != b.N && N != 1 if (M == b.M) { // N != b.N == 1 T* t = b.X; T* u = X; T* v = X + N; do { while (u < v && *u < *t) u++; if (u < v) return false; t++; u += L - N; } while ((v += L) <= &X[M*L]); return true; }; if (b.M == 1) { // N != b.N == 1 && M != b.M T* u = X; T* v = X + N; do { while (u < v && *u < *b.X) u++; if (u < v) return false; u += L - N; } while ((v += L) <= &X[M*L]); return true; }; }; throw whereDimMismatch( "boolean Array::operator<(const Array&)",M,N,b.M,b.N); } // less than or equal tcT boolean Array::operator <= (const T& b) const { if ( !X ) throw nullArray("boolean Array::operator<=(const T&)"); T* u = X; T* v = X + N; do { while (u < v && *u <= b) u++; if (u < v) return false; u += L - N; } while ((v += L) <= &X[M*L]); return true; } tcT /*inline*/ boolean operator <= (const T& a, const Array& b) { return b >= a; } tcT boolean Array::operator <= (const Array& b) const { if ( !X || !b.X ) throw nullArray("boolean Array::operator<=(const Array&)"); if (N == b.N) { if (M == b.M) { // N == b.N T* t = X; T* u = b.X; T* v = b.X + b.N; do { while (u < v && *t <= *u) { t++; u++; }; if (u < v) return false; t += L - N; u += b.L - b.N; } while ((v += b.L) <= &b.X[b.M*b.L]); return true; }; if (M == 1) { // N == b.N && M != b.M T* u = b.X; T* v = b.X + b.N; do { T* t = X; while (u < v && *t <= *u) { t++; u++; }; if (u < v) return false; u += b.L - b.N; } while ((v += b.L) <= &b.X[b.M*b.L]); return true; }; if (b.M == 1) { // N == b.N && M != b.M && M != 1 T* u = X; T* v = X + N; do { T* t = b.X; while (u < v && *u <= *t) { t++; u++; }; if (u < v) return false; t += b.L - b.N; u += L - N; } while ((v += L) <= &X[M*L]); return true; }; }; if (N == 1) { // b.N != N if (M == b.M) { // b.N != N == 1 T* t = X; T* u = b.X; T* v = b.X + b.N; do { while (u < v && *t <= *u) u++; if (u < v) return false; t++; u += b.L - b.N; } while ((v += b.L) <= &b.X[b.M*b.L]); return true; }; if (M == 1) { // b.N != N == 1 && M != b.M T* u = b.X; T* v = b.X + b.N; do { while (u < v && *X <= *u) u++; if (u < v) return false; u += b.L - b.N; } while ((v += b.L) <= &b.X[b.M*b.L]); return true; }; }; if (b.N == 1) { // N != b.N && N != 1 if (M == b.M) { // N != b.N == 1 T* t = b.X; T* u = X; T* v = X + N; do { while (u < v && *u <= *t) u++; if (u < v) return false; t++; u += L - N; } while ((v += L) <= &X[M*L]); return true; }; if (b.M == 1) { // N != b.N == 1 && M != b.M T* u = X; T* v = X + N; do { while (u < v && *u <= *b.X) u++; if (u < v) return false; u += L - N; } while ((v += L) <= &X[M*L]); return true; }; }; throw whereDimMismatch( "boolean Array::operator<=(const Array&)",M,N,b.M,b.N); } // greater than tcT boolean Array::operator > (const T& b) const { if ( !X ) throw nullArray("boolean Array::operator>(const T&)"); T* u = X; T* v = X + N; do { while (u < v && *u > b) u++; if (u < v) return false; u += L - N; } while ((v += L) <= &X[M*L]); return true; } tcT /*inline*/ boolean operator > (const T& a, const Array& b) { return b < a; } tcT boolean Array::operator > (const Array& b) const { if ( !X || !b.X ) throw nullArray("boolean Array::operator>(const Array&)"); if (N == b.N) { if (M == b.M) { // N == b.N T* t = X; T* u = b.X; T* v = b.X + b.N; do { while (u < v && *t > *u) { t++; u++; }; if (u < v) return false; t += L - N; u += b.L - b.N; } while ((v += b.L) <= &b.X[b.M*b.L]); return true; }; if (M == 1) { // N == b.N && M != b.M T* u = b.X; T* v = b.X + b.N; do { T* t = X; while (u < v && *t > *u) { t++; u++; }; if (u < v) return false; u += b.L - b.N; } while ((v += b.L) <= &b.X[b.M*b.L]); return true; }; if (b.M == 1) { // N == b.N && M != b.M && M != 1 T* u = X; T* v = X + N; do { T* t = b.X; while (u < v && *u > *t) { t++; u++; }; if (u < v) return false; t += b.L - b.N; u += L - N; } while ((v += L) <= &X[M*L]); return true; }; }; if (N == 1) { // b.N != N if (M == b.M) { // b.N != N == 1 T* t = X; T* u = b.X; T* v = b.X + b.N; do { while (u < v && *t > *u) u++; if (u < v) return false; t++; u += b.L - b.N; } while ((v += b.L) <= &b.X[b.M*b.L]); return true; }; if (M == 1) { // b.N != N == 1 && M != b.M T* u = b.X; T* v = b.X + b.N; do { while (u < v && *X > *u) u++; if (u < v) return false; u += b.L - b.N; } while ((v += b.L) <= &b.X[b.M*b.L]); return true; }; }; if (b.N == 1) { // N != b.N && N != 1 if (M == b.M) { // N != b.N == 1 T* t = b.X; T* u = X; T* v = X + N; do { while (u < v && *u > *t) u++; if (u < v) return false; t++; u += L - N; } while ((v += L) <= &X[M*L]); return true; }; if (b.M == 1) { // N != b.N == 1 && M != b.M T* u = X; T* v = X + N; do { while (u < v && *u > *b.X) u++; if (u < v) return false; u += L - N; } while ((v += L) <= &X[M*L]); return true; }; }; throw whereDimMismatch( "boolean Array::operator>(const Array&)",M,N,b.M,b.N); } // greater than or equal tcT boolean Array::operator >= (const T& b) const { if ( !X ) throw nullArray("boolean Array::operator>=(const T&)"); T* u = X; T* v = X + N; do { while (u < v && *u >= b) u++; if (u < v) return false; u += L - N; } while ((v += L) <= &X[M*L]); return true; } tcT /*inline*/ boolean operator >= (const T& a, const Array& b) { return b <= a; } tcT boolean Array::operator >= (const Array& b) const { if ( !X || !b.X ) throw nullArray("boolean Array::operator>=(const Array&)"); if (N == b.N) { if (M == b.M) { // N == b.N T* t = X; T* u = b.X; T* v = b.X + b.N; do { while (u < v && *t >= *u) { t++; u++; }; if (u < v) return false; t += L - N; u += b.L - b.N; } while ((v += b.L) <= &b.X[b.M*b.L]); return true; }; if (M == 1) { // N == b.N && M != b.M T* u = b.X; T* v = b.X + b.N; do { T* t = X; while (u < v && *t >= *u) { t++; u++; }; if (u < v) return false; u += b.L - b.N; } while ((v += b.L) <= &b.X[b.M*b.L]); return true; }; if (b.M == 1) { // N == b.N && M != b.M && M != 1 T* u = X; T* v = X + N; do { T* t = b.X; while (u < v && *u >= *t) { t++; u++; }; if (u < v) return false; t += b.L - b.N; u += L - N; } while ((v += L) <= &X[M*L]); return true; }; }; if (N == 1) { // b.N != N if (M == b.M) { // b.N != N == 1 T* t = X; T* u = b.X; T* v = b.X + b.N; do { while (u < v && *t >= *u) u++; if (u < v) return false; t++; u += b.L - b.N; } while ((v += b.L) <= &b.X[b.M*b.L]); return true; }; if (M == 1) { // b.N != N == 1 && M != b.M T* u = b.X; T* v = b.X + b.N; do { while (u < v && *X >= *u) u++; if (u < v) return false; u += b.L - b.N; } while ((v += b.L) <= &b.X[b.M*b.L]); return true; }; }; if (b.N == 1) { // N != b.N && N != 1 if (M == b.M) { // N != b.N == 1 T* t = b.X; T* u = X; T* v = X + N; do { while (u < v && *u >= *t) u++; if (u < v) return false; t++; u += L - N; } while ((v += L) <= &X[M*L]); return true; }; if (b.M == 1) { // N != b.N == 1 && M != b.M T* u = X; T* v = X + N; do { while (u < v && *u >= *b.X) u++; if (u < v) return false; u += L - N; } while ((v += L) <= &X[M*L]); return true; }; }; throw whereDimMismatch( "boolean Array::operator>=(const Array&)",M,N,b.M,b.N); } // equal tcT boolean Array::operator == (const T& b) const { if ( !X ) throw nullArray("boolean Array::operator==(const T&)"); T* u = X; T* v = X + N; do { while (u < v && *u == b) u++; if (u < v) return false; u += L - N; } while ((v += L) <= &X[M*L]); return true; } tcT /*inline*/ boolean operator == (const T& a, const Array& b) { return (b == a); } tcT boolean Array::operator == (const Array& b) const { if ( !X || !b.X ) throw nullArray("boolean Array::operator==(const Array&)"); if (N == b.N) { if (M == b.M) { // N == b.N T* t = X; T* u = b.X; T* v = b.X + b.N; do { while (u < v && *t == *u) { t++; u++; }; if (u < v) return false; t += L - N; u += b.L - b.N; } while ((v += b.L) <= &b.X[b.M*b.L]); return true; }; if (M == 1) { // N == b.N && M != b.M T* u = b.X; T* v = b.X + b.N; do { T* t = X; while (u < v && *t == *u) { t++; u++; }; if (u < v) return false; u += b.L - b.N; } while ((v += b.L) <= &b.X[b.M*b.L]); return true; }; if (b.M == 1) { // N == b.N && M != b.M && M != 1 T* u = X; T* v = X + N; do { T* t = b.X; while (u < v && *u == *t) { t++; u++; }; if (u < v) return false; t += b.L - b.N; u += L - N; } while ((v += L) <= &X[M*L]); return true; }; }; if (N == 1) { // b.N != N if (M == b.M) { // b.N != N == 1 T* t = X; T* u = b.X; T* v = b.X + b.N; do { while (u < v && *t == *u) u++; if (u < v) return false; t++; u += b.L - b.N; } while ((v += b.L) <= &b.X[b.M*b.L]); return true; }; if (M == 1) { // b.N != N == 1 && M != b.M T* u = b.X; T* v = b.X + b.N; do { while (u < v && *X == *u) u++; if (u < v) return false; u += b.L - b.N; } while ((v += b.L) <= &b.X[b.M*b.L]); return true; }; }; if (b.N == 1) { // N != b.N && N != 1 if (M == b.M) { // N != b.N == 1 T* t = b.X; T* u = X; T* v = X + N; do { while (u < v && *u == *t) u++; if (u < v) return false; t++; u += L - N; } while ((v += L) <= &X[M*L]); return true; }; if (b.M == 1) { // N != b.N == 1 && M != b.M T* u = X; T* v = X + N; do { while (u < v && *u == *b.X) u++; if (u < v) return false; u += L - N; } while ((v += L) <= &X[M*L]); return true; }; }; throw whereDimMismatch( "boolean Array::operator==(const Array&)",M,N,b.M,b.N); } // not equal tcT /*inline*/ boolean Array::operator != (const T& b) const { return !(*this == b); } tcT /*inline*/ boolean operator != (const T& a, const Array& b) { return !(b == a); } tcT /*inline*/ boolean Array::operator != (const Array& b) const { return !(*this == b); } /////////////////// // I/O of arrays // /////////////////// //#ifndef _array_output_vars //#define _array_output_vars /*extern*/ int Array_number = 40; /*extern*/ int Array_coefWidth = 6; /*extern*/ int Array_coefPrecision = -1; /*extern*/ boolean Array_charAsChar = false; /*extern*/ boolean Array_longLines = false; //#endif // output tcT ostream& operator << (ostream& s, const Array& b) { if ( !b.X ) throw nullArray( "ostream& operator<<(ostream& s, const Array&)"); T* t = b.X; T* u = b.X; T* v = b.X + b.N; if ( sizeof(T)==1 && !Array_charAsChar ) // T is a char { do { do { t += (0 < Array_number)? Array_number: b.N; t = (t < v)? t: v; while (u < t) if (!(s << coefWidth(Array_coefWidth) << int(*u) && ((++u < t) ? s << "" : s << "\n"))) return s; } while (t < v); t = u += b.L - b.N; } while ((v += b.L) <= &b.X[b.M*b.L]); } else { do { do { t += (0 < Array_number)? Array_number: b.N; t = (t < v)? t: v; while (u < t) { if ( Array_coefPrecision>=0 || *u!=int(*u) ) s << setw(Array_coefWidth) << *u; else s << setw(Array_coefWidth-1+Array_coefPrecision) << int(*u) << setw(-Array_coefPrecision+1) << setfill(' ') << ""; if ( !( s && ((++u < t) ? s << "" : s << "\n"))) return s; } } while (t < v); t = u += b.L - b.N; } while ((v += b.L) <= &b.X[b.M*b.L]); } return s; } // output manipulators ostream& coefWidth(ostream& s, int w) { Array_coefWidth = w; s.width (Array_coefWidth); return s; } OMANIP (int) coefWidth(int w) { return OMANIP (int) (coefWidth, w); } ostream& coefPrecision(ostream& s, int p) { Array_coefPrecision = p; s.unsetf(ios::scientific); s.setf(ios::fixed); if ( p ) s.setf(ios::showpoint); else s.unsetf(ios::showpoint); s.precision(p); return s; } OMANIP (int) coefPrecision(int p) { return OMANIP (int) (coefPrecision, p); } ostream& charAsCharacter (ostream& s) { Array_charAsChar=true; return s; } ostream& charAsNumber (ostream& s) { Array_charAsChar=false; return s; } ostream& coefPerLine (ostream& s, int nbC) { Array_number=nbC; return s; } OMANIP (int) coefPerLine(int nbC) { return OMANIP (int) (coefPerLine, nbC); } // input tcT istream& operator >> (istream& s, const Array& b) { if ( !b.X ) throw nullArray( "istream& operator>>(istream& s, const Array&)"); T* u = b.X; T* v = b.X + b.N; if ( sizeof(T) == 1 && !Array_charAsChar ) // T is a char { int r; do { while (u < v && s >> r) { *u = (T)r; u++; } if (u < v) return s; u += b.L - b.N; } while ((v += b.L) <= &b.X[b.M*b.L]); } else if ( sizeof(T) == 1 && Array_charAsChar ) // T is a char { char r; do { while (u < v && s.get(r)) if (r>31) {*u = r; u++;} if (u < v) return s; u += b.L - b.N; } while ((v += b.L) <= &b.X[b.M*b.L]); } else { do { while (u < v && s >> *u) u++; if (u < v) return s; u += b.L - b.N; } while ((v += b.L) <= &b.X[b.M*b.L]); } return s; } istream& charAsCharacter (istream& s) { Array_charAsChar=true; return s; } istream& charAsNumber (istream& s) { Array_charAsChar=false; return s; } ///////////////////////// // Operators on Arrays // ///////////////////////// // unary plus tcT /*inline*/ Matrix Array::operator + () const { return Matrix(*this); } // unary minus tcT Matrix Array::operator - () const { if ( !X ) throw nullArray("Matrix Array::operator-() const"); T* newX = new T [M*N]; T* t = newX; T* u = X; T* v = X; do while (u < v + N) *t++ = -(*u++); while ((u = v += L) < &X[M*L]); return Matrix(Array(M, N, newX)); } // add (plus) tcT Matrix Array::operator + (const T& b) const { if ( !X ) throw nullArray("Matrix Array::operator+(const T&) const"); T* newx = new T [M*N]; T* t = newx; T* u = X; T* v = X + N; do { while (u < v) *t++ = *u++ + b; u += L - N; } while ((v += L) <= &X[M*L]); return Matrix(Array(M, N, newx)); } tcT Matrix operator + (const T& a, const Array& b) { if ( !b.X ) throw nullArray("Matrix operator+(const T&, const Array&)"); T* newx = new T [b.M*b.N]; T* t = newx; T* u = b.X; T* v = b.X + b.N; do { while (u < v) *t++ = a + *u++; u += b.L - b.N; } while ((v += b.L) <= &b.X[b.M*b.L]); return Matrix(Array(b.M, b.N, newx)); } tcT Matrix Array::operator + (const Array& b) const { if ( !X || !b.X ) throw nullArray("Matrix Array::operator+(const Array&) const"); if (N == b.N) { if (M == b.M) { // N == b.N T* newx = new T [b.M*b.N]; T* s = newx; T* t = newx; T* u = X; T* v = b.X; while ((t += b.N) <= &newx[b.M*b.N]) { while (s < t) *s++ = *u++ + *v++; u += L - N; v += b.L - b.N; }; return Matrix(Array(b.M, b.N, newx)); }; if (M == 1) { // N == b.N && M != b.M T* newx = new T [b.M*b.N]; T* s = newx; T* t = newx; T* v = b.X; while ((t += b.N) <= &newx[b.M*b.N]) { T* u = X; while (s < t) *s++ = *u++ + *v++; v += b.L - b.N; }; return Matrix(Array(b.M, b.N, newx)); }; if (b.M == 1) { // N == b.N && M != b.M && M != 1 T* newx = new T [M*N]; T* s = newx; T* t = newx; T* u = X; while ((t += N) <= &newx[M*N]) { T* v = b.X; while (s < t) *s++ = *u++ + *v++; u += L - N; }; return Matrix(Array(M, N, newx)); }; }; if (N == 1) { // b.N != N if (M == b.M) { // b.N != N == 1 T* newx = new T [b.M*b.N]; T* s = newx; T* t = newx; T* u = X; T* v = b.X; while ((t += b.N) <= &newx[b.M*b.N]) { while (s < t) *s++ = *u + *v++; u++; v += b.L - b.N; }; return Matrix(Array(b.M, b.N, newx)); }; if (M == 1) { // b.N != N == 1 && M != b.M T* newx = new T [b.M*b.N]; T* s = newx; T* t = newx; T* u = X; T* v = b.X; while ((t += b.N) <= &newx[b.M*b.N]) { while (s < t) *s++ = *u + *v++; v += b.L - b.N; }; return Matrix(Array(b.M, b.N, newx)); }; }; if (b.N == 1) { // N != b.N && N != 1 if (M == b.M) { // N != b.N == 1 T* newx = new T [M*N]; T* s = newx; T* t = newx; T* u = X; T* v = b.X; while ((t += N) <= &newx[M*N]) { while (s < t) *s++ = *u++ + *v; u += L - N; v++; }; return Matrix(Array(M, N, newx)); }; if (b.M == 1) { // N != b.N == 1 && M != b.M T* newx = new T [M*N]; T* s = newx; T* t = newx; T* u = X; T* v = b.X; while ((t += N) <= &newx[M*N]) { while (s < t) *s++ = *u++ + *v; u += L - N; }; return Matrix(Array(M, N, newx)); }; }; throw whereDimMismatch( "Matrix Array::operator+(const Array&) const", M,N,b.M,b.N); } // subtract (minus) tcT Matrix Array::operator - (const T& b) const { if ( !X ) throw nullArray("Matrix Array::operator-(const T&) const"); T* newx = new T [M*N]; T* t = newx; T* u = X; T* v = X + N; do { while (u < v) *t++ = *u++ - b; u += L - N; } while ((v += L) <= &X[M*L]); return Matrix(Array(M, N, newx)); } tcT Matrix operator - (const T& a, const Array& b) { if ( !b.X ) throw nullArray("Matrix operator-(const T&, const Array&)"); T* newx = new T [b.M*b.N]; T* t = newx; T* u = b.X; T* v = b.X + b.N; do { while (u < v) *t++ = a - *u++; u += b.L - b.N; } while ((v += b.L) <= &b.X[b.M*b.L]); return Matrix(Array(b.M, b.N, newx)); } tcT Matrix Array::operator - (const Array& b) const { if ( !X || !b.X ) throw nullArray("Matrix Array::operator-(const Array&) const"); if (N == b.N) { if (M == b.M) { // N == b.N T* newx = new T [b.M*b.N]; T* s = newx; T* t = newx; T* u = X; T* v = b.X; while ((t += b.N) <= &newx[b.M*b.N]) { while (s < t) *s++ = *u++ - *v++; u += L - N; v += b.L - b.N; }; return Matrix(Array(b.M, b.N, newx)); }; if (M == 1) { // N == b.N && M != b.M T* newx = new T [b.M*b.N]; T* s = newx; T* t = newx; T* v = b.X; while ((t += b.N) <= &newx[b.M*b.N]) { T* u = X; while (s < t) *s++ = *u++ - *v++; v += b.L - b.N; }; return Matrix(Array(b.M, b.N, newx)); }; if (b.M == 1) { // N == b.N && M != b.M && M != 1 T* newx = new T [M*N]; T* s = newx; T* t = newx; T* u = X; while ((t += N) <= &newx[M*N]) { T* v = b.X; while (s < t) *s++ = *u++ - *v++; u += L - N; }; return Matrix(Array(M, N, newx)); }; }; if (N == 1) { // b.N != N if (M == b.M) { // b.N != N == 1 T* newx = new T [b.M*b.N]; T* s = newx; T* t = newx; T* u = X; T* v = b.X; while ((t += b.N) <= &newx[b.M*b.N]) { while (s < t) *s++ = *u - *v++; u++; v += b.L - b.N; }; return Matrix(Array(b.M, b.N, newx)); }; if (M == 1) { // b.N != N == 1 && M != b.M T* newx = new T [b.M*b.N]; T* s = newx; T* t = newx; T* u = X; T* v = b.X; while ((t += b.N) <= &newx[b.M*b.N]) { while (s < t) *s++ = *u - *v++; v += b.L - b.N; }; return Matrix(Array(b.M, b.N, newx)); }; }; if (b.N == 1) { // N != b.N && N != 1 if (M == b.M) { // N != b.N == 1 T* newx = new T [M*N]; T* s = newx; T* t = newx; T* u = X; T* v = b.X; while ((t += N) <= &newx[M*N]) { while (s < t) *s++ = *u++ - *v; u += L - N; v++; }; return Matrix(Array(M, N, newx)); }; if (b.M == 1) { // N != b.N == 1 && M != b.M T* newx = new T [M*N]; T* s = newx; T* t = newx; T* u = X; T* v = b.X; while ((t += N) <= &newx[M*N]) { while (s < t) *s++ = *u++ - *v; u += L - N; }; return Matrix(Array(M, N, newx)); }; }; throw whereDimMismatch( "Matrix Array::operator-(const Array&) const", M,N,b.M,b.N); } // multiply tcT Matrix Array::operator * (const T& b) const { if ( !X ) throw nullArray("Matrix Array::operator*(const T&) const"); T* newx = new T [M*N]; T* t = newx; T* u = X; T* v = X + N; do { while (u < v) *t++ = *u++ * b; u += L - N; } while ((v += L) <= &X[M*L]); return Matrix(Array(M, N, newx)); } tcT Matrix operator * (const T& a, const Array& b) { if ( !b.X ) throw nullArray("Matrix operator*(const T&, const Array&)"); T* newx = new T [b.M*b.N]; T* t = newx; T* u = b.X; T* v = b.X + b.N; do { while (u < v) *t++ = a * *u++; u += b.L - b.N; } while ((v += b.L) <= &b.X[b.M*b.L]); return Matrix(Array(b.M, b.N, newx)); } tcT Matrix Array::operator * (const Array& b) const { if ( !X || !b.X ) throw nullArray("Matrix Array::operator*(const Array&) const"); if (N == b.N) { if (M == b.M) { // N == b.N T* newx = new T [b.M*b.N]; T* s = newx; T* t = newx; T* u = X; T* v = b.X; while ((t += b.N) <= &newx[b.M*b.N]) { while (s < t) *s++ = *u++ * *v++; u += L - N; v += b.L - b.N; }; return Matrix(Array(b.M, b.N, newx)); }; if (M == 1) { // N == b.N && M != b.M T* newx = new T [b.M*b.N]; T* s = newx; T* t = newx; T* v = b.X; while ((t += b.N) <= &newx[b.M*b.N]) { T* u = X; while (s < t) *s++ = *u++ * *v++; v += b.L - b.N; }; return Matrix(Array(b.M, b.N, newx)); }; if (b.M == 1) { // N == b.N && M != b.M && M != 1 T* newx = new T [M*N]; T* s = newx; T* t = newx; T* u = X; while ((t += N) <= &newx[M*N]) { T* v = b.X; while (s < t) *s++ = *u++ * *v++; u += L - N; }; return Matrix(Array(M, N, newx)); }; }; if (N == 1) { // b.N != N if (M == b.M) { // b.N != N == 1 T* newx = new T [b.M*b.N]; T* s = newx; T* t = newx; T* u = X; T* v = b.X; while ((t += b.N) <= &newx[b.M*b.N]) { while (s < t) *s++ = *u * *v++; u++; v += b.L - b.N; }; return Matrix(Array(b.M, b.N, newx)); }; if (M == 1) { // b.N != N == 1 && M != b.M T* newx = new T [b.M*b.N]; T* s = newx; T* t = newx; T* u = X; T* v = b.X; while ((t += b.N) <= &newx[b.M*b.N]) { while (s < t) *s++ = *u * *v++; v += b.L - b.N; }; return Matrix(Array(b.M, b.N, newx)); }; }; if (b.N == 1) { // N != b.N && N != 1 if (M == b.M) { // N != b.N == 1 T* newx = new T [M*N]; T* s = newx; T* t = newx; T* u = X; T* v = b.X; while ((t += N) <= &newx[M*N]) { while (s < t) *s++ = *u++ * *v; u += L - N; v++; }; return Matrix(Array(M, N, newx)); }; if (b.M == 1) { // N != b.N == 1 && M != b.M T* newx = new T [M*N]; T* s = newx; T* t = newx; T* u = X; T* v = b.X; while ((t += N) <= &newx[M*N]) { while (s < t) *s++ = *u++ * *v; u += L - N; }; return Matrix(Array(M, N, newx)); }; }; throw whereDimMismatch( "Matrix Array::operator*(const Array&) const", M,N,b.M,b.N); } // divide tcT Matrix Array::operator / (const T& b) const { if ( !X ) throw nullArray("Matrix Array::operator/(const T&) const"); T* newx = new T [M*N]; T* t = newx; T* u = X; T* v = X + N; do { while (u < v) *t++ = *u++ / b; u += L - N; } while ((v += L) <= &X[M*L]); return Matrix(Array(M, N, newx)); } tcT Matrix operator / (const T& a, const Array& b) { if ( !b.X ) throw nullArray("Matrix operator/(const T&, const Array&)"); T* newx = new T [b.M*b.N]; T* t = newx; T* u = b.X; T* v = b.X + b.N; do { while (u < v) *t++ = a / *u++; u += b.L - b.N; } while ((v += b.L) <= &b.X[b.M*b.L]); return Matrix(Array(b.M, b.N, newx)); } tcT Matrix Array::operator / (const Array& b) const { if ( !X || !b.X ) throw nullArray("Matrix Array::operator/(const Array&) const"); if (N == b.N) { if (M == b.M) { // N == b.N T* newx = new T [b.M*b.N]; T* s = newx; T* t = newx; T* u = X; T* v = b.X; while ((t += b.N) <= &newx[b.M*b.N]) { while (s < t) *s++ = *u++ / *v++; u += L - N; v += b.L - b.N; }; return Matrix(Array(b.M, b.N, newx)); }; if (M == 1) { // N == b.N && M != b.M T* newx = new T [b.M*b.N]; T* s = newx; T* t = newx; T* v = b.X; while ((t += b.N) <= &newx[b.M*b.N]) { T* u = X; while (s < t) *s++ = *u++ / *v++; v += b.L - b.N; }; return Matrix(Array(b.M, b.N, newx)); }; if (b.M == 1) { // N == b.N && M != b.M && M != 1 T* newx = new T [M*N]; T* s = newx; T* t = newx; T* u = X; while ((t += N) <= &newx[M*N]) { T* v = b.X; while (s < t) *s++ = *u++ / *v++; u += L - N; }; return Matrix(Array(M, N, newx)); }; }; if (N == 1) { // b.N != N if (M == b.M) { // b.N != N == 1 T* newx = new T [b.M*b.N]; T* s = newx; T* t = newx; T* u = X; T* v = b.X; while ((t += b.N) <= &newx[b.M*b.N]) { while (s < t) *s++ = *u / *v++; u++; v += b.L - b.N; }; return Matrix(Array(b.M, b.N, newx)); }; if (M == 1) { // b.N != N == 1 && M != b.M T* newx = new T [b.M*b.N]; T* s = newx; T* t = newx; T* u = X; T* v = b.X; while ((t += b.N) <= &newx[b.M*b.N]) { while (s < t) *s++ = *u / *v++; v += b.L - b.N; }; return Matrix(Array(b.M, b.N, newx)); }; }; if (b.N == 1) { // N != b.N && N != 1 if (M == b.M) { // N != b.N == 1 T* newx = new T [M*N]; T* s = newx; T* t = newx; T* u = X; T* v = b.X; while ((t += N) <= &newx[M*N]) { while (s < t) *s++ = *u++ / *v; u += L - N; v++; }; return Matrix(Array(M, N, newx)); }; if (b.M == 1) { // N != b.N == 1 && M != b.M T* newx = new T [M*N]; T* s = newx; T* t = newx; T* u = X; T* v = b.X; while ((t += N) <= &newx[M*N]) { while (s < t) *s++ = *u++ / *v; u += L - N; }; return Matrix(Array(M, N, newx)); }; }; throw whereDimMismatch( "Matrix Array::operator/(const Array&) const", M,N,b.M,b.N); } // Matrix multiplication(inner product) tcT Matrix Array::operator % (const Array& b) const { if ( !X || !b.X ) throw nullArray("Matrix Array::operator%(const Array&) const"); if (N != b.N) throw whereDimMismatch( "Matrix Array::operator%(const Array&) const", M,N,b.M,b.N); T* newx = new T [M*b.M]; T* s = newx; T* t = X; T* u; T* v; do { v = b.X; while (v < &b.X[b.M*b.L]) { u = t; *s = *u++ * *v++; while (u < t + N) *s += *u++ * *v++; ++s; v += b.L - b.N; }; } while((t += L) < &X[M*L]); return Matrix(Array(M, b.M, newx)); } // Matrix multiplication (outer product) tcT Matrix Array::operator & (const Array& b) const { if ( !X || !b.X ) throw nullArray("Matrix Array::operator&(const Array&) const"); if (M != b.M) throw whereDimMismatch( "Matrix Array::operator&(const Array&) const", M,N,b.M,b.N); T* newx = new T [N*b.N]; T* p = newx; T* s = X; T* t; T* u; T* v; do { t = b.X; do { u = s; v = t; *p = *u * *v; while ((u += L), (v += b.L) < &b.X[b.M*b.L]) *p += *u * *v; p++; } while (++t < &b.X[b.N]); } while (++s < &X[N]); return Matrix(Array(N, b.N, newx)); } // shift left tcT Matrix Array::operator << (int n) const { if ( !X ) throw nullArray("Matrix Array::operator<<(int) const"); T* newx = new T [M*this->N]; T* t = X + n; T* u = newx; T* v = newx; while ((v += this->N) <= &newx[M*this->N]) { while (u < v - n) *u++ = *t++; while (u < v) *u++ = (T)0; t += L - this->N + n; }; return Matrix(Array(M, this->N, newx)); } // shift right tcT Matrix Array::operator >> (int n) const { if ( !X ) throw nullArray("Matrix Array::operator>>(int) const"); T* newx = new T [M*this->N]; T* t = &X[M*L]; T* u = &newx[M*this->N]; T* v = &newx[M*this->N]; while (newx <= (u -= this->N)) { t -= L - this->N + n; while (v > u + n) *--v = *--t; while (v > u) *--v = *t; }; return Matrix(Array(M, this->N, newx)); } // stack tcT Matrix Array::operator ^ (const Array& b) const { if ( !X || !b.X) throw nullArray("Matrix Array::operator^(const Array&) const"); if (N != b.N) throw whereDimMismatch( "Matrix Array::operator^(const Array&) const", M,N,b.M,b.N); T* newx = new T [(M+b.M)*N]; T* t = newx; T* u = X; T* v = X + N; do { while (u < v) *t++ = *u++; u += L - N; } while((v += L) <= &X[M*L]); u = b.X; v = b.X + b.N; do { while (u < v) *t++ = *u++; u += b.L - b.N; } while((v += b.L) <= &b.X[b.M*b.L]); return Matrix(Array(M+b.M, N, newx)); } // adjoin tcT Matrix Array