include/realroot/bernstein_eenv_multiply.hpp File Reference

Go to the source code of this file.

Functions

Function Documentation

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 
) [inline, static]

operations sous forme monomiale

Definition at line 90 of file bernstein_eenv_multiply.hpp.

References assert, and oaddress().

Referenced by mcrossp().

00097 {
00098   /* calcul des addresses de sorties            */
00099   unsigned * loasup, * lobsup; loasup = lobsup = 0;
00100   /* calcul des addresses de sorties de a sur o */
00101   if ( !oasup ) { if ( !asupp ) nas = aenv->data_size(); loasup = new unsigned[ nas ]; oasup = loasup; oaddress( oenv, oasup, aenv, asupp, nas ); };
00102   /* calcul des addresses de sorties de b sur o */
00103   if ( !obsup ) { if ( !bsupp ) nbs = benv->data_size(); lobsup = new unsigned[ nbs ]; obsup = lobsup; oaddress( oenv, obsup, benv, bsupp, nas ); };
00104 
00105   unsigned ia, ib, out; real_t   ca;
00106   switch( (asupp!=0)<<1|(bsupp!=0) )
00107     {
00108     case 0:
00109       assert(nas==aenv->data_size() && nbs == benv->data_size()); 
00110       if ( op == '+' )
00111         {
00112           for ( ia = 0; ia < nas; ia ++ )
00113             for ( ib = 0; ib < nbs;  ib ++ )
00114               dst[oasup[ia]+obsup[ib]] += sra[ia]*srb[ib];
00115           break;
00116         }
00117       if ( op == '-' )
00118         {
00119           for ( ia = 0; ia < nas; ia ++ )
00120             for ( ib = 0; ib < nbs; out += obsup[ib], ib ++ )
00121               dst[oasup[ia]+obsup[ib]] -= sra[ia]*srb[ib];
00122           break;
00123         };
00124       break;
00125     case 1:
00126       assert(nas == aenv->data_size() && bsupp != 0);
00127       if ( op == '+' )
00128         {
00129           assert(nas==aenv->data_size() && nbs == benv->data_size());
00130           for ( ia = 0; ia < nas; ia ++ )
00131             for (  ib = 0; ib < nbs; ib ++ )
00132               dst[oasup[ia]+obsup[ib]] += sra[ia]*srb[bsupp[ib]];
00133           break;
00134         }
00135       if ( op == '-' )
00136         {
00137           for ( ia = 0; ia < nas; ia ++ )
00138             for ( ib = 0; ib < nbs; ib ++ )
00139               dst[oasup[ia]+obsup[ib]] -= sra[ia]*srb[bsupp[ib]];
00140           break;
00141         };
00142       break;
00143     case 2:
00144       assert(nbs == benv->data_size() && asupp != 0);
00145       if ( op == '+' )
00146         {
00147           for ( ia = 0; ia < nas; ia ++ )
00148             for ( ib = 0; ib < nbs; ib ++ )
00149               dst[oasup[ia]+obsup[ib]] += sra[asupp[ia]]*srb[ib];
00150           break;
00151         }
00152       if ( op == '-' )
00153         {
00154           for ( ia = 0; ia < nas; ia ++ )
00155             for ( ib = 0; ib < nbs;  ib ++ )
00156               dst[oasup[ia]+obsup[ib]] -= sra[asupp[ia]]*srb[ib];
00157           break;
00158         };
00159       break;
00160     case 3:
00161       assert(asupp != 0 && bsupp != 0);
00162       if ( op == '+' )
00163         {
00164           for ( ia = 0; ia < nas; ia ++ )
00165             for ( ib = 0; ib < nbs;  ib ++ )
00166               dst[oasup[ia]+obsup[ib]] += sra[asupp[ia]]*srb[bsupp[ib]];
00167           break;
00168         }
00169       if ( op == '-' )
00170         {
00171           assert(nas==aenv->data_size() && nbs == benv->data_size());
00172           for ( ia = 0; ia < nas; ia ++ )
00173             for (  ib = 0; ib < nbs;  ib ++ )
00174               dst[oasup[ia]+obsup[ib]] -= sra[asupp[ia]]*srb[asupp[ib]];
00175           break;
00176         };
00177       break;
00178     };
00179   /* liberation des addresses de sorties   */
00180   if ( loasup ) delete[] loasup; if ( lobsup ) delete[] lobsup;
00181 };

real_t density ( real_t *  src,
const real_t  prec = 1e-12,
unsigned  nsmp = 0 
) [inline]

Definition at line 8 of file bernstein_eenv_multiply.hpp.

References mmx::abs().

00009 {
00010   unsigned c = 0;
00011   if ( nsmp == 0 )
00012     {
00013       nsmp = data_size();
00014       for ( unsigned i = 0; i < nsmp; i ++ )
00015         if ( std::abs(src[i]) > prec ) c++;
00016     }
00017   else
00018     {
00019       for ( unsigned i = 0; i < nsmp; i ++ )
00020         if ( std::abs(src[rand()%data_size()] > prec) ) c++;
00021     };
00022   return ((real_t)c)/nsmp;
00023 };

static void mcrossp ( eenv *  res,
eenv *  a,
eenv *  b,
real_t **  dst,
real_t **  src_a,
real_t **  src_b 
) [inline, static]

Definition at line 313 of file bernstein_eenv_multiply.hpp.

References mcrossp().

00314 { mcrossp( res, dst, a, src_a, b, src_b ); };

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 
) [inline, static]

Definition at line 265 of file bernstein_eenv_multiply.hpp.

References _mvrcvloop_(), mmx::tensor::clear(), mmx::vct::fill(), and oaddress().

Referenced by mcrossp(), and monoms_crossprod().

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

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 
) [inline, static]

Definition at line 328 of file bernstein_eenv_multiply.hpp.

References mmx::tensor::clear(), mmx::vct::fill(), and oaddress().

Referenced by monoms_dotprod().

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

static 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 
) [inline, static]

Definition at line 231 of file bernstein_eenv_multiply.hpp.

References mmx::tensor::clear(), and pmmul().

Referenced by spmmul().

00237 { pmmul( oenv, &dst, aenv, &src_a, benv, &src_b, 1, clear, asupp, nas, oasupp, bsupp, nbs, obsupp ); };

static void monoms_crossprod ( eenv *  res,
eenv *  a,
eenv *  b,
real_t **  dst,
real_t **  src_a,
real_t **  src_b 
) [inline, static]

Definition at line 317 of file bernstein_eenv_multiply.hpp.

References mcrossp().

00320 {
00321   res->op_mul(a,b);
00322   for ( int d = 0; d < 3; d ++ ) dst[d] = new real_t[ res->data_size() ];
00323   mcrossp( res, dst, a, src_a, b, src_b );
00324 };

static void monoms_dotprod ( eenv *  res,
eenv *  a,
eenv *  b,
real_t *&  dst,
real_t **  src_a,
real_t **  src_b,
int  n 
) [inline, static]

Definition at line 349 of file bernstein_eenv_multiply.hpp.

References mdotp().

00350 {
00351   res->op_mul(a,b);
00352   dst = new real_t[res->data_size()];
00353   mdotp( res, dst, a, src_a, b, src_b, n );
00354 };

static void monoms_multiply ( bernstein::eenv *  res,
bernstein::eenv *  a,
bernstein::eenv *  b,
real_t *  dst,
real_t *  src_a,
real_t *  src_b 
) [inline, static]

Definition at line 206 of file bernstein_eenv_multiply.hpp.

References monoms_multiply_loop().

00208 { res->op_mul(a,b); monoms_multiply_loop( res, a, b, dst, src_a, src_b  ); };

static void monoms_multiply_loop ( eenv *  o,
eenv *  a,
eenv *  b,
real_t *  dst,
real_t *  src_a,
real_t *  src_b 
) [inline, static]

Definition at line 202 of file bernstein_eenv_multiply.hpp.

Referenced by monoms_multiply().

00203 { _mvrcvloop_<real_t,'+'>( o, dst, a, src_a, b, src_b ); };

static void oaddress ( bernstein::eenv *  oenv,
unsigned *  osupp,
bernstein::eenv *  ienv,
unsigned *  isupp = 0,
unsigned  nsp = 0 
) [static]

Definition at line 66 of file bernstein_eenv_multiply.hpp.

References vmap().

Referenced by _mvrcvloop_(), mcrossp(), mdotp(), msimplify(), pmmul(), and vrestrict().

00068 {
00069   int v;
00070   unsigned c;
00071   sz_t vmap_[ ienv->m_nvr ];
00072   unsigned addi,addo;
00073   vmap( vmap_, oenv, ienv );
00074   if ( isupp )
00075     for ( c = 0; c < nsp; osupp[c] = addo, c ++ )
00076       for ( addi = isupp[c], addo = 0, v = ienv->m_nvr-1; addi; addi /= ienv->m_szs[v], v-- )
00077         addo += (addi%ienv->m_szs[v])*oenv->m_str[vmap_[v]];
00078   else
00079     {
00080       nsp = ienv->data_size();
00081     for ( c = 0; c < nsp; osupp[c] = addo, c ++ )
00082       for ( addi = c, addo = 0, v = ienv->m_nvr-1; addi; addi /= ienv->m_szs[v], v -- )
00083         addo += (addi%ienv->m_szs[v])*oenv->m_str[vmap_[v]];
00084     };
00085 };

void op_mul ( eenv *  a,
eenv *  b 
)

Definition at line 37 of file bernstein_eenv_multiply.hpp.

References mmx::sparse::copy(), mmx::vctops::max(), and mmx::vctops::set_conversion().

00038 {
00039   if ( nvars() ) this->~eenv();
00040   _alloc_( a->nvars()+b->nvars() );
00041   std::copy( a->m_vrs, a->m_vrs+a->nvars(), m_vrs );
00042   std::copy( b->m_vrs, b->m_vrs+b->nvars(), m_vrs + a->nvars() );
00043   m_nvr = vctops::set_conversion( m_vrs, m_nvr );
00044   for ( sz_t v = 0; v < m_nvr; v ++ )
00045     {
00046       int av = a->indexofvar(var(v));
00047       int bv = b->indexofvar(var(v));
00048       m_szs[ v ] = ((av!=-1)?(a->sz(av)-1):0) + ((bv!=-1)?(b->sz(bv)-1):0) + 1;
00049     };
00050   m_mxvr = std::max(a->m_mxvr,b->m_mxvr);
00051   compute_strides();
00052 };

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 
) [inline, static]

Definition at line 212 of file bernstein_eenv_multiply.hpp.

References mmx::tensor::clear(), mmx::vct::fill(), and oaddress().

Referenced by mmul().

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

static void rvbinoms ( bernstein::eenv *  ienv,
real_t *  bcff,
unsigned *  isupp = 0,
unsigned  nsp = 0,
bernstein::bzenv< real_t > *  bznv = bzenv<real_t>::_default_ 
) [inline, static]

operations sous forme de bernstein

Definition at line 360 of file bernstein_eenv_multiply.hpp.

00362 {
00363   real_t cff;
00364   unsigned c,addi;
00365   int v;
00366   if ( isupp )
00367     for ( c = 0; c < nsp; bcff[c] = cff, c ++ )
00368       for ( cff = (real_t)1.0, addi = isupp[c],  v = ienv->m_nvr-1; addi; addi /= ienv->m_szs[v], v-- )
00369         cff /= bznv->binomial(ienv->m_szs[v]-1,addi % ienv->m_szs[v]);
00370   else
00371     for ( c = 0; c < ienv->data_size(); bcff[c] = cff,  c ++ )
00372       for ( addi = c, v = ienv->m_nvr-1; addi; addi /= ienv->m_szs[v], v -- )
00373         cff /= bznv->binomial(ienv->m_szs[v]-1,addi % ienv->m_szs[v]);
00374 };

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 
) [inline, static]

Definition at line 258 of file bernstein_eenv_multiply.hpp.

References mmx::tensor::clear(), and spmmul().

00260 {
00261   spmmul( oenv, &dst, aenv, &src_a, benv, &src_b, 1, clear, stats );
00262 };

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 
) [inline, static]

Definition at line 240 of file bernstein_eenv_multiply.hpp.

References assert, and mmul().

Referenced by spmmul().

00242 {
00243   unsigned nas,nbs;
00244   unsigned * assup = new unsigned[ aenv->data_size() ];
00245   unsigned * bssup = new unsigned[ benv->data_size() ];
00246   assert( assup && bssup );
00247   for ( int i = 0; i < n; i ++ )
00248     {
00249       nas = aenv->support(assup,src_a[i]);
00250       nbs = benv->support(bssup,src_b[i]);
00251       mmul( oenv, dst[i], aenv, src_a[i], benv, src_b[i], true, assup, nas, 0, bssup, nbs, 0 );
00252     };
00253   delete[] assup;
00254   delete[] bssup;
00255 };

unsigned support ( unsigned *  supp,
real_t *  src,
const real_t &  prec = 1e-6 
) [inline]

Definition at line 27 of file bernstein_eenv_multiply.hpp.

References mmx::abs().

00028 {
00029 //  supp = new unsigned[ data_size() ];
00030   unsigned      c = 0;
00031   for ( unsigned i = 0; i < data_size(); i ++ )
00032     if ( std::abs(src[i]) > prec ) supp[c++] = i;
00033   return c;
00034 };

static void vmap ( sz_t *  vmap,
bernstein::eenv *  o,
bernstein::eenv *  i 
) [static]

Definition at line 55 of file bernstein_eenv_multiply.hpp.

References assert.

Referenced by oaddress().

00056 {
00057   for ( sz_t v = 0; v < i->m_nvr; v ++ )
00058     {
00059       vmap[v] = -1;
00060       for ( sz_t k = 0; k < o->m_nvr; k ++ )
00061         if ( o->m_vrs[k] == i->m_vrs[v] ) vmap[v] = k;
00062       assert(vmap[v]!=-1);
00063     };
00064 };

Generated on 6 Dec 2012 for realroot by  doxygen 1.6.1