include/realroot/bernstein_eenv_multiply.hpp

Go to the documentation of this file.

/********************************************************************
 *   This file is part of the source code of the realroot library.
 *   Author(s): J.P. Pavone, GALAAD, INRIA
 *   $Id: eenv-multiply.hpp,v 1.1 2005/07/11 10:03:55 jppavone Exp $
 ********************************************************************/
/* densite exacte / statistique des coefficients non-nuls */  
template< typename real_t >
real_t density( real_t * src, const real_t prec = 1e-12, unsigned nsmp = 0 )
{
  unsigned c = 0;
  if ( nsmp == 0 )
    {
      nsmp = data_size();
      for ( unsigned i = 0; i < nsmp; i ++ )
        if ( std::abs(src[i]) > prec ) c++;
    }
  else
    {
      for ( unsigned i = 0; i < nsmp; i ++ )
        if ( std::abs(src[rand()%data_size()] > prec) ) c++;
    };
  return ((real_t)c)/nsmp;
};
  
/* calcul le support des coefficients non-nuls dans src */
template<class real_t>
unsigned support( unsigned * supp, real_t * src, const real_t& prec = 1e-6 )
{
//  supp = new unsigned[ data_size() ];
  unsigned      c = 0;
  for ( unsigned i = 0; i < data_size(); i ++ )
    if ( std::abs(src[i]) > prec ) supp[c++] = i;
  return c;
};

/* transforme l'environement en l'environement de a*b */
void op_mul( eenv * a, eenv * b )
{
  if ( nvars() ) this->~eenv();
  _alloc_( a->nvars()+b->nvars() );
  std::copy( a->m_vrs, a->m_vrs+a->nvars(), m_vrs );
  std::copy( b->m_vrs, b->m_vrs+b->nvars(), m_vrs + a->nvars() );
  m_nvr = vctops::set_conversion( m_vrs, m_nvr );
  for ( sz_t v = 0; v < m_nvr; v ++ )
    {
      int av = a->indexofvar(var(v));
      int bv = b->indexofvar(var(v));
      m_szs[ v ] = ((av!=-1)?(a->sz(av)-1):0) + ((bv!=-1)?(b->sz(bv)-1):0) + 1;
    };
  m_mxvr = std::max(a->m_mxvr,b->m_mxvr);
  compute_strides();
};

/* mapping des noms i->o */
static  void vmap( sz_t * vmap, bernstein::eenv * o, bernstein::eenv * i )
{
  for ( sz_t v = 0; v < i->m_nvr; v ++ )
    {
      vmap[v] = -1;
      for ( sz_t k = 0; k < o->m_nvr; k ++ )
        if ( o->m_vrs[k] == i->m_vrs[v] ) vmap[v] = k;
      assert(vmap[v]!=-1);
    };
};
  
static void oaddress( bernstein::eenv * oenv, unsigned * osupp, 
                      bernstein::eenv * ienv, unsigned * isupp = 0, unsigned nsp = 0 )
{
  int v;
  unsigned c;
  sz_t vmap_[ ienv->m_nvr ];
  unsigned addi,addo;
  vmap( vmap_, oenv, ienv );
  if ( isupp )
    for ( c = 0; c < nsp; osupp[c] = addo, c ++ )
      for ( addi = isupp[c], addo = 0, v = ienv->m_nvr-1; addi; addi /= ienv->m_szs[v], v-- )
        addo += (addi%ienv->m_szs[v])*oenv->m_str[vmap_[v]];
  else
    {
      nsp = ienv->data_size();
    for ( c = 0; c < nsp; osupp[c] = addo, c ++ )
      for ( addi = c, addo = 0, v = ienv->m_nvr-1; addi; addi /= ienv->m_szs[v], v -- )
        addo += (addi%ienv->m_szs[v])*oenv->m_str[vmap_[v]];
    };
};

template< class real_t, char op > // multiplication multivariee dans la base monomiale
inline static void _mvrcvloop_( bernstein::eenv * oenv, real_t * dst, 
                                bernstein::eenv * aenv, real_t * sra,
                                bernstein::eenv * benv, real_t * srb,
                                unsigned * asupp = 0, unsigned nas = 0, 
                                unsigned * oasup = 0, 
                                unsigned * bsupp = 0, unsigned nbs = 0, 
                                unsigned * obsup = 0 )
{
  /* calcul des addresses de sorties            */
  unsigned * loasup, * lobsup; loasup = lobsup = 0;
  /* calcul des addresses de sorties de a sur o */
  if ( !oasup ) { if ( !asupp ) nas = aenv->data_size(); loasup = new unsigned[ nas ]; oasup = loasup; oaddress( oenv, oasup, aenv, asupp, nas ); };
  /* calcul des addresses de sorties de b sur o */
  if ( !obsup ) { if ( !bsupp ) nbs = benv->data_size(); lobsup = new unsigned[ nbs ]; obsup = lobsup; oaddress( oenv, obsup, benv, bsupp, nas ); };

  unsigned ia, ib, out; real_t   ca;
  switch( (asupp!=0)<<1|(bsupp!=0) )
    {
    case 0:
      assert(nas==aenv->data_size() && nbs == benv->data_size()); 
      if ( op == '+' )
        {
          for ( ia = 0; ia < nas; ia ++ )
            for ( ib = 0; ib < nbs;  ib ++ )
              dst[oasup[ia]+obsup[ib]] += sra[ia]*srb[ib];
          break;
        }
      if ( op == '-' )
        {
          for ( ia = 0; ia < nas; ia ++ )
            for ( ib = 0; ib < nbs; out += obsup[ib], ib ++ )
              dst[oasup[ia]+obsup[ib]] -= sra[ia]*srb[ib];
          break;
        };
      break;
    case 1:
      assert(nas == aenv->data_size() && bsupp != 0);
      if ( op == '+' )
        {
          assert(nas==aenv->data_size() && nbs == benv->data_size());
          for ( ia = 0; ia < nas; ia ++ )
            for (  ib = 0; ib < nbs; ib ++ )
              dst[oasup[ia]+obsup[ib]] += sra[ia]*srb[bsupp[ib]];
          break;
        }
      if ( op == '-' )
        {
          for ( ia = 0; ia < nas; ia ++ )
            for ( ib = 0; ib < nbs; ib ++ )
              dst[oasup[ia]+obsup[ib]] -= sra[ia]*srb[bsupp[ib]];
          break;
        };
      break;
    case 2:
      assert(nbs == benv->data_size() && asupp != 0);
      if ( op == '+' )
        {
          for ( ia = 0; ia < nas; ia ++ )
            for ( ib = 0; ib < nbs; ib ++ )
              dst[oasup[ia]+obsup[ib]] += sra[asupp[ia]]*srb[ib];
          break;
        }
      if ( op == '-' )
        {
          for ( ia = 0; ia < nas; ia ++ )
            for ( ib = 0; ib < nbs;  ib ++ )
              dst[oasup[ia]+obsup[ib]] -= sra[asupp[ia]]*srb[ib];
          break;
        };
      break;
    case 3:
      assert(asupp != 0 && bsupp != 0);
      if ( op == '+' )
        {
          for ( ia = 0; ia < nas; ia ++ )
            for ( ib = 0; ib < nbs;  ib ++ )
              dst[oasup[ia]+obsup[ib]] += sra[asupp[ia]]*srb[bsupp[ib]];
          break;
        }
      if ( op == '-' )
        {
          assert(nas==aenv->data_size() && nbs == benv->data_size());
          for ( ia = 0; ia < nas; ia ++ )
            for (  ib = 0; ib < nbs;  ib ++ )
              dst[oasup[ia]+obsup[ib]] -= sra[asupp[ia]]*srb[asupp[ib]];
          break;
        };
      break;
    };
  /* liberation des addresses de sorties   */
  if ( loasup ) delete[] loasup; if ( lobsup ) delete[] lobsup;
};


// template< typename real_t >
// static void sparse_multiply( bernstein::eenv * renv, real_t * rdata,
//                           bernstein::eenv * aenv, unsigned  nasupp, unsigned * asupp, real_t * adata,
//                           bernstein::eenv * benv, unsigned  nbsupp, unsigned * bsupp, real_t * bdata )
// {
  
//   unsigned * osuppa = new unsigned[ nasupp ];
//   unsigned * osuppb = new unsigned[ nbsupp ];
//   this-osupp( renv, osuppa, aenv, asupp, nasupp );
//   this->osupp( renv, osuppb, benv, bsupp, nbsupp );
//   for ( unsigned i = 0; i < nasupp; i ++ )
//     for ( unsigned j = 0; j < nbsupp; j ++ )
//       rdata[osuppa[i]+osuppb[j]] += adata[asupp[i]] * bdata[bsupp[j]];
//   delete[] osuppa;
//   delete[] osuppb;
// }

template<typename real_t>
static void monoms_multiply_loop( eenv * o, eenv * a, eenv *  b, real_t * dst, real_t * src_a, real_t * src_b  )
{ _mvrcvloop_<real_t,'+'>( o, dst, a, src_a, b, src_b ); };

template<typename real_t>
static void monoms_multiply( bernstein::eenv * res, bernstein::eenv * a, bernstein::eenv  *b, 
                             real_t * dst, real_t * src_a, real_t * src_b )
{ res->op_mul(a,b); monoms_multiply_loop( res, a, b, dst, src_a, src_b  ); };


template< typename real_t > // multiplication vectorielle
static void pmmul( eenv * oenv, real_t ** dst,  eenv * aenv, real_t ** src_a, eenv * benv, real_t ** src_b, int n, 
                   bool clear = true, 
                   unsigned * asupp  = 0, unsigned nas = 0,
                   unsigned * oasup  = 0, 
                   unsigned * bsupp  = 0, unsigned nbs = 0, 
                   unsigned * obsup  = 0 )
{
  if ( clear ) for (int d=0;d<n;std::fill(dst[d],dst[d]+oenv->data_size(),(real_t)0.0),d++);
  unsigned * loasup, * lobsup; loasup = lobsup = 0;
  /* calcul des addresses de sorties de a sur o */
  if ( !oasup ) { if ( !asupp ) nas = aenv->data_size(); loasup = new unsigned[ nas ]; oasup = loasup; oaddress( oenv, oasup, aenv, asupp, nas ); };
  /* calcul des addresses de sorties de b sur o */
  if ( !obsup ) { if ( !bsupp ) nbs = benv->data_size(); lobsup = new unsigned[ nbs ]; obsup = lobsup; oaddress( oenv, obsup, benv, bsupp, nbs ); };
  for ( unsigned i = 0; i < n; i ++ )
    _mvrcvloop_<real_t,'+'>(oenv,dst[i],aenv,src_a[i],benv,src_b[i],asupp,nas,oasup,bsupp,nbs,obsup);
  if ( loasup ) delete[] loasup; if ( lobsup ) delete[] lobsup;
};

template< typename real_t > // multiplication simple 
static inline void mmul( eenv * oenv, real_t * dst,  eenv * aenv, real_t * src_a, eenv * benv, real_t * src_b,
                         bool clear = true, 
                         unsigned * asupp  = 0, unsigned nas = 0,
                         unsigned * oasupp = 0, 
                         unsigned * bsupp  = 0, unsigned nbs = 0, 
                         unsigned * obsupp = 0 )
{ pmmul( oenv, &dst, aenv, &src_a, benv, &src_b, 1, clear, asupp, nas, oasupp, bsupp, nbs, obsupp ); };

template< typename real_t > // multiplication vectorielle creuse (calcul du support + multiplication)
static void spmmul( eenv * oenv, real_t ** dst,  eenv * aenv, real_t ** src_a, eenv * benv, real_t ** src_b, int n, 
                          bool clear = true, bool stats = false )
{
  unsigned nas,nbs;
  unsigned * assup = new unsigned[ aenv->data_size() ];
  unsigned * bssup = new unsigned[ benv->data_size() ];
  assert( assup && bssup );
  for ( int i = 0; i < n; i ++ )
    {
      nas = aenv->support(assup,src_a[i]);
      nbs = benv->support(bssup,src_b[i]);
      mmul( oenv, dst[i], aenv, src_a[i], benv, src_b[i], true, assup, nas, 0, bssup, nbs, 0 );
    };
  delete[] assup;
  delete[] bssup;
};

template< typename real_t > // multiplication simple creuse (calcul du support + multiplication)
static void spmmul( eenv * oenv, real_t * dst,  eenv * aenv, real_t * src_a, eenv * benv, real_t * src_b, 
                          bool clear = true, bool stats = false )
{
  spmmul( oenv, &dst, aenv, &src_a, benv, &src_b, 1, clear, stats );
};

template<typename real_t>
static void mcrossp( eenv * oenv, real_t ** dst, eenv * aenv, real_t ** src_a, eenv * benv, real_t ** src_b,
                     bool clear = true,
                     unsigned * asupp = 0, unsigned nas = 0, 
                     unsigned * oasup = 0, 
                     unsigned * bsupp = 0, unsigned nbs = 0, 
                     unsigned * obsup = 0 )


{
  if ( clear ) for (int d=0;d<3;std::fill(dst[d],dst[d]+oenv->data_size(),(real_t)0.0),d++);
  unsigned * loasup, * lobsup; loasup = lobsup = 0;
  /* calcul des addresses de sorties de a sur o */
  if ( !oasup ) { if ( !asupp ) nas = aenv->data_size(); loasup = new unsigned[ nas ]; oasup = loasup; oaddress( oenv, oasup, aenv, asupp, nas ); };
  /* calcul des addresses de sorties de b sur o */
  if ( !obsup ) { if ( !bsupp ) nbs = benv->data_size(); lobsup = new unsigned[ nbs ]; obsup = lobsup; oaddress( oenv, obsup, benv, bsupp, nbs ); };
  _mvrcvloop_<real_t,'+'>(oenv,dst[2],aenv,src_a[0],benv,src_b[1],asupp,nas,oasup,bsupp,nbs,obsup);
  _mvrcvloop_<real_t,'-'>(oenv,dst[2],aenv,src_a[1],benv,src_b[0],asupp,nas,oasup,bsupp,nbs,obsup);
  _mvrcvloop_<real_t,'+'>(oenv,dst[0],aenv,src_a[1],benv,src_b[2],asupp,nas,oasup,bsupp,nbs,obsup);
  _mvrcvloop_<real_t,'-'>(oenv,dst[0],aenv,src_a[2],benv,src_b[1],asupp,nas,oasup,bsupp,nbs,obsup);
  _mvrcvloop_<real_t,'+'>(oenv,dst[1],aenv,src_a[2],benv,src_b[0],asupp,nas,oasup,bsupp,nbs,obsup);
  _mvrcvloop_<real_t,'-'>(oenv,dst[1],aenv,src_a[0],benv,src_b[2],asupp,nas,oasup,bsupp,nbs,obsup);
  if ( loasup ) delete[] loasup; if ( lobsup ) delete[] lobsup;
};
 

// template<typename real_t>
// static void spmcrossp( eenv * oenv, real_t ** dst, eenv * aenv, real_t ** src_a, eenv * benv, real_t ** src_b,
//                     bool clear = true, bool stats = false )
// {
//   if ( clear ) for (int d=0;d<3;std::fill(dst[d],dst[d]+oenv->data_size(),(real_t)0.0),d++);
//   unsigned * asupp = new unsigned[ a->data_size() ];
//   unsigned * bsupp = new unsigned [b->data_size()];
//   unsigned * loasup, * lobsup; loasup = lobsup = 0;
//   /* calcul des addresses de sorties de a sur o */
//   if ( !oasup ) { if ( !asupp ) nas = aenv->data_size(); loasup = new unsigned[ nas ]; oasup = loasup; oaddress( oenv, oasup, aenv, asupp, nas ); };
//   /* calcul des addresses de sorties de b sur o */
//   if ( !obsup ) { if ( !bsupp ) nbs = benv->data_size(); lobsup = new unsigned[ nbs ]; obsup = lobsup; oaddress( oenv, obsup, benv, bsupp, nbs ); };
//   _mvrcvloop_<real_t,'+'>(oenv,dst[2],aenv,src_a[0],benv,src_b[1],asupp,nas,oasup,bsupp,nbs,obsup);
//   _mvrcvloop_<real_t,'-'>(oenv,dst[1],aenv,src_a[0],benv,src_b[2],asupp,nas,oasup,bsupp,nbs,obsup);
//   _mvrcvloop_<real_t,'-'>(oenv,dst[2],aenv,src_a[1],benv,src_b[0],asupp,nas,oasup,bsupp,nbs,obsup);
//   _mvrcvloop_<real_t,'+'>(oenv,dst[0],aenv,src_a[1],benv,src_b[2],asupp,nas,oasup,bsupp,nbs,obsup);
//   _mvrcvloop_<real_t,'-'>(oenv,dst[0],aenv,src_a[2],benv,src_b[1],asupp,nas,oasup,bsupp,nbs,obsup);
//   _mvrcvloop_<real_t,'+'>(oenv,dst[1],aenv,src_a[2],benv,src_b[0],asupp,nas,oasup,bsupp,nbs,obsup);
//   if ( loasup ) delete[] loasup; if ( lobsup ) delete[] lobsup;
// };


template<typename real_t> // compatibilite ancien code
static void mcrossp( eenv * res, eenv * a, eenv * b, real_t ** dst, real_t ** src_a, real_t ** src_b )
{ mcrossp( res, dst, a, src_a, b, src_b ); };

template<typename real_t>// compatibilite ancien code
static void monoms_crossprod( eenv * res, eenv * a, eenv * b, 
                              real_t ** dst, 
                              real_t ** src_a, real_t ** src_b )
{
  res->op_mul(a,b);
  for ( int d = 0; d < 3; d ++ ) dst[d] = new real_t[ res->data_size() ];
  mcrossp( res, dst, a, src_a, b, src_b );
};


template<typename real_t>
static void mdotp( eenv * oenv, real_t * dst, eenv * aenv, real_t ** src_a, eenv * benv, real_t ** src_b, int n,
                   bool clear = true,
                   unsigned * asupp = 0, unsigned nas = 0, 
                   unsigned * oasup = 0, 
                   unsigned * bsupp = 0, unsigned nbs = 0, 
                   unsigned * obsup = 0 )


{
  if ( clear ) std::fill(dst,dst+oenv->data_size(),(real_t)0.0);
  unsigned * loasup, * lobsup; loasup = lobsup = 0;
  /* calcul des addresses de sorties de a sur o */
  if ( !oasup ) { if ( !asupp ) nas = aenv->data_size(); loasup = new unsigned[ nas ]; oasup = loasup; oaddress( oenv, oasup, aenv, asupp, nas ); };
  /* calcul des addresses de sorties de b sur o */
  if ( !obsup ) { if ( !bsupp ) nbs = benv->data_size(); lobsup = new unsigned[ nbs ]; obsup = lobsup; oaddress( oenv, obsup, benv, bsupp, nas ); };
  for ( int i = 0; i < n; i ++ )
    _mvrcvloop_<real_t,'+'>(oenv,dst,aenv,src_a[i],benv,src_b[i],asupp,nas,oasup,bsupp,nbs,obsup);
  if ( loasup ) delete[] loasup; if ( lobsup ) delete[] lobsup;
};

template<typename real_t> // compatibilite ancien code
static void monoms_dotprod( eenv * res, eenv * a, eenv * b, real_t *& dst, real_t ** src_a, real_t ** src_b, int n )
{
  res->op_mul(a,b);
  dst = new real_t[res->data_size()];
  mdotp( res, dst, a, src_a, b, src_b, n );
};


template< class real_t >
static void rvbinoms(  bernstein::eenv * ienv, real_t * bcff, 
                       unsigned * isupp = 0, unsigned nsp = 0, bernstein::bzenv<real_t> * bznv = bzenv<real_t>::_default_ )
{
  real_t cff;
  unsigned c,addi;
  int v;
  if ( isupp )
    for ( c = 0; c < nsp; bcff[c] = cff, c ++ )
      for ( cff = (real_t)1.0, addi = isupp[c],  v = ienv->m_nvr-1; addi; addi /= ienv->m_szs[v], v-- )
        cff /= bznv->binomial(ienv->m_szs[v]-1,addi % ienv->m_szs[v]);
  else
    for ( c = 0; c < ienv->data_size(); bcff[c] = cff,  c ++ )
      for ( addi = c, v = ienv->m_nvr-1; addi; addi /= ienv->m_szs[v], v -- )
        cff /= bznv->binomial(ienv->m_szs[v]-1,addi % ienv->m_szs[v]);
};

// template< class real_t, char op > // multiplication multivariee dans la base de bernstein 
// inline static void _bvrcvloop_( bernstein::eenv * oenv, real_t * dst, 
//                              bernstein::eenv * aenv, real_t * sra,
//                              bernstein::eenv * benv, real_t * srb,
//                              bernstein::bzenv< real_t > * bznv = bernstein::bzenv<real_t>::_default_,
//                              unsigned * asupp = 0, unsigned nas = 0, 
//                              unsigned * oasup = 0, 
//                              unsigned * bsupp = 0, unsigned nbs = 0, 
//                              unsigned * obsup = 0,
//                              real_t   * ibzco = 0, 
//                              real_t   * ibzca = 0, real_t * ibzcb = 0 )
// {
//   /* calcul des addresses de sorties            */
//   unsigned * loasup, * lobsup; loasup = lobsup = 0;
//   /* calcul des addresses de sorties de a sur o */
//   if ( !oasup ) { if ( !asupp ) nas = aenv->data_size(); loasup = new unsigned[ nas ]; oasup = loasup; oaddress( oenv, oasup, aenv, asupp, nas ); };
//   /* calcul des addresses de sorties de b sur o */
//   if ( !obsup ) { if ( !bsupp ) nbs = benv->data_size(); lobsup = new unsigned[ nbs ]; obsup = lobsup; oaddress( oenv, obsup, benv, bsupp, nas ); };
  
//   /* calcul des coefficients binomiaux (inverses) */
//   real_t   * libzca, * libzcb, * libzco; libzca = libzcb = libzco = 0;
//   /* calcul des coefficients binomiaux pour aenv  */
//   if ( !ibzca ) { libzca = new real_t[ nas ]; ibzca  = libzca; rvbinoms( aenv, ibzca, asupp, nas, bznv ); };
//   /* calcul des coefficients binomiaux pour benv  */
//   if ( !ibzcb ) { libzcb = new real_t[ nbs ]; ibzcb  = libzcb; rvbinoms( benv, ibzcb, bsupp, nbs, bznv ); };
//   /* calcul des coefficients binomiaux pour oenv  */
//   if ( !ibzco ) { libzco = new real_t[ renv->data_size() ]; ibzco  = libzcb; rvbinoms( renv, ibzco, 0, 0, bznv ); };
//   unsigned ia, ib, oa; real_t   ca;
//   switch( (asupp!=0)<<1|(bsupp!=0) )
//     {
//     case 0:
//       assert(nas==aenv->data_size() && nbs == benv->data_size()); 
//       if ( op == '+' )
//      {
//        for ( ia = 0; ia < nas; ia ++ )
//          for ( oa = oasup[ia], ca = ibzca[ia], ib = 0; ib < nbs; ib ++ )
//            dst[oa+obsup[ib]] += (ibzco[oa+obsup[ib]]/(ca*ibzcb[ib]))*sra[ia]*srb[ib];
//        break;
//      }
//       if ( op == '-' )
//      {
//        for ( ia = 0; ia < nas; ia ++ )
//          for ( oa = oasup[ia], ca = ibzca[ia], ib = 0; ib < nbs; ib ++ )
//            dst[oa+obsup[ib]] -= (ibzco[oa+obsup[ib]]/(ca*ibzcb[ib]))*sra[ia]*srb[ib];
//        break;
//      };
//       break;
//     case 1:
//       assert(nas == aenv->data_size() && bsupp != 0);
//       if ( op == '+' )
//      {
//        assert(nas==aenv->data_size() && nbs == benv->data_size());
//        for ( ia = 0; ia < nas; ia ++ )
//          for ( oa = oasup[ia], ca = ibzca[ia], ib = 0; ib < nbs; ib ++ )
//            dst[oa+obsup[ib]] += (ibzco[oa+obsup[ib]]/(ca*ibzcb[ib]))*sra[ia]*srb[bsupp[ib]];
//        break;
//      }
//       if ( op == '-' )
//      {
//        for ( ia = 0; ia < nas; ia ++ )
//          for ( oa = oasup[ia], ca = ibzca[ia], ib = 0; ib < nbs; ib ++ )
//            dst[oa+obsup[ib]] -= (ibzco[oa+obsup[ib]]/(ca*ibzcb[ib]))*sra[ia]*srb[bsupp[ib]];
//        break;
//      };
//       break;
//     case 2:
//       assert(nbs == benv->data_size() && asupp != 0);
//       if ( op == '+' )
//      {
//        for ( ia = 0; ia < nas; ia ++ )
//          for ( oa = oasup[ia], ca = ibzca[ia], ib = 0; ib < nbs; ib ++ )
//            dst[oa+obsup[ib]] += (ibzco[oa+obsup[ib]]/(ca*ibzcb[ib]))*sra[asupp[ia]]*srb[ib];
//        break;
//      }
//       if ( op == '-' )
//      {
//        for ( ia = 0; ia < nas; ia ++ )
//          for ( oa = oasup[ia], ca = ibzca[ia], ib = 0; ib < nbs; ib ++ )
//            dst[oa+obsup[ib]] -= (ibzco[oa+obsup[ib]]/(ca*ibzcb[ib]))*sra[asupp[ia]]*srb[ib];
//        break;
//      };
//       break;
//     case 3:
//       assert(asupp != 0 && bsupp != 0);
//       if ( op == '+' )
//      {
//        for ( ia = 0; ia < nas; ia ++ )
//          for ( oa = oasup[ia], ca = ibzca[ia], ib = 0; ib < nbs; ib ++ )
//            dst[oa+obsup[ib]] += (ibzco[oa+obsup[ib]]/(ca*ibzcb[ib]))*sra[asupp[ia]]*srb[bsupp[ib]];
//        break;
//      }
//       if ( op == '-' )
//      {
//        assert(nas==aenv->data_size() && nbs == benv->data_size());
//        for ( ia = 0; ia < nas; ia ++ )
//          for ( oa = oasup[ia], ca = ibzca[ia], ib = 0; ib < nbs; ib ++ )
//            dst[oa+obsup[ib]] -= (ibzco[oa+obsup[ib]]/(ca*ibzcb[ib]))*sra[asupp[ia]]*srb[bsupp[ib]];
//        break;
//      };
//       break;
//     };
  
//   /* liberation des addresses de sorties   */
//   if ( loasup ) delete[] loasup; if ( lobsup ) delete[] lobsup;
//   /* liberation des coefficients binomiaux */
//   if ( libzca ) delete[] libzca; if ( libzcb ) delete[] libzcb; if ( libzco ) delete[] libzco;
// };

// template<typename real_t>
// static void bernstein_multiply_loop( eenv * o, eenv * a, eenv *  b, real_t * dst, real_t * src_a, real_t * src_b  )
// { _bvrcvloop_<real_t,'+'>( o, dst, a, src_a, b, src_b ); };

// template<typename real_t>
// static void bcrossp( eenv * oenv, real_t ** dst, eenv * aenv, real_t ** src_a, eenv * benv, real_t ** src_b,
//                   bernstein::bzenv< real_t > * bznv = bernstein::bzenv<real_t>::_default_,
//                   bool clear = true,
//                   unsigned * asupp = 0, unsigned nas = 0, 
//                   unsigned * oasup = 0, 
//                   unsigned * bsupp = 0, unsigned nbs = 0, 
//                   unsigned * obsup = 0 )
// {
//   if ( clear ) for (int d=0;d<3;std::fill(dst[d],dst[d]+oenv->data_size(),(real_t)0.0),d++);
//   unsigned * loasup, * lobsup; loasup = lobsup = 0;
//   /* calcul des addresses de sorties de a sur o */
//   if ( !oasup ) { if ( !asupp ) nas = aenv->data_size(); loasup = new unsigned[ nas ]; oasup = loasup; oaddress( oenv, oasup, aenv, asupp, nas ); };
//   /* calcul des addresses de sorties de b sur o */
//   if ( !obsup ) { if ( !bsupp ) nbs = benv->data_size(); lobsup = new unsigned[ nbs ]; obsup = lobsup; oaddress( oenv, obsup, benv, bsupp, nas ); };
//   /* calcul des coefficients binomiaux (inverses) */
//   real_t   * libzca, * libzcb, * libzco; libzca = libzcb = libzco = 0;
//   /* calcul des coefficients binomiaux pour aenv  */
//   if ( !ibzca ) { libzca = new real_t[ nas ]; ibzca  = libzca; rvbinoms( aenv, ibzca, asupp, nas, bznv ); };
//   /* calcul des coefficients binomiaux pour benv  */
//   if ( !ibzcb ) { libzcb = new real_t[ nbs ]; ibzcb  = libzcb; rvbinoms( benv, ibzcb, bsupp, nbs, bznv ); };
//   /* calcul des coefficients binomiaux pour oenv  */
//   if ( !ibzco ) { libzco = new real_t[ renv->data_size() ]; ibzco  = libzcb; rvbinoms( renv, ibzco, 0, 0, bznv ); };
  
//   _bvrcvloop_<real_t,'-'>(oenv,dst[2],aenv,src_a[1],benv,src_b[0],asupp,nas,obsup,bsupp,nbs,obsup);
//   _bvrcvloop_<real_t,'+'>(oenv,dst[2],aenv,src_a[0],benv,src_b[1],asupp,nas,obsup,bsupp,nbs,obsup);
//   _bvrcvloop_<real_t,'-'>(oenv,dst[0],aenv,src_a[2],benv,src_b[1],asupp,nas,obsup,bsupp,nbs,obsup);
//   _bvrcvloop_<real_t,'+'>(oenv,dst[0],aenv,src_a[1],benv,src_b[2],asupp,nas,obsup,bsupp,nbs,obsup);
//   _bvrcvloop_<real_t,'-'>(oenv,dst[1],aenv,src_a[0],benv,src_b[2],asupp,nas,obsup,bsupp,nbs,obsup);
//   _bvrcvloop_<real_t,'+'>(oenv,dst[1],aenv,src_a[2],benv,src_b[0],asupp,nas,obsup,bsupp,nbs,obsup);
  
//   /* liberation des addresses de sorties   */
//   if ( loasup ) delete[] loasup; if ( lobsup ) delete[] lobsup;
//   /* liberation des coefficients binomiaux */
//   if ( libzca ) delete[] libzca; if ( libzcb ) delete[] libzcb; if ( libzco ) delete[] libzco;
// };
/********************************************************************/

---