array.hpp

Go to the documentation of this file.

/*********************************************************************
*      This file is part of the source code of realroot kernel.        *
*   Author(s): B. Mourrain, GALAAD, INRIA                            *
**********************************************************************/
#ifndef realroot_array_hpp
#define realroot_array_hpp
//------------------------------------------------------------------------
#include <cassert>
#include <realroot/assign.hpp>
//======================================================================
namespace mmx {
//======================================================================

namespace array {
  //----------------------------------------------------------------------
  template <class OS, class R>
  inline OS & print(OS & os, const R & v) {
    os <<"[";
    typename R::const_iterator i=v.begin();
    if (v.size()){
      os<<*i; ++i;
      for(; i != v.end(); ++i) {os <<", "; os<<*i;}
    }
    os <<"]";
    return os;
  }
  //----------------------------------------------------------------------
  template <class OS, class R>
  inline OS & print(OS & os, const R & v, unsigned sz) {
    os <<"[";
    unsigned i=0;
    if (sz){
      os<<v[i]; ++i;
      for(; i <sz; ++i) os <<","<<v[i];
    }
    os <<"]";
    return os;
  }
  //----------------------------------------------------------------------
  template <class IS, class R>
  inline IS & read(IS & is, R & V) {
    typedef typename R::size_type size_type;
    size_type s;
    is >> s;
    V.rep().resize(s);
    for(size_type i=0; i< s; ++i) is >> V[i];
    return(is);
  }
  //----------------------------------------------------------------------
  
  template <class V1, class V2> inline
  bool equal_n( const V1& v1, const V2 & v2, int n )
  {
    int i=0;
    while(i < n && v1[i] == v2[i]) i++;
    return i == n;
  };
  
  template <class V1, class V2> inline
  bool equal(const V1 & v1, const V2 & v2) 
  { return (v1.size() == v2.size()) && equal_n(v1,v2,v1.size()); }

  
  //---------------------------------------------------------------------- 
  template <class V, class C> inline
  void init_n(V & a, const C & v, int n ) { 
    for ( int i = 0; i < n; a[i] = v, i ++ ) {} };

  template <class V, class C>
  void set_cst(V & a, const C & v ) { init_n(a,v,(int)a.size()); };


  template <class V, class W>
  void set(V & a, const W & b ) {for ( int i = 0; i < (int)a.size(); a[i] = b[i], i ++ ) {} };
  
  //---------------------------------------------------------------------- 
  template <class V>
  void reverse_n( V & v, int a, int n )
  {
    int k = n - a;
    for ( int i = 0; i < k/2; std::swap(v[a+i],v[k-1-i]), i ++ ) {}
  };
  
  template <class V>
  void reverse( V & v, int n ) { reverse_n( v, 0, n ); };
  
  template <class V,class I>
  void reverse(V & v, I d) { 
    for(I i=0;i<d/2; std::swap(v[i],v[d-1-i]),++i) {}
  }
  
  //----------------------------------------------------------------------
  template<class V, class W> inline
  void neg_range( V & r, const W & w, int a, int b ) { 
    for ( int i = a; i < b; r[i] = -w[i], i++ ) {}
  }
  template<class V, class W> inline
  void neg_range( V & r, int a, int b )
  { for ( int i = a; i < b; r[i] = -r[i], i++ ){} };
  
  template<class VI> inline
  void neg_i( VI a, VI eda ) { for ( ; a != eda; *a = -(*a), a++ ){} };
  template<class VI> inline
  void neg_i( VI a, VI eda, VI b ) { for ( ; a != eda; *a = -(*b), a++, b++  ){} };
  template<class V> inline
  void neg( V & v ) { neg_i(v.begin(),v.end()); };
  template<class V> inline
  void neg( V & a, const V & b )
  { 
    a.resize( b.size() );
    neg_i( a.begin(), a.end(), b.begin() ); 
  };
  
  template<class V, class W> inline
  void copy_range( V & v, const W & w, int a, int b ) { for ( int i = a; i < b; v[i] = w[i], i ++ ){} }
  template< class V, class W >
  void add_n(V & a, const W & b, int n ) {for ( int i = 0; i < n; a[i] += b[i], i ++ ){} }
  template<class VI, class WI>
  void add_i( VI a, VI ea, WI b ) { for ( ; a != ea; *a += *b, ++a, ++b ){} };
  //----------------------------------------------------------------------
  template <class V, class W>
  void add(V & a, const W & b)
  {
    int asz, bsz;
    asz = a.size(), bsz = b.size();
    if ( asz < bsz  ) 
      { 
    a.resize(bsz);
    add_i(a.begin(),a.end(),b.begin()); 
    copy_range(a,b,asz,bsz);
      }
    else 
      {
    //  PRINT_DEBUG("add");
    add_n(a,b,bsz);
    //  add_i(a.begin(),a.begin()+bsz,b.begin());
      }
  }

  template< class V, class W, class X > inline
  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 ++ ){} };

  //--------------------------------------------------------------------
  template <class V, class W,class X>
  void add_g(V & r, const W & a, const X & b)
  {  
    int asz, bsz, rsz;
    asz = a.size(), bsz = b.size(), rsz = r.size();
    if ( asz > bsz ) 
      {
    if ( asz > rsz ) r.resize(asz);
    add_n(r,a,b,bsz);
    copy_range(r,a,bsz,asz);
      }
    else
      { 
    if ( bsz > rsz ) r.resize(bsz);
    add_n(r,a,b,asz); 
    copy_range(r,b,asz,bsz); 
      };
  }
  //--------------------------------------------------------------------
  template <class V, class W,class X> inline
  void add(V & r, const W & a, const X & b)
  {
    add_g(r,a,b);
  }
  //--------------------------------------------------------------------
  template <class V> inline
  void add(V & r, const V & a, const V & b)
  {
    if(&r==&a) { add(r,b); return; };
    if(&r==&b) { add(r,a); return; };
    add_g(r,a,b);
  }
  //====================================================================
  template< class V, class W >
  void sub_n(V & a, const W & b, int n ) { for ( int i = 0; i < n; a[i] -= b[i], i ++ ){} };

  //----------------------------------------------------------------------
  template <class V, class W>
  void sub(V & a, const W & b)
  {
    int asz, bsz;
    asz = a.size(), bsz = b.size();
    if ( asz < bsz  ) 
      { 
    a.resize(bsz);
    sub_n(a,b,asz);
    neg_range(a,b,asz,bsz); 
      }
    else  
      sub_n(a,b,bsz);
  }
  //--------------------------------------------------------------------
  template< class V, class W, class X > inline
  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 ++ ){} };
  
  //--------------------------------------------------------------------
  template <class V, class W,class X>
  void sub_g(V & r, const W & a, const X & b)
  { 
    int asz, bsz, rsz;
    asz = a.size(), bsz = b.size(), rsz = r.size();
    
    if ( asz > bsz ) 
      {
    if ( asz > rsz ) r.resize(asz);
    sub_n(r,a,b,bsz);
    copy_range(r,a,bsz,asz);
      }
    else
      { 
    if ( bsz > rsz ) r.resize(bsz);
    sub_n(r,a,b,asz);
    neg_range(r,b,asz,bsz);
      };
  }
  template <class V, class W,class X>
  void sub(V & r, const W & a, const X & b)
  {
    sub_g(r,a,b);
  }
  //--------------------------------------------------------------------
  template <class V>
  void sub(V & r, const V & a, const V & b)
  {
    if(&r==&a) { sub(r,b); return; };
    //    if(&r==&b) { sub(r,a); mul_ext(r,-1); return; };
    sub_g(r,a,b);
  }
  //----------------------------------------------------------------------
  template<class V, class W> inline
  void mul_ext_n( V & a, const W & c, int n )
  { for ( int i = 0; i < n; a[i] = a[i]*c, i ++ ){} };
  template<class V, class W, class S> inline
  void mul_ext_n( V & a, const W & b, const S& c, int n ) 
  { for ( int i = 0; i < n; a[i] = b[i]*c, i++ ){} };
  template<class VI, class W> inline
  void mul_ext_i( VI bga, VI eda, const W & c )
    { for ( VI i = bga; i != eda; *i *= c, ++i  ){} };
  template<class VI, class VII, class W> inline
  void mul_ext_i( VI bga, VI eda, VII bgb, const W & c )
  { for ( VI a = bga; a != eda; *a = *bgb * c, ++bgb, ++a ){} };

  template <class V,class W>
  void mul_ext(V & a,  const V & b, const W & c)
  {
    if(&a != &b) {
      a.resize( b.size() );
      mul_ext_i(a.begin(),a.end(),b.begin(),c);
    }
    else 
      mul_ext_i(a.begin(),a.end(),c);
  }
  template <class V,class W>
  void mul_ext(V & a,  const W & c) {
    mul_ext_i(a.begin(),a.end(),c);
  }
  //----------------------------------------------------------------------
  template<class V, class SC> inline
  void div_ext_n( V & a, const SC & sc, int n )
  { for ( int i = 0; i < n; a[i] /= sc, i++){} }
  template<class V1, class V2, class SC>
  void div_ext_n( V1 & a, const V2 & b, const SC& sc, int n )
  { for ( int i = 0; i < n; a[i] = b[i]/sc, i ++ ){} };
  template<class VI, class SC> inline
  void div_ext_i( VI  bga, VI eda, const SC & sc )
  { for ( ; bga < eda; *bga /= sc, bga++){} }
  template<class VI, class SC> inline
  void div_ext_i( VI a, VI eda, VI b, const SC & sc )
  { for ( ; a < eda; *a = *b / sc, a++, b ++ ){} };
  template <class V, class SC>
  void div_ext( V & a,  const V & b, const SC & sc)
  {
    if(&a != &b) {
      a.resize( b.size() );
      div_ext_i(a.begin(),a.end(),b.begin(),sc);
    }
    else         div_ext_i(a.begin(),a.end(),sc);
  }
  template <class V, class W>
  void div_ext(V & a,  const W & c)
  { div_ext_i(a.begin(),a.end(),c); }
  //----------------------------------------------------------------------
  template <class V, class W>
  void div(V & a, const W & c)
  {
    div_ext(a,c);
    //W tmp(c); for(typename V::iterator it=a.begin(); it !=a.end(); ++it) (*it)/=tmp;
  }
  
  //----------------------------------------------------------------------
  template <class R, class S, class C> inline
  void apply_mult(R & result, const S& a, const C & m)
  { 
    typedef typename S::const_iterator  const_iterator;
    typedef typename R::iterator        iterator;
    result.resize(a.size());
    const_iterator b=a.begin(); iterator i=result.begin();
    for(; b!=a.end(); ++i, ++b) *i = (*b) * m;
  }
  //----------------------------------------------------------------------
  template <class C, class R>
  C  norm(const R & v)
  {
    C r=0;
    for(typename R::const_iterator it=v.begin(); it !=v.end(); ++it) r+=(*it)*(*it);
    return sqrt(r);
  }
  //----------------------------------------------------------------------
  template <class C,class R>
  C norm(const R & v, int p)
  {
    C r=0;
    for(typename R::const_iterator it=v.begin(); it !=v.end(); ++it) r+=pow(*it,p);
    return pow(r,1/p);
  }
  //----------------------------------------------------------------------
  template <class R, class S>
  typename R::value_type innerprod_n( const S & v, const R & w, unsigned n )
  {
    typename R::value_type r = 0;
    for ( unsigned i = 0; i < n; r += v[i]*w[i], i++ ) {}
    return r;
  };
  
  template <class R, class S>
  typename R::value_type innerprod(const S & v, const R & w)
  {
    typename R::value_type r =0;
    for(unsigned int i=0; i <v.size(); i++){r += (v[i]*w[i]);}
    return r;
  }
  //--------------------------------------------------------------------
  template<class R>
  typename R::value_type max_abs(const R& v)
  {
    typename R::value_type r=0;
    for(unsigned i=0; i <(unsigned)v.size(); i++)
      {
    if(r>v[i]) r=v[i];
    else if(r>-v[i]) r=-v[i];
      }
    return r;
  }

  //--------------------------------------------------------------------
  template<class U, class R>
  void lcm_denominator(U& r, const R& v)
  {
    for(unsigned i=0; i <v.size(); i++)
      {
    r = lcm(r,denominator(v[i]));
      }
  }
  //----------------------------------------------------------------------
  /* Replace R by its reverse */
  template < typename V> inline 
  void  reverse(V & v)
  { 
    unsigned sz = v.size();
    V temp(v);
    for (unsigned i = 0; i < sz; ++i)
      v[i] = temp[sz-1-i];
}

  template<class Ia, class Ib> inline
  void assign_i( Ia a, Ia eda, Ib b )
  {
    for ( ; a != eda; ++a, ++b  ) let::assign(*a,*b);
  };
  
  template<class Ra, class Rb> inline
  void assign( Ra & a, const Rb & b )
  {
    a.resize( b.size() );
    assign_i( a.begin(), a.end(), b.begin() );
  };

  template<class R, class C> inline
  void vaddsc( R & r, const C & x )
  { for ( typename R::iterator i = r.begin(); i != r.end(); *i += x, i ++ ){} };
  
} //CLOSE


//======================================================================
} // namespace mmx
//======================================================================
#endif // realroot_LINALG_VECTOR_H