tensor_vctops.hpp

Go to the documentation of this file.

#ifndef realroot_MPOLDSE_VCTOPS_H
#define realroot_MPOLDSE_VCTOPS_H
#include <iostream>
#include <realroot/texp_sup.hpp>
#include <realroot/texp_ringof.hpp>

namespace mmx {

namespace vct
{
  
  template<typename A, typename B> inline
  void padd( A * a, const B * b, unsigned sz, int sta = 1, int stb = 1  )
  { 
    for ( A * ea = a; ea != a+sz*sta; *ea += *b , ea += sta, b += stb ) {}
  };
  
  template<typename A, typename B, typename C> inline
  void padd( A * a, const B * b, const C * c, unsigned sz, 
         int sta = 1, int stb = 1, int stc = 1  )
  { 
    for ( A * ea = a; ea != a+sz*sta; *ea = *b+*c, ea += sta, b += stb, c += stc ) ;
  };
  
  template<typename A, typename B> inline
  void psub( A * a, const B * b, unsigned sz, int sta = 1, int stb = 1  )
  { 
    for ( A * ea = a; ea != a+sz*sta; *ea -= *b, ea += sta, b += stb ) ; 
  };
  
  template<typename A, typename B, typename C> inline
  void psub( A * a, const B * b, const C * c, unsigned sz, int sta = 1, int stb = 1, int stc = 1  )
  { 
    for ( A * ea = a; ea != a+sz*sta; *ea = *b-*c, ea += sta, b += stb, c += stc ) ; 
  };
  
  template<typename A, typename B> inline
  void pmul( A * a, const B * b, int sz, int sta = 1, int stb = 1  )
  { 
    for ( A * ea = a; ea != a+sz*sta; *ea *= *b, ea += sta, b += stb ) ; 
  };
  
  template<typename A, typename B, typename C> inline
  void pmul( A * a, const B * b, const C * c, unsigned sz, int sta = 1, int stb = 1, int stc = 1  )
  { 
    for ( A * ea = a; ea != a+sz*sta; *ea = *b**c, ea += sta, b += stb, c += stc ) ; 
  };
  
  template<typename A, typename B> inline
  void pdiv( A * a, const B * b, unsigned sz, int sta = 1, int stb = 1  )
  { 
    for ( A * ea = a; ea != a+sz*sta; *ea /= *b, ea += sta, b += stb ) ; 
  };
  
  template<typename A, typename B, typename C> inline
  void pdiv( A * a, const B * b, const C * c, unsigned sz, int sta = 1, int stb = 1, int stc = 1  )
  { 
    for ( A * ea = a; ea != a+sz*sta; *ea = *b/(*c), ea += sta, b += stb, c += stc ) ; 
  };
  
  template<typename A, typename B> inline
  void scadd( A * a, const B& b, unsigned n, int s = 1 )
  { 
    for ( A * ea = a; ea != a + n*s; *a += b, ea += s ) ;
  };
  
  template<typename A, typename B, typename C> inline
  void scadd( A * a, const B * b, const C& c, unsigned n, int sa = 1, int sb = 1 )
  {
    for ( A * ea = a; ea != a + n*sa; *ea = *b + c, ea += sa, b += sb ) ;
  };
  
  template<typename A, typename B> inline
  void scsub( A * a, const B&  b, unsigned n, int s  = 1 )
  { 
    for ( A * ea = a; ea != a + n*s; *ea -= b, ea += s ) ;
  };
  
  template<typename A, typename B> inline
  void scsub( const B& b, A * a, unsigned n, int s = 1 )
  {
    for ( A * ea = a; ea != a+n*s;  *ea = b-*ea, ea += s ) ;
  }
  
  template<typename A, typename B, typename C> inline
  void scsub( A * a, const B * b, const C& c, unsigned n, int sa = 1, int sb = 1 )
  { 
    for ( A * ea = a; ea != a+n*sa; *ea = *b-c, ea += sa, b += sb ) ;
  }
  
  template<typename A, typename B, typename C> inline
  void scsub( A * a, const B& b, const C * c, unsigned n, int sa = 1, int sb = 1 )
  {
    for ( A * ea = a; ea != a+n*sa; *ea = c-*b, ea += sa, b += sb ) ;
  };
  
  template<typename A, typename B> inline
  void scmul( A * a, const B& b, unsigned n, int s = 1 )
  { 
    for ( A * ea = a; ea != a+n*s; *ea *= b, ea += s ) ;
  }
  
  template<typename A, typename B, typename C> inline
  void scmul( A * a, const B * b, const C& c, unsigned n, int sa = 1, int sb = 1 )
  {
    for ( A * ea = a; ea != a + n*sa; *ea = *b*c, ea += sa, b += sb ) ;
  }
  
  template<typename A, typename B, typename C> inline
  void scmul( A * a, const A * _a, const B * b, const C& c )
  {
    for ( A * ea = a; ea != _a; *ea++ = *b++*c ) ;
  }
  
  template<typename A, typename B> inline
  void scdiv( A * a, const B& b, unsigned n, int s = 1 )
  { 
    for ( A * ea = a; ea != a+n*s; *ea /= b, ea += s ) ;
  }
  
  template<typename A, typename B> inline
  void scdiv( const B& b, A * a, unsigned n, int s = 1 )
  {
    for ( A * ea = a; ea != a+n*s; *ea = b/(*ea), ea += s ) ;
  };
  
  template<typename A, typename B, typename C> inline
  void scdiv( A * a, const B * b, const C& c, unsigned n, int sa = 1, int sb = 1 )
  {
    for ( A * ea = a; ea != a + n*sa; *ea = *b/c, ea += sa, b += sb ) ;
  }
  
  template<typename A, typename B, typename C> inline
  void scdiv( A * a, const C& c, const B * b , unsigned n, int sa = 1, int sb = 1 )
  {
    for ( A * ea = a; ea != a + n*sa; *ea = c/(*b), ea += sa, b += sb ) ;
  }
  
  template<typename A, typename B> inline
  void copy ( A * a, const B * b, unsigned n )
  {
    std::copy(b,b+n,a) ;
  };
  
  template<typename A, typename B> inline 
  void copy ( A * a, const B * b, unsigned n, int sa, int sb )
  {
    for ( A * ea = a; ea != a + n*sa; *ea = *b, ea += sa, b += sb ) ;
  };
  
  template<class C >
  void fill( C * a, unsigned n, int s, const C& x )
  {
    for ( C * pa = a; pa != a+n*s; *pa = x, pa += s ) ;
  };

  template<typename A, typename B> inline
  void icopy( A * a, unsigned * aadd, unsigned nadd, const B * b )
  {
    for ( unsigned * ea = aadd; ea != aadd + nadd; a[*ea++] = *b++ ) ;
  };
  
  template<typename A, typename B> inline
  void ipadd( A * a, unsigned * aadd, unsigned nadd, const B * b )
  {
    for ( unsigned * ea = aadd; ea != aadd + nadd; a[*ea++] += *b++ ) ;
  };
  
  template<typename A, typename B> inline
  void ipsub( A * a, unsigned * aadd, unsigned nadd, const B * b )
  {
    for ( unsigned * ea = aadd; ea != aadd + nadd; a[*ea++] -= *b++ ) ;
  };
  
  /* min - max */
  template<typename A, typename B> inline 
  void accmin( A& mn, B const * v, unsigned n, int s = 1 )
  {
    const B * p;
    for ( p = v; p != v + n*s; p += s ) 
      if ( mn > *p ) mn = *p;
  };
  
  template<typename A, typename B> inline
  void min( A& mn, B const * v, unsigned n, int s = 1 )
  {
    mn = v[0];
    accmin(mn,v+s,n-1,s) ;
  };
  
  template<typename A, typename B> inline 
  void accmax( A& mx, const B  * v, unsigned n, int s = 1 )
  {
    const B * p;
    for (  p = v; p != v+n*s; p += s )
      if ( mx < *p ) mx = *p;
  };
  
  template<typename A, typename B> inline 
  void max( A& mx, B const * v, unsigned n, int s = 1 )
  {
    mx = v[0];
    accmax(mx,v,n,s) ;
  };
  template<typename A> inline
  void print( const A * p, unsigned n, int st = 1, 
          std::ostream& o = std::cout )
  {
    const A * i;
    o << "[ ";
    for ( i = p; i != p + (n-1)*st; i += st )
      o << *i << ", ";
    o << *i << " ]";
  };

  template<class C>
  void conv( C * r, 
         unsigned nr, int str, 
         const 
         C * a, 
         unsigned na, int sta, 
         const 
         C * b, 
         unsigned nb, int stb )
  {
    if ( na < nb ) { conv(r,nr,str,b,nb,stb,a,na,sta); return; };
    int k, l;
    C * er;
    const C * ea, * eb, * eeb, * eea;
    er = r; 
    eb = b;
    
    for ( er = r, eb = b, k = 0; k < std::min(nb,nr); k ++, er += str, eb += stb )
      for ( eeb = eb, ea = a, *er = 0, l = 0; 
        l <= k; 
        *er += *ea**eeb, l ++, eeb -= stb, ea += sta ) ;
    
    for ( ; k < std::min(na,nr) ; k ++,   er += str )
      for ( *er = 0, 
          l = k-nb+1, eeb = b +(nb-1)*stb , ea = a + l*sta; 
        l <= k; 
        *er+= *ea**eeb,
          l ++, eeb -= stb, ea += sta ) ;
    
    for ( ; k < nr; k ++,  er += str )
      for ( *er = 0, 
          l = k-nb+1, eeb = b+(nb-1)*stb,  ea = a + l*sta; 
        l < na; 
        *er += *ea**eeb, 
          l ++, eeb -= stb, ea += sta ) ;
  };
  
  template<class C>
  void accconv( C * r, int str, 
        const C * a, unsigned na, int sta, 
        const C * b, unsigned nb, int stb )
  {
    C * er;
    const C * ea, * eb;
    for ( ea = a; ea != a + na*sta; ea += sta, r += str )
      for ( er = r, eb = b; eb != b + nb*stb; eb += stb, er += str )
    *er += *ea**eb;
  };
  
  template<class C> 
  void conv( C * r, int str, 
         const C * a, unsigned na, int sta, 
         const C * b, unsigned nb, int stb )
  {
    C * er;
    for ( er = r; er != r + (na+nb-1)*str; *er++ = 0 ) ;
    accconv(r,str,a,na,sta,b,nb,stb);
  };

  template<class C> inline
  void minmax( C& m, C& M, const C * a, unsigned n, unsigned s = 1 )
  {
    for ( const C * e = a + n*s-s, m = M = *e; a < e; a += 2*s )
      {
    if ( *a > *(a+s) )
      {
        if ( *a > M ) M = *a; 
        if ( *(a+s) < m ) m = *(a+s) ;
      }
    else 
      {
        if ( *a > M ) M = *a;
        if ( *(a+s) < m ) m = *(a+s) ;
      };
      };
  };

  template<class C> inline
  unsigned mmerge( C * r, const C* a, const C* b, 
           unsigned na, unsigned nb  )
  {
    C * er(r);
    const C * ea(a), * eb(b);
    a += na;
    b += nb;
    for (;;)
      {
    while( ea < a && *ea < *eb ) *er++ = *ea++;
    if (  ea == a ) { std::copy(eb,b,er) ; break; };
    if ( *ea == *eb ) { *er++ = *ea++; eb++;};
    while( eb < b && *eb < *ea ) *er++ = *eb++;
    if (  eb == b ) { std::copy(ea,a,er); break; };
    if ( *eb == *ea ) { *er++ = *eb++; ea++;};
      };
    return er-r;
  };
  
  template<class C> inline
  void mmerge( C * r, int str, 
           const C * a, unsigned na, int sta, 
           const C * b, unsigned nb, int stb )
  {
    C * er(r);
    const C * ea(a), * eb(b);
    
    a += na*sta;
    b += nb*stb;
    
    for (;;)
      {
    for (;ea != a && *ea < *eb; *er = *ea, er += str, ea += sta ) ;
    
    if (  ea == a ) 
      { 
        for ( ; eb != b; *er = *eb, er += str, eb += stb ) ;
        break;
      };
    
    if ( *ea == *eb ) 
      { 
        *er = *ea;
        er += str;
        ea += sta;
        eb += stb;
      };
    
    for (;eb != b && *eb < *ea; *er = *eb, er += str, eb += stb ) ;

    if (  eb == b ) 
      { 
        for ( ; ea != a; *er = *ea, er += str, ea += sta ) ;
        break;
      };
    
    if ( *eb == *ea ) 
      { 
        *er = *eb;
        er += str;
        eb += stb;
        ea += sta;
      };
      };
    return (er - r)/str;
  };
  
  

  template < class C > static
  void shift (C * p, const C & c, unsigned n, int is = 1)
  {
    int j, k, s;
    for (s = n, j = 0; j <= s - 2; j++)
      for (k = s - 2; k >= j; p[k * is] += c * p[k * is + is], k--) ;
  };

  template < class C > static
    void scale (C * p, const C & s, unsigned n, int is = 1)
  {
    C pw = s;
    for (C * cp = p + is; cp != p + n * is; *cp *= pw, pw *= s, cp += is) ;
  };
  
template<class O, class C, class I> inline
void horner( O& o, C const * const c, unsigned n, const I& i, int s = 1  )
{ 
  o = O(0);
  for ( int p = (n-1)*s; p != 0; o +=  as<O>(c[p]), o *= i, p -= s ) ; 
  o+= as<O>(c[0]); 
};

template<typename C, typename I> inline typename texp::sup<C,I>::T 
horner( C const * const c, unsigned n, const I& i, int s = 1  )
{ 
  typename texp::ringof<C,I>::T o(0);
  horner(o,c,n,i,s); 
  return o; 
};

template<class O, class C, class I> inline
void hhorner( O& o, C const * const c, unsigned n, const I& i, const I&i0, int s = 1  )
{ 
  o = O(c[(n-1)*s]);
  for ( int p = (n-2)*s; p >= 0; o *= i, o +=  O(c[p])*i0, p -= s ) ; 
};

  
template<typename C, typename I> inline typename texp::sup<C,I>::T 
hhorner( C const * const c, unsigned n, const I& i,  const I&i0, int s = 1  )
{ 
  typename texp::ringof<C,I>::T o(0);
  hhorner(o,c,n,i,i0,s) ; 
  return o; 
};


template<typename C> inline
void diff( C * dst, C const * const src, unsigned sz, int st = 1 )
{ 
  for ( unsigned i = 0; i < sz-1; dst[i*st] = (i+1)*src[(i+1)*st], i++ ) ; 
};
  
template<typename C, typename I>
void decasteljau( C * c, unsigned n, const I& i, int s = 1 )
{
  C *ec, *p;
  for ( ec = c + (n-1)*s; ec != c; ec -= s )
    for ( p = c; p != ec; p += s )
      *p = (1.0-i)**p+i**(p+s) ;
};

template<class O, class C, class I> inline
void dhorner( O& p, O& dp, C const * const c, unsigned n, const I& t, int s = 1 )
{ 
  p = c[n-1], dp  = 0.0; 
  for ( int j = (n-2)*s; j >=0; dp=dp*t+p, p=p*t+c[j], j -= s ) ; 
};


template<class A>
void inverses( A * a, A * ea ) { for ( ;a != ea; *a++ = 1/(*a) ) ; };


    template<typename real_t>
    real_t min( real_t const * src, int sz, int st = 1 )
    {
      real_t r = src[0];
      for ( int p = st; p != sz*st; p += st ) if ( r > src[p] ) r = src[p];    
      return r;
    };
    
    template<typename real_t>
    real_t max( real_t const * src, int sz, int st = 1 )
    {
      real_t r = src[0];
      for ( int p = st; p != sz*st; p += st ) if ( r < src[p] ) r = src[p];    
      return r;
    };
    



};

} //namespace mmx

#endif