tensor_eenv_loops.hpp

Go to the documentation of this file.

#ifndef realroot_tensor_eenvops_hpp
#define realroot_tensor_eenvops_hpp

#include <realroot/texp_sup.hpp>
#include <realroot/tensor_vctops.hpp>
#include <realroot/binomials.hpp>

namespace mmx {

namespace tensor
{
  template < class R, class C, class eenv, class P >
  void levalm (R & result, const C * data, const eenv & env, const P & p)
  {
    typedef R result_type;

    // query for the environment information

    const int esz = env.sz ();  // number of coefficients
    const int nvr = env.nvr (); // number of variables 
    const int *szs = env.szs ();    // sizes in each variables
    const int *str = env.str ();    // corresponding strides.

    //result_type tmp[env.sz ()], r;    // temporary array used by the inplace evaluation loop
    result_type* tmp = new result_type[env.sz()];
    result_type r;
    std::copy (data, data + env.sz (), tmp);
    for (int v = nvr - 1; v >= 0; v--)
      for (int i = 0; i < esz; i += str[v - 1]) {
    vct::horner(r,tmp + i, szs[v], p[v], str[v]);
    tmp[i] =  r;
      }
//    let::assign(result,tmp[0]);
      result= tmp[0];
  };

  template < class R, class C, class eenv, class P >
  void hevalm (R & result, const C * data, const eenv & env, const P & p)
  {
    typedef R result_type;

    // query for the environment information

    const int esz = env.sz ();  // number of coefficients
    const int nvr = env.nvr (); // number of variables 
    const int *szs = env.szs ();    // sizes in each variables
    const int *str = env.str ();    // corresponding strides.

    //result_type tmp[env.sz ()],r;   // temporary array used by the inplace evaluation loop
    result_type* tmp=new result_type[env.sz ()], r;
    std::copy (data, data + env.sz (), tmp);
    for (int v = nvr - 1; v >= 0; v--)
      for (int i = 0; i < esz; i += str[v - 1]) {
    vct::hhorner(r,tmp + i, szs[v], p[v+1],p[0],str[v]);
    tmp[i] =  r;
      }

    let::assign(result,tmp[0]);

  };

  template<typename real_t, typename value_t>
  void decasteljau( real_t * r, unsigned sz, const value_t& t, int str = 1 )
  {
    real_t *er, *p;
      for ( er = r + (sz-1)*str; er != r; er -= str )
    for ( p = r; p != er; p += str )
      *p = real_t(value_t(1)-t)* *p+t**(p+str);
  };  

  template < class R, class C, class eenv, class P >
  void levalb (R & result, const C * data, const eenv & env, const P & p)
  {
    //    typename sup<typename texp::value_type<P>::result_t, C >::T tmp[env.sz()];
    //XXX    typename texp::ringof_ < typename texp::value_type<P>::T,C >::T

    int szTemp =  env.sz();

    R tmp[szTemp];

    const int esz  (env.sz ());
    const int nvr  (env.nvr ());
    const int *szs (env.szs ());
    const int *str (env.str ());

    std::copy (data, data + env.sz (), tmp);

    for (int v = nvr - 1; v >= 0; v--)
      for (int i = 0; i < esz; i += str[v - 1])
    decasteljau (tmp + i, szs[v], p[v], str[v]);

    let::assign(result,tmp[0]);
  };

  

  template < class OutputIterator, class C, class eenv  >
  void maxs (OutputIterator _maxs_, const C * data, const eenv & env, int v)
  {
    const int *eszs = env.szs ();
    const int *estr = env.str ();
    //    const int envr = env.nvr ();
    int k, i; //, j;
    for ( k = 0; k < eszs[v]; k++, data += estr[v], _maxs_++)
      {
    *_maxs_ = *data;
    for (i = 0; i < estr[-1]; i += estr[v - 1])
      for (const C * edata = data + i; edata < data + i + estr[v]; edata++)
        if (*edata > *_maxs_)
          *_maxs_ = *edata;
      };
  };
  
  template < class OutputIterator, class C, class eenv >
  void mins (OutputIterator _mins_, const C * data, const eenv & env, int v)
  {
    const int *eszs = env.szs ();
    const int *estr = env.str ();
    //    const int envr = env.nvr ();
    int p, k, i, j;
    for (p = 0, k = 0; k < eszs[v]; k++, p += estr[v], _mins_++)
      {
    *_mins_ = data[p];
    for (i = 0; i < estr[-1]; i += estr[v - 1])
      for (j = i; j < i + estr[v]; j++)
        if (data[j + p] < *_mins_)
          *_mins_ = data[j + p];
      };
  };

  

  template < class C, class eenv >
  void vswap( C * data, const eenv & env, int *perm, int *supp = 0, int nsupp = 0  )
  {
    const int *szs = env.szs ();
    const int *str = env.str ();
    int ostr[env.nvr ()];
    int i;
    for (i = env.nvr () - 2, str[env.nvr () - 1] = 1; i >= -1;
     ostr[i] = szs[perm[i + 1]] * str[i + 1], i--) ;
    if (supp == 0)
      {
    int a, o, v;
    for (i = 0; i < env.esz (); i++)
      {
        for (a = i, v = env.nvr () - 1; a; a /= szs[v], v--)
          o += (a % szs[v]) * ostr[perm[v]];
        std::swap (data[i], data[o]) ;
      };
      }
    else
      {
    int a, o, v;
    for (i = 0; i < nsupp; supp[i++] = o)
      {
        for (a = supp[i], v = env.nvr () - 1; a; a /= szs[v], v--)
          o += (a % szs[v]) * ostr[perm[v]];
        std::swap (data[i], data[o]);
      };
      };
    env.vswap (perm);
  };
  

  template < class C, class eenv  >
  void binoms (C * data, const eenv & env, binomials < C > & binms)
  {
    const int *estr = env.str ();
    const int *eszs = env.szs ();
    const int envr = env.nvr ();
    const int esz = env.sz ();
    const C *bins = binms[eszs[envr - 1] - 1];
    for (C * pdata = data; pdata < data + esz; pdata += estr[envr - 2])
      vct::copy (pdata, bins, eszs[envr - 1]);
    for (int v = 0; v < envr - 1; v++)
      {
    bins = binms[eszs[v] - 1];
    //  int k = 0;
    for (int i = 0; i < esz; i += estr[v - 1])
      for (int j = 0; j < estr[v]; j++)
        vct::pmul (data + i + j, bins, eszs[v], estr[v], 1);
      };
  };
  
  template < class C >
  void ibinoms (C * data, const eenv & env, binomials < C > & binms)
  {
    binoms(data,env,binms);
    vct::inverses(data,data+env.sz());
  };
  
  template < class C, class eenv >
  void scale (C * data, const eenv & env, binomials < C > &binoms)
  {
    const int *estr = env.str ();
    const int *eszs = env.szs ();
    const int  envr = env.nvr ();
    const int  esz  = env.sz ();
    const C *  bins = binoms[eszs[envr - 1] - 1];
    
    for (C *  pdata = data; pdata < data + esz; vct::pmul (pdata, bins, eszs[envr - 1]), pdata += estr[envr - 2]) ;
    
    for (int v = 0; v < envr - 1; v++)
      {
    bins = binoms[eszs[v] - 1];
    //  int k = 0;
    for (int i = 0; i < esz; i += estr[v - 1])
      for (int j = 0; j < estr[v]; j++)
        vct::pmul (data + i + j, bins, eszs[v], estr[v], 1);
      };
  };
  
  template < class C, class eenv >
  void uscale (C * data, const eenv & env, binomials < C > &binoms)
  {
    const int *estr = env.str ();
    const int *eszs = env.szs ();
    const int envr = env.nvr ();
    const int esz = env.sz ();
    const C *bins = binoms[eszs[envr - 1] - 1];
    for (C * pdata = data; pdata < data + esz; pdata += estr[envr - 2])
      vct::pdiv (pdata, bins, eszs[envr - 1]);
    for (int v = 0; v < envr - 1; v++)
      {
    bins = binoms[eszs[v] - 1];
    //  int k = 0;
    for (int i = 0; i < esz; i += estr[v - 1])
      for (int j = 0; j < estr[v]; j++)
        vct::pdiv (data + i + j, bins, eszs[v], estr[v], 1);
      };
  };

  template<class C> inline
  void convertb2m( C * bzc, unsigned sz, int st, binomials<C>& binom ) {
    int i,k;
    for ( i = st; i != st*int(sz); i += st )
    for ( k = (sz-1)*st; k != i-st; k -= st )
      bzc[k] -= bzc[k-st];
    vct::pmul(bzc,binom[sz-1],sz,st,1);
  };
  

  template<class C> inline
  void convertm2b( C * bzc, unsigned sz, int st, binomials<C>& binoms ) {
    C *tmp = new C[sz];
    const C * bin;
    unsigned i,k,p,l;
    
    for ( p = 0, i = 0; i < sz; i++, p += st ) 
      { tmp[i] = bzc[p]; bzc[p] = C(0); };
    for ( p = 0, i = 0; i < sz; i++, p += st )
      for ( bin = binoms[sz-i-1], l = p, k = 0; k < sz-i; k++, l += st )
    bzc[l] += tmp[i]*bin[k];
    bin = binoms[sz-1];
    for ( i = 0, l = 0; i < sz; i ++, l += st  )
      bzc[l] /= bin[i];
  };
  
  template < class C, class eenv >
  void convertm2b (C * data, const eenv & e, binomials < C > &bins, int v)
  {
    const int *esz = e.szs ();
    const int *est = e.str ();
    int i, j;
    for (i = 0; i < e.sz (); i += est[v - 1])
      for (j = i; j < i + est[v]; convertm2b (data + j, esz[v], est[v], bins), j++) ;
  };
  
  template < class C, class eenv >
  void convertm2b (C * data, const eenv & env, binomials < C > &bins)
  {
    for (int v = 0; v < env.nvr (); convertm2b (data, env, bins, v++)) ;
  };
  
  template < class C, class eenv >
  void convertb2m (C * data, const eenv & e, binomials < C > &bins, int v)
  {
    const int *esz = e.szs ();
    const int *est = e.str ();
    int i, j;
    for (i = 0; i < e.sz (); i += est[v - 1])
      for (j = i; j < i + est[v];
       convertb2m (data + j, esz[v], est[v], bins), j++) ;
  };
  
  template < class C, class eenv  >
  void convertb2m (C * data, const eenv & env, binomials < C > &bins)
  {
    for (int v = 0; v < env.nvr (); convertb2m (data, env, bins, v++)) ;
  };





  template<typename C> inline 
  void diff( C * b, const C * a, unsigned sz, int sb, int sa )
  {
    assert(sa==sb);
    for ( unsigned i = 1; i < sz; i ++ ) { b[(i-1)*sb] = i*a[i]; };
    //    for ( unsigned i = 1, a += sa; i < sz; *b = *a*i++, b += sb, a += sa ) ; 
  };
  
  template<typename C> inline /* differentiation */
  void diff( C * b, C * e, C const * a, int s )
  { 
    unsigned i;
    for ( i = 1, a += s; b != e; *b = *a*i++, b += s, a += s ) ;
  };

  template<typename C> inline /* differentiation */
  void diff( C * b, C * e, C const * a )
  { 
    for ( int i = 1; b != e; *b++ = *(++a)*i++ ) ;
  };


  template<typename C>
  void m_diff( C * op, const eenv& o, C const * ap, const eenv& a, int v )
  {
    typedef int sz_t;
    sz_t is = 0;
    for ( sz_t i = 0; i < o.sz(); i += o.st(v-1), is += a.st(v-1) )
      for ( sz_t j = 0; j < a.st(v); j ++ )
    diff(op+i+j, ap+is+j,a.sz(v), o.st(v), a.st(v)) ;
  };

  template<typename C>
  void mdiff( C * dst, const eenv& a, C const * const src, const eenv& esrc, int v )
  {
    const int   ssz  = esrc.sz ();
    //    const int * sszs = a.szs();
    const int * sstr = esrc.str();
    int rst = sstr[v-1]-sstr[v];
    unsigned i;
    C * ed, * fd;
    const C * es, * fs;
    for ( ed = dst, es = src; es < src+ssz; ed += rst, es += sstr[v-1] )
      for ( i = 1, fd = ed, fs = es+sstr[v]; fs != es+sstr[v-1]; fd += sstr[v], fs += sstr[v], i ++  )
    vct::scmul(fd,fd+sstr[v],fs,C(i));
  }

  
  
  template<typename C>
  static void bdiff( C * dst, const eenv& a, C const * const src, int v )
  {
    const int   ssz  = a.sz ();
    const int * sszs = a.szs();
    const int * sstr = a.str();
    int rst = sstr[v-1]-sstr[v];
    //    unsigned i;
    C * ed; const C * es;
    for ( ed = dst, es = src; es < src+ssz; ed += rst, es += sstr[v-1] )
      {
    C * fd; const C * fs;
    for ( fs = es, fd = ed; fs < es + sstr[v]; fs++, fd ++ )
      {
        C * dp; const C * p;
        for ( p = fs, dp = fd; p < fs+rst; p += sstr[v], dp += sstr[v] )
          *dp = (*(p+sstr[v])-*p)*(sszs[v]-1);
      };
      };
  }
  
  template<class C, class T> inline 
  void bsplit_uloop( C * r, C * rlast, C * l, const T& t, int st )
  {
    C * er, * p;
    for ( er = rlast - st ; er != r; er -= st, l += st )
      {
    *l = *r;
    for ( p = r; p != er; p += st )
      *p = (1.0-t)**p+t**(p+st);
      };
    *l = *r;
  };
  
  template < class C, class eenv, class T> inline
  void bsplit (C * l, C * r, const eenv & env, const T & t, int v )
  {
    //  const int *sszs = env.szs ();
    const int *sstr = env.str ();
    //  int i, j, k;
    C *er,*el;
    for ( el = l, er = r; el != l + sstr[-1]; el += sstr[v-1], er += sstr[v-1] )
      {
    C *fr,*fl;
    for ( fl = el, fr = er; fl != el + sstr[v]; fl++, fr++ )
      bsplit_uloop(fr,fr+sstr[v-1],fl,t,sstr[v]) ;
      };
  };
  
  
  template<class C> inline 
  void bsplit2_uloop( C * r, C * rlast, C * l,  int st )
  {
    C * er, * p;
    for ( er = rlast - st ; er != r; er -= st, l += st )
      {
    *l = *r;
    for ( p = r; p != er; p += st )
      *p = (*p+*(p+st))/C(2);
      };
    *l = *r;
  };

  template < class C, class eenv > 
  void bsplit2 (C * l, C * r, const eenv & env, int v )
  {
    //    const int *sszs = env.szs ();
    const int *sstr = env.str ();
    C *er,*el;
    for ( el = l, er = r; el != l + sstr[-1]; el += sstr[v-1], er += sstr[v-1] )
      {
    C *fr,*fl;
    for ( fl = el, fr = er; fl != el + sstr[v]; fl++, fr++ )
      bsplit2_uloop(fr,fr+sstr[v-1],fl,sstr[v]);
      };
  };

  
  template<class C,class T> inline 
  void brestrict_uloop( C * r, C * rlast, int st, const T& a, const T& b  )
  {
    C * er, * p;
    for ( er = rlast - st ; er != r; er -= st )
      for ( p = r; p != er; p += st )
    *p = (a**p+b**(p+st));
  };
  


  template<class C, class T> inline
  void brestrictLR  (C * s, const eenv & env, int v, const T& a, const T& b, const T& c, const T& d  )
  {
    //    std::cout << a << " " << b << std::endl;
    //    std::cout << c << " " << d << std::endl;
    //    const int *sszs = env.szs ();
    const int *sstr = env.str ();
    //    int i, j, k;
    for ( C * es = s; es != s + sstr[-1]; es += sstr[v-1] )
      for ( C * fs = es; fs != es + sstr[v]; fs ++ )
    {
      brestrict_uloop(fs,fs+sstr[v-1],sstr[v],a,b);
      brestrict_uloop(fs+sstr[v-1]-sstr[v],fs-sstr[v],-sstr[v],c,d);
    };
  };

  template<class C, class T> inline
  void brestrictRL  (C * s, const eenv & env, int v, const T& a, const T& b, const T& c, const T& d  )
  {
    //    std::cout << a << " " << b << std::endl;
    //    std::cout << c << " " << d << std::endl;
   
    

    //    const int *sszs = env.szs ();
    const int *sstr = env.str ();
    //    int i, j, k;
    for ( C * es = s; es != s + sstr[-1]; es += sstr[v-1] )
      for ( C * fs = es; fs != es + sstr[v]; fs ++ )
    {
      brestrict_uloop(fs+sstr[v-1]-sstr[v],fs-sstr[v],-sstr[v],c,d);    
      brestrict_uloop(fs,fs+sstr[v-1],sstr[v],a,b);
    };
  };

  template<class C, class eenv, class T > inline
  void  brestrict( C * src, const eenv& env, int v, const T& a, const T& b )
  {
    brestrictLR   ( src, env, v, T(1-a), a, T((b-a)/(T(1)-a)), T((T(1)-b)/(T(1)-a)) );
  };


  template<class C>
  void lface( C * dst, const eenv& a, const C * src, int v )
  {
    const int   ssz  = a.sz ();
    //    const int * sszs = a.szs();
    const int * sstr = a.str();

    C * ed; const C * es;
    const C * fs;
    for ( ed = dst, es = src; es < src+ssz; es += sstr[v-1] )
      for ( fs = es; fs < es + sstr[v]; *ed++ = *fs++ ) ;
  };

  template<class C>
  void rface( C * dst, const eenv& a, const C * src, int v )
  {
    const int   ssz  = a.sz ();
    //    const int * sszs = a.szs();
    const int * sstr = a.str();

    C * ed; const C * es;
    const C * fs;
    es = src +  sstr[v-1]-sstr[v];
    for ( ed = dst; es < src+ssz; es += sstr[v-1] )
      {
    for ( fs = es; fs < es + sstr[v]; *ed++ = *fs++ ) ;
      };
  };

  template<class C>
  void slice( C * dst, const eenv& a, const C * src, int v, int n)
  {
    const int   ssz  = a.sz ();
    const int * sstr = a.str();

    C * ed; const C * es;
    const C * fs;
    es = src +  sstr[v-1]-sstr[v];
    for ( ed = dst; es < src+ssz; es += sstr[v-1] )
      {
    //  fs = es + sstr[v-1] - sstr[v];
    for ( fs = es; fs < es + sstr[v]; *ed++ = *fs++ ) ;
      };
  };

  template<class C, class X, class O, class eenv >
  void mpolfill( C * data, const sparse::monomial_seq<X,O>& mpol, const eenv& env )
  {
    const int *vr = env.vrs ();
    const int *st = env.str ();
    const int  nvr = env.nvr();
    typename sparse::monomial_seq<X,O>::const_iterator it;
    C c;
    for (it =mpol.begin (); it != mpol.end (); it++)
      {
    int a = 0;
    for (int v = 0; v < nvr; a += (*it)[vr[v]] * st[v], v++) ;
    let::assign(c,(*it).coeff ());
    data[a] += c;
      };
  }
}

} //namespace mmx
#endif