#ifndef _NTPM_LINALG_H
#define _NTPM_LINALG_H

#include <complex>
#include <numeric>
#include <string>
#include <functional>



namespace NTPM {
    using namespace std;

// The STL minus functional is restrictive in its argument types.
template <class T1, class T2, class To>
struct general_minus : public std::binary_function<T1,T2,To> {
    To operator()(const T1 &t1, const T2 &t2) { return t1-t2; }
};

// The STL multiplies functional is restrictive in its argument types.
template <class _Arg1, class _Arg2, class _Result>
struct general_multiplies : public std::binary_function<_Arg1,_Arg2,_Result> {
  _Result operator()(const _Arg1& __x, const _Arg2& __y) const {
      return _Result(__x * __y); }
};



template <class V1it, class V2T>
inline void cross(V1it ait, V2T b1, V2T b2, V2T b3)
{
    V1it atit = ait;
    iterator_traits<V1it>::value_type a1=*atit, a2=*++atit, a3=*++atit;
    *ait   = a2*b3 - a3*b2;
    *++ait = a3*b1 - a1*b3;
    *++ait = a1*b2 - a2*b1;
}

template <class V1it, class V2it>
inline void cross(V1it ait, const V2it bit)
{
    iterator_traits<V2it>::value_type b1=*bit, b2=*++bit, b3=*++bit;
    cross(ait,b1,b2,b3);
}


template <class T_M_IT, class T_ACC>
typename iterator_traits<T_M_IT>::value_type
inline mean(T_M_IT it, T_M_IT end,
     const T_ACC &t = iterator_traits<T_M_IT>::value_type())
{
    int nelms;
    T_ACC sum = 0;
    for(nelms=0; it != end; ++it, ++nelms)
	sum += *it;
    return typename iterator_traits<T_M_IT>::value_type(sum / nelms);
}


template <class T_M, class T_ACC>
typename T_M::value_type
inline mean(const T_M &m, const T_ACC &t = T_M::value_type())
{
    return mean(m.begin(), m.end(), T_ACC());
}

template <class T_M>
typename T_M::value_type
inline mean(const T_M &m)
{
    return mean(m, typename T_M::value_type());
}

template <class T_M_IT, class T_ACC>
inline T_ACC
product(T_M_IT it, T_M_IT end,
       const T_ACC &t = iterator_traits<T_M_IT>::value_type())
{
    T_ACC acc = accumulate(it, end, T_ACC(1),
                           general_multiplies<T_ACC,
			                    iterator_traits<T_M_IT>::value_type,
					    T_ACC>());  // from valarray
    return acc;
}

template <class T_M, class T_ACC>
inline T_ACC
product(const T_M &m, const T_ACC &t = T_M::value_type())
{
    return product(m.begin(), m.end(), T_ACC());
}

template <class T_M>
inline typename T_M::value_type
product(const T_M &m)
{
    return product(m, typename T_M::value_type());
}


template <class T_M_IT, class T_ACC>
inline T_ACC
sum(T_M_IT it, T_M_IT end,
     const T_ACC &t = iterator_traits<T_M_IT>::value_type())
{
    T_ACC acc = accumulate(it, end, T_ACC(0));  // from valarray
    return acc;
}

template <class T_M, class T_ACC>
inline T_ACC
sum(const T_M &m, const T_ACC &t = T_M::value_type())
{
    return sum(m.begin(), m.end(), T_ACC());
}

template <class T_M>
inline typename T_M::value_type
sum(const T_M &m)
{
    return sum(m, typename T_M::value_type());
}


template <class T_M_IT, class T_ACC>
inline T_ACC
norm(T_M_IT begin, T_M_IT end,
	  const T_ACC &t = typename T_M_IT::value_type())
{
    T_ACC pow = 0;
    for(; begin != end; ++begin)
	pow += (*begin)*(*begin);
    return pow;
}

template <class CONTAINER, class T_ACC>
inline T_ACC
norm(const CONTAINER &c, const T_ACC &t = typename CONTAINER::value_type())
{
    return norm(c.begin(), c.end(), T_ACC());
}



template <class V1, class V2, class T_ACC>
inline T_ACC
dot(const V1 &v1, const V2 &v2, const T_ACC &)
{
    // use stl inner product in <numeric>.
    return inner_product(v1.begin(), v1.end(), v2.begin(), T_ACC());
}

template <class V1, class V2>
inline typename V1::value_type
dot(const V1 &v1, const V2 &v2)
{
    // use stl inner product in <numeric>.
    return dot(v1,v2,typename V1::value_type());
}


template <class V1, class V2>
void copy(V1 &v1, const V2 &v2)
{
    copy(v2.begin(), v2.end(), v1.begin());
}



}

#endif // _NTPM_LINALG_H