00001 
00002 
00003 
00004 
00005 
00006 
00007 
00008 
00009 
00010 
00011 
00012 
00013 #ifndef __MMX_VECTOR_HPP
00014 #define __MMX_VECTOR_HPP
00015 #include <basix/iterator.hpp>
00016 #include <basix/vector_naive.hpp>
00017 
00019 
00020 namespace mmx {
00021 #define Vector_variant(C) vector_variant_helper<C>::VV
00022 #define TMPL_DEF template<typename C, typename V>
00023 
00024 #define TMPL template<typename C, typename V>
00025 #define TMPLK template<typename C, typename V, typename K>
00026 #define Format format<C>
00027 #define Vector vector<C,V>
00028 #define Vector_int vector<int,V>
00029 #define Vector_rep vector_rep<C,V>
00030 #define Abs_vector vector<Abs_type(C),V>
00031 #define Real_vector vector<Real_type(C),V>
00032 #define Center_vector vector<Center_type(C),V>
00033 #define Radius_vector vector<Radius_type(C),V>
00034 #define Lifted_vector vector<Lift_type(C)>
00035 #define Projected_vector vector<Project_type(C)>
00036 #define Reconstructed_vector vector<Reconstruct_type(C)>
00037 #define Truncated_vector vector<Truncate_type(C)>
00038 #define Completed_vector vector<Complete_type(C)>
00039 #define Evaluated_vector vector<Evaluate_type(C,K)>
00040 TMPL class vector;
00041 TMPL inline nat N (const Vector& v);
00042 TMPL inline C* seg (Vector& v);
00043 TMPL inline const C* seg (const Vector& v);
00044 TMPL inline C* seg_inside (const Vector& v);
00045 TMPL inline bool is_a_scalar   (const Vector& v);
00046 TMPL inline bool is_non_scalar (const Vector& v);
00047 
00048 
00049 
00050 
00051 
00052 TMPL_DEF
00053 class vector_rep REP_STRUCT_1(C) {
00054 public:
00055   typedef implementation<vector_linear,V> Vec;
00056 private:
00057   C*   a;           
00058   nat  n;           
00059   nat  l;           
00060   bool scalar_flag; 
00061 public:
00062   inline vector_rep (C* a2, nat n2, bool flag, const Format& fm):
00063     Format (fm), a (a2), n (n2), l (n2), scalar_flag (flag) {}
00064   inline vector_rep (C* a2, nat n2, nat l2, bool flag, const Format& fm):
00065     Format (fm), a (a2), n (n2), l (l2), scalar_flag (flag) {}
00066   inline ~vector_rep () { mmx_delete<C> (a, l); }
00067   void resize (nat n2) {
00068     nat l2;
00069     if (n2 > l) l2= aligned_size<C,V> (max (n2, l<<1));
00070     else if (n2 < (l >> 1)) l2= aligned_size<C,V> (n2);
00071     else { n= n2; return; }
00072     C* b= mmx_new<C> (l2);
00073     Vec::copy (b, a, min (n, n2));
00074     mmx_delete<C> (a, l);
00075     a= b;
00076     n= n2;
00077     l= l2;
00078   }
00079   friend class Vector;
00080   friend nat N LESSGTR (const Vector& v);
00081   friend C* seg LESSGTR (Vector& v);
00082   friend const C* seg LESSGTR (const Vector& v);
00083   friend C* seg_inside LESSGTR (const Vector& v);
00084   friend bool is_a_scalar   LESSGTR (const Vector& v);
00085   friend bool is_non_scalar LESSGTR (const Vector& v);
00086 };
00087 
00088 TMPL_DEF
00089 class vector {
00090 INDIRECT_PROTO_2 (vector, vector_rep, C, V)
00091 public:
00092   typedef implementation<vector_linear,V> Vec;
00093 public:
00094   vector () {
00095     nat n= Vec::def_len;
00096     nat l= aligned_size<C,V> (n);
00097     C* a= mmx_formatted_new<C> (l, Format (no_format ()));
00098     rep= new Vector_rep (a, n, l, false, Format (no_format ()));
00099   }
00100   vector (const Format& fm) {
00101     nat n= Vec::def_len;
00102     nat l= aligned_size<C,V> (n);
00103     C* a= mmx_formatted_new<C> (l, fm);
00104     rep= new Vector_rep (a, n, l, false, fm);
00105   }
00106   template<typename T> 
00107   vector (const T& c) {
00108     Format fm= as<Format > (get_format (c));
00109     nat n= Vec::init_len;
00110     nat l= aligned_size<C,V> (n);
00111     C* a= mmx_formatted_new<C> (l, fm);
00112     Vec::set (a, as<C> (c), n);
00113     rep= new Vector_rep (a, n, l, true, fm);
00114   }
00115   template<typename T>
00116   vector (const T& c, const Format& fm) {
00117     nat n= Vec::init_len;
00118     nat l= aligned_size<C,V> (n);
00119     C* a= mmx_formatted_new<C> (l, fm);
00120     Vec::set_as (a, c, n);
00121     rep= new Vector_rep (a, n, l, true, fm);
00122   }
00123   vector (const C& x) {
00124     Format fm= as<Format > (get_format (x));
00125     nat n= Vec::init_len;
00126     nat l= aligned_size<C,V> (n);
00127     C* a= mmx_formatted_new<C> (l, fm);
00128     Vec::set (a, x, n);
00129     rep= new Vector_rep (a, n, l, true, fm);
00130   }
00131   template<typename T, typename W>
00132   vector (const vector<T,W>& v) {
00133     Format fm= as<Format > (CF(v));
00134     nat l= aligned_size<C,V> (N(v));
00135     C* a= mmx_formatted_new<C> (l, fm);
00136     Vec::cast (a, seg (v), N(v));
00137     rep= new Vector_rep (a, N(v), l, is_a_scalar (v), fm);
00138   }
00139   template<typename T, typename W>
00140   vector (const vector<T,W>& v, const Format& fm) {
00141     nat l= aligned_size<C,V> (N(v));
00142     C* a= mmx_formatted_new<C> (l, fm);
00143     Vec::cast (a, seg (v), N(v));
00144     rep= new Vector_rep (a, N(v), l, is_a_scalar (v), fm);
00145   }
00146   template<typename T> 
00147   vector (const T& c, nat n) {
00148     Format fm= as<Format > (get_format (c));
00149     nat l= aligned_size<C,V> (n);
00150     C* a= mmx_formatted_new<C> (l, fm);
00151     Vec::set (a, as<C> (c), n);
00152     rep= new Vector_rep (a, n, l, true, fm);
00153   }
00154   vector (const C& x, nat n) {
00155     nat l= aligned_size<C,V> (n);
00156     C* a= mmx_formatted_new<C> (l, get_format (x));
00157     Vec::set (a, x, n);
00158     rep= new Vector_rep (a, n, l, false, get_format (x));
00159   }
00160   inline vector (C* a, nat n, const Format& fm) {
00161     rep= new Vector_rep (a, n, false, fm);
00162   }
00163   inline vector (C* a, nat n, nat l, const Format& fm) {
00164     rep= new Vector_rep (a, n, l, false, fm);
00165   }
00166   inline vector (C* a, nat n, nat l, bool flag, const Format& fm) {
00167     rep= new Vector_rep (a, n, l, flag, fm);
00168   }
00169   vector (const iterator<C>& it) {
00170     Format fm= CF(it);
00171     nat i, l= Vec::def_len;
00172     if (l==0) l= aligned_size<C,V> ((nat) 1);
00173     C* a= mmx_formatted_new<C> (l, fm);
00174     for (i=0; busy (it); i++, ++it) {
00175       if (i >= l) {
00176         C* b= mmx_formatted_new<C> (aligned_size<C,V> (l<<1), fm);
00177         Vec::copy (b, a, l);
00178         mmx_delete (a, l);
00179         a= b;
00180         l= l<<1;
00181       }
00182       a[i]= *it;
00183     }
00184     rep= new Vector_rep (a, i, l, false, fm);
00185   }
00186 
00187   inline const C& operator [] (nat i) const {
00188     VERIFY (is_non_scalar (*this), "non-scalar vector expected");
00189     VERIFY (i<N(*this), "index out of range");
00190     return rep->a[i]; }
00191   inline C& operator [] (nat i) {
00192     VERIFY (is_non_scalar (*this), "non-scalar vector expected");
00193     VERIFY (i<N(*this), "index out of range");
00194     secure(); return rep->a[i]; }
00195   inline C scalar () const {
00196     VERIFY (is_a_scalar (*this), "scalar vector expected");
00197     return rep->a[0]; }
00198   inline C& scalar () {
00199     VERIFY (is_a_scalar (*this), "scalar vector expected");
00200     return rep->a[0]; }
00201   inline nat size() const {
00202     return rep->n; }
00203 
00204   Vector& operator << (const C& append);
00205   Vector& operator << (const Vector& w);
00206 };
00207 INDIRECT_IMPL_2 (vector, vector_rep, typename C, C, typename V, V)
00208 DEFINE_UNARY_FORMAT_2 (vector)
00209 
00210 TMPL inline nat N (const Vector& v) { return v->n; }
00211 TMPL inline Format CF (const Vector& v) { return v->tfm (); }
00212 TMPL inline bool is_a_scalar (const Vector& v) { return v->scalar_flag; }
00213 TMPL inline bool is_non_scalar (const Vector& v) { return !v->scalar_flag; }
00214 TMPL inline const C& read (const Vector& v, nat i) { return v[i]; }
00215 TMPL inline const C* seg (const Vector& v) { return v->a; }
00216 TMPL inline C* seg (Vector& v) {
00217   VERIFY (is_non_scalar (v), "non-scalar vector expected");
00218   v.secure(); return v->a; }
00219 TMPL inline C* seg_inside (const Vector& v) { return inside (v)->a; }
00220 template<typename C, typename V, typename C2, typename V2> inline Vector
00221 extend (const Vector& v, const vector<C2,V2>& w) {
00222   VERIFY (is_a_scalar (v), "scalar vector expected");
00223   VERIFY (is_non_scalar (w), "non-scalar vector expected");
00224   return Vector (v.scalar(), N(w)); }
00225 TMPL inline bool is_nil (const Vector& v) {
00226   return N(v) == 0; }
00227 TMPL inline bool is_atom (const Vector& v) {
00228   return is_non_scalar (v) && N(v) == 1; }
00229 
00230 STYPE_TO_TYPE(TMPL,scalar_type,Vector,C);
00231 UNARY_RETURN_TYPE(TMPL,abs_op,Vector,Abs_vector);
00232 UNARY_RETURN_TYPE(TMPL,Re_op,Vector,Real_vector);
00233 UNARY_RETURN_TYPE(TMPL,center_op,Vector,Center_vector);
00234 UNARY_RETURN_TYPE(TMPL,radius_op,Vector,Radius_vector);
00235 UNARY_RETURN_TYPE(TMPL,lift_op,Vector,Lifted_vector);
00236 UNARY_RETURN_TYPE(TMPL,project_op,Vector,Projected_vector);
00237 UNARY_RETURN_TYPE(TMPL,reconstruct_op,Vector,Reconstructed_vector);
00238 BINARY_RETURN_TYPE(TMPL,truncate_op,Vector,nat,Truncated_vector);
00239 UNARY_RETURN_TYPE(TMPL,complete_op,Vector,Completed_vector);
00240 BINARY_RETURN_TYPE(TMPLK,evaluate_op,Vector,K,Evaluated_vector);
00241 
00242 TMPL struct rounding_helper<Vector >: public rounding_helper<C> {};
00243 TMPL struct species_type_information<Vector > {
00244   static const nat id= SPECIES_VECTOR; };
00245 
00246 TMPL
00247 struct fast_helper<Vector > {
00248   typedef Fast_type(C) FC;
00249   typedef vector<FC> fast_type;
00250   typedef typename Vector_variant(FC) FV;
00251   static inline fast_type dd (const Vector& v) {
00252     format<FC> fm= fast (CF(v));
00253     nat l= aligned_size<FC, FV> (N(v));
00254     FC* r= mmx_formatted_new<FC> (l, fm);
00255     for (nat i=0; i<N(v); i++) r[i]= fast (v[i]);
00256     return fast_type (r, N(v), l, is_a_scalar (v), fm); }
00257   static inline Vector uu (const fast_type& v) {
00258     Format fm= slow<C> (CF(v));
00259     nat l= aligned_size<C,V> ((N(v)));
00260     C* r= mmx_formatted_new<C> (l, fm);
00261     for (nat i=0; i<N(v); i++) r[i]= slow<C> (v[i]);
00262     return Vector (r, N(v), l, is_a_scalar (v), fm); }
00263 };
00264 
00265 template<typename T, typename F, typename TV, typename FV>
00266 struct as_helper<vector<T,TV>, vector<F,FV> > {
00267   static inline vector<T,TV>
00268   cv (const vector<F,FV>& v) {
00269     format<T> fm= as<format<T> > (CF(v));
00270     nat l= aligned_size<T,TV> (N(v));
00271     T* c= mmx_formatted_new<T> (l, fm);
00272     for (nat i=0; i<N(v); i++) c[i]= as<T> (v[i]);
00273     return vector<T,TV> (c, N(v), l, is_a_scalar (v), fm); }
00274 };
00275 
00276 template<typename T, typename F, typename TV, typename FV> inline void
00277 set_as (vector<T,TV>& r, const vector<F,FV>& v) {
00278   r= vector<T,TV> (v, CF(r));
00279 }
00280 
00281 template<typename C, typename V, typename T> inline void
00282 set_as (Vector& r, const T& x) {
00283   r= Vector (x, CF(r));
00284 }
00285 
00286 
00287 
00288 
00289 
00290 TMPL
00291 struct vector_as_helper<Vector > {
00292   static inline generic abstract (const Vector& x) {
00293     return as<generic> (as<vector<generic> > (x)); }
00294   static inline Vector concrete (const generic& x) {
00295     return as<Vector > (as<vector<generic> > (x)); }
00296 };
00297 
00298 inline vector<generic>
00299 abstract_vector (const generic& conc) {
00300   return as<vector<generic> > (conc->make_abstract_vector ());
00301 }
00302 
00303 inline generic
00304 concrete_vector (const vector<generic>& abst, const generic& conc) {
00305   return conc->make_concrete_vector (as<generic> (abst));
00306 }
00307 
00308 
00309 
00310 
00311 
00312 template<typename V, typename S>
00313 struct vector_fixed: public V {
00314   typedef typename V::Naive Naive;
00315   typedef vector_fixed<typename V::Aligned,S> Aligned;
00316   typedef vector_fixed<typename V::No_simd,S> No_simd;
00317   typedef vector_fixed<typename V::No_thread,S> No_thread;
00318 };
00319 
00320 template<typename V, typename W, typename S>
00321 struct implementation<vector_defaults,V,vector_fixed<W,S> > {
00322   static const nat def_len = S::val;
00323   static const nat init_len= S::val;
00324 }; 
00325 
00326 template<typename F, typename V, typename W, typename S>
00327 struct implementation<F,V,vector_fixed<W,S> >:
00328   public implementation<F,V,W> {};
00329 
00330 #define FVector     vector    <C,vector_fixed<V,S> >
00331 #define FVector_rep vector_rep<C,vector_fixed<V,S> >
00332 
00333 template<typename C, typename V, typename S>
00334 class FVector_rep REP_STRUCT_1(C) {
00335   C* a;
00336 public:
00337   
00338 
00339 
00340 
00341 
00342 
00343 
00344 
00345 
00346 
00347 
00348   inline vector_rep (C* a2, nat n, bool dummy, const Format& fm):
00349     Format (fm), a (a2) {
00350       VERIFY (n == S::val, "wrong size");
00351       (void) n; (void) dummy; }
00352   inline vector_rep (C* a2, nat n, nat l, bool dummy, const Format& fm):
00353     Format (fm), a (a2) {
00354       VERIFY (n == S::val, "wrong size");
00355       VERIFY ((l == aligned_size<C,V> (S::val)), "wrong allocated size");
00356       (void) n; (void) l; (void) dummy; }
00357   inline virtual ~vector_rep () {
00358     mmx_delete<C> (a, aligned_size<C,V> (S::val)); }
00359   friend class FVector;
00360   friend nat N LESSGTR (const FVector& v);
00361   friend C*  seg LESSGTR (FVector& v);
00362   friend const C* seg LESSGTR (const FVector& v);
00363   friend bool is_a_scalar   LESSGTR (const Vector& v);
00364   friend bool is_non_scalar LESSGTR (const Vector& v);
00365 };
00366 
00367 template<typename C, typename V, typename S> inline nat
00368 N (const FVector& v) { return S::val; }
00369 template<typename C, typename V, typename S> inline bool
00370 is_a_scalar (const FVector& v) { (void) v; return false; }
00371 template<typename C, typename V, typename S> inline bool
00372 is_non_scalar (const FVector& v) { (void) v; return true; }
00373 
00374 #undef FVector
00375 #undef FVector_rep
00376 
00377 
00378 
00379 
00380 
00381 template<typename C> vector<C>
00382 vec () {
00383   Format fm;
00384   nat l= default_aligned_size<C> ((nat) 0);
00385   C* a= mmx_formatted_new<C> (l, fm);
00386   return vector<C> (a, 0, l, fm);
00387 }
00388 
00389 template<typename C> vector<C>
00390 vec (const Format& fm) {
00391   nat l= default_aligned_size<C> ((nat) 0);
00392   C* a= mmx_formatted_new<C> (l, fm);
00393   return vector<C> (a, 0, l, fm);
00394 }
00395 
00396 template<typename C> vector<C>
00397 vec (const C& c1) {
00398   Format fm= get_format (c1);
00399   nat l= default_aligned_size<C> ((nat) 1);
00400   C* a= mmx_formatted_new<C> (l, fm);
00401   a[0]= c1;
00402   return vector<C> (a, 1, l, fm);
00403 }
00404 
00405 template<typename C> vector<C>
00406 vec (const C& c1, const C& c2) {
00407   Format fm= get_format (c1);
00408   nat l= default_aligned_size<C> ((nat) 2);
00409   C* a= mmx_formatted_new<C> (l, fm);
00410   a[0]= c1; a[1]= c2;
00411   return vector<C> (a, 2, l, fm);
00412 }
00413 
00414 template<typename C> vector<C>
00415 vec (const C& c1, const C& c2, const C& c3) {
00416   Format fm= get_format (c1);
00417   nat l= default_aligned_size<C> ((nat) 3);
00418   C* a= mmx_formatted_new<C> (l, fm);
00419   a[0]= c1; a[1]= c2; a[2]= c3;
00420   return vector<C> (a, 3, l, fm);
00421 }
00422 
00423 template<typename C> vector<C>
00424 vec (const C& c1, const C& c2, const C& c3, const C& c4) {
00425   Format fm= get_format (c1);
00426   nat l= default_aligned_size<C> ((nat) 4);
00427   C* a= mmx_formatted_new<C> (l, fm);
00428   a[0]= c1; a[1]= c2; a[2]= c3; a[3]= c4;
00429   return vector<C> (a, 4, l, fm);
00430 }
00431 
00432 template<typename C> vector<C>
00433 vec (const C& c1, const C& c2, const C& c3, const C& c4, const C& c5) {
00434   Format fm= get_format (c1);
00435   nat l= default_aligned_size<C> ((nat) 5);
00436   C* a= mmx_formatted_new<C> (l, fm);
00437   a[0]= c1; a[1]= c2; a[2]= c3; a[3]= c4; a[4]= c5;
00438   return vector<C> (a, 5, l, fm);
00439 }
00440 
00441 template<typename C> vector<C>
00442 vec (const C& c1, const C& c2, const C& c3,
00443      const C& c4, const C& c5, const C& c6) {
00444   Format fm= get_format (c1);
00445   nat l= default_aligned_size<C> ((nat) 6);
00446   C* a= mmx_formatted_new<C> (l, fm);
00447   a[0]= c1; a[1]= c2; a[2]= c3; a[3]= c4; a[4]= c5; a[5]= c6;
00448   return vector<C> (a, 6, l, fm);
00449 }
00450 
00451 template<typename C> vector<C>
00452 vec (const C& c1, const C& c2, const C& c3, const C& c4,
00453      const C& c5, const C& c6, const C& c7) {
00454   Format fm= get_format (c1);
00455   nat l= default_aligned_size<C> ((nat) 7);
00456   C* a= mmx_formatted_new<C> (l, fm);
00457   a[0]= c1; a[1]= c2; a[2]= c3; a[3]= c4; a[4]= c5; a[5]= c6; a[6]= c7;
00458   return vector<C> (a, 7, l, fm);
00459 }
00460 
00461 template<typename C> vector<C>
00462 vec (const C& c1, const C& c2, const C& c3, const C& c4,
00463      const C& c5, const C& c6, const C& c7, const C& c8) {
00464   Format fm= get_format (c1);
00465   nat l= default_aligned_size<C> ((nat) 8);
00466   C* a= mmx_formatted_new<C> (l, fm);
00467   a[0]= c1; a[1]= c2; a[2]= c3; a[3]= c4;
00468   a[4]= c5; a[5]= c6; a[6]= c7; a[7]= c8;
00469   return vector<C> (a, 8, l, fm);
00470 }
00471 
00472 template<typename C> inline vector<C>
00473 fill (nat n) {
00474   Format fm;
00475   nat l= default_aligned_size<C> (n);
00476   C* a= mmx_formatted_new<C> (l, fm);
00477   return vector<C> (a, n, l, fm);
00478 }
00479 
00480 template<typename C> inline vector<C>
00481 fill (nat n, const Format& fm) {
00482   nat l= default_aligned_size<C> (n);
00483   C* a= mmx_formatted_new<C> (l, fm);
00484   return vector<C> (a, n, l, fm);
00485 }
00486 
00487 template<typename C> inline vector<C>
00488 fill (const C& c, nat n) {
00489   Format fm= get_format (c);
00490   nat l= default_aligned_size<C> (n);
00491   C* a= mmx_formatted_new<C> (l, fm);
00492   for (nat i=0; i<n; i++) a[i]= c;
00493   return vector<C> (a, n, l, fm);
00494 }
00495 
00496 
00497 
00498 
00499 
00500 TMPL_DEF
00501 class vector_iterator_rep: public iterator_rep<C> {
00502   const Vector v;  
00503   nat i;           
00504 public:
00505   vector_iterator_rep (const Vector& v2, nat i2= 0):
00506     iterator_rep<C> (CF (v2)), v (v2), i (0) { (void) i2; }
00507 protected:
00508   bool is_busy () { return i < N(v); }
00509   void advance () { i++; }
00510   C current () { return v[i]; }
00511   iterator_rep<C>* clone () { return new vector_iterator_rep (v, i); }
00512 };
00513 
00514 TMPL iterator<C>
00515 iterate (const Vector& v) {
00516   return iterator<C> (new vector_iterator_rep<C,V> (v));
00517 }
00518 
00519 
00520 
00521 
00522 
00523 TMPL syntactic
00524 flatten (const Vector& v) {
00525   if (is_a_scalar (v)) return flatten (v.scalar());
00526   nat i, n= N(v);
00527   vector<syntactic> r= fill <syntactic> (n);
00528   for (i=0; i<n; i++)
00529     r[i]= flatten (v[i]);
00530   return apply (GEN_SQTUPLE, r);
00531 }
00532 
00533 TMPL
00534 struct binary_helper<Vector > {
00535   static inline string short_type_name () {
00536     return "V" * Short_type_name (C); }
00537   static inline generic full_type_name () {
00538     return gen ("Vector", Full_type_name (C)); }
00539   static nat size (const Vector& v) {
00540     if (is_a_scalar (v)) return 0;
00541     else return N(v); }
00542   static generic access (const Vector& v, nat i) {
00543     ASSERT (!is_a_scalar (v), "non-scalar vector expected");
00544     return as<generic> (v[i]); }
00545   static generic disassemble (const Vector& v) {
00546     if (is_a_scalar (v)) return as<generic> (v.scalar ());
00547     return as<generic> (as<vector<generic> > (v)); }
00548   static Vector assemble (const generic& v) {
00549     
00550     if (!is<vector<generic> > (v)) return Vector (as<C> (v));
00551     else return as<Vector > (as<vector<generic> > (v)); }
00552   static void write (const port& out, const Vector& v) {
00553     binary_write<Format > (out, CF (v));
00554     binary_write<bool> (out, is_a_scalar (v));
00555     if (is_a_scalar (v)) binary_write<C> (out, v.scalar ());
00556     else {
00557       binary_write<nat> (out, N(v));
00558       for (nat i=0; i<N(v); i++)
00559         binary_write<C> (out, v[i]); } }
00560   static Vector read (const port& in) {
00561     Format fm= binary_read<Format > (in);
00562     bool sc= binary_read<bool> (in);
00563     if (sc) return Vector (binary_read<C> (in));
00564     else {
00565       nat n= binary_read<nat> (in);
00566       nat l= aligned_size<C,V> (n);
00567       C*  r= mmx_new<C> (l);
00568       for (nat i=0; i<n; i++)
00569         r[i]= binary_read<C> (in);
00570       return Vector (r, n, l, fm); } }
00571 };
00572 
00573 
00574 
00575 
00576 
00577 TMPL Vector
00578 range (const Vector& v, nat start, nat end) {
00579   typedef implementation<vector_linear,V> Vec;
00580   VERIFY(end >= start, "bad range");
00581   nat n= end - start;
00582   if (is_a_scalar (v)) return Vector (v.scalar(), n);
00583   nat l= aligned_size<C,V> (n);
00584   C* a= mmx_formatted_new<C> (l, CF(v));
00585   Vec::copy (a, seg (v) + start, n);
00586   return Vector (a, n, l, CF(v));
00587 }
00588 
00589 TMPL Vector
00590 extract_mod (const Vector& v, nat k, nat p) {
00591   nat n= (N(v) - k + p - 1) / p;
00592   nat l= aligned_size<C,V> (n);
00593   C* r= mmx_formatted_new<C> (l, CF(v));
00594   for (nat i=0; i<n; i++) r[i]= v[i*p+k];
00595   return Vector (r, n, l, CF(v));
00596 }
00597 
00598 TMPL Vector
00599 extract (const Vector& v, vector<nat> e) {
00600   nat l= aligned_size<C,V> (N(e));
00601   C* r= mmx_formatted_new<C> (l, CF(v));
00602   for (nat i=0; i<N(e); i++)
00603     if (e[i] < N(v)) r[i]= v[e[i]];
00604   return Vector (r, N(e), l, CF(v));
00605 }
00606 
00607 TMPL Vector
00608 remove (const Vector& v, vector<nat> e) {
00609   nat n= N(v);
00610   vector<nat> f ((nat) 0, n);
00611   for (nat i= 0; i < n; i++) f[i]= i;
00612   for (nat i= 0; i < N(e); i++) if (e[i] < n) f[e[i]]= n;
00613   Vector r;
00614   for (nat i= 0; i < n; i++) {
00615     if (f[i] == n) continue;
00616     r << v[i];
00617   }
00618   return r;
00619 }
00620 
00621 TMPL inline C car (const Vector& v) { return v[0]; }
00622 TMPL inline Vector cdr (const Vector& v) { return range (v, 1, N(v)); }
00623 
00624 TMPL Vector
00625 append (const Vector& v, const C& c) {
00626   typedef implementation<vector_linear,V> Vec;
00627   ASSERT (is_non_scalar (v), "non-scalar vector expected");
00628   nat n= N(v);
00629   nat l= aligned_size<C,V> (n+1);
00630   C* a= mmx_formatted_new<C> (l, CF(v));
00631   a[n]= c;
00632   Vec::copy (a, seg (v), n);
00633   return Vector (a, n+1, l, CF(v));
00634 }
00635 
00636 TMPL Vector
00637 append (const Vector& v, const Vector& w) {
00638   typedef implementation<vector_linear,V> Vec;
00639   ASSERT (is_non_scalar (v) && is_non_scalar (w),
00640           "non-scalar vectors expected");
00641   nat n= N(v), p= N(w);
00642   nat l= aligned_size<C,V> (n+p);
00643   C* a= mmx_formatted_new<C> (l, CF(v));
00644   Vec::copy (a, seg (v), n);
00645   Vec::copy (a+n, seg (w), p);
00646   return Vector (a, n+p, l, CF(v));
00647 }
00648 
00649 TMPL Vector
00650 cons (const C& c, const Vector& v) {
00651   typedef implementation<vector_linear,V> Vec;
00652   ASSERT (is_non_scalar (v), "non-scalar vector expected");
00653   nat n= N(v);
00654   nat l= aligned_size<C,V> (n+1);
00655   C* a= mmx_formatted_new<C> (l, CF(v));
00656   a[0]= c;
00657   Vec::copy (a+1, seg (v), n);
00658   return Vector (a, n+1, l, CF(v));
00659 }
00660 
00661 TMPL Vector&
00662 Vector::operator << (const C& append) {
00663   ASSERT (is_non_scalar (*this), "non-scalar vector expected");
00664   secure ();
00665   nat n= rep->n;
00666   rep->resize (n+1);
00667   rep->a[n]= append;
00668   return *this;
00669 }
00670 
00671 TMPL Vector&
00672 Vector::operator << (const Vector& w) {
00673   typedef implementation<vector_linear,V> Vec;
00674   ASSERT (is_non_scalar (*this) && is_non_scalar (w),
00675           "non-scalar vectors expected");
00676   secure ();
00677   nat n= rep->n, p= N(w);
00678   rep->resize (n + p);
00679   Vec::copy (rep->a + n, seg (w), p);
00680   return *this;
00681 }
00682 
00683 TMPL void
00684 inside_append (const Vector& v, const Vector& w) {
00685   typedef implementation<vector_linear,V> Vec;
00686   ASSERT (is_non_scalar (v) && is_non_scalar (w),
00687           "non-scalar vectors expected");
00688   nat n= N(v), p= N(w);
00689   inside (v)->resize (n + p);
00690   Vec::copy (const_cast<C*> (seg (v)) + n, seg (w), p);
00691 }
00692 
00693 TMPL Vector
00694 reverse (const Vector& v) {
00695   typedef implementation<vector_linear,V> Vec;
00696   if (is_a_scalar (v)) return v;
00697   nat n= N(v);
00698   nat l= aligned_size<C,V> (n);
00699   C* a= mmx_formatted_new<C> (l, CF(v));
00700   Vec::vec_reverse (a, seg (v), n);
00701   return Vector (a, n, l, CF(v));
00702 }
00703 
00704 TMPL int
00705 find (const Vector& v, const C& x) {
00706   nat i, n=N(v);
00707   for (i=0; i<n; i++)
00708     if (v[i] == x) return (int) i;
00709   return (int) i;
00710 }
00711 
00712 TMPL bool
00713 contains (const Vector& v, const C& x) {
00714   return find (v, x) < ((int) N(v));
00715 }
00716 
00717 TMPL Vector
00718 insert (const Vector& v, const C& x) {
00719   if (contains (v, x)) return v;
00720   else return append (v, x);
00721 }
00722 
00723 TMPL Vector
00724 insert (const Vector& v, const C& x, nat pos) {
00725   
00726   ASSERT (pos <= N(v), "position out of range");
00727   nat n= N(v) + 1, i; nat l= aligned_size<C,V> (n);
00728   C* a= mmx_formatted_new<C> (l, CF(v));
00729   const C* p= seg (v);
00730   for (i= 0; i < pos; i++, p++) a[i]= *p;
00731   a[pos]= x; i++;
00732   for (;i < n; i++, p++) a[i]= *p;
00733   return Vector (a, n, l, CF(v));
00734 }
00735 
00736 template<typename C,typename V,typename W> Vector
00737 insert (const Vector& v, const Vector& x, const vector<nat,W>& pos) {
00738   
00739   ASSERT (N(x) == N(pos), "wrong sizes");
00740   nat n= N(v) + N(x), i= 0; nat l= aligned_size<C,V> (n);  
00741   C* a= mmx_formatted_new<C> (l, CF(v));
00742   const C* p= seg (v);
00743   for (nat k= 0; k < N(pos); k++) {
00744     ASSERT (pos[k] < n, "position out of range");
00745     for (; i < pos[k]; i++, p++) a[i]= *p;
00746     a[pos[k]]= x[k]; i++;
00747   }
00748   for (;i < n; i++, p++) a[i]= *p;
00749   return Vector (a, n, l, CF(v));
00750 }
00751 
00752 
00753 
00754 
00755 
00756 template<typename Op, typename C, typename V> nat
00757 unary_hash (const vector<C,V>& v) {
00758   register nat i, h= 12345, n= N(v);
00759   for (i=0; i<n; i++)
00760     h= (h<<1) ^ (h<<5) ^ (h>>27) ^ Op::op (v[i]);
00761   return h;
00762 }
00763 
00764 template<typename Op, typename C, typename V>
00765 vector<Unary_return_type(Op,C),V>
00766 unary_map (const vector<C,V>& v) {
00767   typedef implementation<vector_linear,V> Vec;
00768   typedef Unary_return_type(Op,C) T;
00769   format<T> fm= unary_map<Op> (CF(v));
00770   if (is_a_scalar (v)) return vector<T,V> (Op::op (v.scalar()));
00771   nat n= N(v);
00772   nat l= aligned_size<T,V> (n);
00773   T* r= mmx_formatted_new<T> (l, fm);
00774   Vec::template vec_unary<Op> (r, seg (v), n);
00775   return vector<T,V> (r, n, l, fm);
00776 }
00777 
00778 template<typename Op, typename C1, typename C2, typename V>
00779 vector<Binary_return_type(Op,C1,C2),V>
00780 binary_map (const vector<C1,V>& v, const vector<C2,V>& w) {
00781   typedef implementation<vector_linear,V> Vec;
00782   typedef Binary_return_type(Op,C1,C2) T;
00783   format<T> fm= binary_map<Op> (CF(v), CF(w));
00784   if (is_a_scalar (v) || is_a_scalar (w)) {
00785     if (is_non_scalar (v)) return binary_map<Op> (v, extend (w, v));
00786     if (is_non_scalar (w)) return binary_map<Op> (extend (v, w), w);
00787     return vector<T,V> (Op::op (v.scalar(), w.scalar()));
00788   }
00789   nat n= N(v);
00790   ASSERT (N(w) == n, "lengths don't match");
00791   nat l= aligned_size<T,V> (n);
00792   T* r= mmx_formatted_new<T> (l, fm);
00793   Vec::template vec_binary<Op> (r, seg (v), seg (w), n);
00794   return vector<T,V> (r, n, l, fm);
00795 }
00796 
00797 template<typename Op, typename C, typename V, typename X>
00798 vector<Binary_return_type(Op,C,X),V>
00799 binary_map_scalar (const vector<C,V>& v, const X& x) {
00800   typedef implementation<vector_linear,V> Vec;
00801   typedef Binary_return_type(Op,C,X) T;
00802   format<T> fm= binary_map_scalar<C> (CF(v), x);
00803   if (is_a_scalar (v)) return vector<T,V> (Op::op (v.scalar(), x));
00804   nat n= N(v);
00805   nat l= aligned_size<T,V> (n);
00806   T* r= mmx_formatted_new<T> (l, fm);
00807   Vec::template vec_binary_scalar<Op> (r, seg (v), x, n);
00808   return vector<T,V> (r, n, l, fm);
00809 }
00810 
00811 template<typename Op, typename C, typename V> vector<C,V>&
00812 nullary_set (vector<C,V>& v) {
00813   typedef implementation<vector_linear,V> Vec;
00814   nat n= N(v);
00815   Vec::template vec_nullary<Op> (seg (v), n);
00816   return v;
00817 }
00818 
00819 template<typename Op, typename T, typename C, typename V> vector<T,V>&
00820 unary_set (vector<T,V>& v, const vector<C,V>& w) {
00821   typedef implementation<vector_linear,V> Vec;
00822   if (is_a_scalar (v) || is_a_scalar (w)) {
00823     if (is_non_scalar (v))
00824       return unary_set<Op> (v, extend (w, v));
00825     else if (is_non_scalar (w))
00826       v= extend (v, w);
00827     else {
00828       Op::set_op (v.scalar(), w.scalar());
00829       return v;
00830     }
00831   }
00832   nat n= N(v);
00833   ASSERT (N(w) == n, "lengths don't match");
00834   Vec::template vec_unary<Op> (seg (v), seg (w), n);
00835   return v;
00836 }
00837 
00838 template<typename Op, typename T, typename V, typename X> vector<T,V>&
00839 unary_set_scalar (vector<T,V>& v, const X& x) {
00840   typedef implementation<vector_linear,V> Vec;
00841   if (is_a_scalar (v)) {
00842     Op::set_op (v.scalar(), x);
00843     return v;
00844   }
00845   nat n= N(v);
00846   Vec::template vec_unary_scalar<Op> (seg (v), x, n);
00847   return v;
00848 }
00849 
00850 template<typename Op, typename C1, typename C2, typename V> bool
00851 binary_test (const vector<C1,V>& v, const vector<C2,V>& w) {
00852   typedef implementation<vector_linear,V> Vec;
00853   if (is_a_scalar (v) || is_a_scalar (w)) {
00854     if (is_non_scalar (v)) return binary_test<Op> (v, extend (w, v));
00855     if (is_non_scalar (w)) return binary_test<Op> (extend (v, w), w);
00856     return Op::op (v.scalar(), w.scalar());
00857   }
00858   nat n= N(v);
00859   if (N(w) != n) return false;
00860   return Vec::template vec_binary_test<Op> (seg (v), seg (w), n);
00861 }
00862 
00863 template<typename Op, typename C, typename V, typename X> bool
00864 binary_test_scalar (const vector<C,V>& v, const X& c) {
00865   typedef implementation<vector_linear,V> Vec;
00866   if (is_a_scalar (v)) return Op::op (v.scalar(), c);
00867   return Vec::template vec_binary_test_scalar<Op> (seg (v), c, N(v));
00868 }
00869 
00870 template<typename Op, typename C, typename V> Unary_return_type(Op,C)
00871 big (const vector<C,V>& v) {
00872   typedef implementation<vector_linear,V> Vec;
00873   ASSERT (is_non_scalar (v), "non scalar vector expected");
00874   return Vec::template vec_unary_big<Op> (seg (v), N(v), CF(v));
00875 }
00876 
00877 template<typename Op, typename C, typename V> Unary_return_type(Op,C)
00878 big_dicho (const vector<C,V>& v) {
00879   typedef implementation<vector_linear,V> Vec;
00880   ASSERT (is_non_scalar (v), "non scalar vector expected");
00881   return Vec::template vec_unary_big_dicho<Op> (seg (v), N(v), CF(v));
00882 }
00883 
00884 
00885 
00886 
00887 
00888 template<typename S1, typename D> vector<D>
00889 map (const function_1<D,Argument(S1) >& fun, const vector<S1>& v1) {
00890   format<D> fm= map (fun, CF(v1));
00891   ASSERT (is_non_scalar (v1), "non-scalar vector expected");
00892   nat n= N(v1);
00893   nat l= default_aligned_size<D> (n);
00894   D* r= mmx_formatted_new<D> (l, fm);
00895   for (nat i=0; i<n; i++) r[i]= fun (v1[i]);
00896   return vector<D> (r, n, l, fm);
00897 }
00898 
00899 template<typename S1, typename D> vector<D>
00900 map (D (*fun) (const S1&), const vector<S1>& v1) {
00901   format<D> fm= map (fun, CF(v1));
00902   ASSERT (is_non_scalar (v1), "non-scalar vector expected");
00903   nat n= N(v1);
00904   nat l= default_aligned_size<D> (n);
00905   D* r= mmx_formatted_new<D> (l, fm);
00906   for (nat i=0; i<n; i++) r[i]= fun (v1[i]);
00907   return vector<D> (r, n, l, fm);
00908 }
00909 
00910 template<typename S1, typename D, typename Fun> vector<D>
00911 map (const Fun& fun, const vector<S1>& v1, const format<D>& fm) {
00912   ASSERT (is_non_scalar (v1), "non-scalar vector expected");
00913   nat n= N(v1);
00914   nat l= default_aligned_size<D> (n);
00915   D* r= mmx_formatted_new<D> (l, fm);
00916   for (nat i=0; i<n; i++) r[i]= fun (v1[i]);
00917   return vector<D> (r, n, l, fm);
00918 }
00919 
00920 template<typename S1, typename S2, typename D> vector<D>
00921 map (const function_2<D,Argument(S1),Argument(S2) >& fun,
00922      const vector<S1>& v1, const vector<S2>& v2)
00923 {
00924   format<D> fm; 
00925   ASSERT (is_non_scalar (v1), "non-scalar vector expected");
00926   ASSERT (is_non_scalar (v2), "non-scalar vector expected");
00927   ASSERT (N(v1) == N(v2), "lengths don't match");
00928   nat n= N(v1);
00929   nat l= default_aligned_size<D> (n);
00930   D* r= mmx_formatted_new<D> (l, fm);
00931   for (nat i=0; i<n; i++) r[i]= fun (v1[i], v2[i]);
00932   return vector<D> (r, n, l, fm);
00933 }
00934 
00935 template<typename S1, typename S2, typename D> vector<D>
00936 map (D (*fun) (const S1&, const S2&),
00937      const vector<S1>& v1, const vector<S2>& v2)
00938 {
00939   format<D> fm; 
00940   ASSERT (is_non_scalar (v1), "non-scalar vector expected");
00941   ASSERT (is_non_scalar (v2), "non-scalar vector expected");
00942   ASSERT (N(v1) == N(v2), "lengths don't match");
00943   nat n= N(v1);
00944   nat l= default_aligned_size<D> (n);
00945   D* r= mmx_formatted_new<D> (l, fm);
00946   for (nat i=0; i<n; i++) r[i]= fun (v1[i], v2[i]);
00947   return vector<D> (r, n, l, fm);
00948 }
00949 
00950 template<typename S1, typename S2, typename D, typename Fun> vector<D>
00951 map (const Fun& fun,
00952      const vector<S1>& v1, const vector<S2>& v2, const format<D>& fm)
00953 {
00954   ASSERT (is_non_scalar (v1), "non-scalar vector expected");
00955   ASSERT (is_non_scalar (v2), "non-scalar vector expected");
00956   ASSERT (N(v1) == N(v2), "lengths don't match");
00957   nat n= N(v1);
00958   nat l= default_aligned_size<D> (n);
00959   D* r= mmx_formatted_new<D> (l, fm);
00960   for (nat i=0; i<n; i++) r[i]= fun (v1[i], v2[i]);
00961   return vector<D> (r, n, l, fm);
00962 }
00963 
00964 template<typename S1, typename S2, typename S3, typename D> vector<D>
00965 map (const function_3<D,Argument(S1),Argument(S2),Argument(S3) >& fun,
00966      const vector<S1>& v1, const vector<S2>& v2, const vector<S3>& v3)
00967 {
00968   format<D> fm; 
00969   ASSERT (is_non_scalar (v1), "non-scalar vector expected");
00970   ASSERT (is_non_scalar (v2), "non-scalar vector expected");
00971   ASSERT (is_non_scalar (v3), "non-scalar vector expected");
00972   ASSERT (N(v1) == N(v2) && N(v2) == N(v3), "lengths don't match");
00973   nat n= N(v1);
00974   nat l= default_aligned_size<D> (n);
00975   D* r= mmx_formatted_new<D> (l, fm);
00976   for (nat i=0; i<n; i++) r[i]= fun (v1[i], v2[i], v3[i]);
00977   return vector<D> (r, n, l, fm);
00978 }
00979 
00980 template<typename S1, typename S2, typename S3,
00981          typename D, typename Fun> vector<D>
00982 map (const Fun& fun,
00983      const vector<S1>& v1, const vector<S2>& v2, const vector<S3>& v3,
00984      const format<D>& fm)
00985 {
00986   ASSERT (is_non_scalar (v1), "non-scalar vector expected");
00987   ASSERT (is_non_scalar (v2), "non-scalar vector expected");
00988   ASSERT (is_non_scalar (v3), "non-scalar vector expected");
00989   ASSERT (N(v1) == N(v2) && N(v2) == N(v3), "lengths don't match");
00990   nat n= N(v1);
00991   nat l= default_aligned_size<D> (n);
00992   D* r= mmx_formatted_new<D> (l, fm);
00993   for (nat i=0; i<n; i++) r[i]= fun (v1[i], v2[i], v3[i]);
00994   return vector<D> (r, n, l, fm);
00995 }
00996 
00997 template<typename S1, typename Fun> void
00998 foreach (const Fun& fun, const vector<S1>& v1) {
00999   ASSERT (is_non_scalar (v1), "non-scalar vector expected");
01000   for (nat i=0; i<N(v1); i++) fun (v1[i]);
01001 }
01002 
01003 template<typename S1, typename S2, typename Fun> void
01004 foreach (const Fun& fun, const vector<S1>& v1, const vector<S2>& v2) {
01005   ASSERT (is_non_scalar (v1) && is_non_scalar (v2),
01006           "non-scalar vectors expected");
01007   ASSERT (N(v1) == N(v2), "vectors of unequal lengths");
01008   for (nat i=0; i<N(v1); i++) fun (v1[i], v2[i]);
01009 }
01010 
01011 template<typename S1, typename S2, typename S3, typename Fun> void
01012 foreach (const Fun& fun,
01013          const vector<S1>& v1, const vector<S2>& v2, const vector<S3>& v3) {
01014   ASSERT (is_non_scalar (v1) && is_non_scalar (v2) && is_non_scalar (v3),
01015           "non-scalar vectors expected");
01016   ASSERT (N(v1) == N(v2) && N(v2) == N(v3), "vectors of unequal lengths");
01017   for (nat i=0; i<N(v1); i++) fun (v1[i], v2[i], v3[i]);
01018 }
01019 
01020 
01021 
01022 
01023 
01024 TMPL void set_default (Vector& v) { nullary_set<default_as_op> (v); }
01025 TMPL void set_pi (Vector& v) { nullary_set<pi_as_op> (v); }
01026 TMPL void set_log2 (Vector& v) { nullary_set<log2_as_op> (v); }
01027 TMPL void set_euler (Vector& v) { nullary_set<euler_as_op> (v); }
01028 TMPL void set_catalan (Vector& v) { nullary_set<catalan_as_op> (v); }
01029 TMPL void set_imaginary (Vector& v) { nullary_set<imaginary_as_op> (v); }
01030 TMPL void set_nan (Vector& v) { nullary_set<nan_as_op> (v); }
01031 TMPL void set_fuzz (Vector& v) { nullary_set<fuzz_as_op> (v); }
01032 TMPL void set_smallest (Vector& v) { nullary_set<smallest_as_op> (v); }
01033 TMPL void set_largest (Vector& v) { nullary_set<largest_as_op> (v); }
01034 TMPL void set_accuracy (Vector& v) { nullary_set<accuracy_as_op> (v); }
01035 TMPL void set_infinity (Vector& v) { nullary_set<infinity_as_op> (v); }
01036 TMPL void set_maximal (Vector& v) { nullary_set<maximal_as_op> (v); }
01037 TMPL void set_minimal (Vector& v) { nullary_set<minimal_as_op> (v); }
01038 
01039 
01040 
01041 
01042 
01043 TMPL Vector copy (const Vector& v) {
01044   return unary_map<id_op> (v); }
01045 TMPL Vector duplicate (const Vector& v) {
01046   return unary_map<duplicate_op> (v); }
01047 
01048 TMPL Vector operator - (const Vector& v) {
01049   return unary_map<neg_op> (v); }
01050 TMPL Vector operator + (const Vector& v, const Vector& w) {
01051   return binary_map<add_op> (v, w); }
01052 TMPL Vector operator + (const C& v, const Vector& w) {
01053   return binary_map<add_op> (Vector (v), w); }
01054 template<typename V> Vector_int 
01055 operator + (const int& v, const Vector_int& w) {
01056   return binary_map<add_op> (Vector_int (v), w); }
01057 TMPL Vector operator + (const Vector& w, const C& v) {
01058   return binary_map<add_op> (w, Vector (v)); }
01059 template<typename V> Vector_int
01060 operator + (const Vector_int& w, const int& v) {
01061   return binary_map<add_op> (w, Vector_int (v)); }
01062 TMPL Vector operator - (const Vector& v, const Vector& w) {
01063   return binary_map<sub_op> (v, w); }
01064 TMPL Vector operator - (const Vector& v, const C& w) {
01065   return binary_map<sub_op> (v, Vector (w)); }
01066 template<typename V> Vector_int
01067 operator - (const Vector_int& v, const int& w) {
01068   return binary_map<sub_op> (v, Vector_int (w)); }
01069 TMPL Vector operator - (const C& v, const Vector& w) {
01070   return binary_map<sub_op> (Vector (v), w); }
01071 template<typename V> Vector_int
01072 operator - (const int& v, const Vector_int& w) {
01073   return binary_map<sub_op> (Vector_int (v), w); }
01074 TMPL Vector operator * (const Vector& v, const Vector& w) {
01075   return binary_map<mul_op> (v, w); }
01076 TMPL Vector operator / (const Vector& v, const Vector& w) {
01077   return binary_map<div_op> (v, w); }
01078 TMPL Vector quo (const Vector& v, const Vector& w) {
01079   return binary_map<quo_op> (v, w); }
01080 TMPL Vector rem (const Vector& v, const Vector& w) {
01081   return binary_map<rem_op> (v, w); }
01082 TMPL Vector operator * (const Vector& v, const C& c) {
01083   return binary_map_scalar<rmul_op> (v, c); }
01084 template<typename V> Vector_int
01085 operator * (const Vector_int& v, const int& c) {
01086   return binary_map_scalar<rmul_op> (v, c); }
01087 TMPL Vector operator * (const C& c, const Vector& v) {
01088   return binary_map_scalar<lmul_op> (v, c); }
01089 template<typename V> Vector_int
01090 operator * (const int& c, const Vector_int& v) {
01091   return binary_map_scalar<lmul_op> (v, c); }
01092 TMPL Vector operator / (const Vector& v, const C& c) {
01093   return binary_map_scalar<rdiv_op> (v, c); }
01094 template<typename V> Vector_int
01095 operator / (const Vector_int& v, const int& c) {
01096   return binary_map_scalar<rdiv_op> (v, c); }
01097 TMPL Vector operator / (const C& c, const Vector& v) {
01098   return binary_map_scalar<ldiv_op> (v, c); }
01099 template<typename V> Vector_int
01100 operator / (const int& c, const Vector_int& v) {
01101   return binary_map_scalar<ldiv_op> (v, c); }
01102 TMPL Vector quo (const Vector& v, const C& c) {
01103   return binary_map_scalar<rquo_op> (v, c); }
01104 TMPL Vector rem (const Vector& v, const C& c) {
01105   return binary_map_scalar<rrem_op> (v, c); }
01106 TMPL Vector operator | (const Vector& v, const Vector& w) {
01107   return binary_map<or_op> (v, w); }
01108 TMPL Vector operator & (const Vector& v, const Vector& w) {
01109   return binary_map<and_op> (v, w); }
01110 
01111 TMPL Vector& operator += (Vector& v, const Vector& w) {
01112   return unary_set<add_op> (v, w); }
01113 TMPL Vector& operator -= (Vector& v, const Vector& w) {
01114   return unary_set<sub_op> (v, w); }
01115 TMPL Vector& operator *= (Vector& v, const Vector& w) {
01116   return unary_set<mul_op> (v, w); }
01117 TMPL Vector& operator /= (Vector& v, const Vector& w) {
01118   return unary_set<div_op> (v, w); }
01119 TMPL Vector& operator *= (Vector& v, const C& c) {
01120   return unary_set_scalar<rmul_op> (v, c); }
01121 TMPL Vector& operator /= (Vector& v, const C& c) {
01122   return unary_set_scalar<rdiv_op> (v, c); }
01123 
01124 TMPL Vector inf (const Vector& v, const Vector& w) {
01125   return binary_map<inf_op> (v, w); }
01126 TMPL Vector sup (const Vector& v, const Vector& w) {
01127   return binary_map<sup_op> (v, w); }
01128 TMPL bool operator <= (const Vector& v, const Vector& w) {
01129   return binary_test<lesseq_op> (v, w); }
01130 TMPL bool operator >= (const Vector& v, const Vector& w) {
01131   return binary_test<gtreq_op> (v, w); }
01132 TMPL bool operator == (const Vector& v, const C& c) {
01133   return binary_test_scalar<equal_op> (v, c); }
01134 TMPL bool operator != (const Vector& v, const C& c) {
01135   return !binary_test_scalar<equal_op> (v, c); }
01136 TMPL bool operator <= (const Vector& v, const C& c) {
01137   return binary_test_scalar<lesseq_op> (v, c); }
01138 TMPL inline bool operator < (const Vector& v, const C& c) {
01139   return v <= c && v != c; }
01140 TMPL bool operator >= (const Vector& v, const C& c) {
01141   return binary_test_scalar<gtreq_op> (v, c); }
01142 TMPL inline bool operator > (const Vector& v, const C& c) {
01143   return v >= c && v != c; }
01144 TMPL inline bool operator == (const C& c, const Vector& v) { return v == c; }
01145 TMPL inline bool operator != (const C& c, const Vector& v) { return v != c; }
01146 TMPL inline bool operator <  (const C& c, const Vector& v) { return v >  c; }
01147 TMPL inline bool operator >  (const C& c, const Vector& v) { return v <  c; }
01148 TMPL inline bool operator <= (const C& c, const Vector& v) { return v >= c; }
01149 TMPL inline bool operator >= (const C& c, const Vector& v) { return v <= c; }
01150 
01151 TRUE_IDENTITY_OP_SUGAR(TMPL,Vector)
01152 EXACT_IDENTITY_OP_SUGAR(TMPL,Vector)
01153 STRICT_COMPARE_SUGAR(TMPL,Vector)
01154 ARITH_SCALAR_INT_SUGAR(TMPL,Vector);
01155 EQUAL_INT_SUGAR(TMPL,Vector);
01156 COMPARE_INT_SUGAR(TMPL,Vector);
01157 
01158 template<typename V> bool
01159 operator == (const Vector_int& v, const int& c) {
01160   return binary_test_scalar<equal_op> (v, c); }
01161 template<typename V> bool
01162 operator != (const Vector_int& v, const int& c) {
01163   return !binary_test_scalar<equal_op> (v, c); }
01164 template<typename V> bool
01165 operator <= (const Vector_int& v, const int& c) {
01166   return binary_test_scalar<lesseq_op> (v, c); }
01167 template<typename V> inline bool
01168 operator < (const Vector_int& v, const int& c) {
01169   return v <= c && v != c; }
01170 template<typename V> bool
01171 operator >= (const Vector_int& v, const int& c) {
01172   return binary_test_scalar<gtreq_op> (v, c); }
01173 template<typename V> inline bool
01174 operator > (const Vector_int& v, const int& c) {
01175   return v >= c && v != c; }
01176 
01177 TMPL Vector sqrt (const Vector& v) { return unary_map<sqrt_op> (v); }
01178 TMPL Vector exp (const Vector& v) { return unary_map<exp_op> (v); }
01179 TMPL Vector log (const Vector& v) { return unary_map<log_op> (v); }
01180 TMPL Vector cos (const Vector& v) { return unary_map<cos_op> (v); }
01181 TMPL Vector sin (const Vector& v) { return unary_map<sin_op> (v); }
01182 TMPL Vector tan (const Vector& v) { return unary_map<tan_op> (v); }
01183 TMPL Vector acos (const Vector& v) { return unary_map<acos_op> (v); }
01184 TMPL Vector asin (const Vector& v) { return unary_map<asin_op> (v); }
01185 TMPL Vector atan (const Vector& v) { return unary_map<atan_op> (v); }
01186 
01187 TMPL Vector pow (const Vector& v, const Vector& w) {
01188   return binary_map<pow_op> (v, w); }
01189 TMPL Vector pow (const Vector& v, const C& w) {
01190   return binary_map_scalar<rpow_op> (v, w); }
01191 TMPL Vector pow (const C& w, const Vector& v) {
01192   return binary_map_scalar<lpow_op> (v, w); }
01193 TMPL Vector pow (const Vector& v, const int& w) {
01194   return binary_map_scalar<rpow_op> (v, w); }
01195 TMPL Vector pow (const int& w, const Vector& v) {
01196   return binary_map_scalar<lpow_op> (v, w); }
01197 
01198 TMPL Vector derive (const Vector& v) {
01199   return unary_map<derive_op> (v); }
01200 TMPL Vector integrate (const Vector& v) {
01201   return unary_map<integrate_op> (v); }
01202 template<typename C, typename V, typename K, typename W> Vector
01203 integrate_init (const Vector& v, const vector<K,W>& w) {
01204   return binary_map<integrate_init_op> (v, as<vector<K,V> > (w)); }
01205 template<typename C,typename V,typename X> Vector
01206 derive (const Vector& v, const X& x) {
01207   return binary_map_scalar<derive_op> (v, x); }
01208 template<typename C,typename V,typename X> Vector
01209 integrate (const Vector& v, const X& x) {
01210   return binary_map_scalar<integrate_op> (v, x); }
01211 
01212 TMPL inline C big_add (const Vector& v) { return big<add_op> (v); }
01213 TMPL inline C big_mul (const Vector& v) { return big<mul_op> (v); }
01214 TMPL inline C big_or  (const Vector& v) { return big< or_op> (v); }
01215 TMPL inline C big_and (const Vector& v) { return big<and_op> (v); }
01216 TMPL inline C big_min (const Vector& v) { return big<min_op> (v); }
01217 TMPL inline C big_max (const Vector& v) { return big<max_op> (v); }
01218 TMPL inline C big_inf (const Vector& v) { return big<inf_op> (v); }
01219 TMPL inline C big_sup (const Vector& v) { return big<sup_op> (v); }
01220 
01221 
01222 
01223 
01224 
01225 TMPL
01226 struct lift_helper<Vector> {
01227   template<typename T,typename W> static void
01228   set_op (vector<T,W>& w, const Vector& v) {
01229     typedef implementation<vector_linear,V> Vec;
01230     format<T> fm= CF(w);
01231     if (is_a_scalar (v)) {
01232       T a= get_sample (fm);
01233       set_lift (a, v.scalar());
01234       w= vector<T,W> (a, fm);
01235     } else {
01236       nat n= N(v);
01237       nat l= aligned_size<T,W> (n);
01238       T* r= mmx_formatted_new<T> (l, fm);
01239       Vec::template vec_unary<lift_op> (r, seg (v), n);
01240       w= vector<T,W> (r, n, l, fm); } }
01241   static inline Lifted_vector
01242   op (const Vector& v) {
01243     Lifted_vector w (lift (CF(v)));
01244     set_op (w, v);
01245     return w; }
01246 };
01247 
01248 TMPL
01249 struct project_helper<Vector> {
01250   template<typename T,typename W> static void
01251   set_op (vector<T,W>& w, const Vector& v) {
01252     typedef implementation<vector_linear,V> Vec;
01253     format<T> fm= CF(w);
01254     if (is_a_scalar (v)) {
01255       T a= get_sample (fm);
01256       set_project (a, v.scalar());
01257       w= vector<T,W> (a, fm);
01258     } else {
01259       nat n= N(v);
01260       nat l= aligned_size<T,W> (n);
01261       T* r= mmx_formatted_new<T> (l, fm);
01262       Vec::template vec_unary<project_op> (r, seg (v), n);
01263       w= vector<T,W> (r, n, l, fm); } }
01264   static inline Projected_vector
01265   op (const Vector& v) {
01266     Projected_vector w (project (CF(v)));
01267     set_op (w, v);
01268     return w; }
01269 };
01270 
01271 TMPL Reconstructed_vector reconstruct (const Vector& v) {
01272   return as<Reconstructed_vector> (unary_map<reconstruct_op> (v)); }
01273 
01274 template<typename C,typename V,typename W> bool
01275 is_reconstructible (const Vector& v, vector<Reconstruct_type(C),W>& w) {
01276   w= fill (N(v), CF(w));
01277   for (nat i= 0; i < N(v); i++)
01278     if (!is_reconstructible (v[i], w[i])) return false;
01279   return true; }
01280 
01281 
01282 
01283 
01284 
01285 TMPLK Evaluated_vector evaluate (const Vector& v, const K& x) {
01286   return as<Evaluated_vector> (binary_map_scalar<evaluate_op> (v, x)); }
01287 
01288 template<typename C,typename V,typename K,typename W> bool
01289 is_evaluable (const Vector& v, const K& a,
01290               vector<Evaluate_type(C,K),W>& w) {
01291   w= fill (N(v), CF(w));
01292   for (nat i= 0; i < N(v); i++)
01293     if (!is_evaluable (v[i], a, w[i])) return false;
01294   return true; }
01295 
01296 
01297 
01298 
01299 
01300 TMPL Truncated_vector truncate (const Vector& v, nat n) {
01301   return as<Truncated_vector> (binary_map_scalar<truncate_op> (v, n)); }
01302 
01303 TMPL Completed_vector complete (const Vector& v) {
01304   return as<Completed_vector> (unary_map<complete_op> (v)); }
01305 
01306 
01307 
01308 
01309 
01310 TMPL inline bool is_finite (const Vector& v) {
01311   if (is_a_scalar (v)) return is_finite (v.scalar());
01312   return big<and_is_finite_op> (v); }
01313 TMPL inline bool is_nan (const Vector& v) {
01314   if (is_a_scalar (v)) return is_nan (v.scalar());
01315   return big<or_is_nan_op> (v); }
01316 TMPL inline bool is_infinite (const Vector& v) {
01317   if (is_a_scalar (v)) return is_infinite (v.scalar());
01318   return !is_nan (v) && big<or_is_infinite_op> (v); }
01319 TMPL inline bool is_fuzz (const Vector& v) {
01320   if (is_a_scalar (v)) return is_fuzz (v.scalar());
01321   return !is_nan (v) && big<or_is_fuzz_op> (v); }
01322 TMPL inline bool is_reliable (const Vector& v) {
01323   if (is_a_scalar (v)) return is_reliable (v.scalar());
01324   return is_reliable (C (0)); }
01325 
01326 TMPL inline Vector change_precision (const Vector& v, xnat p) {
01327   return binary_map_scalar<change_precision_op> (v, p); }
01328 TMPL inline xnat precision (const Vector& v) {
01329   return big<min_precision_op> (v); }
01330 
01331 TMPL inline xint exponent (const Vector& v) {
01332   return big<max_exponent_op> (v); }
01333 TMPL inline double magnitude (const Vector& v) {
01334   return big<max_magnitude_op> (v); }
01335 
01336 
01337 
01338 
01339 
01340 TMPL Abs_vector abs (const Vector& v) { return unary_map<abs_op> (v); }
01341 TMPL Real_vector Re (const Vector& v) { return unary_map<Re_op> (v); }
01342 TMPL Real_vector Im (const Vector& v) { return unary_map<Im_op> (v); }
01343 TMPL Vector conj (const Vector& v) { return unary_map<conj_op> (v); }
01344 
01345 TMPL Center_vector center (const Vector& v) {
01346   return unary_map<center_op> (v); }
01347 TMPL Radius_vector radius (const Vector& v) {
01348   return unary_map<radius_op> (v); }
01349 TMPL Center_vector lower (const Vector& v) {
01350   return unary_map<lower_op> (v); }
01351 TMPL Center_vector upper (const Vector& v) {
01352   return unary_map<upper_op> (v); }
01353 TMPL inline Vector sharpen (const Vector& v) {
01354   return unary_map<sharpen_op> (v); }
01355 TMPLK inline Vector blur (const Vector& v, const K& x) {
01356   return binary_map_scalar<blur_op> (v, x); }
01357 template<typename C, typename V, typename C2, typename V2> inline Vector
01358 blur (const Vector& v, const vector<C2,V2>& w) {
01359   return binary_map<blur_op> (v, w); }
01360 TMPL inline bool included (const Vector& v, const Vector& w) {
01361   return binary_test<included_op> (v, w); }
01362 
01363 
01364 
01365 
01366 
01367 template<typename C> inline C
01368 trig_op::diff_op (const C& me, const C& x) {
01369   return from_vector (vec (-derive (access (x, 0)) * access (me, 1),
01370                            derive (access (x, 0)) * access (me, 0)));
01371 }
01372 
01373 template<typename C> inline C
01374 tan_op::def (const C& me, const C& f) {
01375   (void) me; const vector<C> v= cos_sin (f); return v[1] / v[0];
01376 }
01377 
01378 template<typename C> inline C
01379 solve_vector_lde_op::diff_op (const C& me, const C& x) {
01380   typedef As_vector_type (C) V;
01381   typedef Scalar_type (V) K;
01382   const V w= as_vector (me);
01383   const V v= as_vector (x);
01384   nat i, n= N(v);
01385   K sum (0);
01386   V d (K (0), n);
01387   for (i=0; i<n; i++) {
01388     if (i==0) sum= v[i] * w[i];
01389     else sum= sum + v[i] * w[i];
01390     if (i<n-1) d[i]= w[i+1];
01391     else d[i]= sum;
01392   }
01393   return from_vector (d);
01394 }
01395 
01396 
01397 
01398 
01399 
01400 TMPL C dot (const Vector& v, const Vector& w) {
01401   typedef implementation<vector_linear,V> Vec;
01402   ASSERT (is_non_scalar (v) && is_non_scalar (w),
01403           "non-scalar vectors expected");
01404   ASSERT (N(v) == N(w), "lengths don't match");
01405   return Vec::inn_prod (seg (v), seg (w), N(v), CF (v)); }
01406 
01407 template<typename C> vector<C> cos_sin (const C& x) {
01408   return vec (cos (x), sin (x)); }
01409 template<typename C> vector<C> trig (const vector<C>& v) {
01410   return vec (cos (v[0]), sin (v[0])); }
01411 
01412 #undef TMPL_DEF
01413 #undef TMPL
01414 #undef TMPLK
01415 #undef Format
01416 #undef Vector
01417 #undef Vector_int
01418 #undef Vector_rep
01419 #undef Abs_vector
01420 #undef Real_vector
01421 #undef Center_vector
01422 #undef Radius_vector
01423 #undef Lifted_vector
01424 #undef Projected_vector
01425 #undef Reconstructed_vector
01426 #undef Truncated_vector
01427 #undef Completed_vector
01428 #undef Evaluate_vector
01429 } 
01430 #endif // __MMX_VECTOR_HPP