00001 #ifndef realroot_SBDSLV_EENV_C
00002 #define realroot_SBDSLV_EENV_C
00003 #include <realroot/tensor_vctops.hpp>
00004 namespace mmx {
00005 namespace tensor
00006 {
00007 inline void eenv::swap (eenv & e) { std::swap (data, e.data); };
00008
00009 inline int eenv::nvr () const { return data[1]; };
00010
00011 inline const int *eenv::szs () const { return data + 2; };
00012 inline int *eenv::szs () { return cdup ().data + 2; };
00013 inline int *eenv::szs_data () { return data + 2; };
00014
00015 inline int eenv::sz (int v) const { return szs ()[v]; };
00016 inline int eenv::sz () const { return str ()[-1]; };
00017
00018 inline int *eenv::vrs () { return szs () + nvr (); };
00019 inline const int *eenv::vrs () const { return szs () + nvr (); };
00020
00021 inline int *eenv::str () { return vrs () + nvr () + 1; };
00022 inline const int *eenv::str () const { return vrs () + nvr () + 1; };
00023
00024 inline int eenv::st (int v) const { return str ()[v]; };
00025
00026 inline int eenv::vr (int v) const { return vrs ()[v]; };
00027
00028 inline int eenv::msz (int nvr_) { return 1 + 1 + nvr_ + nvr_ + nvr_ + 1 + 2 * nvr_; };
00029 inline int eenv::msz () const { return eenv::msz (nvr ()); };
00030
00031 inline bool eenv::operator< (const eenv & b) const { return sz () < b.sz (); };
00032
00033 inline int *eenv::alloc (int nvr_) {
00034
00035
00036 new (data)int[eenv::msz (nvr_)];
00037 return data;
00038 }
00039
00040 inline void eenv::construct (int nvr_, int const *szs_, int const *vrs_,
00041 int *mem)
00042 {
00043 #ifdef EENV_DEBUG
00044 std::cout << "eenv::construct begin\n";
00045 std::cout <<nvr_<<" "<<mem<<" "<<msz (nvr_)<<std::endl;
00046 #endif
00047 data = (mem==0)?alloc(nvr_):mem;
00048 data[0] = 1;
00049 data[1] = nvr_;
00050 if (vrs_)
00051 {
00052 std::pair<int,int> *tmp = new std::pair<int, int>[nvr_];
00053
00054 for ( int i = 0; i < nvr_; i ++ ) {
00055 tmp[i].first = vrs_[i];
00056 tmp[i].second = szs_[i];
00057 }
00058
00059 std::sort(tmp,tmp+nvr_);
00060 for ( int i = 0; i < nvr(); i ++ ) {
00061 data[2+nvr()+i] = tmp[i].first;
00062 data[2+i] = tmp[i].second;
00063 };
00064 cstr ();
00065
00066 delete tmp;
00067 }
00068 else
00069 if (szs_)
00070 {
00071 std::copy (szs_, szs_ + nvr_, data+2);
00072 cstr ();
00073
00074 for (int v = 0; v < nvr () ; data[2+nvr()+v] = v, v++ ) ;
00075 };
00076 };
00077
00078 inline bool eenv::equal (const eenv & a, const eenv & b)
00079 {
00080 if (a.data == b.data) return true;
00081
00082 for (int i = 1; i < 2 + 2 * a.nvr (); i++)
00083 if (a.data[i] != b.data[i]) return false;
00084 return true;
00085 };
00086
00087 inline void eenv::dealloc () { if (data && (--data[0] == 0) ) delete[] data; };
00088
00089 inline eenv & eenv::copy (eenv & e, const eenv & b)
00090 {
00091
00092 e.dealloc ();
00093 e.data = b.data;
00094 if (e.data)
00095 e.data[0]++;
00096 return e;
00097 };
00098
00099 inline eenv & eenv::cdup ()
00100 {
00101 if (data[0] == 1)
00102 return *this;
00103
00104 data[0]--;
00105 int *sv = data;
00106 int _sz = msz();
00107 data = alloc (nvr ());
00108 std::copy (sv, sv + _sz, data);
00109 return *this;
00110 };
00111
00112 inline void eenv::new_copy (eenv & e, const eenv & b)
00113 {
00114 e.dealloc ();
00115 unsigned sz = b.msz();
00116 e.data = new int[sz];
00117 for(unsigned i=0;i<sz;i++) e.data[i]=b.data[i];
00118 e.data[0]=1;
00119 }
00120
00121 inline void eenv::cstr ()
00122 {
00123 #ifdef EENV_DEBUG
00124 std::cout << "eenv::cstr()" << std::endl;
00125 #endif
00126 int i, s = 2+2*nvr()+1;
00127 const int *_sz(szs_data());
00128 for ( i = nvr()-2, data[s+nvr()-1] = 1; i >= -1; data[s+i]=_sz[i+1]*data[s+i+1], i-- ) ;
00129 #ifdef EENV_DEBUG
00130 std::cout << (*this) << std::endl;
00131 #endif
00132
00133 };
00134
00135 inline bool eenv::hasvar (int &lv, int gv) const
00136 {
00137 int g, d, m;
00138 const int * _vr = vrs();
00139 if ((gv < _vr[0]) || (gv > _vr[(d=nvr())-1])) return false;
00140 if ( gv == _vr[g = 0] ) { lv = 0; return true; };
00141 if ( gv == _vr[d=nvr()-1]) { lv = nvr()-1; return true;};
00142 do
00143 {
00144 m = (g + d) / 2;
00145 if ( _vr[m] == gv ) { lv = m; return true; };
00146 if ( _vr[m] < gv ) g = m; else d = m;
00147 }
00148 while (d - g > 1);
00149 return false;
00150 };
00151
00152
00153 inline eenv::~eenv () { dealloc (); };
00154
00155 inline eenv::eenv (int nvr_, int const *szs_, int const *vrs_, int *mem)
00156 :data(new int[msz(nvr_)])
00157 {
00158 construct (nvr_, szs_, vrs_, mem);
00159 };
00160
00161 inline eenv::eenv (const eenv & e ):data(e.data) {
00162
00163 if ( data ) data[0]++;
00164 };
00165 inline eenv & eenv::operator= (const eenv & e) { return copy (*this, e); };
00166 inline eenv eenv::mrgvrs (const eenv & a, const eenv & b)
00167 {
00168
00169 unsigned c = 0;
00170 int va, vb;
00171 int prv = (int) (-1);
00172 int anvr = a.nvr ();
00173 int bnvr = b.nvr ();
00174 const int * avr = a.vrs();
00175 const int * bvr = b.vrs() ;
00176
00177
00178
00179 int tmp[anvr + bnvr];
00180
00181 c = va = vb = 0;
00182
00183 do
00184 {
00185 if ( avr[va] <= bvr[vb] )
00186 {
00187 if ( avr[va] != prv ) tmp[c++] = prv = avr[va];
00188 va++;
00189 continue;
00190 };
00191
00192 if ( bvr[vb] < avr[va] )
00193 {
00194 if ( bvr[vb] != prv ) tmp[c++] = prv = bvr[vb];
00195 vb++;
00196 continue;
00197 };
00198 }
00199 while (va < anvr && vb < bnvr);
00200
00201 if (va < anvr) for (; va < anvr; va++)
00202 if ( avr[va] != prv ) tmp[c++] = avr[va];
00203
00204 if (vb < bnvr) for (; vb < bnvr; vb++)
00205 if ( bvr[vb] != prv ) tmp[c++] = bvr[vb];
00206
00207 eenv oenv (c);
00208
00209
00210 std::copy(tmp,tmp+oenv.nvr(),oenv.vrs());
00211 return oenv;
00212 };
00213
00214
00215 inline eenv eenv::common (const eenv & a, const eenv & b)
00216 {
00217 if ( a == b ) return eenv(a);
00218 int v, lva, lvb;
00219 eenv oenv (mrgvrs (a, b));
00220
00221 const int *ovr = oenv.vrs ();
00222 const int *aszs = a.szs ();
00223 const int *bszs = b.szs ();
00224 int *oszs = oenv.szs ();
00225
00226 for (v = 0; v < oenv.nvr (); v++)
00227 {
00228 if (a.hasvar (lva, ovr[v]))
00229 {
00230 if (b.hasvar (lvb, ovr[v]))
00231 {
00232 oszs[v] = std::max(aszs[lva], bszs[lvb]);
00233 }
00234 else
00235 oszs[v] = aszs[lva];
00236 }
00237 else
00238 {
00239 if (b.hasvar (lvb, ovr[v]))
00240 oszs[v] = bszs[lvb];
00241 else
00242 std::cout << "erreur dans common\n";
00243 };
00244 };
00245
00246 oenv.cstr ();
00247 return oenv;
00248 };
00249
00250 inline eenv eenv::mul (const eenv & a, const eenv & b)
00251 {
00252
00253 int v, lva, lvb;
00254 eenv oenv (mrgvrs (a, b));
00255
00256 int *oszs = oenv.szs ();
00257 const int *ovrs = oenv.vrs ();
00258 const int *aszs = a.szs ();
00259 const int *bszs = b.szs ();
00260
00261 for (v = 0; v < oenv.nvr (); v++)
00262 {
00263 oszs[v] = 0;
00264 if (a.hasvar (lva, ovrs[v]))
00265 oszs[v] += aszs[lva] - 1;
00266 if (b.hasvar (lvb, ovrs[v]))
00267 oszs[v] += bszs[lvb] - 1;
00268 oszs[v]++;
00269 };
00270 oenv.cstr ();
00271 return oenv;
00272 };
00273
00274 inline eenv eenv::diff (const eenv & b, int lv)
00275 {
00276
00277 eenv oenv; eenv::new_copy(oenv,b);
00278 int *oszs = oenv.szs();
00279 oszs[lv] = std::max (oszs[lv] - 1, 1);
00280 oenv.cstr ();
00281 return oenv;
00282 };
00283
00284 inline bool eenv::vmap (int *vmap, const eenv & o, const eenv & i)
00285 {
00286 for (int v = 0; v < i.nvr (); v++)
00287 {
00288 int lv;
00289 if (!o.hasvar (lv, i.vr (v)))
00290 {
00291 return false;
00292 };
00293 vmap[v] = lv;
00294 };
00295 return true;
00296 };
00297
00298
00299 inline bool eenv::oaddress (const eenv & oenv, unsigned *osupp,
00300 const eenv & ienv, unsigned *isupp, unsigned nsp)
00301 {
00302
00303
00304 int v;
00305 unsigned c;
00306 int vmap_[ienv.nvr ()];
00307 unsigned addi, addo;
00308 if (!vmap (vmap_, oenv, ienv))
00309 return false;
00310 if (isupp)
00311 for (c = 0; c < nsp; osupp[c] = addo, c++)
00312 for (addi = isupp[c], addo = 0, v = ienv.nvr () - 1; addi;
00313 addi /= ienv.sz (v), v--)
00314 addo += (addi % ienv.sz (v)) * oenv.st (vmap_[v]);
00315 else
00316 {
00317 const int *istr = ienv.str ();
00318 const int *ostr = oenv.str ();
00319 osupp[0] = 0;
00320 for (v = 0; v < ienv.nvr (); v++)
00321 for (int j = 0; j < istr[-1]; j += istr[v - 1])
00322 for (int i = j + istr[v]; i < j + istr[v - 1]; i += istr[v])
00323 osupp[i] = osupp[i - istr[v]] + ostr[vmap_[v]];
00324 };
00325 return true;
00326 };
00327
00328 inline
00329 std::ostream & operator<< (std::ostream & out, const eenv & env)
00330 {
00331 out << "*eenv*\n";
00332 out << "\tpos: " << env.data << "\n";
00333 out << "\tref: " << env.data[0] << "\n";
00334 out << "\tnvr: " << env.nvr() << "\n";
00335 out << "\tesz: " << env.sz () << "\n";
00336 out << "\tszs: ";
00337 vct::print (env.szs (), env.nvr (), 1, out);
00338 out << "\n";
00339 out << "\tvrs: ";
00340 vct::print (env.vrs (), env.nvr (), 1, out);
00341 out << "\n";
00342 out << "\tstr: ";
00343 vct::print (env.str (), env.nvr (), 1, out);
00344 return out;
00345 };
00346
00347 inline
00348 bool eenv::subset (const eenv & a, const eenv & b)
00349 {
00350 const int *avrs = a.vrs ();
00351 const int *aszs = a.szs ();
00352 const int *bszs = b.szs ();
00353 int lv;
00354 for (int i = 0; i < a.nvr (); i++)
00355 if ((!b.hasvar (lv, avrs[i])) || aszs[i] > bszs[lv])
00356 return false;
00357 return true;
00358 };
00359
00360 inline bool eenv::operator== (const eenv & e) const
00361 {
00362 return eenv::equal (*this, e);
00363 };
00364
00365 inline bool eenv::operator!= (const eenv & e) const
00366 {
00367 return !eenv::equal (*this, e);
00368 };
00369
00370 inline eenv eenv::elevation (const eenv & o, const eenv & i)
00371 {
00372 int nszs[o.nvr ()];
00373 const int * ovrs = o.vrs();
00374 const int * oszs = o.szs ();
00375 const int * iszs = i.szs ();
00376 int lv;
00377 for (int v = 0; v < o.nvr (); v++)
00378 if (i.hasvar (lv, ovrs[v]))
00379 nszs[v] = oszs[v] - iszs[lv] + 1;
00380 else
00381 nszs[v] = oszs[v];
00382 eenv tmp(o.nvr (), nszs, o.vrs ());
00383
00384 return tmp;
00385 };
00386
00387
00388
00389
00390 struct vd
00391 {
00392 int n, d;
00393 vd ():n (-1), d (0)
00394 {
00395 };
00396 };
00397
00398 inline std::ostream & operator<< (std::ostream & o, const vd & vd_)
00399 {
00400 o << "vd(" << vd_.n << ", " << vd_.d << ")";
00401 return o;
00402 };
00403
00404 template < class MONOM > void mescan (int &last, vd * vdeg, const MONOM & m)
00405 {
00406 for (unsigned i = 0; i < m.size (); i++)
00407 {
00408 if (m[i] != 0)
00409 {
00410 if (vdeg[i].n != -1)
00411 vdeg[i].d = std::max ((int) m[i], vdeg[i].d);
00412 else
00413 {
00414 if (last != -1)
00415 {
00416 vdeg[i].n = vdeg[last].n;
00417 vdeg[i].d = m[i];
00418 vdeg[last].n = i;
00419 last = i;
00420 }
00421 else
00422 {
00423 last = i;
00424 vdeg[i].n = i;
00425 vdeg[i].d = m[i];
00426 };
00427 };
00428 };
00429 };
00430 }
00431
00432
00433 struct oulala {
00434 int initial; int name;
00435 bool operator<( const oulala& b ) const { return name < b.name; };
00436 };
00437
00438 template < class C, class O >
00439 inline eenv::eenv (const sparse::monomial_seq< C, O > &mpol)
00440 {
00441 int nvr = nbvar (mpol);
00442 vd *lenv = new vd[nvr];
00443 int last = -1;
00444 typename sparse::monomial_seq < C, O >::const_iterator it;
00445 for (it = mpol.begin (); it != mpol.end (); mescan (last, lenv, *it++)) ;
00446 if ( last == -1 )
00447 {
00448 nvr = 1; int szs = 1;
00449 new (this) eenv(1,&szs,0);
00450 return;
00451 };
00452 int c = last;
00453 nvr = 0;
00454 do
00455 {
00456 c = lenv[c].n;
00457 nvr++;
00458 }
00459 while (c != last);
00460 int vrs[nvr];
00461 int szs[nvr];
00462 c = last;
00463
00464 int k = 0;
00465 do
00466 {
00467 vrs[k] = c;
00468 szs[k] = lenv[c].d + 1;
00469 c = lenv[c].n;
00470 lenv[vrs[k]].n = -1;
00471 lenv[vrs[k]].d = 0;
00472 k++;
00473 }
00474 while (c != last);
00475 new (this) eenv (nvr, szs, vrs );
00476
00477 delete lenv;
00478 };
00479
00480 inline bool eenv::equiv (const eenv & a, const eenv & b)
00481 {
00482 return equal(a,b);
00483 };
00484
00485
00486 inline eenv eenv::rewrite( const eenv& o, const eenv& i )
00487 {
00488 int vmap_[ o.nvr() ];
00489 if ( ! vmap( vmap_, o, i ) ) { std::cerr << "vmap error(rewrite)\n"; exit(1); };
00490 eenv tmp(o);
00491 for ( int v = 0; v < o.nvr(); v ++ ) tmp.szs()[v] = 1;
00492 for ( int v = 0; v < i.nvr(); v ++ ) tmp.szs()[vmap_[v]] = i.szs()[v];
00493 tmp.cstr();
00494 return tmp;
00495 };
00496
00497 }
00498
00499 }
00500 #endif