tensor_convert.hpp

Go to the documentation of this file.

#ifndef realroot_MPOLDSE_MODULE_H
#define realroot_MPOLDSE_MODULE_H
//====================================================================
#include <realroot/tensor_eenv.hpp> 
#include <realroot/tensor_vctops.hpp>
#include <realroot/binomials.hpp>
//====================================================================
namespace mmx {

namespace tensor {

  typedef tensor::eenv eenv;
  
  template<class MPLBASE, class C >
  void scadd (MPLBASE & mpl, const C & c)
  {
    vct::scadd (mpl.begin (), c, mpl.esz (), 1);
  };
  
  template<class MPLBASE, class MPLBASE1, class C >
  void scadd ( MPLBASE& r, const MPLBASE1 &a, const C & c)
  {
    if (r.env != a.env) 
      scadd (r = a, c);
    else 
      vct::scadd (r.begin (), a.begin (), c, r.esz ());
  };
  
  template<class MPLBASE0,class MPLBASE1>
    void conv( MPLBASE0 &r, const MPLBASE1 &a );

  template<class MPLBASE0, class MPLBASE1, class MPLBASE2 >
  void conv( MPLBASE0& r, const MPLBASE1& a, const MPLBASE2& b )
  {
    if (&a == &r) { conv (r, b); return; };
    if (&b == &r) { conv (r, a); return; };
    clear (r, eenv::mul (a.env, b.env));

    unsigned *oa = new unsigned[a.esz () + b.esz ()];
    unsigned *ob = oa + a.esz ();

    eenv::oaddress (r.env, oa, a.env);
    eenv::oaddress (r.env, ob, b.env);
    
    typedef typename MPLBASE0::value_type rreal;
    typedef typename MPLBASE1::value_type areal;
    typedef typename MPLBASE2::value_type breal;
    
    unsigned *ea;
    unsigned *eb;
    rreal * pr = r.begin();
    const areal * pa;
    const breal * pb;
    
    for (ea = oa, pa = a.begin (); ea != oa + a.esz (); ea++, pa++)
      for (eb = ob, pb = b.begin (); eb != ob + b.esz ();
       pr[*ea + *eb++] += (*pa) * (*pb++) ) ;
    delete[]oa;
  };
  
  template<class MPLBASE0,class MPLBASE1> void 
  conv( MPLBASE0 &r, const MPLBASE1 &a ) { 
    MPLBASE0 tmp;
    r.swap (tmp);
    conv (r, tmp, a);
  };
  
 
  template<class UPOL,class MPOL>
  UPOL convert(const MPOL & p, unsigned v)
  {
    unsigned N = degree(p,v)+1;
    UPOL r(N, AsSize());
    for(unsigned i =0; i<N;i++) { slice(r[i].rep(),p.rep(),v,i); }
    return r;
  }

} // namespace dense
//======================================================================
} //namespace mmx
#endif