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