#ifndef Foundation_Math_pe_INCLUDED
#define Foundation_Math_pe_INCLUDED

namespace Tn {

namespace math {

template <int I, typename T>
struct pe {

    static void zero(T* v) {
        pe<I-1,T>::zero(v);
        v[I] = T(0);
    }

    static void neg(T* d, const T* s) {
        pe<I-1,T>::neg(d, s);
        d[I] = -s[I];
    }

    static void copy(T* d, const T* s) {
        pe<I-1,T>::copy(d, s);
        d[I] = s[I];
    }

    template <typename S>
    static void cast(T* d, const S* s) {
        pe<I-1,T>::cast(d, s);
        d[I] = T(s[I]);
    }

    static T dot(const T* a, const T* b) {
        return a[I] * b[I] + pe<I-1,T>::dot(a,b);
    }

    static void add(T* r, const T* a, const T* b) {
        pe<I-1,T>::add(r,a,b);
        r[I] = a[I] + b[I];
    }

    static void add(T* r, const T* a, const T* b, const T* c) {
        pe<I-1,T>::add(r,a,b,c);
        r[I] = a[I] + b[I] + c[I];
    }

    static void sub(T* r, const T* a, const T* b) {
        pe<I-1,T>::sub(r,a,b);
        r[I] = a[I] - b[I];
    }

    static void scale(T* r, const T* a, T s) {
        pe<I-1,T>::scale(r,a,s);
        r[I] = a[I] * s;
    }

    static void scale(T* r, T s, const T* a) {
        pe<I-1,T>::scale(r,a,s);
        r[I] = s * a[I];
    }

    static bool isZero(T const * v) {
        return pe<I-1,T>::isZero(v) && (v[I]==T(0));
    }
    
    static bool isEqual(T const * a, T const * b) {
        return pe<I-1,T>::isEqual(a, b) && (a[I]==b[I]);
    }
};


template <typename T>
struct pe<0,T> {

    static void zero(T* v) {
        v[0] = T(0);
    }

    static void neg(T* d, const T* s) {
        d[0] = -s[0];
    }

    static void copy(T* d, const T* s) {
        d[0] = s[0];
    }

    template <typename S>
    static void cast(T* d, const S* s) {
        d[0] = T(s[0]);
    }

    static T dot(const T* a, const T* b) {
        return a[0] * b[0];
    }

    static void add(T* r, const T* a, const T* b) {
        r[0] = a[0] + b[0];
    }

    static void add(T* r, const T* a, const T* b, const T* c) {
        r[0] = a[0] + b[0] + c[0];
    }

    static void sub(T* r, const T* a, const T* b) {
        r[0] = a[0] - b[0];
    }

    static void scale(T* r, const T* a, T s) {
        r[0] = a[0] * s;
    }

    static void scale(T* r, T s, const T* a) {
        r[0] = s * a[0];
    }

    static bool isZero(T const * v) {
        return (v[0]==T(0));
    }
    
    static bool isEqual(T const * a, T const * b) {
        return (a[0]==b[0]);
    }
    
};

} // end namespace math

} // end namespace Tn

#endif

