realroot: include/realroot/array.hpp Source File

00001 /*********************************************************************
00002 *      This file is part of the source code of realroot kernel.        *
00003 *   Author(s): B. Mourrain, GALAAD, INRIA                            *
00004 **********************************************************************/
00006 #ifndef realroot_array_hpp
00007 #define realroot_array_hpp
00008 //------------------------------------------------------------------------
00009 #include <cassert>
00010 #include <realroot/assign.hpp>
00011 //======================================================================
00012 namespace mmx {
00013 //======================================================================
00014 
00035 namespace array {
00036   //----------------------------------------------------------------------
00038   template <class OS, class R>
00039   inline OS & print(OS & os, const R & v) {
00040     os <<"[";
00041     typename R::const_iterator i=v.begin();
00042     if (v.size()){
00043       os<<*i; ++i;
00044       for(; i != v.end(); ++i) {os <<", "; os<<*i;}
00045     }
00046     os <<"]";
00047     return os;
00048   }
00049   //----------------------------------------------------------------------
00051   template <class OS, class R>
00052   inline OS & print(OS & os, const R & v, unsigned sz) {
00053     os <<"[";
00054     unsigned i=0;
00055     if (sz){
00056       os<<v[i]; ++i;
00057       for(; i <sz; ++i) os <<","<<v[i];
00058     }
00059     os <<"]";
00060     return os;
00061   }
00062   //----------------------------------------------------------------------
00067   template <class IS, class R>
00068   inline IS & read(IS & is, R & V) {
00069     typedef typename R::size_type size_type;
00070     size_type s;
00071     is >> s;
00072     V.rep().resize(s);
00073     for(size_type i=0; i< s; ++i) is >> V[i];
00074     return(is);
00075   }
00076   //----------------------------------------------------------------------
00077   
00078   template <class V1, class V2> inline
00079   bool equal_n( const V1& v1, const V2 & v2, int n )
00080   {
00081     int i=0;
00082     while(i < n && v1[i] == v2[i]) i++;
00083     return i == n;
00084   };
00085   
00087   template <class V1, class V2> inline
00088   bool equal(const V1 & v1, const V2 & v2) 
00089   { return (v1.size() == v2.size()) && equal_n(v1,v2,v1.size()); }
00090 
00091   
00092   //---------------------------------------------------------------------- 
00093   template <class V, class C> inline
00094   void init_n(V & a, const C & v, int n ) { 
00095     for ( int i = 0; i < n; a[i] = v, i ++ ) {} };
00096 
00098   template <class V, class C>
00099   void set_cst(V & a, const C & v ) { init_n(a,v,(int)a.size()); };
00100 
00101 
00102   template <class V, class W>
00103   void set(V & a, const W & b ) {for ( int i = 0; i < (int)a.size(); a[i] = b[i], i ++ ) {} };
00104   
00105   //---------------------------------------------------------------------- 
00107   template <class V>
00108   void reverse_n( V & v, int a, int n )
00109   {
00110     int k = n - a;
00111     for ( int i = 0; i < k/2; std::swap(v[a+i],v[k-1-i]), i ++ ) {}
00112   };
00113   
00114   template <class V>
00115   void reverse( V & v, int n ) { reverse_n( v, 0, n ); };
00116   
00118   template <class V,class I>
00119   void reverse(V & v, I d) { 
00120     for(I i=0;i<d/2; std::swap(v[i],v[d-1-i]),++i) {}
00121   }
00122   
00123   //----------------------------------------------------------------------
00124   template<class V, class W> inline
00125   void neg_range( V & r, const W & w, int a, int b ) { 
00126     for ( int i = a; i < b; r[i] = -w[i], i++ ) {}
00127   }
00128   template<class V, class W> inline
00129   void neg_range( V & r, int a, int b )
00130   { for ( int i = a; i < b; r[i] = -r[i], i++ ){} };
00131   
00132   template<class VI> inline
00133   void neg_i( VI a, VI eda ) { for ( ; a != eda; *a = -(*a), a++ ){} };
00134   template<class VI> inline
00135   void neg_i( VI a, VI eda, VI b ) { for ( ; a != eda; *a = -(*b), a++, b++  ){} };
00136   template<class V> inline
00137   void neg( V & v ) { neg_i(v.begin(),v.end()); };
00138   template<class V> inline
00139   void neg( V & a, const V & b )
00140   { 
00141     a.resize( b.size() );
00142     neg_i( a.begin(), a.end(), b.begin() ); 
00143   };
00144   
00145   template<class V, class W> inline
00146   void copy_range( V & v, const W & w, int a, int b ) { for ( int i = a; i < b; v[i] = w[i], i ++ ){} }
00147   template< class V, class W >
00148   void add_n(V & a, const W & b, int n ) {for ( int i = 0; i < n; a[i] += b[i], i ++ ){} }
00149   template<class VI, class WI>
00150   void add_i( VI a, VI ea, WI b ) { for ( ; a != ea; *a += *b, ++a, ++b ){} };
00151   //----------------------------------------------------------------------
00153   template <class V, class W>
00154   void add(V & a, const W & b)
00155   {
00156     int asz, bsz;
00157     asz = a.size(), bsz = b.size();
00158     if ( asz < bsz  ) 
00159       { 
00160         a.resize(bsz);
00161         add_i(a.begin(),a.end(),b.begin()); 
00162         copy_range(a,b,asz,bsz);
00163       }
00164     else 
00165       {
00166         //      PRINT_DEBUG("add");
00167         add_n(a,b,bsz);
00168         //      add_i(a.begin(),a.begin()+bsz,b.begin());
00169       }
00170   }
00171 
00172   template< class V, class W, class X > inline
00173   void add_n(V & a, const W & b, const X& c, int n ) { for ( int i = 0; i < n; a[i] = b[i]+c[i], i ++ ){} };
00174 
00175   //--------------------------------------------------------------------
00178   template <class V, class W,class X>
00179   void add_g(V & r, const W & a, const X & b)
00180   {  
00181     int asz, bsz, rsz;
00182     asz = a.size(), bsz = b.size(), rsz = r.size();
00183     if ( asz > bsz ) 
00184       {
00185         if ( asz > rsz ) r.resize(asz);
00186         add_n(r,a,b,bsz);
00187         copy_range(r,a,bsz,asz);
00188       }
00189     else
00190       { 
00191         if ( bsz > rsz ) r.resize(bsz);
00192         add_n(r,a,b,asz); 
00193         copy_range(r,b,asz,bsz); 
00194       };
00195   }
00196   //--------------------------------------------------------------------
00199   template <class V, class W,class X> inline
00200   void add(V & r, const W & a, const X & b)
00201   {
00202     add_g(r,a,b);
00203   }
00204   //--------------------------------------------------------------------
00207   template <class V> inline
00208   void add(V & r, const V & a, const V & b)
00209   {
00210     if(&r==&a) { add(r,b); return; };
00211     if(&r==&b) { add(r,a); return; };
00212     add_g(r,a,b);
00213   }
00214   //====================================================================
00215   template< class V, class W >
00216   void sub_n(V & a, const W & b, int n ) { for ( int i = 0; i < n; a[i] -= b[i], i ++ ){} };
00217 
00218   //----------------------------------------------------------------------
00220   template <class V, class W>
00221   void sub(V & a, const W & b)
00222   {
00223     int asz, bsz;
00224     asz = a.size(), bsz = b.size();
00225     if ( asz < bsz  ) 
00226       { 
00227         a.resize(bsz);
00228         sub_n(a,b,asz);
00229         neg_range(a,b,asz,bsz); 
00230       }
00231     else  
00232       sub_n(a,b,bsz);
00233   }
00234   //--------------------------------------------------------------------
00235   template< class V, class W, class X > inline
00236   void sub_n(V & a, const W & b, const X& c, int n ) { for ( int i = 0; i < n; a[i] = b[i]-c[i], i ++ ){} };
00237   
00238   //--------------------------------------------------------------------
00239   template <class V, class W,class X>
00240   void sub_g(V & r, const W & a, const X & b)
00241   { 
00242     int asz, bsz, rsz;
00243     asz = a.size(), bsz = b.size(), rsz = r.size();
00244     
00245     if ( asz > bsz ) 
00246       {
00247         if ( asz > rsz ) r.resize(asz);
00248         sub_n(r,a,b,bsz);
00249         copy_range(r,a,bsz,asz);
00250       }
00251     else
00252       { 
00253         if ( bsz > rsz ) r.resize(bsz);
00254         sub_n(r,a,b,asz);
00255         neg_range(r,b,asz,bsz);
00256       };
00257   }
00260   template <class V, class W,class X>
00261   void sub(V & r, const W & a, const X & b)
00262   {
00263     sub_g(r,a,b);
00264   }
00265   //--------------------------------------------------------------------
00268   template <class V>
00269   void sub(V & r, const V & a, const V & b)
00270   {
00271     if(&r==&a) { sub(r,b); return; };
00272     //    if(&r==&b) { sub(r,a); mul_ext(r,-1); return; };
00273     sub_g(r,a,b);
00274   }
00275   //----------------------------------------------------------------------
00276   template<class V, class W> inline
00277   void mul_ext_n( V & a, const W & c, int n )
00278   { for ( int i = 0; i < n; a[i] = a[i]*c, i ++ ){} };
00279   template<class V, class W, class S> inline
00280   void mul_ext_n( V & a, const W & b, const S& c, int n ) 
00281   { for ( int i = 0; i < n; a[i] = b[i]*c, i++ ){} };
00282   template<class VI, class W> inline
00283   void mul_ext_i( VI bga, VI eda, const W & c )
00284     { for ( VI i = bga; i != eda; *i *= c, ++i  ){} };
00285   template<class VI, class VII, class W> inline
00286   void mul_ext_i( VI bga, VI eda, VII bgb, const W & c )
00287   { for ( VI a = bga; a != eda; *a = *bgb * c, ++bgb, ++a ){} };
00288 
00290   template <class V,class W>
00291   void mul_ext(V & a,  const V & b, const W & c)
00292   {
00293     if(&a != &b) {
00294       a.resize( b.size() );
00295       mul_ext_i(a.begin(),a.end(),b.begin(),c);
00296     }
00297     else 
00298       mul_ext_i(a.begin(),a.end(),c);
00299   }
00300   template <class V,class W>
00301   void mul_ext(V & a,  const W & c) {
00302     mul_ext_i(a.begin(),a.end(),c);
00303   }
00304   //----------------------------------------------------------------------
00305   template<class V, class SC> inline
00306   void div_ext_n( V & a, const SC & sc, int n )
00307   { for ( int i = 0; i < n; a[i] /= sc, i++){} }
00308   template<class V1, class V2, class SC>
00309   void div_ext_n( V1 & a, const V2 & b, const SC& sc, int n )
00310   { for ( int i = 0; i < n; a[i] = b[i]/sc, i ++ ){} };
00311   template<class VI, class SC> inline
00312   void div_ext_i( VI  bga, VI eda, const SC & sc )
00313   { for ( ; bga < eda; *bga /= sc, bga++){} }
00314   template<class VI, class SC> inline
00315   void div_ext_i( VI a, VI eda, VI b, const SC & sc )
00316   { for ( ; a < eda; *a = *b / sc, a++, b ++ ){} };
00318   template <class V, class SC>
00319   void div_ext( V & a,  const V & b, const SC & sc)
00320   {
00321     if(&a != &b) {
00322       a.resize( b.size() );
00323       div_ext_i(a.begin(),a.end(),b.begin(),sc);
00324     }
00325     else         div_ext_i(a.begin(),a.end(),sc);
00326   }
00328   template <class V, class W>
00329   void div_ext(V & a,  const W & c)
00330   { div_ext_i(a.begin(),a.end(),c); }
00331   //----------------------------------------------------------------------
00333   template <class V, class W>
00334   void div(V & a, const W & c)
00335   {
00336     div_ext(a,c);
00337     //W tmp(c); for(typename V::iterator it=a.begin(); it !=a.end(); ++it) (*it)/=tmp;
00338   }
00339   
00340   //----------------------------------------------------------------------
00341   template <class R, class S, class C> inline
00342   void apply_mult(R & result, const S& a, const C & m)
00343   { 
00344     typedef typename S::const_iterator  const_iterator;
00345     typedef typename R::iterator        iterator;
00346     result.resize(a.size());
00347     const_iterator b=a.begin(); iterator i=result.begin();
00348     for(; b!=a.end(); ++i, ++b) *i = (*b) * m;
00349   }
00350   //----------------------------------------------------------------------
00352   template <class C, class R>
00353   C  norm(const R & v)
00354   {
00355     C r=0;
00356     for(typename R::const_iterator it=v.begin(); it !=v.end(); ++it) r+=(*it)*(*it);
00357     return sqrt(r);
00358   }
00359   //----------------------------------------------------------------------
00361   template <class C,class R>
00362   C norm(const R & v, int p)
00363   {
00364     C r=0;
00365     for(typename R::const_iterator it=v.begin(); it !=v.end(); ++it) r+=pow(*it,p);
00366     return pow(r,1/p);
00367   }
00368   //----------------------------------------------------------------------
00370   template <class R, class S>
00371   typename R::value_type innerprod_n( const S & v, const R & w, unsigned n )
00372   {
00373     typename R::value_type r = 0;
00374     for ( unsigned i = 0; i < n; r += v[i]*w[i], i++ ) {}
00375     return r;
00376   };
00377   
00378   template <class R, class S>
00379   typename R::value_type innerprod(const S & v, const R & w)
00380   {
00381     typename R::value_type r =0;
00382     for(unsigned int i=0; i <v.size(); i++){r += (v[i]*w[i]);}
00383     return r;
00384   }
00385   //--------------------------------------------------------------------
00386   template<class R>
00387   typename R::value_type max_abs(const R& v)
00388   {
00389     typename R::value_type r=0;
00390     for(unsigned i=0; i <(unsigned)v.size(); i++)
00391       {
00392         if(r>v[i]) r=v[i];
00393         else if(r>-v[i]) r=-v[i];
00394       }
00395     return r;
00396   }
00397 
00398   //--------------------------------------------------------------------
00399   template<class U, class R>
00400   void lcm_denominator(U& r, const R& v)
00401   {
00402     for(unsigned i=0; i <v.size(); i++)
00403       {
00404         r = lcm(r,denominator(v[i]));
00405       }
00406   }
00407   //----------------------------------------------------------------------
00408   /* Replace R by its reverse */
00409   template < typename V> inline 
00410   void  reverse(V & v)
00411   { 
00412     unsigned sz = v.size();
00413     V temp(v);
00414     for (unsigned i = 0; i < sz; ++i)
00415       v[i] = temp[sz-1-i];
00416 }
00417 
00418   template<class Ia, class Ib> inline
00419   void assign_i( Ia a, Ia eda, Ib b )
00420   {
00421     for ( ; a != eda; ++a, ++b  ) let::assign(*a,*b);
00422   };
00423   
00424   template<class Ra, class Rb> inline
00425   void assign( Ra & a, const Rb & b )
00426   {
00427     a.resize( b.size() );
00428     assign_i( a.begin(), a.end(), b.begin() );
00429   };
00430 
00431   template<class R, class C> inline
00432   void vaddsc( R & r, const C & x )
00433   { for ( typename R::iterator i = r.begin(); i != r.end(); *i += x, i ++ ){} };
00434   
00435 } //CLOSE
00436 
00437 
00438 //======================================================================
00439 } // namespace mmx
00440 //======================================================================
00441 #endif // realroot_LINALG_VECTOR_H
00442