/Users/mourrain/Devel/mmx/shape/include/shape/fxv.hpp

Go to the documentation of this file.

#ifndef SYNAPS_VECTOR_FXV_H
#define SYNAPS_VECTOR_FXV_H
#include <shape/ssi_base_array_ops.hpp>
#include <realroot/assign.hpp>
#include <realroot/texp_expression.hpp>
#include <realroot/texp_operators.hpp>
#include <vector>

namespace mmx {

struct fxv_empty {};
inline std::ostream & operator<<( std::ostream & o, const fxv_empty & ) { return o; };

template<class C, unsigned N, class H> struct fxv;
template<class C, unsigned N>
struct fxv_dlink
{ 
  fxv< C, N, fxv_dlink > * l0, * l1; 
  fxv_dlink() : l0(0), l1(0){};
};

template<class C, unsigned N>
std::ostream& operator<<( std::ostream & o, const fxv_dlink<C,N>& a )
{
  o << a.l0 << ", " << a.l1 << std::endl; return o;
}; 

template< typename C, unsigned N, class H = fxv_empty >
struct fxv : H 
{
  typedef fxv<C,N,H> self_t;
  static const int _dimension_ = N;
  typedef C value_type;
  C data[N];
  C& operator[]( unsigned i ) { return data[i]; };
  const C& operator[]( unsigned i ) const  { return data[i]; };
  template<typename K,class H2>
  fxv& operator=( const fxv<K,N,H2>& v );
  self_t & operator=(const C& x );
  template<class S>
  self_t & operator=(const texp::template_expression<S> & M)
  { 
    using namespace let;
    assign(*this,M);
    return *this; 
  };
  bool operator==(int n) const;
  bool operator!=(int n) const {return !(*this==n);}
  bool operator==(const self_t & p) const;
  bool operator!=(const self_t & p) const {return !(*this==p);}
  template<class X> self_t& operator+=(const X& x) { add(*this,x); return *this; };
  template<class X> self_t& operator-=(const X& x) { sub(*this,x); return *this; };
  template<class X> self_t& operator*=(const X& x) { mul(*this,x); return *this; };
  template<class X> self_t& operator/=(const X& x) { div(*this,x); return *this; };
  template<class X> self_t& operator%=(const X& x) { mod(*this,x); return *this; };
};


template<typename C, unsigned N, class H> inline
void add( fxv<C,N,H>& a, const fxv<C,N,H>& b ) { add(a.data,b.data); };
template<typename C, unsigned N, class H> inline
void add( fxv<C,N,H>& a, const fxv<C,N,H>& b, const fxv<C,N,H>& c ) { add(a.data,b.data,c.data); };
template<typename C, unsigned N, class H> inline
void sub( fxv<C,N,H>& a, const fxv<C,N,H>& b ) { sub(a.data,b.data); };
template<typename C, unsigned N, class H> inline
void sub( fxv<C,N,H>& a, const fxv<C,N,H>& b, const fxv<C,N,H>& c) { sub(a.data,b.data,c.data); };
template<typename C, unsigned N, class H> inline
void scmul( fxv<C,N,H>& a, const C& c ) { scmul(a.data,c); };
template<typename C, unsigned N, class H> inline
void scmul( fxv<C,N,H>& a, const fxv<C,N,H>& b, const C& c ) { scmul(a.data,b.data,c); };
template<typename A, typename B, unsigned N> inline
void scdiv( fxv<A,N>& a, const B& s ) { scdiv(a.data,s); };
template<typename C, unsigned N, class H> inline
void scdiv( fxv<C,N,H>& a, const fxv<C,N,H>& b, const C& s ) { scdiv(a.data,b.data,s); };
template<typename C, unsigned N, class H> inline 
C dotprod( const fxv<C,N,H>& a, const fxv<C,N,H>& b ) { return dotprod(a.data,b.data); };
template<class C, unsigned N,class H> inline
C norm( const fxv<C,N,H>& v, int p ) { return norm(v.data,p); };        
template<class A, class B, unsigned N> inline
void init( const fxv<A,N>& v, const B& k ) { init(v.data,k); };
template<typename C, typename H> inline
void crossprod( fxv<C,3,H>& res, const fxv<C,3,H>& a, const fxv<C,3,H>& b )
{ res[0]=a[1]*b[2]-b[1]*a[2]; res[1]=a[2]*b[0]-b[2]*a[0]; res[2]=a[0]*b[1]-b[0]*a[1]; };
template<typename C, unsigned N, class H> inline
void add( fxv<C,N,H>& a, const C& c ) { scadd(a,c); };

template<typename C, unsigned N, class H>
void add( fxv<C,N,H>& a, const fxv<C,N,H>& b, const C& c ) { scadd(a,b,c); };

template<typename C, unsigned N, class H> inline
void sub( fxv<C,N,H>& a, const C& c ) { scsub(a,c); };

template<typename C, unsigned N, class H> inline
void mul( fxv<C,N,H>& a, const C& c ) { scmul(a,c); };

template<typename C, unsigned N, class H> inline
void mul( fxv<C,N,H>& a, const fxv<C,N,H>& b, const C& c )  { scmul(a.data,b.data,c); };

template<typename C, unsigned N, class H> inline
void mul( fxv<C,N,H>& a, const C& c, const fxv<C,N,H>& b )  { scmul(a.data,b.data,c); };

template<typename C, unsigned N, class H> inline
void mul( C& r, const fxv<C,N,H>& a, const fxv<C,N,H>& b ) { r = dotprod(a,b); };

template<typename C, unsigned N, class H> inline
void div( fxv<C,N,H>& a, const C& s ) { scdiv(a.data,s); };

template<typename C, unsigned N, class H> inline
void div( fxv<C,N,H>& a, const fxv<C,N,H>& b, const C& s ) { scdiv(a.data,b.data,s); };

template<typename C, unsigned N, class H> inline
void mod( fxv<C,N,H>& a, const fxv<C,N,H>& b, const fxv<C,N,H>& c ) { a = b-((b*c)/(c*c))*c; };

template<typename C, unsigned N, class R > inline
C norm2( const fxv<C,N,R>& v ) { return dotprod(v,v); };

template<typename C, unsigned N, class H > inline
C norm( const fxv<C,N,H>& v  ) { return sqrt(norm2(v)); };

template<typename C, unsigned N, class H > inline
C max( const fxv<C,N,H>& v ) { return max(v.data); };

template<typename C, unsigned N, class H > inline 
C min( const fxv<C,N,H>& v ) { return min(v); };

template<typename C, unsigned N, class H> inline
void urand( fxv<C,N,H>&  v, const C& a = (C)0.0, const C& b = (C)1.0 ) { urand(v.data,a,b); };

template<typename C, unsigned N, class H> inline
std::ostream& operator<<( std::ostream& o, const fxv<C,N,H>& v ) { 
  o << *((const H*)&v);
  o << "[ data = ";
  print(o,v.data);
  o << "]";
  return o;
};

template<typename C, unsigned N, class H> inline
void init( fxv<C,N,H>& v, const C& k ) { init(v.rep(),k); };


namespace texp {
  
  template<typename Ca,typename Cb,unsigned N,class H>
    struct ringof< fxv< Ca, N, H >, fxv< Cb, N, H > > { typedef fxv< typename ringof<Ca,Cb>::T,N,H> T; };
  template<typename Ca,typename Cb,unsigned N,class H>
  struct binary_operator_prototype< _mul_, fxv<Ca,N,H>, fxv<Cb,N,H> > 
  {
    typedef typename ringof< fxv<Ca,N,H>, fxv<Cb,N,H> >::T U; 
    typedef U V;
    typedef typename U::value_type  F;
  };
  // scalar * vector prototype.
  template < typename Ca, typename Cb, unsigned N, class H >
  struct binary_operator_prototype< _mul_, fxv<Ca,N,H>, Cb >
  {
    typedef typename ringof<Ca,Cb>::T               C;
    typedef C V;
    typedef fxv<C,N,H>                             F;
    typedef F                                       U;
  };
  
  template < typename Ca, typename Cb, unsigned N, class H >
  struct binary_operator_prototype< _div_, fxv<Ca,N,H>, Cb >
  {
    typedef typename fieldof<Ca,Cb>::T               C;
    typedef C V;
    typedef fxv<C,N,H>                             F;
    typedef F                                       U;
  };
  
  template < typename Ca, typename Cb, unsigned N, class H >
  struct binary_operator_prototype< _div_, Cb, fxv<Ca,N,H> >
  {
    typedef binary_operator_prototype< _div_, fxv<Ca,N,H >, Cb > X;
    typedef typename X::V U;
    typedef typename X::U V;
    typedef typename X::F F;
  };
  
  template < typename Ca, typename Cb, unsigned N, class H >
  struct binary_operator_prototype< _mul_, Cb, fxv<Ca,N,H> >
  {
    typedef binary_operator_prototype< _mul_, fxv<Ca,N,H >, Cb > X;
    typedef typename X::V U;
    typedef typename X::U V;
    typedef typename X::F F;
  };
};

template<class C, unsigned N, class H>
bool fxv<C,N,H>::operator==(const fxv<C,N,H>& v ) const 
{
  return eqxual(this->data,v.data); // type::equal
};

template<class C, unsigned N, class H>
template<typename C2,class H2>
fxv<C,N,H>& fxv<C,N,H>::operator=( const fxv<C2,N,H2>& v ) { copy(this->data,v.data); };

template<class C, unsigned N, class H> 
C  max_abs( const fxv<C,N,H> & v )
{
  using std::abs;
  C s(abs(v[0]));
  for ( unsigned i = 1; i < N; i ++ )
    s = std::max(s,abs(v[i]));
  return s;
};

template<class C, unsigned N, class H> inline
const C & distance( const fxv<C,N,H> & a, const fxv<C,N,H> & b )
{
  return distance(a.data,b.data);
};

namespace texp
{
  template<class C, int N, class H>
  struct structureof< fxv<C,N,H> > { typedef structure::vector T ; };
};
#define head template<typename Ca,typename Cb,unsigned N,class H>
#define parm0 fxv<Ca,N,H>
#define parm1 fxv<Cb,N,H>
declare_binary_operator(head,parm0,parm1,texp::_add_,operator+);
declare_binary_operator(head,parm0,parm1,texp::_add_,operator-);
declare_binary_operator(head,parm0,parm1,texp::_mul_,operator*);
#undef  head
#define head template<class C,unsigned N, class H>
#undef  parm0
#define parm0 fxv<C,N,H>
#undef  parm1
#define parm1 typename fxv<C,N,H>::value_type
declare_binary_operator(head,parm0,parm1,texp::_mul_,operator*);
declare_binary_operator(head,parm1,parm0,texp::_mul_,operator*);
declare_binary_operator(head,parm0,parm1,texp::_div_,operator/);
#undef  parm1
#undef  head
#define head template<class C, class H, class K, unsigned N>
#define parm1 typename K::integer
declare_binary_operator(head,parm0,parm1,texp::_mul_,operator*);
declare_binary_operator(head,parm1,parm0,texp::_mul_,operator*);
declare_binary_operator(head,parm0,parm1,texp::_div_,operator/);
#undef  parm1
#define parm1 typename K::rational
declare_binary_operator(head,parm0,parm1,texp::_mul_,operator*);
declare_binary_operator(head,parm1,parm0,texp::_mul_,operator*);
declare_binary_operator(head,parm0,parm1,texp::_div_,operator/);
#undef  parm1
#define parm1 typename K::floating
declare_binary_operator(head,parm0,parm1,texp::_mul_,operator*);
declare_binary_operator(head,parm1,parm0,texp::_mul_,operator*);
declare_binary_operator(head,parm0,parm1,texp::_div_,operator/);
#undef head
#undef parm0
#undef parm1


} //namespace mmx

#endif