//#####################################################################
// Copyright 2006-2007, Geoffrey Irving, Craig Schroeder, Andrew Selle, Tamar Shinar.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
// Class VECTOR
//#####################################################################
#ifndef __VECTOR__
#define __VECTOR__

#include <Core/Arrays/ARRAY_BASE.h>
#include <Core/Data_Structures/ELEMENT_ID.h>
#include <Core/Math_Tools/Componentwise_Min_Max.h>
#include <Core/Math_Tools/int_div.h>
#include <Core/Math_Tools/wrap.h>
#include <cmath>
#include <numbers>
#include <algorithm>
#ifdef DIFFERENCE // Windows workaround.
#undef DIFFERENCE
#endif
namespace PhysBAM{

using ::std::exp;
using ::std::sin;
using ::std::cos;
using ::std::abs;
using ::std::floor;
using ::std::rint;
using ::std::ceil;
using ::std::sqrt;
using ::std::log;
using ::std::numbers::pi;
using ::std::numbers::pi_v;

template<class TV> class RANGE;
template<class T_ARRAY,class T_INDICES> class INDIRECT_ARRAY;
template<class T,int d> struct IS_ARRAY<VECTOR<T,d> > {static const bool value=true;};
template<class T,int d> struct HAS_POINTER<VECTOR<T,d> > {static const bool value=true;};

#define MAKE_VECTOR_DATA                                                \
    VECTOR_DATA():array{}{}                                             \
    template<class ...Args> VECTOR_DATA(const T& a,Args &&... args):array{a,std::forward<Args>(args)...}{} \
    template<int... i,class T2,class A> VECTOR_DATA(const ARRAY_BASE<T2,A>& v,std::integer_sequence<int,i...>):array{((T)v(i))...}{} \
    template<int... i,class A> VECTOR_DATA(const ARRAY_BASE<T,A>& v,std::integer_sequence<int,i...>):array{(v(i))...}{} \
    template<int... i> VECTOR_DATA(VECTOR_DATA&& v,std::integer_sequence<int,i...>):array{(std::move(v.array[i]))...}{}

#define COMMA ,
template<class T,int d> struct VECTOR_DATA
{
    T array[d+!d];
    MAKE_VECTOR_DATA;
};
template<class T> struct VECTOR_DATA<T,1>
{
    union{struct{T x;};T array[1];};
    MAKE_VECTOR_DATA;
    ~VECTOR_DATA(){x.~T();}
};
template<class T> struct VECTOR_DATA<T,2>
{
    union{struct{T x COMMA y;};T array[2];};
    MAKE_VECTOR_DATA;
    ~VECTOR_DATA(){x.~T();y.~T();}
};
template<class T> struct VECTOR_DATA<T,3>
{
    union{struct{T x COMMA y COMMA z;};T array[3];};
    MAKE_VECTOR_DATA;
    ~VECTOR_DATA(){x.~T();y.~T();z.~T();}
};
#undef COMMA

template<class T,int d>
class VECTOR:public VECTOR_DATA<T,d>,public ARRAY_BASE<T,VECTOR<T,d> >
{
    static_assert(d>=0);
    struct UNUSABLE{};
public:
    typedef VECTOR_DATA<T,d> DATA;
    static constexpr int m=d;
    static constexpr int fixed_m=m;
    static constexpr int type=vec_gen;
    typedef ARRAY_BASE<T,VECTOR<T,m>> BASE;
    typedef SCALAR_TYPE<T> SCALAR;
    typedef T ELEMENT;
    typedef std::conditional_t<(d>3||d==0),UNUSABLE,VECTOR<T,(d-(d<3))>> SPIN;
    typedef int INDEX;
    typedef T value_type; // for stl
    typedef T* iterator; // for stl
    typedef const T* const_iterator; // for stl
    typedef T& RESULT_TYPE;
    typedef const T& CONST_RESULT_TYPE;
    using BASE::Dot;using BASE::Fill_Uniform;
    using DATA::array;

    VECTOR() = default;

    explicit VECTOR(INITIAL_SIZE n)
    {
        assert(n==INITIAL_SIZE(d));
    }

    template<int m_=0,class R=std::enable_if_t<d==1+m_>>
    VECTOR(const T& a) :DATA{a}{}

    template<int m_=0,class R=std::enable_if_t<d==2+m_>>
    VECTOR(const T& a,const T& b) :DATA{a,b}{}

    template<int m_=0,class R=std::enable_if_t<d==3+m_>>
    VECTOR(const T& a,const T& b,const T& c) :DATA{a,b,c}{}

    template<class ...Args,class R=std::enable_if_t<d==4+sizeof...(Args)>>
    VECTOR(const T& a,const T& b,const T& c,const T& e,Args &&... args)
        :DATA{a,b,c,e,(const T&)args...}
    {}

    template<class T_VECTOR,class R=std::enable_if_t<compat(m,T_VECTOR::fixed_m)>>
    explicit VECTOR(const ARRAY_BASE<T,T_VECTOR>& v)
        :DATA(v,std::make_integer_sequence<int,d>{})
    {
        if constexpr(T_VECTOR::fixed_m<0){assert(v.Size()==Size());}
    }

    template<class T2,int d2,class R=std::enable_if_t<d2<=d>>
    explicit VECTOR(const VECTOR<T2,d2>& v)
        :DATA(v,std::make_integer_sequence<int,d2>{})
    {
    }

    VECTOR(const VECTOR& v)
        :DATA(v,std::make_integer_sequence<int,d>{})
    {
        static_assert(sizeof(VECTOR)==(d+!d)*sizeof(T));
    }

    VECTOR(VECTOR&& v)
        :DATA(std::move(v),std::make_integer_sequence<int,d>{})
    {}

    template<int n>
    VECTOR(const VECTOR<T,n>& v1,const VECTOR<T,d-n>& v2)
        :DATA(v1,std::make_integer_sequence<int,n>{})
    {
        std::copy(v2.begin(),v2.end(),begin()+n);
    }

    VECTOR& operator=(const VECTOR& v)
    {
        []<int...i>(auto&a,const auto&b,std::integer_sequence<int,i...>)
        {((a(i)=b(i)),...);}(*this,v,std::make_integer_sequence<int,d>{});
        return *this;
    }
    VECTOR& operator=(VECTOR&& v)
    {
        []<int...i>(auto&a,auto&&b,std::integer_sequence<int,i...>)
        {((a(i)=std::move(b(i))),...);}(*this,v,std::make_integer_sequence<int,d>{});
        return *this;
    }

    template<class T_VECTOR,class R=std::enable_if_t<compat(m,T_VECTOR::fixed_m)>>
    VECTOR& operator=(const ARRAY_BASE<T,T_VECTOR>& v)
    {
        if constexpr(T_VECTOR::fixed_m<0){assert(v.Size()==Size());}
        []<int...i>(auto&a,const auto&b,std::integer_sequence<int,i...>)
        {((a(i)=b(i)),...);}(*this,v,std::make_integer_sequence<int,d>{});
        return *this;
    }

    ~VECTOR() = default;

    int Size() const
    {return m;}

    const T& operator[](const int i) const
    {assert((unsigned)i<d);return array[i];}

    T& operator[](const int i)
    {assert((unsigned)i<d);return array[i];}

    const T& operator()(const int i) const
    {assert((unsigned)i<d);return array[i];}

    T& operator()(const int i)
    {assert((unsigned)i<d);return array[i];}

    bool operator==(const VECTOR& v) const
    {
        return []<int...i>(const auto&a,const auto&b,std::integer_sequence<int,i...>)
        {return ((a(i)==b(i))&&...);}(*this,v,std::make_integer_sequence<int,d>{});
    }

    bool operator!=(const VECTOR& v) const
    {return !((*this)==v);}

    FIXED_NUMBER<T,0> Dot(const ZERO_VECTOR<T,m>&) const
    {return FIXED_NUMBER<T,0>();}

    template<int p>
    auto Cross(const VECTOR<T,p>& v) const
    {
        if constexpr(m==0&&p==1) return VECTOR<T,1>();
        else if constexpr(m==1&&p==0) return VECTOR<T,1>();
        else if constexpr(m==1&&p==1) return VECTOR<T,0>();
        else if constexpr(m==1&&p==2) return VECTOR<T,2>(-array[0]*v(1),array[0]*v(0));
        else if constexpr(m==2&&p==1) return VECTOR<T,2>(array[1]*v(0),-array[0]*v(0));
        else if constexpr(m==2&&p==2) return VECTOR<T,1>(array[0]*v(1)-array[1]*v(0));
        else
        {
            static_assert(m==3&&p==3);
            return VECTOR<T,3>(
                array[1]*v(2)-array[2]*v(1),
                array[2]*v(0)-array[0]*v(2),
                array[0]*v(1)-array[1]*v(0));
        }
    }

    template<int p> auto
    Cross(const ZERO_VECTOR<T,p>& v) const
        -> ZERO_VECTOR<T,decltype(Cross(VECTOR<T,p>()))::m>
    {return {};}

    T Angle_Between(const VECTOR& v) const // 0 .. pi
    {
        if constexpr(m==0) return 0;
        else if constexpr(m==1) return (array[0]>0)!=(v(0)>0)?pi_v<T>:0;
        else if constexpr(m==2 || m==3)
        {
            T s=Cross(v).Magnitude(),c=Dot(v);
            return atan2(s,c);
        }
        else return BASE::Angle_Between(v);
    }

    T Oriented_Angle_Between(const VECTOR& v)
    {
        static_assert(d==2);
        T s=Cross(v).x,c=Dot(v);
        return atan2(s,c);
    }

    VECTOR Normalized() const
    {VECTOR r(*this);r.Normalize();return r;}

    bool Elements_Equal() const
    {bool equal=true;
    for(int i=1;i<d;i++) equal&=(array[i]==array[0]);
    return equal;}

    VECTOR Componentwise_Min(const VECTOR& v) const
    {VECTOR r;
    for(int i=0;i<d;i++) r.array[i]=min(array[i],v.array[i]);
    return r;}

    VECTOR Componentwise_Max(const VECTOR& v) const
    {VECTOR r;
    for(int i=0;i<d;i++) r.array[i]=max(array[i],v.array[i]);
    return r;}

    VECTOR Projected_On_Unit_Direction(const VECTOR& direction) const
    {return Dot(direction)*direction;}

    VECTOR Projected(const VECTOR& direction) const // un-normalized direction
    {return Dot(direction)/direction.Magnitude_Squared()*direction;}

    void Project_On_Unit_Direction(const VECTOR& direction)
    {*this=Dot(direction)*direction;}

    void Project(const VECTOR& direction) // un-normalized direction
    {*this=Dot(direction)/direction.Magnitude_Squared()*direction;}

    VECTOR Projected_Orthogonal_To_Unit_Direction(const VECTOR& direction) const
    {return *this-Dot(direction)*direction;}

    void Project_Orthogonal_To_Unit_Direction(const VECTOR& direction)
    {*this-=Dot(direction)*direction;}

    int Number_True() const
    {static_assert((is_same<T,bool>::value));int count=0;
    for(int i=0;i<d;i++)if(array[i]) count++;
    return count;}

    static VECTOR Axis_Vector(const int axis)
    {VECTOR r;r(axis)=(T)1;return r;}

    VECTOR Add_Axis(int axis,T value) const
    {VECTOR vec=*this;vec(axis)+=value;return vec;}

    static VECTOR Constant_Vector(const T& constant)
    {VECTOR r;
    for(int i=0;i<d;i++) r.array[i]=constant;
    return r;}

    static VECTOR All_Ones_Vector()
    {return Constant_Vector(1);}

    void Fill(const T& constant)
    {for(int i=0;i<d;i++) array[i]=constant;}

    template<class ...Args,class R=std::enable_if_t<d==sizeof...(Args)>>
    void Get(Args&&... args) const
    {
        int k=0;
        ((args=array[k++]),...);
    }

    template<class ...Args,class R=std::enable_if_t<d==sizeof...(Args)>>
    void Set(Args&&... args)
    {
        int k=0;
        ((array[k++]=args),...);
    }

    int Find(const T& element) const
    {for(int i=0;i<d;i++) if(array[i]==element) return i;
    return -1;}

    bool Contains(const T& element) const
    {for(int i=0;i<d;i++) if(array[i]==element) return true;
    return false;}

    template<class T_ARRAY>
    bool Contains_Any(const T_ARRAY& elements) const
    {static_assert((is_same<typename T_ARRAY::ELEMENT,T>::value));
    for(int i=0;i<elements.Size();i++) if(Contains(elements(i))) return true;
    return false;}

    VECTOR<T,d-1> Remove_Index(const int index) const
    {assert((unsigned)index<d);VECTOR<T,d-1> r;
    for(int i=0;i<d-1;i++) r[i]=(*this)[i+(i>=index)];
    return r;}

    VECTOR<T,d+1> Insert(const T& element,const int index) const
    {VECTOR<T,d+1> r;r[index]=element;
    for(int i=0;i<d;i++) r[i+(i>=index)]=(*this)[i];
    return r;}

    VECTOR<T,d+1> Prepend(const T& element) const
    {VECTOR<T,d+1> r;r[0]=element;
    for(int i=0;i<d;i++) r[i+1]=(*this)[i];
    return r;}

    VECTOR<T,d+1> Append(const T& element) const
    {VECTOR<T,d+1> r;
    for(int i=0;i<d;i++) r[i]=(*this)[i];
    r[d]=element;
    return r;}

    template<int d2> VECTOR<T,d+d2> Append(const VECTOR<T,d2>& elements) const
    {VECTOR<T,d+d2> r;
    for(int i=0;i<d;i++) r[i]=(*this)[i];
    for(int i=0;i<d2;i++) r[i+d]=elements[i];
    return r;}

    VECTOR Sorted() const
    {VECTOR r(*this);r.Sort();return r;}

    VECTOR Reversed() const
    {VECTOR r;
    for(int i=0;i<d;i++) r.array[d-i-1]=array[i];
    return r;}

    template<int n> void Split(VECTOR<T,n>& v1,VECTOR<T,d-n>& v2) const
    {for(int i=0;i<n;i++) v1(i)=(*this)(i);
    for(int i=n;i<d;i++) v2(i-n)=(*this)(i);}

    VECTOR Orthogonal_Vector() const
    {
        if constexpr(d<=1){PHYSBAM_FATAL_ERROR();}
        else if constexpr(d==2) return VECTOR(-array[1],array[0]);
        else if constexpr(d==3)
        {
            T abs_x=abs(array[0]),abs_y=abs(array[1]),abs_z=abs(array[2]);
            if(abs_x<abs_y) return abs_x<abs_z?VECTOR(0,array[2],-array[1]):VECTOR(array[1],-array[0],0);
            else return abs_y<abs_z?VECTOR(-array[2],0,array[0]):VECTOR(array[1],-array[0],0);
        }
        else
        {
            int a=Dominant_Axis(),b=!a;
            VECTOR r;
            r(a)=(*this)(b);
            r(b)=-(*this)(a);
            return r;
        }
    }

    VECTOR Unit_Orthogonal_Vector() const // roughly 6 mults, 2 adds, 1 div, 1 sqrt
    {return Orthogonal_Vector().Normalized();}
    
    VECTOR Rotate_Clockwise_90() const
    {static_assert(d==2);return VECTOR(array[1],-array[0]);}

    VECTOR Rotate_Counterclockwise_90() const
    {static_assert(d==2);return VECTOR(-array[1],array[0]);}

    VECTOR Perpendicular() const
    {static_assert(d==2);return VECTOR(-array[1],array[0]);}

    template<class T_VECTOR>
    void Set_Subvector(const int istart,const T_VECTOR& v)
    {for(int i=0;i<v.Size();i++) (*this)(istart+i)=v(i);}

    template<class T_VECTOR>
    void Add_Subvector(const int istart,const T_VECTOR& v)
    {for(int i=0;i<v.Size();i++) (*this)(istart+i)+=v(i);}
    
    template<class T_VECTOR>
    void Get_Subvector(const int istart,T_VECTOR& v) const
    {for(int i=0;i<v.Size();i++) v(i)=(*this)(istart+i);}

    template<class F>
    auto Map(F f) const -> VECTOR<std::remove_const_t<std::remove_reference_t<decltype(f(T()))>>,m>
    {
        VECTOR<std::remove_const_t<std::remove_reference_t<decltype(f(T()))>>,m> r;
        for(int i=0;i<m;i++) r(i)=f(array[i]);
        return r;
    }
    
    T* Get_Array_Pointer()
    {return array;}

    const T* Get_Array_Pointer() const
    {return array;}

    T* begin() // for stl
    {return array;}

    const T* begin() const // for stl
    {return array;}

    T* end() // for stl
    {return array+d;}

    const T* end() const // for stl
    {return array+d;}

    template<class RW> void Read(UNTYPED_ISTREAM<RW>& input)
    {Read_Binary_Array(input,array,d);}

    template<class RW> void Write(UNTYPED_OSTREAM<RW>& output) const
    {Write_Binary_Array(output,array,d);}

    void Fill_Uniform(RANDOM_NUMBERS<T>& rand,const RANGE<VECTOR>& box)
    {for(int i=0;i<m;i++) rand.Fill_Uniform(array[i],box.min_corner(i),box.max_corner(i));}

//#####################################################################
    int Dominant_Axis() const
    {return this->Arg_Maximum_Magnitude();}
//#####################################################################
};

template<class T> inline T
Triple_Product(const VECTOR<T,3>& u,const VECTOR<T,3>& v,const VECTOR<T,3>& w)
{
    return u.Dot(v.Cross(w));
}

template<class T,int d> inline VECTOR<T,d>
Inverse(const VECTOR<T,d>& v)
{
    VECTOR<T,d> r;
    for(int i=0;i<d;i++) r[i]=1/v[i];
    return r;
}
//#####################################################################
// Function abs
//#####################################################################
template<class T,int d> inline VECTOR<T,d>
abs(const VECTOR<T,d>& v)
{
    VECTOR<T,d> r;
    for(int i=0;i<d;i++) r[i]=abs(v[i]);
    return r;
}
//#####################################################################
// Function floor
//#####################################################################
template<class T,int d> inline VECTOR<T,d>
floor(const VECTOR<T,d>& v)
{
    VECTOR<T,d> r;
    for(int i=0;i<d;i++) r[i]=floor(v[i]);
    return r;
}
//#####################################################################
// Function ceil
//#####################################################################
template<class T,int d> inline VECTOR<T,d>
ceil(const VECTOR<T,d>& v)
{
    VECTOR<T,d> r;
    for(int i=0;i<d;i++) r[i]=ceil(v[i]);
    return r;
}
//#####################################################################
// Function rint
//#####################################################################
template<class T,int d> inline VECTOR<T,d>
rint(const VECTOR<T,d>& v)
{
    VECTOR<T,d> r;
    for(int i=0;i<d;i++) r[i]=rint(v[i]);
    return r;
}
//#####################################################################
// Function fdiv
//#####################################################################
template<int d> inline VECTOR<int,d>
fdiv(const VECTOR<int,d>& a,const VECTOR<int,d>& b)
{
    VECTOR<int,d> r;
    for(int i=0;i<d;i++) r[i]=fdiv(a[i],b[i]);
    return r;
}
//#####################################################################
// Function fdiv
//#####################################################################
template<int d> inline VECTOR<int,d>
fdiv(const VECTOR<int,d>& a,int b)
{
    VECTOR<int,d> r;
    for(int i=0;i<d;i++) r[i]=fdiv(a[i],b);
    return r;
}
//#####################################################################
// Function cdiv
//#####################################################################
template<int d> inline VECTOR<int,d>
cdiv(const VECTOR<int,d>& a,const VECTOR<int,d>& b)
{
    VECTOR<int,d> r;
    for(int i=0;i<d;i++) r[i]=cdiv(a[i],b[i]);
    return r;
}
//#####################################################################
// Function cdiv
//#####################################################################
template<int d> inline VECTOR<int,d>
cdiv(const VECTOR<int,d>& a,int b)
{
    VECTOR<int,d> r;
    for(int i=0;i<d;i++) r[i]=cdiv(a[i],b);
    return r;
}
//#####################################################################
// Function exp
//#####################################################################
template<class T,int d> inline VECTOR<T,d>
exp(const VECTOR<T,d>& v)
{
    VECTOR<T,d> r;
    for(int i=0;i<d;i++) r[i]=exp(v[i]);
    return r;
}
//#####################################################################
// Function log
//#####################################################################
template<class T,int d> inline VECTOR<T,d>
log(const VECTOR<T,d>& v)
{
    VECTOR<T,d> r;
    for(int i=0;i<d;i++) r[i]=log(v[i]);
    return r;
}
//#####################################################################
// Function sin
//#####################################################################
template<class T,int d> inline VECTOR<T,d>
sin(const VECTOR<T,d>& v)
{
    VECTOR<T,d> r;
    for(int i=0;i<d;i++) r[i]=sin(v[i]);
    return r;
}
//#####################################################################
// Function cos
//#####################################################################
template<class T,int d> inline VECTOR<T,d>
cos(const VECTOR<T,d>& v)
{
    VECTOR<T,d> r;
    for(int i=0;i<d;i++) r[i]=cos(v[i]);
    return r;
}
//#####################################################################
// Function sqrt
//#####################################################################
template<class T,int d> inline VECTOR<T,d>
sqrt(const VECTOR<T,d>& v)
{
    VECTOR<T,d> r;
    for(int i=0;i<d;i++) r[i]=sqrt(v[i]);
    return r;
}
template<class T> inline std::enable_if_t<is_fundamental<T>::value,T>
clamp(const T& v,const T& a,const T& b)
{return std::clamp(v,a,b);}
//#####################################################################
// Function clamp
//#####################################################################
template<class T,int d> inline VECTOR<T,d>
clamp(const VECTOR<T,d>& v,const VECTOR<T,d>& vmin,const VECTOR<T,d>& vmax)
{
    VECTOR<T,d> r;
    for(int i=0;i<d;i++) r[i]=clamp(v[i],vmin[i],vmax[i]);
    return r;
}
//#####################################################################
// Function clamp
//#####################################################################
template<class T,int d> inline VECTOR<T,d>
clamp(const VECTOR<T,d>& v,const ND<T>& min,const ND<T>& max)
{
    VECTOR<T,d> r;
    for(int i=0;i<d;i++) r[i]=clamp(v[i],min,max);
    return r;
}
//#####################################################################
// Function Componentwise_Max
//#####################################################################
template<class T,int d> inline VECTOR<T,d>
Componentwise_Max(const VECTOR<T,d>& v,const VECTOR<T,d>& vmin)
{
    VECTOR<T,d> r;
    for(int i=0;i<d;i++) r[i]=Componentwise_Max(v[i],vmin[i]);
    return r;
}
//#####################################################################
// Function Componentwise_Max
//#####################################################################
template<class T,int d> inline VECTOR<T,d>
Componentwise_Max(const VECTOR<T,d>& v,const ND<T>& min)
{
    VECTOR<T,d> r;
    for(int i=0;i<d;i++) r[i]=Componentwise_Max(v[i],min);
    return r;
}
//#####################################################################
// Function Componentwise_Min
//#####################################################################
template<class T,int d> inline VECTOR<T,d>
Componentwise_Min(const VECTOR<T,d>& v,const VECTOR<T,d>& vmax)
{
    VECTOR<T,d> r;
    for(int i=0;i<d;i++) r[i]=Componentwise_Min(v[i],vmax[i]);
    return r;
}
//#####################################################################
// Function Componentwise_Min
//#####################################################################
template<class T,int d> inline VECTOR<T,d>
Componentwise_Min(const VECTOR<T,d>& v,const ND<T>& max)
{
    VECTOR<T,d> r;
    for(int i=0;i<d;i++) r[i]=Componentwise_Min(v[i],max);
    return r;
}
//#####################################################################
// Function wrap
//#####################################################################
template<class T,int d> inline VECTOR<T,d>
wrap(const VECTOR<T,d>& v,const VECTOR<T,d>& vmin,const VECTOR<T,d>& vmax)
{
    VECTOR<T,d> r;
    for(int i=0;i<d;i++) r(i)=wrap(v(i),vmin(i),vmax(i));
    return r;
}
//#####################################################################
// Function Zip
//#####################################################################
template<class T,int d,class F,class... Args> auto
Zip(F f,const VECTOR<T,d>& v,Args&&... args)
{
    VECTOR<std::remove_const_t<std::remove_reference_t<decltype(f(v(0),args(0)...))>>,d> r;
    for(int i=0;i<d;i++) r(i)=f(v(i),args(i)...);
    return r;
}
template<int d> VECTOR<bool,d> Complement(const VECTOR<bool,d>& v)
{
    VECTOR<bool,d> r;
    for(int i=0;i<d;i++) r(i)=!v(i);
    return r;
}
//#####################################################################
template<class T> struct HASH_REDUCE<VECTOR<T,0> >
{static int H(const VECTOR<T,0>& key){return missing_element_hash;}};
template<class T> struct HASH_REDUCE<VECTOR<T,1> >
{static int H(const VECTOR<T,1>& key){return int_hash(HASH_REDUCE<T>::H(key.x));}};
template<class T> struct HASH_REDUCE<VECTOR<T,2> >
{static int H(const VECTOR<T,2>& key){return int_hash(HASH_REDUCE<T>::H(key.x),HASH_REDUCE<T>::H(key.y));}};
template<class T> struct HASH_REDUCE<VECTOR<T,3> >
{static int H(const VECTOR<T,3>& key){return int_hash(HASH_REDUCE<T>::H(key.x),HASH_REDUCE<T>::H(key.y),HASH_REDUCE<T>::H(key.z));}};
template<class T> struct HASH_REDUCE<VECTOR<T,4> >
{static int H(const VECTOR<T,4>& key){return int_hash(int_hash(HASH_REDUCE<T>::H(key[0]),HASH_REDUCE<T>::H(key[1])),HASH_REDUCE<T>::H(key[2]),HASH_REDUCE<T>::H(key[3]));}};
template<class T,int d> struct HASH_REDUCE<VECTOR<T,d> >
{static int H(const VECTOR<T,d>& key){return Hash_Reduce_Array_Helper(key);}};
void Ignore(std::istream& input,char c);
template<class T,int d>
inline std::istream& operator>>(std::istream& input,VECTOR<T,d>& v)
{
    Ignore(input,'(');
    Ignore(input,'[');
    for(int i=0;i<d;i++) input>>v(i);
    Ignore(input,']');
    Ignore(input,')');
    return input;
}
}
#endif
