Developer documentation

bernstein_eenv_multiply.hpp
Go to the documentation of this file.
1 /********************************************************************
2  * This file is part of the source code of the realroot library.
3  * Author(s): J.P. Pavone, GALAAD, INRIA
4  * $Id: eenv-multiply.hpp,v 1.1 2005/07/11 10:03:55 jppavone Exp $
5  ********************************************************************/
6 /* densite exacte / statistique des coefficients non-nuls */
7 template< typename real_t >
8 real_t density( real_t * src, const real_t prec = 1e-12, unsigned nsmp = 0 )
9 {
10  unsigned c = 0;
11  if ( nsmp == 0 )
12  {
13  nsmp = data_size();
14  for ( unsigned i = 0; i < nsmp; i ++ )
15  if ( std::abs(src[i]) > prec ) c++;
16  }
17  else
18  {
19  for ( unsigned i = 0; i < nsmp; i ++ )
20  if ( std::abs(src[rand()%data_size()] > prec) ) c++;
21  };
22  return ((real_t)c)/nsmp;
23 };
24 
25 /* calcul le support des coefficients non-nuls dans src */
26 template<class real_t>
27 unsigned support( unsigned * supp, real_t * src, const real_t& prec = 1e-6 )
28 {
29 // supp = new unsigned[ data_size() ];
30  unsigned c = 0;
31  for ( unsigned i = 0; i < data_size(); i ++ )
32  if ( std::abs(src[i]) > prec ) supp[c++] = i;
33  return c;
34 };
35 
36 /* transforme l'environement en l'environement de a*b */
37 void op_mul( eenv * a, eenv * b )
38 {
39  if ( nvars() ) this->~eenv();
40  _alloc_( a->nvars()+b->nvars() );
41  std::copy( a->m_vrs, a->m_vrs+a->nvars(), m_vrs );
42  std::copy( b->m_vrs, b->m_vrs+b->nvars(), m_vrs + a->nvars() );
43  m_nvr = vctops::set_conversion( m_vrs, m_nvr );
44  for ( sz_t v = 0; v < m_nvr; v ++ )
45  {
46  int av = a->indexofvar(var(v));
47  int bv = b->indexofvar(var(v));
48  m_szs[ v ] = ((av!=-1)?(a->sz(av)-1):0) + ((bv!=-1)?(b->sz(bv)-1):0) + 1;
49  };
50  m_mxvr = std::max(a->m_mxvr,b->m_mxvr);
51  compute_strides();
52 };
53 
54 /* mapping des noms i->o */
55 static void vmap( sz_t * vmap, bernstein::eenv * o, bernstein::eenv * i )
56 {
57  for ( sz_t v = 0; v < i->m_nvr; v ++ )
58  {
59  vmap[v] = -1;
60  for ( sz_t k = 0; k < o->m_nvr; k ++ )
61  if ( o->m_vrs[k] == i->m_vrs[v] ) vmap[v] = k;
62  assert(vmap[v]!=-1);
63  };
64 };
65 
66 static void oaddress( bernstein::eenv * oenv, unsigned * osupp,
67  bernstein::eenv * ienv, unsigned * isupp = 0, unsigned nsp = 0 )
68 {
69  int v;
70  unsigned c;
71  sz_t vmap_[ ienv->m_nvr ];
72  unsigned addi,addo;
73  vmap( vmap_, oenv, ienv );
74  if ( isupp )
75  for ( c = 0; c < nsp; osupp[c] = addo, c ++ )
76  for ( addi = isupp[c], addo = 0, v = ienv->m_nvr-1; addi; addi /= ienv->m_szs[v], v-- )
77  addo += (addi%ienv->m_szs[v])*oenv->m_str[vmap_[v]];
78  else
79  {
80  nsp = ienv->data_size();
81  for ( c = 0; c < nsp; osupp[c] = addo, c ++ )
82  for ( addi = c, addo = 0, v = ienv->m_nvr-1; addi; addi /= ienv->m_szs[v], v -- )
83  addo += (addi%ienv->m_szs[v])*oenv->m_str[vmap_[v]];
84  };
85 };
86 
89 template< class real_t, char op > // multiplication multivariee dans la base monomiale
90 inline static void _mvrcvloop_( bernstein::eenv * oenv, real_t * dst,
91  bernstein::eenv * aenv, real_t * sra,
92  bernstein::eenv * benv, real_t * srb,
93  unsigned * asupp = 0, unsigned nas = 0,
94  unsigned * oasup = 0,
95  unsigned * bsupp = 0, unsigned nbs = 0,
96  unsigned * obsup = 0 )
97 {
98  /* calcul des addresses de sorties */
99  unsigned * loasup, * lobsup; loasup = lobsup = 0;
100  /* calcul des addresses de sorties de a sur o */
101  if ( !oasup ) { if ( !asupp ) nas = aenv->data_size(); loasup = new unsigned[ nas ]; oasup = loasup; oaddress( oenv, oasup, aenv, asupp, nas ); };
102  /* calcul des addresses de sorties de b sur o */
103  if ( !obsup ) { if ( !bsupp ) nbs = benv->data_size(); lobsup = new unsigned[ nbs ]; obsup = lobsup; oaddress( oenv, obsup, benv, bsupp, nas ); };
104 
105  unsigned ia, ib, out; real_t ca;
106  switch( (asupp!=0)<<1|(bsupp!=0) )
107  {
108  case 0:
109  assert(nas==aenv->data_size() && nbs == benv->data_size());
110  if ( op == '+' )
111  {
112  for ( ia = 0; ia < nas; ia ++ )
113  for ( ib = 0; ib < nbs; ib ++ )
114  dst[oasup[ia]+obsup[ib]] += sra[ia]*srb[ib];
115  break;
116  }
117  if ( op == '-' )
118  {
119  for ( ia = 0; ia < nas; ia ++ )
120  for ( ib = 0; ib < nbs; out += obsup[ib], ib ++ )
121  dst[oasup[ia]+obsup[ib]] -= sra[ia]*srb[ib];
122  break;
123  };
124  break;
125  case 1:
126  assert(nas == aenv->data_size() && bsupp != 0);
127  if ( op == '+' )
128  {
129  assert(nas==aenv->data_size() && nbs == benv->data_size());
130  for ( ia = 0; ia < nas; ia ++ )
131  for ( ib = 0; ib < nbs; ib ++ )
132  dst[oasup[ia]+obsup[ib]] += sra[ia]*srb[bsupp[ib]];
133  break;
134  }
135  if ( op == '-' )
136  {
137  for ( ia = 0; ia < nas; ia ++ )
138  for ( ib = 0; ib < nbs; ib ++ )
139  dst[oasup[ia]+obsup[ib]] -= sra[ia]*srb[bsupp[ib]];
140  break;
141  };
142  break;
143  case 2:
144  assert(nbs == benv->data_size() && asupp != 0);
145  if ( op == '+' )
146  {
147  for ( ia = 0; ia < nas; ia ++ )
148  for ( ib = 0; ib < nbs; ib ++ )
149  dst[oasup[ia]+obsup[ib]] += sra[asupp[ia]]*srb[ib];
150  break;
151  }
152  if ( op == '-' )
153  {
154  for ( ia = 0; ia < nas; ia ++ )
155  for ( ib = 0; ib < nbs; ib ++ )
156  dst[oasup[ia]+obsup[ib]] -= sra[asupp[ia]]*srb[ib];
157  break;
158  };
159  break;
160  case 3:
161  assert(asupp != 0 && bsupp != 0);
162  if ( op == '+' )
163  {
164  for ( ia = 0; ia < nas; ia ++ )
165  for ( ib = 0; ib < nbs; ib ++ )
166  dst[oasup[ia]+obsup[ib]] += sra[asupp[ia]]*srb[bsupp[ib]];
167  break;
168  }
169  if ( op == '-' )
170  {
171  assert(nas==aenv->data_size() && nbs == benv->data_size());
172  for ( ia = 0; ia < nas; ia ++ )
173  for ( ib = 0; ib < nbs; ib ++ )
174  dst[oasup[ia]+obsup[ib]] -= sra[asupp[ia]]*srb[asupp[ib]];
175  break;
176  };
177  break;
178  };
179  /* liberation des addresses de sorties */
180  if ( loasup ) delete[] loasup; if ( lobsup ) delete[] lobsup;
181 };
182 
183 
184 // template< typename real_t >
185 // static void sparse_multiply( bernstein::eenv * renv, real_t * rdata,
186 // bernstein::eenv * aenv, unsigned nasupp, unsigned * asupp, real_t * adata,
187 // bernstein::eenv * benv, unsigned nbsupp, unsigned * bsupp, real_t * bdata )
188 // {
189 
190 // unsigned * osuppa = new unsigned[ nasupp ];
191 // unsigned * osuppb = new unsigned[ nbsupp ];
192 // this-osupp( renv, osuppa, aenv, asupp, nasupp );
193 // this->osupp( renv, osuppb, benv, bsupp, nbsupp );
194 // for ( unsigned i = 0; i < nasupp; i ++ )
195 // for ( unsigned j = 0; j < nbsupp; j ++ )
196 // rdata[osuppa[i]+osuppb[j]] += adata[asupp[i]] * bdata[bsupp[j]];
197 // delete[] osuppa;
198 // delete[] osuppb;
199 // }
200 
201 template<typename real_t>
202 static void monoms_multiply_loop( eenv * o, eenv * a, eenv * b, real_t * dst, real_t * src_a, real_t * src_b )
203 { _mvrcvloop_<real_t,'+'>( o, dst, a, src_a, b, src_b ); };
204 
205 template<typename real_t>
206 static void monoms_multiply( bernstein::eenv * res, bernstein::eenv * a, bernstein::eenv *b,
207  real_t * dst, real_t * src_a, real_t * src_b )
208 { res->op_mul(a,b); monoms_multiply_loop( res, a, b, dst, src_a, src_b ); };
209 
210 
211 template< typename real_t > // multiplication vectorielle
212 static void pmmul( eenv * oenv, real_t ** dst, eenv * aenv, real_t ** src_a, eenv * benv, real_t ** src_b, int n,
213  bool clear = true,
214  unsigned * asupp = 0, unsigned nas = 0,
215  unsigned * oasup = 0,
216  unsigned * bsupp = 0, unsigned nbs = 0,
217  unsigned * obsup = 0 )
218 {
219  if ( clear ) for (int d=0;d<n;std::fill(dst[d],dst[d]+oenv->data_size(),(real_t)0.0),d++);
220  unsigned * loasup, * lobsup; loasup = lobsup = 0;
221  /* calcul des addresses de sorties de a sur o */
222  if ( !oasup ) { if ( !asupp ) nas = aenv->data_size(); loasup = new unsigned[ nas ]; oasup = loasup; oaddress( oenv, oasup, aenv, asupp, nas ); };
223  /* calcul des addresses de sorties de b sur o */
224  if ( !obsup ) { if ( !bsupp ) nbs = benv->data_size(); lobsup = new unsigned[ nbs ]; obsup = lobsup; oaddress( oenv, obsup, benv, bsupp, nbs ); };
225  for ( unsigned i = 0; i < n; i ++ )
226  _mvrcvloop_<real_t,'+'>(oenv,dst[i],aenv,src_a[i],benv,src_b[i],asupp,nas,oasup,bsupp,nbs,obsup);
227  if ( loasup ) delete[] loasup; if ( lobsup ) delete[] lobsup;
228 };
229 
230 template< typename real_t > // multiplication simple
231 static inline void mmul( eenv * oenv, real_t * dst, eenv * aenv, real_t * src_a, eenv * benv, real_t * src_b,
232  bool clear = true,
233  unsigned * asupp = 0, unsigned nas = 0,
234  unsigned * oasupp = 0,
235  unsigned * bsupp = 0, unsigned nbs = 0,
236  unsigned * obsupp = 0 )
237 { pmmul( oenv, &dst, aenv, &src_a, benv, &src_b, 1, clear, asupp, nas, oasupp, bsupp, nbs, obsupp ); };
238 
239 template< typename real_t > // multiplication vectorielle creuse (calcul du support + multiplication)
240 static void spmmul( eenv * oenv, real_t ** dst, eenv * aenv, real_t ** src_a, eenv * benv, real_t ** src_b, int n,
241  bool clear = true, bool stats = false )
242 {
243  unsigned nas,nbs;
244  unsigned * assup = new unsigned[ aenv->data_size() ];
245  unsigned * bssup = new unsigned[ benv->data_size() ];
246  assert( assup && bssup );
247  for ( int i = 0; i < n; i ++ )
248  {
249  nas = aenv->support(assup,src_a[i]);
250  nbs = benv->support(bssup,src_b[i]);
251  mmul( oenv, dst[i], aenv, src_a[i], benv, src_b[i], true, assup, nas, 0, bssup, nbs, 0 );
252  };
253  delete[] assup;
254  delete[] bssup;
255 };
256 
257 template< typename real_t > // multiplication simple creuse (calcul du support + multiplication)
258 static void spmmul( eenv * oenv, real_t * dst, eenv * aenv, real_t * src_a, eenv * benv, real_t * src_b,
259  bool clear = true, bool stats = false )
260 {
261  spmmul( oenv, &dst, aenv, &src_a, benv, &src_b, 1, clear, stats );
262 };
263 
264 template<typename real_t>
265 static void mcrossp( eenv * oenv, real_t ** dst, eenv * aenv, real_t ** src_a, eenv * benv, real_t ** src_b,
266  bool clear = true,
267  unsigned * asupp = 0, unsigned nas = 0,
268  unsigned * oasup = 0,
269  unsigned * bsupp = 0, unsigned nbs = 0,
270  unsigned * obsup = 0 )
271 
272 
273 {
274  if ( clear ) for (int d=0;d<3;std::fill(dst[d],dst[d]+oenv->data_size(),(real_t)0.0),d++);
275  unsigned * loasup, * lobsup; loasup = lobsup = 0;
276  /* calcul des addresses de sorties de a sur o */
277  if ( !oasup ) { if ( !asupp ) nas = aenv->data_size(); loasup = new unsigned[ nas ]; oasup = loasup; oaddress( oenv, oasup, aenv, asupp, nas ); };
278  /* calcul des addresses de sorties de b sur o */
279  if ( !obsup ) { if ( !bsupp ) nbs = benv->data_size(); lobsup = new unsigned[ nbs ]; obsup = lobsup; oaddress( oenv, obsup, benv, bsupp, nbs ); };
280  _mvrcvloop_<real_t,'+'>(oenv,dst[2],aenv,src_a[0],benv,src_b[1],asupp,nas,oasup,bsupp,nbs,obsup);
281  _mvrcvloop_<real_t,'-'>(oenv,dst[2],aenv,src_a[1],benv,src_b[0],asupp,nas,oasup,bsupp,nbs,obsup);
282  _mvrcvloop_<real_t,'+'>(oenv,dst[0],aenv,src_a[1],benv,src_b[2],asupp,nas,oasup,bsupp,nbs,obsup);
283  _mvrcvloop_<real_t,'-'>(oenv,dst[0],aenv,src_a[2],benv,src_b[1],asupp,nas,oasup,bsupp,nbs,obsup);
284  _mvrcvloop_<real_t,'+'>(oenv,dst[1],aenv,src_a[2],benv,src_b[0],asupp,nas,oasup,bsupp,nbs,obsup);
285  _mvrcvloop_<real_t,'-'>(oenv,dst[1],aenv,src_a[0],benv,src_b[2],asupp,nas,oasup,bsupp,nbs,obsup);
286  if ( loasup ) delete[] loasup; if ( lobsup ) delete[] lobsup;
287 };
288 
289 
290 // template<typename real_t>
291 // static void spmcrossp( eenv * oenv, real_t ** dst, eenv * aenv, real_t ** src_a, eenv * benv, real_t ** src_b,
292 // bool clear = true, bool stats = false )
293 // {
294 // if ( clear ) for (int d=0;d<3;std::fill(dst[d],dst[d]+oenv->data_size(),(real_t)0.0),d++);
295 // unsigned * asupp = new unsigned[ a->data_size() ];
296 // unsigned * bsupp = new unsigned [b->data_size()];
297 // unsigned * loasup, * lobsup; loasup = lobsup = 0;
298 // /* calcul des addresses de sorties de a sur o */
299 // if ( !oasup ) { if ( !asupp ) nas = aenv->data_size(); loasup = new unsigned[ nas ]; oasup = loasup; oaddress( oenv, oasup, aenv, asupp, nas ); };
300 // /* calcul des addresses de sorties de b sur o */
301 // if ( !obsup ) { if ( !bsupp ) nbs = benv->data_size(); lobsup = new unsigned[ nbs ]; obsup = lobsup; oaddress( oenv, obsup, benv, bsupp, nbs ); };
302 // _mvrcvloop_<real_t,'+'>(oenv,dst[2],aenv,src_a[0],benv,src_b[1],asupp,nas,oasup,bsupp,nbs,obsup);
303 // _mvrcvloop_<real_t,'-'>(oenv,dst[1],aenv,src_a[0],benv,src_b[2],asupp,nas,oasup,bsupp,nbs,obsup);
304 // _mvrcvloop_<real_t,'-'>(oenv,dst[2],aenv,src_a[1],benv,src_b[0],asupp,nas,oasup,bsupp,nbs,obsup);
305 // _mvrcvloop_<real_t,'+'>(oenv,dst[0],aenv,src_a[1],benv,src_b[2],asupp,nas,oasup,bsupp,nbs,obsup);
306 // _mvrcvloop_<real_t,'-'>(oenv,dst[0],aenv,src_a[2],benv,src_b[1],asupp,nas,oasup,bsupp,nbs,obsup);
307 // _mvrcvloop_<real_t,'+'>(oenv,dst[1],aenv,src_a[2],benv,src_b[0],asupp,nas,oasup,bsupp,nbs,obsup);
308 // if ( loasup ) delete[] loasup; if ( lobsup ) delete[] lobsup;
309 // };
310 
311 
312 template<typename real_t> // compatibilite ancien code
313 static void mcrossp( eenv * res, eenv * a, eenv * b, real_t ** dst, real_t ** src_a, real_t ** src_b )
314 { mcrossp( res, dst, a, src_a, b, src_b ); };
315 
316 template<typename real_t>// compatibilite ancien code
317 static void monoms_crossprod( eenv * res, eenv * a, eenv * b,
318  real_t ** dst,
319  real_t ** src_a, real_t ** src_b )
320 {
321  res->op_mul(a,b);
322  for ( int d = 0; d < 3; d ++ ) dst[d] = new real_t[ res->data_size() ];
323  mcrossp( res, dst, a, src_a, b, src_b );
324 };
325 
326 
327 template<typename real_t>
328 static void mdotp( eenv * oenv, real_t * dst, eenv * aenv, real_t ** src_a, eenv * benv, real_t ** src_b, int n,
329  bool clear = true,
330  unsigned * asupp = 0, unsigned nas = 0,
331  unsigned * oasup = 0,
332  unsigned * bsupp = 0, unsigned nbs = 0,
333  unsigned * obsup = 0 )
334 
335 
336 {
337  if ( clear ) std::fill(dst,dst+oenv->data_size(),(real_t)0.0);
338  unsigned * loasup, * lobsup; loasup = lobsup = 0;
339  /* calcul des addresses de sorties de a sur o */
340  if ( !oasup ) { if ( !asupp ) nas = aenv->data_size(); loasup = new unsigned[ nas ]; oasup = loasup; oaddress( oenv, oasup, aenv, asupp, nas ); };
341  /* calcul des addresses de sorties de b sur o */
342  if ( !obsup ) { if ( !bsupp ) nbs = benv->data_size(); lobsup = new unsigned[ nbs ]; obsup = lobsup; oaddress( oenv, obsup, benv, bsupp, nas ); };
343  for ( int i = 0; i < n; i ++ )
344  _mvrcvloop_<real_t,'+'>(oenv,dst,aenv,src_a[i],benv,src_b[i],asupp,nas,oasup,bsupp,nbs,obsup);
345  if ( loasup ) delete[] loasup; if ( lobsup ) delete[] lobsup;
346 };
347 
348 template<typename real_t> // compatibilite ancien code
349 static void monoms_dotprod( eenv * res, eenv * a, eenv * b, real_t *& dst, real_t ** src_a, real_t ** src_b, int n )
350 {
351  res->op_mul(a,b);
352  dst = new real_t[res->data_size()];
353  mdotp( res, dst, a, src_a, b, src_b, n );
354 };
355 
356 
359 template< class real_t >
360 static void rvbinoms( bernstein::eenv * ienv, real_t * bcff,
361  unsigned * isupp = 0, unsigned nsp = 0, bernstein::bzenv<real_t> * bznv = bzenv<real_t>::_default_ )
362 {
363  real_t cff;
364  unsigned c,addi;
365  int v;
366  if ( isupp )
367  for ( c = 0; c < nsp; bcff[c] = cff, c ++ )
368  for ( cff = (real_t)1.0, addi = isupp[c], v = ienv->m_nvr-1; addi; addi /= ienv->m_szs[v], v-- )
369  cff /= bznv->binomial(ienv->m_szs[v]-1,addi % ienv->m_szs[v]);
370  else
371  for ( c = 0; c < ienv->data_size(); bcff[c] = cff, c ++ )
372  for ( addi = c, v = ienv->m_nvr-1; addi; addi /= ienv->m_szs[v], v -- )
373  cff /= bznv->binomial(ienv->m_szs[v]-1,addi % ienv->m_szs[v]);
374 };
375 
376 // template< class real_t, char op > // multiplication multivariee dans la base de bernstein
377 // inline static void _bvrcvloop_( bernstein::eenv * oenv, real_t * dst,
378 // bernstein::eenv * aenv, real_t * sra,
379 // bernstein::eenv * benv, real_t * srb,
380 // bernstein::bzenv< real_t > * bznv = bernstein::bzenv<real_t>::_default_,
381 // unsigned * asupp = 0, unsigned nas = 0,
382 // unsigned * oasup = 0,
383 // unsigned * bsupp = 0, unsigned nbs = 0,
384 // unsigned * obsup = 0,
385 // real_t * ibzco = 0,
386 // real_t * ibzca = 0, real_t * ibzcb = 0 )
387 // {
388 // /* calcul des addresses de sorties */
389 // unsigned * loasup, * lobsup; loasup = lobsup = 0;
390 // /* calcul des addresses de sorties de a sur o */
391 // if ( !oasup ) { if ( !asupp ) nas = aenv->data_size(); loasup = new unsigned[ nas ]; oasup = loasup; oaddress( oenv, oasup, aenv, asupp, nas ); };
392 // /* calcul des addresses de sorties de b sur o */
393 // if ( !obsup ) { if ( !bsupp ) nbs = benv->data_size(); lobsup = new unsigned[ nbs ]; obsup = lobsup; oaddress( oenv, obsup, benv, bsupp, nas ); };
394 
395 // /* calcul des coefficients binomiaux (inverses) */
396 // real_t * libzca, * libzcb, * libzco; libzca = libzcb = libzco = 0;
397 // /* calcul des coefficients binomiaux pour aenv */
398 // if ( !ibzca ) { libzca = new real_t[ nas ]; ibzca = libzca; rvbinoms( aenv, ibzca, asupp, nas, bznv ); };
399 // /* calcul des coefficients binomiaux pour benv */
400 // if ( !ibzcb ) { libzcb = new real_t[ nbs ]; ibzcb = libzcb; rvbinoms( benv, ibzcb, bsupp, nbs, bznv ); };
401 // /* calcul des coefficients binomiaux pour oenv */
402 // if ( !ibzco ) { libzco = new real_t[ renv->data_size() ]; ibzco = libzcb; rvbinoms( renv, ibzco, 0, 0, bznv ); };
403 // unsigned ia, ib, oa; real_t ca;
404 // switch( (asupp!=0)<<1|(bsupp!=0) )
405 // {
406 // case 0:
407 // assert(nas==aenv->data_size() && nbs == benv->data_size());
408 // if ( op == '+' )
409 // {
410 // for ( ia = 0; ia < nas; ia ++ )
411 // for ( oa = oasup[ia], ca = ibzca[ia], ib = 0; ib < nbs; ib ++ )
412 // dst[oa+obsup[ib]] += (ibzco[oa+obsup[ib]]/(ca*ibzcb[ib]))*sra[ia]*srb[ib];
413 // break;
414 // }
415 // if ( op == '-' )
416 // {
417 // for ( ia = 0; ia < nas; ia ++ )
418 // for ( oa = oasup[ia], ca = ibzca[ia], ib = 0; ib < nbs; ib ++ )
419 // dst[oa+obsup[ib]] -= (ibzco[oa+obsup[ib]]/(ca*ibzcb[ib]))*sra[ia]*srb[ib];
420 // break;
421 // };
422 // break;
423 // case 1:
424 // assert(nas == aenv->data_size() && bsupp != 0);
425 // if ( op == '+' )
426 // {
427 // assert(nas==aenv->data_size() && nbs == benv->data_size());
428 // for ( ia = 0; ia < nas; ia ++ )
429 // for ( oa = oasup[ia], ca = ibzca[ia], ib = 0; ib < nbs; ib ++ )
430 // dst[oa+obsup[ib]] += (ibzco[oa+obsup[ib]]/(ca*ibzcb[ib]))*sra[ia]*srb[bsupp[ib]];
431 // break;
432 // }
433 // if ( op == '-' )
434 // {
435 // for ( ia = 0; ia < nas; ia ++ )
436 // for ( oa = oasup[ia], ca = ibzca[ia], ib = 0; ib < nbs; ib ++ )
437 // dst[oa+obsup[ib]] -= (ibzco[oa+obsup[ib]]/(ca*ibzcb[ib]))*sra[ia]*srb[bsupp[ib]];
438 // break;
439 // };
440 // break;
441 // case 2:
442 // assert(nbs == benv->data_size() && asupp != 0);
443 // if ( op == '+' )
444 // {
445 // for ( ia = 0; ia < nas; ia ++ )
446 // for ( oa = oasup[ia], ca = ibzca[ia], ib = 0; ib < nbs; ib ++ )
447 // dst[oa+obsup[ib]] += (ibzco[oa+obsup[ib]]/(ca*ibzcb[ib]))*sra[asupp[ia]]*srb[ib];
448 // break;
449 // }
450 // if ( op == '-' )
451 // {
452 // for ( ia = 0; ia < nas; ia ++ )
453 // for ( oa = oasup[ia], ca = ibzca[ia], ib = 0; ib < nbs; ib ++ )
454 // dst[oa+obsup[ib]] -= (ibzco[oa+obsup[ib]]/(ca*ibzcb[ib]))*sra[asupp[ia]]*srb[ib];
455 // break;
456 // };
457 // break;
458 // case 3:
459 // assert(asupp != 0 && bsupp != 0);
460 // if ( op == '+' )
461 // {
462 // for ( ia = 0; ia < nas; ia ++ )
463 // for ( oa = oasup[ia], ca = ibzca[ia], ib = 0; ib < nbs; ib ++ )
464 // dst[oa+obsup[ib]] += (ibzco[oa+obsup[ib]]/(ca*ibzcb[ib]))*sra[asupp[ia]]*srb[bsupp[ib]];
465 // break;
466 // }
467 // if ( op == '-' )
468 // {
469 // assert(nas==aenv->data_size() && nbs == benv->data_size());
470 // for ( ia = 0; ia < nas; ia ++ )
471 // for ( oa = oasup[ia], ca = ibzca[ia], ib = 0; ib < nbs; ib ++ )
472 // dst[oa+obsup[ib]] -= (ibzco[oa+obsup[ib]]/(ca*ibzcb[ib]))*sra[asupp[ia]]*srb[bsupp[ib]];
473 // break;
474 // };
475 // break;
476 // };
477 
478 // /* liberation des addresses de sorties */
479 // if ( loasup ) delete[] loasup; if ( lobsup ) delete[] lobsup;
480 // /* liberation des coefficients binomiaux */
481 // if ( libzca ) delete[] libzca; if ( libzcb ) delete[] libzcb; if ( libzco ) delete[] libzco;
482 // };
483 
484 // template<typename real_t>
485 // static void bernstein_multiply_loop( eenv * o, eenv * a, eenv * b, real_t * dst, real_t * src_a, real_t * src_b )
486 // { _bvrcvloop_<real_t,'+'>( o, dst, a, src_a, b, src_b ); };
487 
488 // template<typename real_t>
489 // static void bcrossp( eenv * oenv, real_t ** dst, eenv * aenv, real_t ** src_a, eenv * benv, real_t ** src_b,
490 // bernstein::bzenv< real_t > * bznv = bernstein::bzenv<real_t>::_default_,
491 // bool clear = true,
492 // unsigned * asupp = 0, unsigned nas = 0,
493 // unsigned * oasup = 0,
494 // unsigned * bsupp = 0, unsigned nbs = 0,
495 // unsigned * obsup = 0 )
496 // {
497 // if ( clear ) for (int d=0;d<3;std::fill(dst[d],dst[d]+oenv->data_size(),(real_t)0.0),d++);
498 // unsigned * loasup, * lobsup; loasup = lobsup = 0;
499 // /* calcul des addresses de sorties de a sur o */
500 // if ( !oasup ) { if ( !asupp ) nas = aenv->data_size(); loasup = new unsigned[ nas ]; oasup = loasup; oaddress( oenv, oasup, aenv, asupp, nas ); };
501 // /* calcul des addresses de sorties de b sur o */
502 // if ( !obsup ) { if ( !bsupp ) nbs = benv->data_size(); lobsup = new unsigned[ nbs ]; obsup = lobsup; oaddress( oenv, obsup, benv, bsupp, nas ); };
503 // /* calcul des coefficients binomiaux (inverses) */
504 // real_t * libzca, * libzcb, * libzco; libzca = libzcb = libzco = 0;
505 // /* calcul des coefficients binomiaux pour aenv */
506 // if ( !ibzca ) { libzca = new real_t[ nas ]; ibzca = libzca; rvbinoms( aenv, ibzca, asupp, nas, bznv ); };
507 // /* calcul des coefficients binomiaux pour benv */
508 // if ( !ibzcb ) { libzcb = new real_t[ nbs ]; ibzcb = libzcb; rvbinoms( benv, ibzcb, bsupp, nbs, bznv ); };
509 // /* calcul des coefficients binomiaux pour oenv */
510 // if ( !ibzco ) { libzco = new real_t[ renv->data_size() ]; ibzco = libzcb; rvbinoms( renv, ibzco, 0, 0, bznv ); };
511 
512 // _bvrcvloop_<real_t,'-'>(oenv,dst[2],aenv,src_a[1],benv,src_b[0],asupp,nas,obsup,bsupp,nbs,obsup);
513 // _bvrcvloop_<real_t,'+'>(oenv,dst[2],aenv,src_a[0],benv,src_b[1],asupp,nas,obsup,bsupp,nbs,obsup);
514 // _bvrcvloop_<real_t,'-'>(oenv,dst[0],aenv,src_a[2],benv,src_b[1],asupp,nas,obsup,bsupp,nbs,obsup);
515 // _bvrcvloop_<real_t,'+'>(oenv,dst[0],aenv,src_a[1],benv,src_b[2],asupp,nas,obsup,bsupp,nbs,obsup);
516 // _bvrcvloop_<real_t,'-'>(oenv,dst[1],aenv,src_a[0],benv,src_b[2],asupp,nas,obsup,bsupp,nbs,obsup);
517 // _bvrcvloop_<real_t,'+'>(oenv,dst[1],aenv,src_a[2],benv,src_b[0],asupp,nas,obsup,bsupp,nbs,obsup);
518 
519 // /* liberation des addresses de sorties */
520 // if ( loasup ) delete[] loasup; if ( lobsup ) delete[] lobsup;
521 // /* liberation des coefficients binomiaux */
522 // if ( libzca ) delete[] libzca; if ( libzcb ) delete[] libzcb; if ( libzco ) delete[] libzco;
523 // };
524 /********************************************************************/
mmx::vctops::set_conversion
unsigned set_conversion(real_t *src, unsigned sz, const real_t &epsilon=(real_t) 0)
Definition: loops_vctops.hpp:206
mmx::vct::fill
void fill(C *a, unsigned n, int s, const C &x)
Definition: tensor_vctops.hpp:152
mmx::b
const C & b
Definition: Interval_glue.hpp:25
mmx::tensor::clear
void clear(monomials< C > &monoms)
Definition: tensor_monomials_fcts.hpp:87
mmx::tensor::eenv
tensor::eenv eenv
Definition: tensor_convert.hpp:12
mmx::sparse::copy
TMPL void copy(Polynomial &r, const Polynomial &a)
Copy of a in r.
Definition: sparse_monomials.hpp:613
op_mul
void op_mul(eenv *a, eenv *b)
Definition: bernstein_eenv_multiply.hpp:37
mmx::abs
void abs(Interval< C, r > &x, const Interval< C, r > &a)
Definition: Interval_fcts.hpp:185
mmx::c
const C & c
Definition: Interval_glue.hpp:45
density
real_t density(real_t *src, const real_t prec=1e-12, unsigned nsmp=0)
Definition: bernstein_eenv_multiply.hpp:8
support
unsigned support(unsigned *supp, real_t *src, const real_t &prec=1e-6)
Definition: bernstein_eenv_multiply.hpp:27
mmx::max
Interval< T, r > max(const Interval< T, r > &a, const Interval< T, r > &b)
Definition: Interval_fcts.hpp:135
assert
#define assert(expr, msg)
Definition: shared_object.hpp:57
Home 
Generated on Sat Apr 4 2015 09:58:00 for realroot- by