/Users/mourrain/Devel/mmx/basix/include/basix/vector.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : vector.hpp
* DESCRIPTION: Dense univariate vectors
* COPYRIGHT  : (C) 2004  Joris van der Hoeven
*******************************************************************************
* This software falls under the GNU general public license and comes WITHOUT
* ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
* If you don't have this file, write to the Free Software Foundation, Inc.,
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
******************************************************************************/

#ifndef __MMX_VECTOR_HPP
#define __MMX_VECTOR_HPP
#include <basix/iterator.hpp>
#include <basix/vector_naive.hpp>


namespace mmx {
#define Vector_variant(C) vector_variant_helper<C>::VV
#define TMPL_DEF template<typename C, typename V>
// NOTE: default argument vector_variant_helper<C>::VV specified in basix.hpp
#define TMPL template<typename C, typename V>
#define TMPLK template<typename C, typename V, typename K>
#define Format format<C>
#define Vector vector<C,V>
#define Vector_int vector<int,V>
#define Vector_rep vector_rep<C,V>
#define Abs_vector vector<Abs_type(C),V>
#define Real_vector vector<Real_type(C),V>
#define Center_vector vector<Center_type(C),V>
#define Radius_vector vector<Radius_type(C),V>
TMPL class vector;
TMPL inline nat N (const Vector& v);
TMPL inline C* seg (Vector& v);
TMPL inline const C* seg (const Vector& v);
TMPL inline bool is_a_scalar   (const Vector& v);
TMPL inline bool is_non_scalar (const Vector& v);

/******************************************************************************
* Dense vectors of arbitrary sizes
******************************************************************************/

TMPL_DEF
class vector_rep REP_STRUCT_1(C) {
public:
  typedef implementation<vector_linear,V> Vec;
private:
  C*   a;           // entries of vector
  nat  n;           // dimension of vector
  nat  l;           // allocated number of entries
  bool scalar_flag; // is the vector a scalar value?
public:
  inline vector_rep (C* a2, nat n2, bool flag, const Format& fm):
    Format (fm), a (a2), n (n2), l (n2), scalar_flag (flag) {}
  inline vector_rep (C* a2, nat n2, nat l2, bool flag, const Format& fm):
    Format (fm), a (a2), n (n2), l (l2), scalar_flag (flag) {}
  inline ~vector_rep () { mmx_delete<C> (a, l); }
  void resize (nat n2) {
    nat l2;
    if (n2 > l) l2= aligned_size<C,V> (max (n2, l<<1));
    else if (n2 < (l >> 1)) l2= aligned_size<C,V> (n2);
    else { n= n2; return; }
    C* b= mmx_new<C> (l2);
    Vec::copy (b, a, min (n, n2));
    mmx_delete<C> (a, l);
    a= b;
    n= n2;
    l= l2;
  }
  friend class Vector;
  friend nat N LESSGTR (const Vector& v);
  friend C* seg LESSGTR (Vector& v);
  friend const C* seg LESSGTR (const Vector& v);
  friend bool is_a_scalar   LESSGTR (const Vector& v);
  friend bool is_non_scalar LESSGTR (const Vector& v);
};

TMPL_DEF
class vector {
INDIRECT_PROTO_2 (vector, vector_rep, C, V)
public:
  typedef implementation<vector_linear,V> Vec;
public:
  vector () {
    nat n= Vec::def_len;
    nat l= aligned_size<C,V> (n);
    C* a= mmx_formatted_new<C> (l, Format (no_format ()));
    rep= new Vector_rep (a, n, l, false, Format (no_format ()));
  }
  vector (const Format& fm) {
    nat n= Vec::def_len;
    nat l= aligned_size<C,V> (n);
    C* a= mmx_formatted_new<C> (l, fm);
    rep= new Vector_rep (a, n, l, false, fm);
  }
  template<typename T> 
  vector (const T& c) {
    Format fm= as<Format > (get_format (c));
    nat n= Vec::init_len;
    nat l= aligned_size<C,V> (n);
    C* a= mmx_formatted_new<C> (l, fm);
    Vec::set (a, as<C> (c), n);
    rep= new Vector_rep (a, n, l, true, fm);
  }
  template<typename T>
  vector (const T& c, const Format& fm) {
    nat n= Vec::init_len;
    nat l= aligned_size<C,V> (n);
    C* a= mmx_formatted_new<C> (l, fm);
    Vec::set_as (a, c, n);
    rep= new Vector_rep (a, n, l, true, fm);
  }
  vector (const C& x) {
    Format fm= as<Format > (get_format (x));
    nat n= Vec::init_len;
    nat l= aligned_size<C,V> (n);
    C* a= mmx_formatted_new<C> (l, fm);
    Vec::set (a, x, n);
    rep= new Vector_rep (a, n, l, true, fm);
  }
  template<typename T, typename W>
  vector (const vector<T,W>& v) {
    Format fm= as<Format > (CF(v));
    nat l= aligned_size<C,V> (N(v));
    C* a= mmx_formatted_new<C> (l, fm);
    Vec::cast (a, seg (v), N(v));
    rep= new Vector_rep (a, N(v), l, is_a_scalar (v), fm);
  }
  template<typename T, typename W>
  vector (const vector<T,W>& v, const Format& fm) {
    nat l= aligned_size<C,V> (N(v));
    C* a= mmx_formatted_new<C> (l, fm);
    Vec::cast (a, seg (v), N(v));
    rep= new Vector_rep (a, N(v), l, is_a_scalar (v), fm);
  }
  template<typename T> 
  vector (const T& c, nat n) {
    Format fm= as<Format > (get_format (c));
    nat l= aligned_size<C,V> (n);
    C* a= mmx_formatted_new<C> (l, fm);
    Vec::set (a, as<C> (c), n);
    rep= new Vector_rep (a, n, l, true, fm);
  }
  vector (const C& x, nat n) {
    nat l= aligned_size<C,V> (n);
    C* a= mmx_formatted_new<C> (l, get_format (x));
    Vec::set (a, x, n);
    rep= new Vector_rep (a, n, l, false, get_format (x));
  }
  inline vector (C* a, nat n, const Format& fm) {
    rep= new Vector_rep (a, n, false, fm);
  }
  inline vector (C* a, nat n, nat l, const Format& fm) {
    rep= new Vector_rep (a, n, l, false, fm);
  }
  inline vector (C* a, nat n, nat l, bool flag, const Format& fm) {
    rep= new Vector_rep (a, n, l, flag, fm);
  }
  vector (const iterator<C>& it) {
    Format fm= CF(it);
    nat i, l= Vec::def_len;
    if (l==0) l= aligned_size<C,V> ((nat) 1);
    C* a= mmx_formatted_new<C> (l, fm);
    for (i=0; busy (it); i++, ++it) {
      if (i >= l) {
        C* b= mmx_formatted_new<C> (aligned_size<C,V> (l<<1), fm);
        Vec::copy (b, a, l);
        mmx_delete (a, l);
        a= b;
        l= l<<1;
      }
      a[i]= *it;
    }
    rep= new Vector_rep (a, i, l, false, fm);
  }

  inline const C& operator [] (nat i) const {
    VERIFY (is_non_scalar (*this), "non-scalar vector expected");
    VERIFY (i<N(*this), "index out of range");
    return rep->a[i]; }
  inline C& operator [] (nat i) {
    VERIFY (is_non_scalar (*this), "non-scalar vector expected");
    VERIFY (i<N(*this), "index out of range");
    secure(); return rep->a[i]; }
  inline C scalar () const {
    VERIFY (is_a_scalar (*this), "scalar vector expected");
    return rep->a[0]; }
  inline C& scalar () {
    VERIFY (is_a_scalar (*this), "scalar vector expected");
    return rep->a[0]; }
  inline nat size() const {
    return rep->n; }

  Vector& operator << (const C& append);
  Vector& operator << (const Vector& w);
};
INDIRECT_IMPL_2 (vector, vector_rep, typename C, C, typename V, V)
DEFINE_UNARY_FORMAT_2 (vector)

TMPL inline nat N (const Vector& v) { return v->n; }
TMPL inline Format CF (const Vector& v) { return v->tfm (); }
TMPL inline bool is_a_scalar (const Vector& v) { return v->scalar_flag; }
TMPL inline bool is_non_scalar (const Vector& v) { return !v->scalar_flag; }
TMPL inline const C& read (const Vector& v, nat i) { return v[i]; }
TMPL inline const C* seg (const Vector& v) { return v->a; }
TMPL inline C* seg (Vector& v) {
  VERIFY (is_non_scalar (v), "non-scalar vector expected");
  v.secure(); return v->a; }
template<typename C, typename V, typename C2, typename V2> inline Vector
extend (const Vector& v, const vector<C2,V2>& w) {
  VERIFY (is_a_scalar (v), "scalar vector expected");
  VERIFY (is_non_scalar (w), "non-scalar vector expected");
  return Vector (v.scalar(), N(w)); }
TMPL inline bool is_nil (const Vector& v) {
  return N(v) == 0; }
TMPL inline bool is_atom (const Vector& v) {
  return is_non_scalar (v) && N(v) == 1; }

STYPE_TO_TYPE(TMPL,scalar_type,Vector,C);
UNARY_RETURN_TYPE(TMPL,abs_op,Vector,Abs_vector);
UNARY_RETURN_TYPE(TMPL,Re_op,Vector,Real_vector);
UNARY_RETURN_TYPE(TMPL,center_op,Vector,Center_vector);
UNARY_RETURN_TYPE(TMPL,radius_op,Vector,Radius_vector);
TMPL struct rounding_helper<Vector >: public rounding_helper<C> {};
TMPL struct species_type_information<Vector > {
  static const nat id= SPECIES_VECTOR; };

TMPL
struct fast_helper<Vector > {
  typedef Fast_type(C) FC;
  typedef vector<FC> fast_type;
  typedef typename Vector_variant(FC) FV;
  static inline fast_type dd (const Vector& v) {
    format<FC> fm= fast (CF(v));
    nat l= aligned_size<FC, FV> (N(v));
    FC* r= mmx_formatted_new<FC> (l, fm);
    for (nat i=0; i<N(v); i++) r[i]= fast (v[i]);
    return fast_type (r, N(v), l, is_a_scalar (v), fm); }
  static inline Vector uu (const fast_type& v) {
    Format fm= slow<C> (CF(v));
    nat l= aligned_size<C,V> ((N(v)));
    C* r= mmx_formatted_new<C> (l, fm);
    for (nat i=0; i<N(v); i++) r[i]= slow<C> (v[i]);
    return Vector (r, N(v), l, is_a_scalar (v), fm); }
};

template<typename T, typename F, typename TV, typename FV>
struct as_helper<vector<T,TV>, vector<F,FV> > {
  static inline vector<T,TV>
  cv (const vector<F,FV>& v) {
    format<T> fm= as<format<T> > (CF(v));
    nat l= aligned_size<T,TV> (N(v));
    T* c= mmx_formatted_new<T> (l, fm);
    for (nat i=0; i<N(v); i++) c[i]= as<T> (v[i]);
    return vector<T,TV> (c, N(v), l, is_a_scalar (v), fm); }
};

template<typename T, typename F, typename TV, typename FV> inline void
set_as (vector<T,TV>& r, const vector<F,FV>& v) {
  r= vector<T,TV> (v, CF(r));
}

template<typename C, typename V, typename T> inline void
set_as (Vector& r, const T& x) {
  r= Vector (x, CF(r));
}

/******************************************************************************
* Abstract vectors
******************************************************************************/

TMPL
struct vector_as_helper<Vector > {
  static inline generic abstract (const Vector& x) {
    return as<generic> (as<vector<generic> > (x)); }
  static inline Vector concrete (const generic& x) {
    return as<Vector > (as<vector<generic> > (x)); }
};

inline vector<generic>
abstract_vector (const generic& conc) {
  return as<vector<generic> > (conc->make_abstract_vector ());
}

inline generic
concrete_vector (const vector<generic>& abst, const generic& conc) {
  return conc->make_concrete_vector (as<generic> (abst));
}

/******************************************************************************
* Fixed size dense vectors
******************************************************************************/

template<typename V, typename S>
struct vector_fixed: public V {
  typedef typename V::Naive Naive;
  typedef vector_fixed<typename V::Aligned,S> Aligned;
  typedef vector_fixed<typename V::No_simd,S> No_simd;
  typedef vector_fixed<typename V::No_thread,S> No_thread;
};

template<typename V, typename W, typename S>
struct implementation<vector_defaults,V,vector_fixed<W,S> > {
  static const nat def_len = S::val;
  static const nat init_len= S::val;
}; // implementation<vector_defaults,V,vector_fixed<W,S> >

template<typename F, typename V, typename W, typename S>
struct implementation<F,V,vector_fixed<W,S> >:
  public implementation<F,V,W> {};

#define FVector     vector    <C,vector_fixed<V,S> >
#define FVector_rep vector_rep<C,vector_fixed<V,S> >

template<typename C, typename V, typename S>
class FVector_rep REP_STRUCT_1(C) {
  C* a;
public:
  /*
  inline vector_rep (C* a2, nat n, bool dummy):
    a (a2) {
      VERIFY (n == S::val, "wrong size");
      (void) n; (void) dummy; }
  inline vector_rep (C* a2, nat n, nat l, bool dummy):
    a (a2) {
      VERIFY (n == S::val, "wrong size");
      VERIFY ((l == aligned_size<C,V> (S::val)), "wrong allocated size");
      (void) n; (void) l; (void) dummy; }
  */
  inline vector_rep (C* a2, nat n, bool dummy, const Format& fm):
    Format (fm), a (a2) {
      VERIFY (n == S::val, "wrong size");
      (void) n; (void) dummy; }
  inline vector_rep (C* a2, nat n, nat l, bool dummy, const Format& fm):
    Format (fm), a (a2) {
      VERIFY (n == S::val, "wrong size");
      VERIFY ((l == aligned_size<C,V> (S::val)), "wrong allocated size");
      (void) n; (void) l; (void) dummy; }
  inline virtual ~vector_rep () {
    mmx_delete<C> (a, aligned_size<C,V> (S::val)); }
  friend class FVector;
  friend nat N LESSGTR (const FVector& v);
  friend C*  seg LESSGTR (FVector& v);
  friend const C* seg LESSGTR (const FVector& v);
  friend bool is_a_scalar   LESSGTR (const Vector& v);
  friend bool is_non_scalar LESSGTR (const Vector& v);
};

template<typename C, typename V, typename S> inline nat
N (const FVector& v) { return S::val; }
template<typename C, typename V, typename S> inline bool
is_a_scalar (const FVector& v) { (void) v; return false; }
template<typename C, typename V, typename S> inline bool
is_non_scalar (const FVector& v) { (void) v; return true; }

#undef FVector
#undef FVector_rep

/******************************************************************************
* Explicit dense vectors
******************************************************************************/

template<typename C> vector<C>
vec () {
  Format fm;
  nat l= default_aligned_size<C> ((nat) 0);
  C* a= mmx_formatted_new<C> (l, fm);
  return vector<C> (a, 0, l, fm);
}

template<typename C> vector<C>
vec (const Format& fm) {
  nat l= default_aligned_size<C> ((nat) 0);
  C* a= mmx_formatted_new<C> (l, fm);
  return vector<C> (a, 0, l, fm);
}

template<typename C> vector<C>
vec (const C& c1) {
  Format fm= get_format (c1);
  nat l= default_aligned_size<C> ((nat) 1);
  C* a= mmx_formatted_new<C> (l, fm);
  a[0]= c1;
  return vector<C> (a, 1, l, fm);
}

template<typename C> vector<C>
vec (const C& c1, const C& c2) {
  Format fm= get_format (c1);
  nat l= default_aligned_size<C> ((nat) 2);
  C* a= mmx_formatted_new<C> (l, fm);
  a[0]= c1; a[1]= c2;
  return vector<C> (a, 2, l, fm);
}

template<typename C> vector<C>
vec (const C& c1, const C& c2, const C& c3) {
  Format fm= get_format (c1);
  nat l= default_aligned_size<C> ((nat) 3);
  C* a= mmx_formatted_new<C> (l, fm);
  a[0]= c1; a[1]= c2; a[2]= c3;
  return vector<C> (a, 3, l, fm);
}

template<typename C> vector<C>
vec (const C& c1, const C& c2, const C& c3, const C& c4) {
  Format fm= get_format (c1);
  nat l= default_aligned_size<C> ((nat) 4);
  C* a= mmx_formatted_new<C> (l, fm);
  a[0]= c1; a[1]= c2; a[2]= c3; a[3]= c4;
  return vector<C> (a, 4, l, fm);
}

template<typename C> vector<C>
vec (const C& c1, const C& c2, const C& c3, const C& c4, const C& c5) {
  Format fm= get_format (c1);
  nat l= default_aligned_size<C> ((nat) 5);
  C* a= mmx_formatted_new<C> (l, fm);
  a[0]= c1; a[1]= c2; a[2]= c3; a[3]= c4; a[4]= c5;
  return vector<C> (a, 5, l, fm);
}

template<typename C> vector<C>
vec (const C& c1, const C& c2, const C& c3,
     const C& c4, const C& c5, const C& c6) {
  Format fm= get_format (c1);
  nat l= default_aligned_size<C> ((nat) 6);
  C* a= mmx_formatted_new<C> (l, fm);
  a[0]= c1; a[1]= c2; a[2]= c3; a[3]= c4; a[4]= c5; a[5]= c6;
  return vector<C> (a, 6, l, fm);
}

template<typename C> vector<C>
vec (const C& c1, const C& c2, const C& c3, const C& c4,
     const C& c5, const C& c6, const C& c7) {
  Format fm= get_format (c1);
  nat l= default_aligned_size<C> ((nat) 7);
  C* a= mmx_formatted_new<C> (l, fm);
  a[0]= c1; a[1]= c2; a[2]= c3; a[3]= c4; a[4]= c5; a[5]= c6; a[6]= c7;
  return vector<C> (a, 7, l, fm);
}

template<typename C> vector<C>
vec (const C& c1, const C& c2, const C& c3, const C& c4,
     const C& c5, const C& c6, const C& c7, const C& c8) {
  Format fm= get_format (c1);
  nat l= default_aligned_size<C> ((nat) 8);
  C* a= mmx_formatted_new<C> (l, fm);
  a[0]= c1; a[1]= c2; a[2]= c3; a[3]= c4;
  a[4]= c5; a[5]= c6; a[6]= c7; a[7]= c8;
  return vector<C> (a, 8, l, fm);
}

template<typename C> inline vector<C>
fill (nat n) {
  Format fm;
  nat l= default_aligned_size<C> (n);
  C* a= mmx_formatted_new<C> (l, fm);
  return vector<C> (a, n, l, fm);
}

template<typename C> inline vector<C>
fill (nat n, const Format& fm) {
  nat l= default_aligned_size<C> (n);
  C* a= mmx_formatted_new<C> (l, fm);
  return vector<C> (a, n, l, fm);
}

template<typename C> inline vector<C>
fill (const C& c, nat n) {
  Format fm= get_format (c);
  nat l= default_aligned_size<C> (n);
  C* a= mmx_formatted_new<C> (l, fm);
  for (nat i=0; i<n; i++) a[i]= c;
  return vector<C> (a, n, l, fm);
}

/******************************************************************************
* Vector iteration
******************************************************************************/

TMPL_DEF
class vector_iterator_rep: public iterator_rep<C> {
  const Vector v;  // the vector to traverse
  nat i;           // current position
public:
  vector_iterator_rep (const Vector& v2, nat i2= 0):
    iterator_rep<C> (CF (v2)), v (v2), i (0) { (void) i2; }
protected:
  bool is_busy () { return i < N(v); }
  void advance () { i++; }
  C current () { return v[i]; }
  iterator_rep<C>* clone () { return new vector_iterator_rep (v, i); }
};

TMPL iterator<C>
iterate (const Vector& v) {
  return iterator<C> (new vector_iterator_rep<C,V> (v));
}

/******************************************************************************
* Hashing and generic output routines for iterators and vectors
******************************************************************************/

TMPL syntactic
flatten (const Vector& v) {
  if (is_a_scalar (v)) return flatten (v.scalar());
  nat i, n= N(v);
  vector<syntactic> r= fill <syntactic> (n);
  for (i=0; i<n; i++)
    r[i]= flatten (v[i]);
  return apply (GEN_SQTUPLE, r);
}

TMPL
struct binary_helper<Vector > {
  static inline string short_type_name () {
    return "V" * Short_type_name (C); }
  static inline generic full_type_name () {
    return gen ("Vector", Full_type_name (C)); }
  static nat size (const Vector& v) {
    if (is_a_scalar (v)) return 0;
    else return N(v); }
  static generic access (const Vector& v, nat i) {
    ASSERT (!is_a_scalar (v), "non-scalar vector expected");
    return as<generic> (v[i]); }
  static generic disassemble (const Vector& v) {
    if (is_a_scalar (v)) return as<generic> (v.scalar ());
    return as<generic> (as<vector<generic> > (v)); }
  static Vector assemble (const generic& v) {
    // FIXME: what to do with the format?
    if (!is<vector<generic> > (v)) return Vector (as<C> (v));
    else return as<Vector > (as<vector<generic> > (v)); }
  static void write (const port& out, const Vector& v) {
    binary_write<Format > (out, CF (v));
    binary_write<bool> (out, is_a_scalar (v));
    if (is_a_scalar (v)) binary_write<C> (out, v.scalar ());
    else {
      binary_write<nat> (out, N(v));
      for (nat i=0; i<N(v); i++)
        binary_write<C> (out, v[i]); } }
  static Vector read (const port& in) {
    Format fm= binary_read<Format > (in);
    bool sc= binary_read<bool> (in);
    if (sc) return Vector (binary_read<C> (in));
    else {
      nat n= binary_read<nat> (in);
      nat l= aligned_size<C,V> (n);
      C*  r= mmx_new<C> (l);
      for (nat i=0; i<n; i++)
        r[i]= binary_read<C> (in);
      return Vector (r, n, l, fm); } }
};

/******************************************************************************
* Extracting ranges of elements and appending new elements
******************************************************************************/

TMPL Vector
range (const Vector& v, nat start, nat end) {
  typedef implementation<vector_linear,V> Vec;
  VERIFY(end >= start, "bad range");
  nat n= end - start;
  if (is_a_scalar (v)) return Vector (v.scalar(), n);
  nat l= aligned_size<C,V> (n);
  C* a= mmx_formatted_new<C> (l, CF(v));
  Vec::copy (a, seg (v) + start, n);
  return Vector (a, n, l, CF(v));
}

TMPL Vector
extract_mod (const Vector& v, nat k, nat p) {
  nat n= (N(v) - k + p - 1) / p;
  nat l= aligned_size<C,V> (n);
  C* r= mmx_formatted_new<C> (l, CF(v));
  for (nat i=0; i<n; i++) r[i]= v[i*p+k];
  return Vector (r, n, l, CF(v));
}

TMPL Vector
extract (const Vector& v, vector<nat> e) {
  nat l= aligned_size<C,V> (N(e));
  C* r= mmx_formatted_new<C> (l, CF(v));
  for (nat i=0; i<N(e); i++)
    if (e[i] < N(v)) r[i]= v[e[i]];
  return Vector (r, N(e), l, CF(v));
}

TMPL Vector
remove (const Vector& v, vector<nat> e) {
  nat n= N(v);
  vector<nat> f ((nat) 0, n);
  for (nat i= 0; i < n; i++) f[i]= i;
  for (nat i= 0; i < N(e); i++) if (e[i] < n) f[e[i]]= n;
  Vector r;
  for (nat i= 0; i < n; i++) {
    if (f[i] == n) continue;
    r << v[i];
  }
  return r;
}

TMPL inline C car (const Vector& v) { return v[0]; }
TMPL inline Vector cdr (const Vector& v) { return range (v, 1, N(v)); }

TMPL Vector
append (const Vector& v, const C& c) {
  typedef implementation<vector_linear,V> Vec;
  ASSERT (is_non_scalar (v), "non-scalar vector expected");
  nat n= N(v);
  nat l= aligned_size<C,V> (n+1);
  C* a= mmx_formatted_new<C> (l, CF(v));
  a[n]= c;
  Vec::copy (a, seg (v), n);
  return Vector (a, n+1, l, CF(v));
}

TMPL Vector
append (const Vector& v, const Vector& w) {
  typedef implementation<vector_linear,V> Vec;
  ASSERT (is_non_scalar (v) && is_non_scalar (w),
          "non-scalar vectors expected");
  nat n= N(v), p= N(w);
  nat l= aligned_size<C,V> (n+p);
  C* a= mmx_formatted_new<C> (l, CF(v));
  Vec::copy (a, seg (v), n);
  Vec::copy (a+n, seg (w), p);
  return Vector (a, n+p, l, CF(v));
}

TMPL Vector
cons (const C& c, const Vector& v) {
  typedef implementation<vector_linear,V> Vec;
  ASSERT (is_non_scalar (v), "non-scalar vector expected");
  nat n= N(v);
  nat l= aligned_size<C,V> (n+1);
  C* a= mmx_formatted_new<C> (l, CF(v));
  a[0]= c;
  Vec::copy (a+1, seg (v), n);
  return Vector (a, n+1, l, CF(v));
}

TMPL Vector&
Vector::operator << (const C& append) {
  ASSERT (is_non_scalar (*this), "non-scalar vector expected");
  secure ();
  nat n= rep->n;
  rep->resize (n+1);
  rep->a[n]= append;
  return *this;
}

TMPL Vector&
Vector::operator << (const Vector& w) {
  typedef implementation<vector_linear,V> Vec;
  ASSERT (is_non_scalar (*this) && is_non_scalar (w),
          "non-scalar vectors expected");
  secure ();
  nat n= rep->n, p= N(w);
  rep->resize (n + p);
  Vec::copy (rep->a + n, seg (w), p);
  return *this;
}

TMPL void
inside_append (const Vector& v, const Vector& w) {
  typedef implementation<vector_linear,V> Vec;
  ASSERT (is_non_scalar (v) && is_non_scalar (w),
          "non-scalar vectors expected");
  nat n= N(v), p= N(w);
  inside (v)->resize (n + p);
  Vec::copy (const_cast<C*> (seg (v)) + n, seg (w), p);
}

TMPL Vector
reverse (const Vector& v) {
  typedef implementation<vector_linear,V> Vec;
  if (is_a_scalar (v)) return v;
  nat n= N(v);
  nat l= aligned_size<C,V> (n);
  C* a= mmx_formatted_new<C> (l, CF(v));
  Vec::vec_reverse (a, seg (v), n);
  return Vector (a, n, l, CF(v));
}

TMPL int
find (const Vector& v, const C& x) {
  nat i, n=N(v);
  for (i=0; i<n; i++)
    if (v[i] == x) return (int) i;
  return (int) i;
}

TMPL bool
contains (const Vector& v, const C& x) {
  return find (v, x) < ((int) N(v));
}

TMPL Vector
insert (const Vector& v, const C& x) {
  if (contains (v, x)) return v;
  else return append (v, x);
}

TMPL Vector
insert (const Vector& v, const C& x, nat pos) {
  // insert x at 'pos' in v
  ASSERT (pos <= N(v), "position out of range");
  nat n= N(v) + 1, i; nat l= aligned_size<C,V> (n);
  C* a= mmx_formatted_new<C> (l, CF(v));
  const C* p= seg (v);
  for (i= 0; i < pos; i++, p++) a[i]= *p;
  a[pos]= x; i++;
  for (;i < n; i++, p++) a[i]= *p;
  return Vector (a, n, l, CF(v));
}

template<typename C,typename V,typename W> Vector
insert (const Vector& v, const Vector& x, const vector<nat,W>& pos) {
  // insert x[0] at pos[0] in v, then x[1] at pos[1], etc
  ASSERT (N(x) == N(pos), "wrong sizes");
  nat n= N(v) + N(x), i= 0; nat l= aligned_size<C,V> (n);  
  C* a= mmx_formatted_new<C> (l, CF(v));
  const C* p= seg (v);
  for (nat k= 0; k < N(pos); k++) {
    ASSERT (pos[k] < n, "position out of range");
    for (; i < pos[k]; i++, p++) a[i]= *p;
    a[pos[k]]= x[k]; i++;
  }
  for (;i < n; i++, p++) a[i]= *p;
  return Vector (a, n, l, CF(v));
}

/******************************************************************************
* Routines on dense vectors
******************************************************************************/

template<typename Op, typename C, typename V> nat
unary_hash (const vector<C,V>& v) {
  register nat i, h= 12345, n= N(v);
  for (i=0; i<n; i++)
    h= (h<<1) ^ (h<<5) ^ (h>>27) ^ Op::op (v[i]);
  return h;
}

template<typename Op, typename C, typename V>
vector<Unary_return_type(Op,C),V>
unary_map (const vector<C,V>& v) {
  typedef implementation<vector_linear,V> Vec;
  typedef Unary_return_type(Op,C) T;
  format<T> fm= unary_map<Op> (CF(v));
  if (is_a_scalar (v)) return vector<T,V> (Op::op (v.scalar()));
  nat n= N(v);
  nat l= aligned_size<T,V> (n);
  T* r= mmx_formatted_new<T> (l, fm);
  Vec::template vec_unary<Op> (r, seg (v), n);
  return vector<T,V> (r, n, l, fm);
}

template<typename Op, typename C1, typename C2, typename V>
vector<Binary_return_type(Op,C1,C2),V>
binary_map (const vector<C1,V>& v, const vector<C2,V>& w) {
  typedef implementation<vector_linear,V> Vec;
  typedef Binary_return_type(Op,C1,C2) T;
  format<T> fm= binary_map<Op> (CF(v), CF(w));
  if (is_a_scalar (v) || is_a_scalar (w)) {
    if (is_non_scalar (v)) return binary_map<Op> (v, extend (w, v));
    if (is_non_scalar (w)) return binary_map<Op> (extend (v, w), w);
    return vector<T,V> (Op::op (v.scalar(), w.scalar()));
  }
  nat n= N(v);
  ASSERT (N(w) == n, "lengths don't match");
  nat l= aligned_size<T,V> (n);
  T* r= mmx_formatted_new<T> (l, fm);
  Vec::template vec_binary<Op> (r, seg (v), seg (w), n);
  return vector<T,V> (r, n, l, fm);
}

template<typename Op, typename C, typename V, typename X>
vector<Binary_return_type(Op,C,X),V>
binary_map_scalar (const vector<C,V>& v, const X& x) {
  typedef implementation<vector_linear,V> Vec;
  typedef Binary_return_type(Op,C,X) T;
  format<T> fm= binary_map_scalar<C> (CF(v), x);
  if (is_a_scalar (v)) return vector<T,V> (Op::op (v.scalar(), x));
  nat n= N(v);
  nat l= aligned_size<T,V> (n);
  T* r= mmx_formatted_new<T> (l, fm);
  Vec::template vec_binary_scalar<Op> (r, seg (v), x, n);
  return vector<T,V> (r, n, l, fm);
}

template<typename Op, typename C, typename V> vector<C,V>&
nullary_set (vector<C,V>& v) {
  typedef implementation<vector_linear,V> Vec;
  nat n= N(v);
  Vec::template vec_nullary<Op> (seg (v), n);
  return v;
}

template<typename Op, typename T, typename C, typename V> vector<T,V>&
unary_set (vector<T,V>& v, const vector<C,V>& w) {
  typedef implementation<vector_linear,V> Vec;
  if (is_a_scalar (v) || is_a_scalar (w)) {
    if (is_non_scalar (v))
      return unary_set<Op> (v, extend (w, v));
    else if (is_non_scalar (w))
      v= extend (v, w);
    else {
      Op::set_op (v.scalar(), w.scalar());
      return v;
    }
  }
  nat n= N(v);
  ASSERT (N(w) == n, "lengths don't match");
  Vec::template vec_unary<Op> (seg (v), seg (w), n);
  return v;
}

template<typename Op, typename T, typename V, typename X> vector<T,V>&
unary_set_scalar (vector<T,V>& v, const X& x) {
  typedef implementation<vector_linear,V> Vec;
  if (is_a_scalar (v)) {
    Op::set_op (v.scalar(), x);
    return v;
  }
  nat n= N(v);
  Vec::template vec_unary_scalar<Op> (seg (v), x, n);
  return v;
}

template<typename Op, typename C1, typename C2, typename V> bool
binary_test (const vector<C1,V>& v, const vector<C2,V>& w) {
  typedef implementation<vector_linear,V> Vec;
  if (is_a_scalar (v) || is_a_scalar (w)) {
    if (is_non_scalar (v)) return binary_test<Op> (v, extend (w, v));
    if (is_non_scalar (w)) return binary_test<Op> (extend (v, w), w);
    return Op::op (v.scalar(), w.scalar());
  }
  nat n= N(v);
  if (N(w) != n) return false;
  return Vec::template vec_binary_test<Op> (seg (v), seg (w), n);
}

template<typename Op, typename C, typename V, typename X> bool
binary_test_scalar (const vector<C,V>& v, const X& c) {
  typedef implementation<vector_linear,V> Vec;
  if (is_a_scalar (v)) return Op::op (v.scalar(), c);
  return Vec::template vec_binary_test_scalar<Op> (seg (v), c, N(v));
}

template<typename Op, typename C, typename V> Unary_return_type(Op,C)
big (const vector<C,V>& v) {
  typedef implementation<vector_linear,V> Vec;
  ASSERT (is_non_scalar (v), "non scalar vector expected");
  return Vec::template vec_unary_big<Op> (seg (v), N(v), CF(v));
}

template<typename Op, typename C, typename V> Unary_return_type(Op,C)
big_dicho (const vector<C,V>& v) {
  typedef implementation<vector_linear,V> Vec;
  ASSERT (is_non_scalar (v), "non scalar vector expected");
  return Vec::template vec_unary_big_dicho<Op> (seg (v), N(v), CF(v));
}

/******************************************************************************
* Dynamic mappers
******************************************************************************/

template<typename S1, typename D> vector<D>
map (const function_1<D,Argument(S1) >& fun, const vector<S1>& v1) {
  format<D> fm= map (fun, CF(v1));
  ASSERT (is_non_scalar (v1), "non-scalar vector expected");
  nat n= N(v1);
  nat l= default_aligned_size<D> (n);
  D* r= mmx_formatted_new<D> (l, fm);
  for (nat i=0; i<n; i++) r[i]= fun (v1[i]);
  return vector<D> (r, n, l, fm);
}

template<typename S1, typename D> vector<D>
map (D (*fun) (const S1&), const vector<S1>& v1) {
  format<D> fm= map (fun, CF(v1));
  ASSERT (is_non_scalar (v1), "non-scalar vector expected");
  nat n= N(v1);
  nat l= default_aligned_size<D> (n);
  D* r= mmx_formatted_new<D> (l, fm);
  for (nat i=0; i<n; i++) r[i]= fun (v1[i]);
  return vector<D> (r, n, l, fm);
}

template<typename S1, typename D> vector<D>
map (const function_1<D,Argument(S1) >& fun,
     const vector<S1>& v1, const format<D>& fm)
{
  ASSERT (is_non_scalar (v1), "non-scalar vector expected");
  nat n= N(v1);
  nat l= default_aligned_size<D> (n);
  D* r= mmx_formatted_new<D> (l, fm);
  for (nat i=0; i<n; i++) r[i]= fun (v1[i]);
  return vector<D> (r, n, l, fm);
}

template<typename S1, typename D> vector<D>
map (D (*fun) (const S1&),
     const vector<S1>& v1, const format<D>& fm) {
  ASSERT (is_non_scalar (v1), "non-scalar vector expected");
  nat n= N(v1);
  nat l= default_aligned_size<D> (n);
  D* r= mmx_formatted_new<D> (l, fm);
  for (nat i=0; i<n; i++) r[i]= fun (v1[i]);
  return vector<D> (r, n, l, fm);
}

template<typename S1, typename S2, typename D> vector<D>
map (const function_2<D,Argument(S1),Argument(S2) >& fun,
     const vector<S1>& v1, const vector<S2>& v2)
{
  format<D> fm; //= map (fun, CF(v1), CF(v2));
  ASSERT (is_non_scalar (v1), "non-scalar vector expected");
  ASSERT (is_non_scalar (v2), "non-scalar vector expected");
  ASSERT (N(v1) == N(v2), "lengths don't match");
  nat n= N(v1);
  nat l= default_aligned_size<D> (n);
  D* r= mmx_formatted_new<D> (l, fm);
  for (nat i=0; i<n; i++) r[i]= fun (v1[i], v2[i]);
  return vector<D> (r, n, l, fm);
}

template<typename S1, typename S2, typename S3, typename D> vector<D>
map (const function_3<D,Argument(S1),Argument(S2),Argument(S3) >& fun,
     const vector<S1>& v1, const vector<S2>& v2, const vector<S3>& v3)
{
  format<D> fm; //= map (fun, CF(v1), CF(v2), CF(v3));
  ASSERT (is_non_scalar (v1), "non-scalar vector expected");
  ASSERT (is_non_scalar (v2), "non-scalar vector expected");
  ASSERT (is_non_scalar (v3), "non-scalar vector expected");
  ASSERT (N(v1) == N(v2) && N(v2) == N(v3), "lengths don't match");
  nat n= N(v1);
  nat l= default_aligned_size<D> (n);
  D* r= mmx_formatted_new<D> (l, fm);
  for (nat i=0; i<n; i++) r[i]= fun (v1[i], v2[i], v3[i]);
  return vector<D> (r, n, l, fm);
}

/******************************************************************************
* Constants
******************************************************************************/

TMPL void set_pi (Vector& v) { nullary_set<pi_as_op> (v); }
TMPL void set_log2 (Vector& v) { nullary_set<log2_as_op> (v); }
TMPL void set_euler (Vector& v) { nullary_set<euler_as_op> (v); }
TMPL void set_catalan (Vector& v) { nullary_set<catalan_as_op> (v); }
TMPL void set_imaginary (Vector& v) { nullary_set<imaginary_as_op> (v); }
TMPL void set_nan (Vector& v) { nullary_set<nan_as_op> (v); }
TMPL void set_fuzz (Vector& v) { nullary_set<fuzz_as_op> (v); }
TMPL void set_smallest (Vector& v) { nullary_set<smallest_as_op> (v); }
TMPL void set_largest (Vector& v) { nullary_set<largest_as_op> (v); }
TMPL void set_accuracy (Vector& v) { nullary_set<accuracy_as_op> (v); }
TMPL void set_infinity (Vector& v) { nullary_set<infinity_as_op> (v); }
TMPL void set_maximal (Vector& v) { nullary_set<maximal_as_op> (v); }
TMPL void set_minimal (Vector& v) { nullary_set<minimal_as_op> (v); }

/******************************************************************************
* Instantiations of the abstract routines
******************************************************************************/

TMPL Vector copy (const Vector& v) {
  return unary_map<id_op> (v); }
TMPL Vector duplicate (const Vector& v) {
  return unary_map<duplicate_op> (v); }

TMPL Vector operator - (const Vector& v) {
  return unary_map<neg_op> (v); }
TMPL Vector operator + (const Vector& v, const Vector& w) {
  return binary_map<add_op> (v, w); }
TMPL Vector operator + (const C& v, const Vector& w) {
  return binary_map<add_op> (Vector (v), w); }
template<typename V> Vector_int 
operator + (const int& v, const Vector_int& w) {
  return binary_map<add_op> (Vector_int (v), w); }
TMPL Vector operator + (const Vector& w, const C& v) {
  return binary_map<add_op> (w, Vector (v)); }
template<typename V> Vector_int
operator + (const Vector_int& w, const int& v) {
  return binary_map<add_op> (w, Vector_int (v)); }
TMPL Vector operator - (const Vector& v, const Vector& w) {
  return binary_map<sub_op> (v, w); }
TMPL Vector operator - (const Vector& v, const C& w) {
  return binary_map<sub_op> (v, Vector (w)); }
template<typename V> Vector_int
operator - (const Vector_int& v, const int& w) {
  return binary_map<sub_op> (v, Vector_int (w)); }
TMPL Vector operator - (const C& v, const Vector& w) {
  return binary_map<sub_op> (Vector (v), w); }
template<typename V> Vector_int
operator - (const int& v, const Vector_int& w) {
  return binary_map<sub_op> (Vector_int (v), w); }
TMPL Vector operator * (const Vector& v, const Vector& w) {
  return binary_map<mul_op> (v, w); }
TMPL Vector operator / (const Vector& v, const Vector& w) {
  return binary_map<div_op> (v, w); }
TMPL Vector quo (const Vector& v, const Vector& w) {
  return binary_map<quo_op> (v, w); }
TMPL Vector rem (const Vector& v, const Vector& w) {
  return binary_map<rem_op> (v, w); }
TMPL Vector operator * (const Vector& v, const C& c) {
  return binary_map_scalar<rmul_op> (v, c); }
template<typename V> Vector_int
operator * (const Vector_int& v, const int& c) {
  return binary_map_scalar<rmul_op> (v, c); }
TMPL Vector operator * (const C& c, const Vector& v) {
  return binary_map_scalar<lmul_op> (v, c); }
template<typename V> Vector_int
operator * (const int& c, const Vector_int& v) {
  return binary_map_scalar<lmul_op> (v, c); }
TMPL Vector operator / (const Vector& v, const C& c) {
  return binary_map_scalar<rdiv_op> (v, c); }
template<typename V> Vector_int
operator / (const Vector_int& v, const int& c) {
  return binary_map_scalar<rdiv_op> (v, c); }
TMPL Vector operator / (const C& c, const Vector& v) {
  return binary_map_scalar<ldiv_op> (v, c); }
template<typename V> Vector_int
operator / (const int& c, const Vector_int& v) {
  return binary_map_scalar<ldiv_op> (v, c); }
TMPL Vector quo (const Vector& v, const C& c) {
  return binary_map_scalar<rquo_op> (v, c); }
TMPL Vector rem (const Vector& v, const C& c) {
  return binary_map_scalar<rrem_op> (v, c); }
TMPL Vector operator | (const Vector& v, const Vector& w) {
  return binary_map<or_op> (v, w); }
TMPL Vector operator & (const Vector& v, const Vector& w) {
  return binary_map<and_op> (v, w); }

TMPL Vector& operator += (Vector& v, const Vector& w) {
  return unary_set<add_op> (v, w); }
TMPL Vector& operator -= (Vector& v, const Vector& w) {
  return unary_set<sub_op> (v, w); }
TMPL Vector& operator *= (Vector& v, const Vector& w) {
  return unary_set<mul_op> (v, w); }
TMPL Vector& operator /= (Vector& v, const Vector& w) {
  return unary_set<div_op> (v, w); }
TMPL Vector& operator *= (Vector& v, const C& c) {
  return unary_set_scalar<rmul_op> (v, c); }
TMPL Vector& operator /= (Vector& v, const C& c) {
  return unary_set_scalar<rdiv_op> (v, c); }

TMPL Vector inf (const Vector& v, const Vector& w) {
  return binary_map<inf_op> (v, w); }
TMPL Vector sup (const Vector& v, const Vector& w) {
  return binary_map<sup_op> (v, w); }
TMPL bool operator <= (const Vector& v, const Vector& w) {
  return binary_test<lesseq_op> (v, w); }
TMPL bool operator >= (const Vector& v, const Vector& w) {
  return binary_test<gtreq_op> (v, w); }
TMPL bool operator == (const Vector& v, const C& c) {
  return binary_test_scalar<equal_op> (v, c); }
template<typename V> bool
operator == (const Vector_int& v, const int& c) {
  return binary_test_scalar<equal_op> (v, c); }
TMPL bool operator <= (const Vector& v, const C& c) {
  return binary_test_scalar<lesseq_op> (v, c); }
template<typename V> bool
operator <= (const Vector_int& v, const int& c) {
  return binary_test_scalar<lesseq_op> (v, c); }
TMPL bool operator >= (const Vector& v, const C& c) {
  return binary_test_scalar<gtreq_op> (v, c); }
template<typename V> bool
operator >= (const Vector_int& v, const int& c) {
  return binary_test_scalar<gtreq_op> (v, c); }

TRUE_IDENTITY_OP_SUGAR(TMPL,Vector)
EXACT_IDENTITY_OP_SUGAR(TMPL,Vector)
STRICT_COMPARE_SUGAR(TMPL,Vector)
ARITH_SCALAR_INT_SUGAR(TMPL,Vector);
EQUAL_INT_SUGAR(TMPL,Vector);
COMPARE_INT_SUGAR(TMPL,Vector);

TMPL Vector sqrt (const Vector& v) { return unary_map<sqrt_op> (v); }
TMPL Vector exp (const Vector& v) { return unary_map<exp_op> (v); }
TMPL Vector log (const Vector& v) { return unary_map<log_op> (v); }
TMPL Vector cos (const Vector& v) { return unary_map<cos_op> (v); }
TMPL Vector sin (const Vector& v) { return unary_map<sin_op> (v); }
TMPL Vector tan (const Vector& v) { return unary_map<tan_op> (v); }
TMPL Vector acos (const Vector& v) { return unary_map<acos_op> (v); }
TMPL Vector asin (const Vector& v) { return unary_map<asin_op> (v); }
TMPL Vector atan (const Vector& v) { return unary_map<atan_op> (v); }

TMPL Vector pow (const Vector& v, const Vector& w) {
  return binary_map<pow_op> (v, w); }
TMPL Vector pow (const Vector& v, const int& w) {
  return binary_map_scalar<rpow_op> (v, w); }
TMPL Vector pow (const int& w, const Vector& v) {
  return binary_map_scalar<lpow_op> (v, w); }

TMPL Vector derive (const Vector& v) {
  return unary_map<derive_op> (v); }
TMPL Vector integrate (const Vector& v) {
  return unary_map<integrate_op> (v); }
template<typename C,typename V,typename X> Vector
derive (const Vector& v, const X& x) {
  return binary_map_scalar<derive_op> (v, x); }
template<typename C,typename V,typename X> Vector
integrate (const Vector& v, const X& x) {
  return binary_map_scalar<integrate_op> (v, x); }

TMPL inline C big_add (const Vector& v) { return big<add_op> (v); }
TMPL inline C big_mul (const Vector& v) { return big<mul_op> (v); }
TMPL inline C big_or  (const Vector& v) { return big< or_op> (v); }
TMPL inline C big_and (const Vector& v) { return big<and_op> (v); }
TMPL inline C big_min (const Vector& v) { return big<min_op> (v); }
TMPL inline C big_max (const Vector& v) { return big<max_op> (v); }
TMPL inline C big_inf (const Vector& v) { return big<inf_op> (v); }
TMPL inline C big_sup (const Vector& v) { return big<sup_op> (v); }

/******************************************************************************
* Floating point related functions
******************************************************************************/

TMPL inline bool is_finite (const Vector& v) {
  if (is_a_scalar (v)) return is_finite (v.scalar());
  return big<and_is_finite_op> (v); }
TMPL inline bool is_nan (const Vector& v) {
  if (is_a_scalar (v)) return is_nan (v.scalar());
  return big<or_is_nan_op> (v); }
TMPL inline bool is_infinite (const Vector& v) {
  if (is_a_scalar (v)) return is_infinite (v.scalar());
  return !is_nan (v) && big<or_is_infinite_op> (v); }
TMPL inline bool is_fuzz (const Vector& v) {
  if (is_a_scalar (v)) return is_fuzz (v.scalar());
  return !is_nan (v) && big<or_is_fuzz_op> (v); }
TMPL inline bool is_reliable (const Vector& v) {
  if (is_a_scalar (v)) return is_reliable (v.scalar());
  return is_reliable (C (0)); }

TMPL inline Vector change_precision (const Vector& v, xnat p) {
  return binary_map_scalar<change_precision_op> (v, p); }
TMPL inline xnat precision (const Vector& v) {
  return big<min_precision_op> (v); }

TMPL inline xint exponent (const Vector& v) {
  return big<max_exponent_op> (v); }
TMPL inline double magnitude (const Vector& v) {
  return big<max_magnitude_op> (v); }

/******************************************************************************
* Complex and ball vectors
******************************************************************************/

TMPL Abs_vector abs (const Vector& v) { return unary_map<abs_op> (v); }
TMPL Real_vector Re (const Vector& v) { return unary_map<Re_op> (v); }
TMPL Real_vector Im (const Vector& v) { return unary_map<Im_op> (v); }
TMPL Vector conj (const Vector& v) { return unary_map<conj_op> (v); }

TMPL Center_vector center (const Vector& v) {
  return unary_map<center_op> (v); }
TMPL Radius_vector radius (const Vector& v) {
  return unary_map<radius_op> (v); }
TMPL Center_vector lower (const Vector& v) {
  return unary_map<lower_op> (v); }
TMPL Center_vector upper (const Vector& v) {
  return unary_map<upper_op> (v); }
TMPL inline Vector sharpen (const Vector& v) {
  return unary_map<sharpen_op> (v); }
TMPLK inline Vector blur (const Vector& v, const K& x) {
  return binary_map_scalar<blur_op> (v, x); }
template<typename C, typename V, typename C2, typename V2> inline Vector
blur (const Vector& v, const vector<C2,V2>& w) {
  return binary_map<blur_op> (v, w); }

/******************************************************************************
* Implementation of some routines on operators
******************************************************************************/

template<typename C> inline C
trig_op::diff_op (const C& me, const C& x) {
  return from_vector (vec (-derive (access (x, 0)) * access (me, 1),
                           derive (access (x, 0)) * access (me, 0)));
}

template<typename C> inline C
tan_op::def (const C& me, const C& f) {
  (void) me; const vector<C> v= cos_sin (f); return v[1] / v[0];
}

template<typename C> inline C
solve_vector_lde_op::diff_op (const C& me, const C& x) {
  typedef As_vector_type (C) V;
  typedef Scalar_type (V) K;
  const V w= as_vector (me);
  const V v= as_vector (x);
  nat i, n= N(v);
  K sum (0);
  V d (K (0), n);
  for (i=0; i<n; i++) {
    if (i==0) sum= v[i] * w[i];
    else sum= sum + v[i] * w[i];
    if (i<n-1) d[i]= w[i+1];
    else d[i]= sum;
  }
  return from_vector (d);
}

/******************************************************************************
* Further operations
******************************************************************************/

TMPL C dot (const Vector& v, const Vector& w) {
  typedef implementation<vector_linear,V> Vec;
  ASSERT (is_non_scalar (v) && is_non_scalar (w),
          "non-scalar vectors expected");
  ASSERT (N(v) == N(w), "lengths don't match");
  return Vec::inn_prod (seg (v), seg (w), N(v), CF (v)); }

template<typename C> vector<C> cos_sin (const C& x) {
  return vec (cos (x), sin (x)); }
template<typename C> vector<C> trig (const vector<C>& v) {
  return vec (cos (v[0]), sin (v[0])); }

#undef TMPL_DEF
#undef TMPL
#undef TMPLK
#undef Format
#undef Vector
#undef Vector_int
#undef Vector_rep
#undef Abs_vector
#undef Real_vector
#undef Center_vector
#undef Radius_vector
} // namespace mmx
#endif // __MMX_VECTOR_HPP