#ifndef __SUNINTERVAL_MATRIX_H
#define __SUNINTERVAL_MATRIX_H

#include <valarray>
#include <iostream.h>
#include "suninterval_vector.h"

namespace SUNW_interval {

template< typename T > class mref;
template< typename T > class slice_iter;
template< typename T > class cslice_iter;

//
// nmatrix
//
template< typename T >
class nmatrix {
    struct mrep;
    mrep *rep;
public:
    nmatrix(size_t dd1, size_t dd2);
    nmatrix(const T& value, size_t dd1, size_t dd2);
    nmatrix(const T* pointer, size_t dd1, size_t dd2);
    nmatrix(const nmatrix& m);
    nmatrix(const std::valarray<T>& v, size_t dd1, size_t dd2);
    ~nmatrix();
        
    nmatrix<T>& operator= (const T& val);
    nmatrix<T>& operator= (const nmatrix<T>& m);
    mref<T> operator() (size_t i, size_t j);
    T operator() (size_t i, size_t j) const;
    slice_iter<T> operator() (size_t i);
    cslice_iter<T> operator() (size_t i) const;
    slice_iter<T> operator[] (size_t i);
    cslice_iter<T> operator[] (size_t i) const;

    slice_iter<T> row(size_t i);
    cslice_iter<T> row(size_t i) const;
    slice_iter<T> column(size_t i);
    cslice_iter<T> column(size_t i) const;

    nmatrix<T> operator+() const;
    nmatrix<T> operator-() const;
    nmatrix<T>& operator+= (const nmatrix<T>& m);
    nmatrix<T>& operator-= (const nmatrix<T>& m);
    nmatrix<T>& operator*= (const nmatrix<T>& m);
    nmatrix<T>& operator/= (const nmatrix<T>& m);
    nmatrix<T>& operator+= (const T& val);
    nmatrix<T>& operator-= (const T& val);
    nmatrix<T>& operator*= (const T& val);
    nmatrix<T>& operator/= (const T& val);

    size_t size() const { return rep->d1 * rep->d2; }
    size_t dim1() const { return rep->d1; }
    size_t dim2() const { return rep->d2; }

    T& ref(size_t i);
    T& direct_ref(size_t i) { return rep->v[i]; };
    std::valarray<T>& array() const { return rep->v; }
};

//
// reference to nmatrix element
//
template< typename T >
class mref {
    nmatrix<T>& m;
    size_t i;
public:
    mref(nmatrix<T>& mm, size_t ii) : m(mm), i(ii) { }
    operator T () const { return m.array()[i]; };
    void operator= (const T& x) {
        m.ref(i) = x;
    }
    void operator= (const mref& r) {
        m.ref(i) = (T)r;
    }
    void operator+= (const T& x) {
        m.ref(i) += x;
    }
    void operator-= (const T& x) {
        m.ref(i) -= x;
    }
    void operator*= (const T& x) {
        m.ref(i) *= x;
    }
    void operator/= (const T& x) {
        m.ref(i) /= x;
    }
};

template< typename T > 
bool operator== (const slice_iter<T>& p, const slice_iter<T>& q);
template< typename T > 
bool operator!= (const slice_iter<T>& p, const slice_iter<T>& q);
template< typename T > 
bool operator<  (const slice_iter<T>& p, const slice_iter<T>& q);
template< typename T > 
bool operator== (const cslice_iter<T>& p, const cslice_iter<T>& q);
template< typename T > 
bool operator!= (const cslice_iter<T>& p, const cslice_iter<T>& q);
template< typename T > 
bool operator<  (const cslice_iter<T>& p, const cslice_iter<T>& q);
	
//
// slice of array
//
template< typename T >
class slice_iter {
    nmatrix<T>& m;
    std::slice s;
    size_t curr;
    size_t index(size_t i) { return s.start() + i*s.stride(); }
public:
    slice_iter(nmatrix<T>& mm, std::slice ss) : 
        m(mm), s(ss), curr(0) { }

    slice_iter end() const {
        slice_iter t = *this;
        t.curr = s.size();
        return t;
    }
    slice_iter& operator++ () { curr++; return *this; }
    slice_iter operator++ (int) { slice_iter t = *this; curr++; return t; }

    mref<T> operator[] (size_t i) { return mref<T>(m, index(i)); }
    mref<T> operator() (size_t i) { return mref<T>(m, index(i)); }
    mref<T> operator* () { return mref<T>(m, index(curr)); }

    friend bool operator== <>(const slice_iter<T>& p, const slice_iter<T>& q);
    friend bool operator!= <>(const slice_iter<T>& p, const slice_iter<T>& q);
    friend bool operator< <>(const slice_iter<T>& p, const slice_iter<T>& q);
};

//
// constant slice of array
//
template< typename T >
class cslice_iter {
    std::valarray<T> *v;
    std::slice s;
    size_t curr;
    const T& ref(size_t i) const { 
        return (*v)[s.start() + i*s.stride()]; 
    }
public:
    cslice_iter(std::valarray<T>* vv, std::slice ss) : 
        v(vv), s(ss), curr(0) { }

    cslice_iter end() const {
        cslice_iter t = *this;
        t.curr = s.size();
        return t;
    }
    cslice_iter& operator++ () { curr++; return *this; }
    cslice_iter operator++ (int) { cslice_iter t = *this; curr++; return t; }

    const T& operator[] (size_t i) const { return ref(i); }
    const T& operator() (size_t i) const { return ref(i); }
    const T& operator* () const { return ref(curr); }
    
    friend bool operator== <>(const cslice_iter<T>& p, const cslice_iter<T>& q);
    friend bool operator!= <>(const cslice_iter<T>& p, const cslice_iter<T>& q);
    friend bool operator< <>(const cslice_iter<T>& p, const cslice_iter<T>& q);
};

template< typename T >
inline bool operator== (const slice_iter<T>& p, const slice_iter<T>& q)
{
    return p.curr == q.curr && p.s.stride() == q.s.stride() &&
       p.s.start() == q.s.start(); 
}

template< typename T >
inline bool operator!= (const slice_iter<T>& p, const slice_iter<T>& q)
{
    return !(p==q);
}

template< typename T >
inline bool operator< (const slice_iter<T>& p, const slice_iter<T>& q)
{
    return p.curr < q.curr && p.s.stride() == q.s.stride() &&
       p.s.start() == q.s.start(); 
}

template< typename T >
inline bool operator== (const cslice_iter<T>& p, const cslice_iter<T>& q)
{
    return p.curr == q.curr && p.s.stride() == q.s.stride() &&
       p.s.start() == q.s.start(); 
}

template< typename T >
inline bool operator!= (const cslice_iter<T>& p, const cslice_iter<T>& q)
{
    return !(p==q);
}

template< typename T >
inline bool operator< (const cslice_iter<T>& p, const cslice_iter<T>& q)
{
    return p.curr < q.curr && p.s.stride() == q.s.stride() &&
       p.s.start() == q.s.start(); 
}

//
// nmatrix representation class
//
template< typename T >
struct nmatrix<T>::mrep {
    std::valarray<T> v;
    size_t d1, d2;
    unsigned n;
    mrep(size_t dd1, size_t dd2);
    mrep(const T& value, size_t dd1, size_t dd2);
    mrep(const T* pointer, size_t dd1, size_t dd2);
    mrep(const std::valarray<T>& vec, size_t dd1, size_t dd2);

    mrep *get_own_copy() {
        if (n == 1) 
            return this;
        n--;
        return new mrep(v, d1, d2);
    }
};

template< typename T >
inline nmatrix<T>::mrep::mrep(size_t dd1, size_t dd2) : v(dd1*dd2) 
{
    n = 1;
    d1 = dd1;
    d2 = dd2;
}

template< typename T >
inline nmatrix<T>::mrep::mrep(const T& value, 
                              size_t dd1, size_t dd2) : v(value, dd1*dd2)
{
    n = 1;
    d1 = dd1;
    d2 = dd2;
}

template< typename T >
inline nmatrix<T>::mrep::mrep(const T* pointer, 
                              size_t dd1, size_t dd2) : v(pointer, dd1*dd2)
{
    n = 1;
    d1 = dd1;
    d2 = dd2;
}

template< typename T >
inline nmatrix<T>::mrep::mrep(const std::valarray<T>& vec, 
                              size_t dd1, size_t dd2) : v(vec)
{
    n = 1;
    d1 = dd1;
    d2 = dd2;
}

template< typename T >
inline T& nmatrix<T>::ref(size_t i)
{
    rep = rep->get_own_copy();
    return rep->v[i];
}
    
template< typename T >
inline nmatrix<T>::nmatrix(size_t dd1, size_t dd2) 
{
    rep = new mrep(dd1, dd2);
}

template< typename T >
inline nmatrix<T>::nmatrix(const T& value, size_t dd1, size_t dd2) 
{
    rep = new mrep(value, dd1, dd2);
}
        
template< typename T >
inline nmatrix<T>::nmatrix(const T* pointer, size_t dd1, size_t dd2) 
{
    rep = new mrep(pointer, dd1, dd2);
}
        
template< typename T >
inline nmatrix<T>::nmatrix(const nmatrix& m) 
{
    mrep *t = m.rep;
    t->n++;
    rep = t;
}
    
template< typename T >
inline nmatrix<T>::nmatrix
(const std::valarray<T>& vec, size_t dd1, size_t dd2)
{
    rep = new mrep(vec, dd1, dd2);
}

template< typename T >
inline nmatrix<T>::~nmatrix() 
{
    if (--rep->n == 0) 
        delete rep;
}

template< typename T >
inline nmatrix<T>& 
nmatrix<T>::operator= (const T& val) 
{
    rep = rep->get_own_copy();
    rep->v = val;
    return *this;
}

template< typename T >
inline nmatrix<T>& 
nmatrix<T>::operator= (const nmatrix<T>& m) 
{
    mrep *t = m.rep;
    t->n++;
    if (--rep->n == 0)
        delete rep;
    rep = t;
    return *this;
}

template< typename T >
inline slice_iter<T> nmatrix<T>::row(size_t i) 
{
    return slice_iter<T>(*this, std::slice(i, rep->d1, rep->d2));
}

template< typename T >
inline cslice_iter<T> nmatrix<T>::row(size_t i) const 
{
    return cslice_iter<T>(&rep->v, std::slice(i, rep->d1, rep->d2));
}

template< typename T >
inline slice_iter<T> nmatrix<T>::column(size_t i) 
{
    return slice_iter<T>(*this, std::slice(i*rep->d2, rep->d2, 1));
}

template< typename T >
inline cslice_iter<T> nmatrix<T>::column(size_t i) const 
{
    return cslice_iter<T>(&rep->v, std::slice(i*rep->d2, rep->d2, 1));
}

template< typename T >
inline mref<T> nmatrix<T>::operator() (size_t i, size_t j) 
{ 
    return row(i)[j]; 
}

template< typename T >
inline T nmatrix<T>::operator() (size_t i, size_t j) const 
{ 
    return row(i)[j]; 
}

template< typename T >
inline slice_iter<T> nmatrix<T>::operator() (size_t i) 
{ 
    return row(i); 
}

template< typename T >
inline cslice_iter<T> nmatrix<T>::operator() (size_t i) const 
{ 
    return row(i); 
}
    
template< typename T >
inline slice_iter<T> nmatrix<T>::operator[] (size_t i) 
{ 
    return row(i); 
}

template< typename T >
inline cslice_iter<T> nmatrix<T>::operator[] (size_t i) const 
{ 
    return row(i); 
}

template< typename T >
inline nmatrix<T> nmatrix<T>::operator+() const 
{ 
    return *this; 
}

template< typename T >
inline nmatrix<T> nmatrix<T>::operator-() const 
{ 
    return nmatrix<T>(-rep->v, dim1(), dim2());
}
    
template< typename T >
inline nmatrix<T>& 
nmatrix<T>::operator+= (const nmatrix<T>& m) 
{
    rep = rep->get_own_copy();
    rep->v += m.array();
    return *this;
}

template< typename T >
inline nmatrix<T>& 
nmatrix<T>::operator-= (const nmatrix<T>& m) 
{
    rep = rep->get_own_copy();
    rep->v -= m.array();
    return *this;
}

template< typename T >
inline nmatrix<T>& 
nmatrix<T>::operator*= (const nmatrix<T>& m) 
{
    rep = rep->get_own_copy();
    rep->v *= m.array();
    return *this;
}

template< typename T >
inline nmatrix<T>& 
nmatrix<T>::operator/= (const nmatrix<T>& m) 
{
    rep = rep->get_own_copy();
    rep->v /= m.array();
    return *this;
}

template< typename T >
inline nmatrix<T>& 
nmatrix<T>::operator+= (const T& val) 
{
    rep = rep->get_own_copy();
    rep->v += val;
    return *this;
}

template< typename T >
inline nmatrix<T>& 
nmatrix<T>::operator-= (const T& val) 
{ 
    rep = rep->get_own_copy();
    rep->v -= val;
    return *this;
}

template< typename T >
inline nmatrix<T>& 
nmatrix<T>::operator*= (const T& val) 
{ 
    rep = rep->get_own_copy();
    rep->v *= val;
    return *this;
}

template< typename T >
inline nmatrix<T>& 
nmatrix<T>::operator/= (const T& val) 
{ 
    rep = rep->get_own_copy();
    rep->v /= val;
    return *this;
}

// nonmember operators and functions with nmatrix

template< typename T >
inline nmatrix<T> operator+ (const nmatrix<T>& a, const nmatrix<T>& b)
{
    if (a.dim1() != b.dim1() || a.dim2() != b.dim2())
        cerr << "incompatible operands shapes" << endl;
    return nmatrix<T>(a.array() + b.array(), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T> operator- (const nmatrix<T>& a, const nmatrix<T>& b)
{
    if (a.dim1() != b.dim1() || a.dim2() != b.dim2())
        cerr << "incompatible operands shapes" << endl;
    return nmatrix<T>(a.array() - b.array(), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T> operator* (const nmatrix<T>& a, const nmatrix<T>& b)
{
    if (a.dim1() != b.dim1() || a.dim2() != b.dim2())
        cerr << "incompatible operands shapes" << endl;
    return nmatrix<T>(a.array() * b.array(), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T> operator/ (const nmatrix<T>& a, const nmatrix<T>& b)
{
    if (a.dim1() != b.dim1() || a.dim2() != b.dim2())
        cerr << "incompatible operands shapes" << endl;
    return nmatrix<T>(a.array() / b.array(), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T> operator+ (const nmatrix<T>& a, const T& b)
{
    return nmatrix<T>(a.array() + b, a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T> operator- (const nmatrix<T>& a, const T& b)
{
    return nmatrix<T>(a.array() - b, a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T> operator* (const nmatrix<T>& a, const T& b)
{
    return nmatrix<T>(a.array() * b, a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T> operator/ (const nmatrix<T>& a, const T& b)
{
    return nmatrix<T>(a.array() / b, a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T> operator+ (const T& a, const nmatrix<T>& b)
{
    return nmatrix<T>(a + b.array(), b.dim1(), b.dim2());
}

template< typename T >
inline nmatrix<T> operator- (const T& a, const nmatrix<T>& b)
{
    return nmatrix<T>(a - b.array(), b.dim1(), b.dim2());
}

template< typename T >
inline nmatrix<T> operator* (const T& a, const nmatrix<T>& b)
{
    return nmatrix<T>(a * b.array(), b.dim1(), b.dim2());
}

template< typename T >
inline nmatrix<T> operator/ (const T& a, const nmatrix<T>& b)
{
    return nmatrix<T>(a / b.array(), b.dim1(), b.dim2());
}

template< typename T >
inline nmatrix<bool>
operator== (const nmatrix<T>& a, const nmatrix<T>& b)
{
    if (a.dim1() != b.dim1() || a.dim2() != b.dim2())
        cerr << "incompatible operands shapes" << endl;
    return nmatrix<bool>(a.array() == b.array(), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<bool>
operator== (const nmatrix<T>& a, const T& b)
{
    return nmatrix<bool>(a.array() == b, a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<bool>
operator== (const T& a, const nmatrix<T>& b)
{
    return nmatrix<bool>(a == b.array(), b.dim1(), b.dim2());
}

template< typename T >
inline nmatrix<bool>
operator!= (const nmatrix<T>& a, const nmatrix<T>& b)
{
    if (a.dim1() != b.dim1() || a.dim2() != b.dim2())
        cerr << "incompatible operands shapes" << endl;
    return nmatrix<bool>(a.array() != b.array(), a.dim1(), a.dim2());
}
  
template< typename T >
inline nmatrix<bool>
operator!= (const nmatrix<T>& a, const T& b)
{
    return nmatrix<bool>(a.array() != b, a.dim1(), a.dim2());
}
  
template< typename T >
inline nmatrix<bool>
operator!= (const T& a, const nmatrix<T>& b)
{
    return nmatrix<bool>(a != b.array(), b.dim1(), b.dim2());
}
  
/* Extraction Functions: infimim, supremum, midpoint, width, magnitude,
mignitude */

template< typename T >
inline nmatrix<T> inf(const nmatrix< interval<T> >& a)
{
    return nmatrix<T>(inf(a.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T> sup(const nmatrix< interval<T> >& a)
{
    return nmatrix<T>(sup(a.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T> mid(const nmatrix< interval<T> >& a)
{
    return nmatrix<T>(mid(a.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T> wid(const nmatrix< interval<T> >& a)
{
    return nmatrix<T>(wid(a.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T> mag(const nmatrix< interval<T> >& a)
{
    return nmatrix<T>(mag(a.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T> mig(const nmatrix< interval<T> >& a)
{
    return nmatrix<T>(mig(a.array()), a.dim1(), a.dim2());
}

/* Intersection Division */

template< typename T >
inline nmatrix<T>
divix(const nmatrix<T>& a, const nmatrix<T>& b, const nmatrix<T>& c)
{
    return nmatrix<T>(divix(a.array(), b.array(), c.array()), a.dim1(), a.dim2());
}

/* Hull Functions: regular hull, intersection */

template< typename T >
inline nmatrix<T>
interval_hull(const nmatrix<T>& a, const nmatrix<T>& b)
{
    return nmatrix<T>(interval_hull(a.array(), b.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T>
intersect(const nmatrix<T>& a, const nmatrix<T>& b)
{
    return nmatrix<T>(intersect(a.array(), b.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T>
interval_hull(const nmatrix<T>& a, const T& b)
{
    return nmatrix<T>(interval_hull(a.array(), b), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T>
intersect(const nmatrix<T>& a, const T& b)
{
    return nmatrix<T>(intersect(a.array(), b), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T>
interval_hull(const T& a, const nmatrix<T>& b)
{
    return nmatrix<T>(interval_hull(a, b.array()), b.dim1(), b.dim2());
}

template< typename T >
inline nmatrix<T>
intersect(const T& a, const nmatrix<T>& b)
{
    return nmatrix<T>(intersect(a, b.array()), b.dim1(), b.dim2());
}

/* Query Functions: is empty, ndigits */

template< typename T >
inline nmatrix<interval_bool> is_empty(const nmatrix<T>& a)
{
    return nmatrix<interval_bool>(is_empty(a.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<int> ndigits(const nmatrix<T>& a)
{
    return nmatrix<int>(ndigits(a.array()), a.dim1(), a.dim2());
}

/* Membership Functions */

template< typename T1, typename T2 >
inline nmatrix<interval_bool>
in(const nmatrix<T1>& a, const nmatrix<T2>& b)
{
    return nmatrix<interval_bool>(in(a.array(), b.array()), a.dim1(), a.dim2());
}

/* Comparison Functions: disjoint, in interior, subset, superset,
proper subset, proper superset */

template< typename T >
inline nmatrix<interval_bool>
disjoint(const nmatrix<T>& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(disjoint(a.array(), b.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
in_interior(const nmatrix<T>& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(in_interior(a.array(), b.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
subset(const nmatrix<T>& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(subset(a.array(), b.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
superset(const nmatrix<T>& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(superset(a.array(), b.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
proper_subset(const nmatrix<T>& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(proper_subset(a.array(), b.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
proper_superset(const nmatrix<T>& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(proper_superset(a.array(), b.array()), a.dim1(), a.dim2());
}

/* Comparison Functions: disjoint, in interior, subset, superset,
proper subset, proper superset */

template< typename T >
inline nmatrix<interval_bool>
disjoint(const nmatrix<T>& a, const T& b)
{
    return nmatrix<interval_bool>(disjoint(a.array(), b), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
in_interior(const nmatrix<T>& a, const T& b)
{
    return nmatrix<interval_bool>(in_interior(a.array(), b), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
subset(const nmatrix<T>& a, const T& b)
{
    return nmatrix<interval_bool>(subset(a.array(), b), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
superset(const nmatrix<T>& a, const T& b)
{
    return nmatrix<interval_bool>(superset(a.array(), b), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
proper_subset(const nmatrix<T>& a, const T& b)
{
    return nmatrix<interval_bool>(proper_subset(a.array(), b), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
proper_superset(const nmatrix<T>& a, const T& b)
{
    return nmatrix<interval_bool>(proper_superset(a.array(), b), a.dim1(), a.dim2());
}

/* Comparison Functions: disjoint, in interior, subset, superset,
proper subset, proper superset */

template< typename T >
inline nmatrix<interval_bool>
disjoint(const T& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(disjoint(a, b.array()), b.dim1(), b.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
in_interior(const T& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(in_interior(a, b.array()), b.dim1(), b.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
subset(const T& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(subset(a, b.array()), b.dim1(), b.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
superset(const T& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(superset(a, b.array()), b.dim1(), b.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
proper_subset(const T& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(proper_subset(a, b.array()), b.dim1(), b.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
proper_superset(const T& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(proper_superset(a, b.array()), b.dim1(), b.dim2());
}

/* Set Relations: equal, not equal, less than, less than or equal,
greater than, greater than or equal */

template< typename T >
inline nmatrix<interval_bool>
seq(const nmatrix<T>& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(seq(a.array(), b.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
sne(const nmatrix<T>& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(sne(a.array(), b.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
slt(const nmatrix<T>& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(slt(a.array(), b.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
sle(const nmatrix<T>& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(sle(a.array(), b.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
sgt(const nmatrix<T>& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(sgt(a.array(), b.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
sge(const nmatrix<T>& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(sge(a.array(), b.array()), a.dim1(), a.dim2());
}

/* Certainly Relations: equal, not equal, less than, less than or
equal, greater than, greater than or equal */

template< typename T >
inline nmatrix<interval_bool>
ceq(const nmatrix<T>& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(ceq(a.array(), b.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
cne(const nmatrix<T>& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(cne(a.array(), b.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
clt(const nmatrix<T>& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(clt(a.array(), b.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
cle(const nmatrix<T>& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(cle(a.array(), b.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
cgt(const nmatrix<T>& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(cgt(a.array(), b.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
cge(const nmatrix<T>& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(cge(a.array(), b.array()), a.dim1(), a.dim2());
}

/* Possibly Relations: equal, not equal, less than, less than or equal,
greater than, greater than or equal */

template< typename T >
inline nmatrix<interval_bool>
peq(const nmatrix<T>& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(peq(a.array(), b.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
pne(const nmatrix<T>& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(pne(a.array(), b.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
plt(const nmatrix<T>& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(plt(a.array(), b.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
ple(const nmatrix<T>& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(ple(a.array(), b.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
pgt(const nmatrix<T>& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(pgt(a.array(), b.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
pge(const nmatrix<T>& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(pge(a.array(), b.array()), a.dim1(), a.dim2());
}

/* Set Relations: equal, not equal, less than, less than or equal,
greater than, greater than or equal */

template< typename T >
inline nmatrix<interval_bool>
seq(const nmatrix<T>& a, const T& b)
{
    return nmatrix<interval_bool>(seq(a.array(), b), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
sne(const nmatrix<T>& a, const T& b)
{
    return nmatrix<interval_bool>(sne(a.array(), b), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
slt(const nmatrix<T>& a, const T& b)
{
    return nmatrix<interval_bool>(slt(a.array(), b), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
sle(const nmatrix<T>& a, const T& b)
{
    return nmatrix<interval_bool>(sle(a.array(), b), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
sgt(const nmatrix<T>& a, const T& b)
{
    return nmatrix<interval_bool>(sgt(a.array(), b), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
sge(const nmatrix<T>& a, const T& b)
{
    return nmatrix<interval_bool>(sge(a.array(), b), a.dim1(), a.dim2());
}

/* Certainly Relations: equal, not equal, less than, less than or
equal, greater than, greater than or equal */

template< typename T >
inline nmatrix<interval_bool>
ceq(const nmatrix<T>& a, const T& b)
{
    return nmatrix<interval_bool>(ceq(a.array(), b), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
cne(const nmatrix<T>& a, const T& b)
{
    return nmatrix<interval_bool>(cne(a.array(), b), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
clt(const nmatrix<T>& a, const T& b)
{
    return nmatrix<interval_bool>(clt(a.array(), b), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
cle(const nmatrix<T>& a, const T& b)
{
    return nmatrix<interval_bool>(cle(a.array(), b), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
cgt(const nmatrix<T>& a, const T& b)
{
    return nmatrix<interval_bool>(cgt(a.array(), b), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
cge(const nmatrix<T>& a, const T& b)
{
    return nmatrix<interval_bool>(cge(a.array(), b), a.dim1(), a.dim2());
}

/* Possibly Relations: equal, not equal, less than, less than or equal,
greater than, greater than or equal */

template< typename T >
inline nmatrix<interval_bool>
peq(const nmatrix<T>& a, const T& b)
{
    return nmatrix<interval_bool>(peq(a.array(), b), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
pne(const nmatrix<T>& a, const T& b)
{
    return nmatrix<interval_bool>(pne(a.array(), b), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
plt(const nmatrix<T>& a, const T& b)
{
    return nmatrix<interval_bool>(plt(a.array(), b), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
ple(const nmatrix<T>& a, const T& b)
{
    return nmatrix<interval_bool>(ple(a.array(), b), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
pgt(const nmatrix<T>& a, const T& b)
{
    return nmatrix<interval_bool>(pgt(a.array(), b), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
pge(const nmatrix<T>& a, const T& b)
{
    return nmatrix<interval_bool>(pge(a.array(), b), a.dim1(), a.dim2());
}

/* Set Relations: equal, not equal, less than, less than or equal,
greater than, greater than or equal */

template< typename T >
inline nmatrix<interval_bool>
seq(const T& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(seq(a, b.array()), b.dim1(), b.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
sne(const T& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(sne(a, b.array()), b.dim1(), b.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
slt(const T& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(slt(a, b.array()), b.dim1(), b.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
sle(const T& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(sle(a, b.array()), b.dim1(), b.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
sgt(const T& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(sgt(a, b.array()), b.dim1(), b.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
sge(const T& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(sge(a, b.array()), b.dim1(), b.dim2());
}

/* Certainly Relations: equal, not equal, less than, less than or
equal, greater than, greater than or equal */

template< typename T >
inline nmatrix<interval_bool>
ceq(const T& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(ceq(a, b.array()), b.dim1(), b.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
cne(const T& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(cne(a, b.array()), b.dim1(), b.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
clt(const T& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(clt(a, b.array()), b.dim1(), b.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
cle(const T& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(cle(a, b.array()), b.dim1(), b.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
cgt(const T& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(cgt(a, b.array()), b.dim1(), b.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
cge(const T& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(cge(a, b.array()), b.dim1(), b.dim2());
}

/* Possibly Relations: equal, not equal, less than, less than or equal,
greater than, greater than or equal */

template< typename T >
inline nmatrix<interval_bool>
peq(const T& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(peq(a, b.array()), b.dim1(), b.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
pne(const T& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(pne(a, b.array()), b.dim1(), b.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
plt(const T& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(plt(a, b.array()), b.dim1(), b.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
ple(const T& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(ple(a, b.array()), b.dim1(), b.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
pgt(const T& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(pgt(a, b.array()), b.dim1(), b.dim2());
}

template< typename T >
inline nmatrix<interval_bool>
pge(const T& a, const nmatrix<T>& b)
{
    return nmatrix<interval_bool>(pge(a, b.array()), b.dim1(), b.dim2());
}

/* Arithmetic Functions: absolute value, min, max, ceil, floor, fmod */
  
template< typename T >
inline nmatrix<T> fabs(const nmatrix<T>& a)
{
    return nmatrix<T>(fabs(a.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T>
minimum(const nmatrix<T>& a, const nmatrix<T>& b)
{
    return nmatrix<T>(minimum(a.array(), b.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T>
maximum(const nmatrix<T>& a, const nmatrix<T>& b)
{
    return nmatrix<T>(maximum(a.array(), b.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T>
minimum(const nmatrix<T>& a, const T& b)
{
    return nmatrix<T>(minimum(a.array(), b), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T>
maximum(const nmatrix<T>& a, const T& b)
{
    return nmatrix<T>(maximum(a.array(), b), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T>
minimum(const T& a, const nmatrix<T>& b)
{
    return nmatrix<T>(minimum(a, b.array()), b.dim1(), b.dim2());
}

template< typename T >
inline nmatrix<T>
maximum(const T& a, const nmatrix<T>& b)
{
    return nmatrix<T>(maximum(a, b.array()), b.dim1(), b.dim2());
}

template< typename T >
inline nmatrix<T> ceil(const nmatrix< interval<T> >& a)
{
    return nmatrix<T>(ceil(a.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T> floor(const nmatrix< interval<T> >& a)
{
    return nmatrix<T>(floor(a.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T>
fmod(const nmatrix<T>& a, const nmatrix<T>& b)
{
    return nmatrix<T>(fmod(a.array(), b.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T>
fmod(const nmatrix<T>& a, const T& b)
{
    return nmatrix<T>(fmod(a.array(), b), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T>
fmod(const T& a, const nmatrix<T>& b)
{
    return nmatrix<T>(fmod(a, b.array()), b.dim1(), b.dim2());
}

/* Power Functions: square, square root, cube root, natural power,
natural log, log base 10, integral power, floating power */

template< typename T >
inline nmatrix<T> sqr(const nmatrix<T>& a)
{
    return nmatrix<T>(sqr(a.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T> sqrt(const nmatrix<T>& a)
{
    return nmatrix<T>(sqrt(a.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T> cbrt(const nmatrix<T>& a)
{
    return nmatrix<T>(cbrt(a.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T> exp(const nmatrix<T>& a)
{
    return nmatrix<T>(exp(a.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T> log(const nmatrix<T>& a)
{
    return nmatrix<T>(log(a.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T> log10(const nmatrix<T>& a)
{
    return nmatrix<T>(log10(a.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T> 
pow(const nmatrix<T>& a, const nmatrix<T>& b)
{
    return nmatrix<T>(pow(a.array(), b.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T> 
pow(const nmatrix<T>& a, const T& b)
{
    return nmatrix<T>(pow(a.array(), b), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T> 
pow(const T& a, const nmatrix<T>& b)
{
    return nmatrix<T>(pow(a, b.array()), b.dim1(), b.dim2());
}

template< typename T >
inline nmatrix<T> 
pow(const nmatrix<T>& a, const nmatrix<int>& b)
{
    return nmatrix<T>(pow(a.array(), b.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T> 
pow(const nmatrix<T>& a, const int& b)
{
    return nmatrix<T>(pow(a.array(), b), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T> 
pow(const T& a, const nmatrix<int>& b)
{
    return nmatrix<T>(pow(a, b.array()), b.dim1(), b.dim2());
}

/* Trigonometric Functions */

template< typename T >
inline nmatrix<T> acos(const nmatrix<T>& a)
{
    return nmatrix<T>(acos(a.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T> asin(const nmatrix<T>& a)
{
    return nmatrix<T>(asin(a.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T> atan(const nmatrix<T>& a)
{
    return nmatrix<T>(atan(a.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T> cos(const nmatrix<T>& a)
{
    return nmatrix<T>(cos(a.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T> cosh(const nmatrix<T>& a)
{
    return nmatrix<T>(cosh(a.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T> sin(const nmatrix<T>& a)
{
    return nmatrix<T>(sin(a.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T> sinh(const nmatrix<T>& a)
{
    return nmatrix<T>(sinh(a.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T> tan(const nmatrix<T>& a)
{
    return nmatrix<T>(tan(a.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T> tanh(const nmatrix<T>& a)
{
    return nmatrix<T>(tanh(a.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T> 
atan2(const nmatrix<T>& a, const nmatrix<T>& b)
{
    return nmatrix<T>(atan2(a.array(), b.array()), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T> 
atan2(const nmatrix<T>& a, const T& b)
{
    return nmatrix<T>(atan2(a.array(), b), a.dim1(), a.dim2());
}

template< typename T >
inline nmatrix<T> 
atan2(const T& a, const nmatrix<T>& b)
{
    return nmatrix<T>(atan2(a, b.array()), b.dim1(), b.dim2());
}

template< typename T >
nmatrix<T> matmul(const nmatrix<T>& a, const nmatrix<T>& b);

template< typename T >
nvector<T> matmul(const nmatrix<T>& m, const nvector<T>& a);

template< typename T >
nvector<T> matmul(const nvector<T>& a, const nmatrix<T>& m);

template< typename T >
nmatrix<T> transpose(const nmatrix<T>& m);
    
template< typename T >
T 
dot_product(const cslice_iter<T>& a, const std::valarray<T>& b);

template< typename T >
T 
dot_product(const cslice_iter<T>& a, const cslice_iter<T>& b);

template< typename T >
ostream& operator<<(ostream& out, const nmatrix<T>& m);

}

#endif
