tensor_bernstein_fcts.hpp

Go to the documentation of this file.

#ifndef realroot_MPOLDSE_BERNSTEIN_C
#define realroot_MPOLDSE_BERNSTEIN_C
//====================================================================
namespace mmx {

namespace tensor
{
  template< class C >  inline void 
  binoms( bernstein<C>& mpl ) {
    binoms( mpl.begin(), mpl.env, binomials<C>::default_ ) ;
  }
  
  template < class C > inline void
  scale (bernstein < C > &mpl)
  {
   scale (mpl.begin (), mpl.env, binomials < C >::default_);
  };
  
  template < class C > inline void
  uscale (bernstein < C > &mpl)
  {
    uscale (mpl.begin (), mpl.env, binomials < C >::default_);
  };

  template < class C > inline void 
  convertb2m ( monomials<C>  & mpl ) { 
    convertb2m (mpl.begin (), mpl.env, binomials < C >::default_); 
  }
 
  template < class C > inline void
  convertm2b (bernstein < C > &mpl) {
    convertm2b (mpl.begin (), mpl.env, binomials < C >::default_);
  };
  
  template < class C , class V> inline void
  convertm2b (bernstein < C > &mpl, const V& bx) {
    convertm2b (mpl);
    for(unsigned i=0;i<bx.size();i+=2)
      restrict(mpl,i/2,bx[i],bx[i+1]);
  };
  
  template < class C > void
  elevate (bernstein < C > &r, const eenv & elev)
  {
    bernstein < C > tmp (elev);
    tensor::binoms (tmp);
    scale (r);
    tensor::conv(r, tmp);
    uscale (r);
  };

  template < class C > void
  add (bernstein < C > &mpl, const C & c)
  {
    vct::scadd (mpl.begin (), c, mpl.esz (), 1);
  };
  
  template < class C > void
  add (bernstein < C > &r, const bernstein < C > &a, const C & c)
  {
    if (r.env != a.env)
      add (r = a, c);
    else
      vct::scadd (r.begin (), a.begin (), c, r.esz (), 1);
  };

  template < class C > void
  add (bernstein < C > &r, const C & c, const bernstein < C > &a) {
    add(r,a,c);
  }

  template<class C> void
  rewrite( bernstein<C>& r, const eenv& newe )
  {
    if ( newe == r.env ) return;
    eenv elenv = eenv::elevation (newe, r.env);
    elevate (r, elenv) ;
  };
  
  template < class C > void
  add (bernstein < C > &r, const bernstein < C > &a)
  {
    eenv cm(eenv::common(r.env,a.env));
    if ( r.env != cm ) { rewrite(r,cm); };
    if ( a.env == cm ) 
      {
    vct::padd (r.begin (), a.begin (), r.esz (), 1, 1); 
      }
    else 
      {
    bernstein<C> tmp(a);
    rewrite(tmp,cm);
    vct::padd (r.begin (), tmp.begin (), r.esz (), 1, 1);
      };
  };


template < class C > void
add (bernstein < C > &r, const bernstein < C > &a, const bernstein < C > &b) {
  add (r=a, b);
}


template < class C > void
sub (bernstein < C > &r, const bernstein < C > &a) {
  eenv cm(eenv::common(r.env,a.env));
  if ( r.env != cm ) { rewrite(r,cm); };
  if ( a.env == cm ) 
    {
      vct::psub (r.begin (), a.begin (), r.esz (), 1, 1); 
    }
  else 
    {
      bernstein<C> tmp(a);
      rewrite(tmp,cm);
      vct::psub (r.begin (), tmp.begin (), r.esz (), 1, 1);
    };
}


template < class C > void
sub (bernstein < C > &r, const bernstein < C > &a, const bernstein < C > &b) {
  sub (r = a, b);
}

template < class C > void
sub (bernstein < C > &mpl, const C & c) {
  vct::scsub (mpl.begin (), c, mpl.esz (), 1);
}

template < class C > void
sub (bernstein < C > &mpl, const C & c, const bernstein < C > &a) {
  add(mpl,a,-c);
}

template < class C > void
sub (bernstein < C > &mpl, const bernstein < C > &a, const C& c) {
  add(mpl,a,-c);
}

template < class C > void
mul (bernstein < C > &r, const bernstein < C > &a, const bernstein < C > &b) {
  bernstein < C > a_ (a);
  scale (a_);
  bernstein < C > b_ (b);
  scale (b_);
  tensor::conv (r, a_, b_);
  uscale (r);
};

template < class C > void
mul (bernstein < C > &r, const bernstein < C > &a) {
  bernstein < C > a_;
  r.swap (a_);
  scale (a_);
  bernstein < C > b_ (a);
  scale (b_);
  tensor::conv (r, a_, b_);
  uscale (r);
};

template < class C > inline void
mul (bernstein < C > &r, const bernstein < C > &a, const C &c) {
  mul(r=a,c);
};

template < class C > inline void
mul (bernstein < C > &r, const C &c, const bernstein < C > &a) {
  mul(r=a,c);
};

template < class C > void
mul (bernstein < C > &r, const C &c) {
  vct::scmul(r.begin (),c, r.esz (), 1);
}

template < class C > inline void
div (bernstein < C > &r, const bernstein < C > &a, const C &c) {
  div(r=a,c);
}

template < class C > bernstein < C >::bernstein (const eenv & e): base_type(e) {};

template < class C > bernstein < C >::bernstein (const eenv & e, const C & c): 
  base_type(e,c){};

template < class C > 
bernstein < C >::bernstein (int nvr, const int *szs, const int *vrs): base_type(nvr,szs,vrs)
{
}

template < class C > bernstein < C >::bernstein (const C & x) : base_type(x) {}

template < class Result, class Coeff, class Parameters > inline 
void eval (Result & result, const bernstein < Coeff > & controls , const Parameters & parameters )
{
  levalb (result, controls.begin (), controls.env,parameters);
};

template < class Coeff> inline 
void eval (Coeff& result, const bernstein < Coeff > & controls , const Coeff & v)
{
  Coeff tmp[controls.env.sz()];
  std::copy (controls.begin(), controls.end(), tmp);
  decasteljau (tmp , controls.env.sz(), v, 1);
  result=tmp[0];
};

template < class Result, class Coeff, class Parameters > inline 
void eval (Result & result, const bernstein < Coeff > & controls , const Parameters & parameters , unsigned n) { 
  //assert(parameters.size()<n);
  levalb( result, controls.begin (), controls.env, parameters ); 
};

template < class Result, class Coeff, class Parameters > inline 
void subs0 (Result *result, const bernstein < Coeff > & controls , const Parameters & parameters )
{ 
  subs0b( result, controls.begin (), controls.env, parameters ); 
};

template < class C > void
casteljau (bernstein < C > &a, bernstein < C > &b, const C & t, int v)
{
  if (a.env != b.env)
    realloc (a, b.env);
  casteljau (a.begin (), b.begin (), t, v);
};

template < class C > void
casteljau (bernstein < C > &left, bernstein < C > &right,
       const bernstein < C > &source, const C & t, int v)
{
  casteljau (left, right = source, t, v);
};


  template < class C > void
  casteljau (bernstein < C > &a, bernstein < C > &b, int v)
  {
    if (a.env != b.env)
      realloc (a, b.env);
    casteljau (a.begin (), b.begin (), v);
  };
  
  template < class C > void
  casteljau (bernstein < C > &left, bernstein < C > &right,
         const bernstein < C > &source, int v)
  {
    casteljau (left, right = source, v);
  };

  template < class C >
  std::ostream & operator<< (std::ostream & o, const bernstein < C > &mpl)
  {
    o << "bernstein object:\n";
    o << "\t(addresses) o: " << &mpl << " e: " << mpl.env.data <<  "d: " << mpl.
      begin () << std::endl;
    o << "\t(eenv)" << mpl.env << std::endl;
    o << "\t(data)";
    array::print (o, mpl.data);
    o << std::endl;
    o << ":object bernstein\n";
    return o;
  };

  template<class C> void
  assign( monomials<C>& monoms, const bernstein<C>& controls )
  {
    monoms = controls;
    tensor::convertb2m(monoms);
  };

 
//   template < class C, class O, class R > void
//   convert (mpoldst < C, O, R > &mpol, const bernstein < C > &mpl)
//   {
//     monomials<C> tmp;
//     assign(tmp,mpl);
//     convert(mpol,tmp);
//   };

  template < class C > void
  cfdump (std::ostream & o, const bernstein < C > &mpl)
  {
    for (const C * a = mpl.begin (); a != mpl.end (); o << " " << *a++){}
  };

  template <class OSTREAM, class C > void
  mprint (OSTREAM & o, const bernstein < C > &mpl)
  {
    monomials < C > mp (0);
    convert (mp, mpl);
    o << mp;
    
  };

  template<class OSTREAM, class C > void 
  print(OSTREAM& o , const bernstein<C>& mpl ) { 
    monomials<C> mp(mpl);
    convertb2m(mp);
    print(o,mp);
  };

  template<class OSTREAM, class C > void inline 
  print(OSTREAM& o , const bernstein<C>& mpl , const variables& v) { 
    monomials<C> mp(mpl);
    convertb2m(mp);
    print(o,mp,v);
  };

  template < class C > bool bernstein < C >::operator== (const bernstein & mpl) const 
  {
    return base_type::operator==(mpl); 
  };

  template<class C> inline
  void diff( bernstein<C>& res, const bernstein<C>& src, int v )
  {
    clear( res, eenv::diff(src.env,v) );
    bdiff( res.begin(), src.env, src.begin(), v );
  };

  template<class C>
  bernstein<C>& bernstein<C>::operator=( const C& x ) 
  {
    this->~bernstein<C>();
    new (this) bernstein<C>(x);
    return *this;
  };

  template<class C,class T> //inline
  void split( bernstein<C>& l, bernstein<C>& r, int v, const T& t   )
  {
    if ( l.env != r.env ) { realloc(l,r.env); };
    bsplit( l.begin(), r.begin(), r.env, t, r.env.localv(v) );
  };


  template<class C> inline
  void split( bernstein<C>& l, bernstein<C>& r, int v  )
  {
    if ( l.env != r.env ) { realloc(l,r.env); };
    bsplit2( l.begin(), r.begin(), r.env, r.env.localv(v) );
  };

  template<class C,class T>
  void restrict( bernstein<C>& rst, int v, const T& a, const T& b ) { 
    v = rst.env.vridx(v);
    if (v < rst.env.nvr())
      brestrict( rst.begin(), rst.env, v, a, b ); 
  };
  
template<class C> inline
void lface ( bernstein<C>& f, const bernstein<C>& src, int v )
{
  lface( f.begin(), src.env, src.begin(), v );
};


template<class C> inline
void rface ( bernstein<C>& f, const bernstein<C>& src, int v )
{
  rface( f.begin(), src.env, src.begin(), v );
};


template<class C> inline
void face ( bernstein<C>& f, const bernstein<C>& src, int v, int n )
{
  assert(n==0||n==1);
  eenv fenv(face_env(src.env,src.env.localv(v)));
  if ( f.env != fenv ) realloc( f, fenv );
  if ( n ) rface( f, src, src.env.localv(v) );
  else     lface( f, src, src.env.localv(v) );
};


}


namespace let {
  template<class C, class U>  void assign(C& p, const U& q);

  template<class C, class U>
  void assign(tensor::bernstein<C>& p,  const tensor::bernstein<U>& q)
  {
    p.env=q.env;
    p.data = typename tensor::bernstein<C>::vector_type(q.data.size());
    for(unsigned i=0;i<q.data.size();i++)  assign(p.data[i],q.data[i]);
  }

  template<class C, class D> inline void 
  assign( tensor::bernstein<C>& b, const tensor::monomials<D>& m)
  {
    b = tensor::bernstein<C>(m.env);
    for(int i=0;i<m.size();i++) b[i] =as<C>(m[i]);
    tensor::convertm2b(b);
  }
  
  template<class C, class D, class O> inline void 
  assign( tensor::bernstein<C>& b, const sparse::monomial_seq<D,O>& m)
  {
    b = tensor::bernstein<C>(m);
  } 
  template<class C, class D, class O, class DOM> inline void 
  assign( tensor::bernstein<C>& b, const sparse::monomial_seq<D,O>& m, const DOM& dom)
  {
    b = tensor::bernstein<C>(m,dom);
  }

}//namespace let

} //namespace mmx
//====================================================================
#include "realroot/tensor_convert.hpp"
//====================================================================
#endif /* realroot_MPOLDSE_BERNSTEIN_C */