1 2 3 4 #ifndef POLYCPP_HEADER 5 #define POLYCPP_HEADER 6 #include <kernel/mod2.h> 7 #include <kernel/IIntvec.h> 8 #include <coeffs/numbers.h> 9 #include <kernel/Number.h> 10 #include <kernel/polys.h> 11 #include <polys/monomials/ring.h> 12 13 14 #include <boost/shared_ptr.hpp> 15 16 #include <vector> 17 #include <exception> 18 using std::exception; 19 20 #define BOOST_DISABLE_THREADS 21 22 class DifferentDomainException: public exception{ 23 24 }; 25 class ExceptionBasedErrorHandler{ 26 public: 27 static const bool handleErrors=true; handleDifferentRing(ring r,ring s)28 static void handleDifferentRing(ring r, ring s){ 29 PrintS("throwing"); 30 throw DifferentDomainException(); 31 } 32 }; 33 34 //PolyImpl is a 08/15 poly wrapper 35 //Poly wraps around PolyImpl with reference counting using boost 36 class TrivialErrorHandler{ 37 public: 38 static const bool handleErrors=false; handleDifferentRing(ring r,ring s)39 static void handleDifferentRing(ring r, ring s){ 40 } 41 }; 42 typedef TrivialErrorHandler MyErrorHandler; 43 class PolyImpl{ 44 template <poly_variant,class,class> friend class PolyBase; 45 //friend class PolyBase<POLY_VARIANT_RING,Poly,TrivialErrorHandler>; 46 //friend class PolyBase<POLY_VARIANT_MODUL,Vector, TrivialErrorHandler>; 47 //friend class PolyBase<POLY_VARIANT_MODUL>; 48 friend class Poly; 49 friend class Vector; 50 //friend class Number; 51 protected: getInternalReference()52 poly getInternalReference() const{ 53 return p; 54 } 55 public: getRing()56 ring getRing() const{ 57 return r.get(); 58 } 59 friend inline bool operator==(const Poly& p1, const Poly& p2); 60 friend inline bool operator==(const Vector& p1, const Vector& p2); 61 friend PolyImpl operator+(const PolyImpl& p1, const PolyImpl& n2); 62 friend PolyImpl operator-(const PolyImpl& p1, const PolyImpl& n2); 63 friend PolyImpl operator/(const PolyImpl& p1, const PolyImpl& n2); 64 friend PolyImpl operator*(const PolyImpl& p1, const PolyImpl& n2); 65 friend bool operator==(const PolyImpl& p1, const PolyImpl& n2); 66 friend PolyImpl operator+(const PolyImpl& p1, int n2); 67 friend PolyImpl operator-(const PolyImpl& p1, int n2); 68 friend PolyImpl operator/(const PolyImpl& p1, int n2); 69 friend PolyImpl operator*(const PolyImpl& p1, int n2); 70 friend bool operator==(const PolyImpl& p1, int n2); leadCoef()71 Number leadCoef(){ 72 return Number(p->coef,r.get()); 73 } 74 PolyImpl& operator=(const PolyImpl& p2){ 75 //durch Reihenfolge Selbstzuweisungen ber�cksichtigt 76 if (this==&p2) return *this; 77 poly pc=p_Copy(p2.p,p2.r.get()); 78 if(r!=NULL) 79 p_Delete(&p,r.get()); 80 r=p2.r; 81 p=pc; 82 return *this; 83 } 84 PolyImpl operator-(){ 85 PolyImpl t(*this); 86 t.p=p_Copy(p,r.get()); 87 t.p=p_Neg(t.p,r.get()); 88 return t; 89 } 90 PolyImpl& operator+=(const PolyImpl & p2){ 91 if (r!=p2.r){ 92 WerrorS("not the same ring"); 93 return *this; 94 } 95 if (this==&p2){ 96 number two=n_Init(2,r.get()); 97 p_Mult_nn(p,two,r.get()); 98 return *this; 99 } 100 p=p_Add_q(p,p_Copy(p2.p,p2.r.get()),r.get()); 101 102 return *this; 103 } 104 PolyImpl& operator*=(const PolyImpl & p2){ 105 if (r!=p2.r){ 106 WerrorS("not the same ring"); 107 return *this; 108 } 109 if (this==&p2){ 110 poly pc=p_Copy(p,r.get()); 111 p=p_Mult_q(p,p2.p,r.get()); 112 return *this; 113 } 114 p=p_Mult_q(p,p_Copy(p2.p,p2.r.get()),r.get()); 115 return *this; 116 } 117 PolyImpl& operator*=(const Number & n){ 118 if (r!=n.r){ 119 WerrorS("not the same ring"); 120 return *this; 121 } 122 123 p=p_Mult_nn(p,n.n,r.get()); 124 return *this; 125 } 126 PolyImpl& operator-=(const PolyImpl & p2){ 127 if (r!=p2.r){ 128 WerrorS("not the same ring"); 129 return *this; 130 } 131 if (this==&p2){ 132 p_Delete(&p,r.get()); 133 p=NULL; 134 return *this; 135 } 136 137 poly pc=p_Copy(p2.p,p2.r.get()); 138 pc=p_Neg(pc,r.get()); 139 p=p_Add_q(p,pc,r.get()); 140 141 142 return *this; 143 } 144 145 146 PolyImpl& operator=(int n){ 147 148 p_Delete(&p,r.get()); 149 p=p_ISet(n,r.get()); 150 return *this; 151 152 } 153 154 PolyImpl()155 PolyImpl(){ 156 r=currRing; 157 p=NULL; 158 } PolyImpl(const PolyImpl & p)159 PolyImpl(const PolyImpl & p){ 160 r=p.r; 161 this->p=p_Copy(p.p,r.get()); 162 } PolyImpl(poly p,intrusive_ptr<ip_sring> r)163 PolyImpl(poly p, intrusive_ptr<ip_sring> r){ 164 this->p=p_Copy(p,r.get()); 165 this->r=r; 166 } PolyImpl(poly p,intrusive_ptr<ip_sring> r,int)167 PolyImpl(poly p, intrusive_ptr<ip_sring> r,int){ 168 this->p=p; 169 this->r=r; 170 } PolyImpl(int n,intrusive_ptr<ip_sring> r)171 PolyImpl(int n, intrusive_ptr<ip_sring> r){ 172 this->p=p_ISet(n,r.get()); 173 this->r=r; 174 } PolyImpl(const Number & n)175 PolyImpl(const Number & n){ 176 177 r=n.r.get(); 178 this->p=p_NSet(n_Copy(n.n,r.get()),r.get()); 179 180 } PolyImpl(int n)181 explicit PolyImpl(int n){ 182 r=currRing; 183 this->p=p_ISet(n,r.get()); 184 } print()185 void print(){ 186 p_Write(p,r.get(),r.get()); 187 } 188 ~PolyImpl()189 virtual ~PolyImpl(){ 190 if (r!=NULL) 191 p_Delete(&p,r.get()); 192 } 193 194 protected: 195 poly p; 196 intrusive_ptr<ip_sring> r; 197 198 }; 199 200 inline PolyImpl operator+(const PolyImpl &p1, const PolyImpl& p2){ 201 PolyImpl erg(p1); 202 erg+=p2; 203 return erg; 204 } 205 inline PolyImpl operator*(const PolyImpl &p1, const PolyImpl& p2){ 206 PolyImpl erg(p1); 207 erg*=p2; 208 return erg; 209 } 210 inline PolyImpl operator-(const PolyImpl &p1, const PolyImpl& p2){ 211 PolyImpl erg(p1); 212 erg-=p2; 213 return erg; 214 } 215 216 217 218 inline PolyImpl operator+(const PolyImpl &p1, int p2){ 219 PolyImpl erg(p1); 220 erg+=PolyImpl(p2,p1.r.get()); 221 return erg; 222 } 223 inline PolyImpl operator*(const PolyImpl &p1, int p2){ 224 PolyImpl erg(p1); 225 erg*=PolyImpl(p2,p1.r.get()); 226 return erg; 227 } 228 inline PolyImpl operator-(const PolyImpl &p1, int p2){ 229 PolyImpl erg(p1); 230 erg-=PolyImpl(p2,p1.r); 231 return erg; 232 } 233 234 235 inline PolyImpl operator+(int p1, const PolyImpl& p2){ 236 PolyImpl erg(p2); 237 return erg+=PolyImpl(p1,p2.getRing()); 238 } 239 240 inline PolyImpl operator*(int p1, const PolyImpl& p2){ 241 PolyImpl erg(p2); 242 return erg*=PolyImpl(p1,p2.getRing()); 243 } 244 245 using namespace boost; 246 247 248 template<class T> class ConstTermReference{ 249 private: 250 ring r; 251 poly t; 252 public: T()253 operator T() const { 254 return T(p_Head(t,r),r); 255 } ConstTermReference(poly p,ring r)256 ConstTermReference(poly p, ring r){ 257 this->t=p; 258 this->r=r; 259 } isConstant()260 bool isConstant() const{ 261 return p_LmIsConstant(t,r); 262 } 263 264 }; 265 266 template<class T> class PolyInputIterator: 267 public std::iterator<std::input_iterator_tag,T,int, shared_ptr<const T>,ConstTermReference<T> > 268 { 269 270 271 private: 272 poly t; 273 ring r; 274 public: 275 bool operator==(const PolyInputIterator& t2){ 276 return t2.t==t; 277 } 278 bool operator!=(const PolyInputIterator& t2){ 279 return t2.t!=t; 280 } PolyInputIterator(poly p,ring r)281 PolyInputIterator(poly p, ring r){ 282 t=p; 283 this->r=r; 284 } PolyInputIterator(const PolyInputIterator & it)285 PolyInputIterator(const PolyInputIterator& it){ 286 t=it.t; 287 r=it.r; 288 } 289 PolyInputIterator& operator++(){ 290 t=t->next; 291 return *this; 292 } 293 PolyInputIterator operator++(int){ 294 PolyInputIterator it(*this); 295 ++(*this); 296 return it; 297 } 298 const ConstTermReference<T> operator*(){ 299 return ConstTermReference<T> (t,r); 300 } 301 shared_ptr<const T> operator->(){ 302 return shared_ptr<const T>(new T(p_Head(t,r),r,0)); 303 } 304 305 }; 306 307 template<poly_variant variant, class create_type_input, class error_handle_traits> class PolyBase{ 308 private: 309 typedef PolyBase<variant,create_type_input,error_handle_traits> ThisType; 310 public: as_poly()311 poly as_poly() const{ 312 return p_Copy(ptr->p,ptr->getRing()); 313 } checkIsSameRing(T & p)314 template<class T> void checkIsSameRing(T& p){ 315 if (error_handle_traits::handleErrors){ 316 if (p.getRing()!=this->getRing()){ 317 error_handle_traits::handleDifferentRing( 318 this->getRing(), 319 p.getRing() 320 ); 321 } 322 } 323 } 324 typedef create_type_input create_type; 325 typedef PolyInputIterator<create_type> iterator; leadExp()326 Intvec leadExp(){ 327 int nvars=rVar(ptr->r); 328 Intvec res(nvars); 329 for(int i=0;i<nvars;i++){ 330 res[i]=p_GetExp(ptr->p,i+1,ptr->getRing()); 331 } 332 return res; 333 } copy_on_write()334 void copy_on_write(){ 335 if (!ptr.unique()){ 336 ptr.reset(new PolyImpl(*ptr)); 337 } 338 } print()339 void print() const { 340 ptr->print(); 341 } 342 //* ressource managed by Singular c_string()343 char* c_string() const{ 344 345 return p_String(ptr->p,ptr->getRing(),ptr->getRing()); 346 } 347 348 PolyBase(ring r=currRing):ptr(new PolyImpl((poly) NULL,r)){ 349 } 350 351 PolyBase(const char* c, ring r=currRing):ptr(new PolyImpl((poly)NULL,r)){ 352 //p_Read takes no const so do 353 char* cp=(char*) omalloc((strlen(c)+1)*sizeof(char)); 354 strcpy(cp,c); 355 p_Read(cp,ptr->p,r); 356 omfree(cp); 357 } PolyBase(const PolyBase & p)358 PolyBase(const PolyBase&p):ptr(p.ptr){ 359 } 360 361 362 363 PolyBase& operator+=(const PolyBase& p2){ 364 checkIsSameRing(p2); 365 copy_on_write(); 366 *ptr += *p2.ptr; 367 368 return *this; 369 } 370 PolyBase& operator*=(const Poly & p2); 371 PolyBase& operator*=(Number n){ 372 copy_on_write(); 373 *ptr *=n; 374 375 return *this; 376 } 377 /* void print(){ 378 StringSetS(""); 379 write(); 380 { char* s = StringEndS(); PrintS(s); omFree(s); } 381 }*/ ~PolyBase()382 virtual ~PolyBase(){} PolyBase(poly p,ring r)383 PolyBase(poly p, ring r):ptr(new PolyImpl(p_Copy(p,r),r)){ 384 } PolyBase(poly p,ring r,int)385 PolyBase(poly p, ring r,int):ptr(new PolyImpl(p,r,0)){ 386 } 387 /*Poly(Poly& p){ 388 ptr=p.ptr; 389 }*/ 390 begin()391 PolyInputIterator<create_type> begin(){ 392 return PolyInputIterator<create_type>(ptr->p,ptr->getRing()); 393 } end()394 PolyInputIterator<create_type> end(){ 395 return PolyInputIterator<create_type>(NULL, ptr->getRing()); 396 } getRing()397 ring getRing() const{ 398 return ptr->getRing(); 399 } lmTotalDegree()400 int lmTotalDegree() const{ 401 return pTotaldegree(ptr->p); 402 } leadCoef()403 Number leadCoef(){ 404 return ptr->leadCoef(); 405 } 406 create_type operator-(){ 407 create_type erg(*this); 408 erg*=Number(-1,ptr->getRing()); 409 return erg; 410 } 411 protected: 412 PolyBase(PolyImpl & impl)413 PolyBase(PolyImpl& impl):ptr(&impl){ 414 415 } getInternalReference()416 poly getInternalReference(){ 417 return ptr->getInternalReference(); 418 } 419 protected: 420 421 shared_ptr<PolyImpl> ptr; 422 423 }; 424 425 class Poly: public PolyBase<POLY_VARIANT_RING, Poly, MyErrorHandler>{ 426 private: 427 typedef PolyBase<POLY_VARIANT_RING, Poly,MyErrorHandler> Base; 428 friend class Vector; 429 friend class PolyBase<POLY_VARIANT_MODUL,Vector,MyErrorHandler>; 430 public: 431 432 Poly(ring r=currRing):Base ((poly)NULL,r,0){ 433 } 434 Poly(int n, ring r=currRing):Base(*(new PolyImpl(n,r))){ 435 436 } Base(c,r)437 Poly(const char* c, ring r=currRing):Base(c,r){ 438 439 } Poly(const Base & p)440 Poly(const Base& p):Base(p){ 441 } 442 Poly(const Number & n)443 Poly(const Number& n):Base(*(new PolyImpl(n))){ 444 445 } Poly(poly p,ring r)446 Poly(poly p, ring r):Base(p,r){ 447 448 } Poly(poly p,ring r,int)449 Poly(poly p, ring r, int):Base(p,r,0){ 450 } 451 Poly(const std::vector<int>& v, ring r=currRing):Base(*(new PolyImpl((poly) NULL,r))){ 452 unsigned int i; 453 int s=v.size(); 454 poly p=p_ISet(1,r); 455 for(i=0;i<v.size();i++){ 456 pSetExp(p,i+1,v[i]); 457 } 458 pSetm(p); 459 ptr.reset(new PolyImpl(p,r)); 460 } 461 /* Poly& operator+=(const Number& n){ 462 Poly p2(n); 463 ((PolyBase<POLY_VARIANT_RING, Poly>&) (*this))+=p2; 464 return *this; 465 }*/ 466 Poly& operator+=(const Poly& p ){ 467 468 ((Base&)*this)+=p; 469 return *this; 470 } 471 Poly& operator+=(const Base& p ){ 472 473 ((Base&)*this)+=p; 474 return *this; 475 } 476 friend inline bool operator==(const Poly& p1, const Poly& p2); 477 478 }; 479 class Vector: public PolyBase<POLY_VARIANT_MODUL, Vector, MyErrorHandler>{ 480 private: 481 typedef PolyBase<POLY_VARIANT_MODUL, Vector, MyErrorHandler> Base; 482 public: 483 484 Vector(ring r=currRing):Base ((poly)NULL,r,0){ 485 } 486 Vector(int n, ring r=currRing):Base(*(new PolyImpl(n,r))){ 487 488 } Base(c,r)489 Vector(const char* c, ring r=currRing):Base(c,r){ 490 491 } Vector(const Base & p)492 Vector(const Base& p):Base(p){ 493 } 494 495 Vector(poly p,ring r)496 Vector(poly p, ring r):Base(p,r){ 497 498 } Vector(poly p,ring r,int)499 Vector(poly p, ring r, int):Base(p,r,0){ 500 } 501 Vector(std::vector<int> v, ring r=currRing):Base(*(new PolyImpl((poly) NULL,r))){ 502 unsigned int i; 503 int s=v.size(); 504 poly p=p_ISet(1,r); 505 for(i=0;i<v.size();i++){ 506 pSetExp(p,i+1,v[i]); 507 } 508 pSetm(p); 509 ptr.reset(new PolyImpl(p,r)); 510 } 511 /* Poly& operator+=(const Number& n){ 512 Poly p2(n); 513 ((PolyBase<POLY_VARIANT_MODUL, Poly>&) (*this))+=p2; 514 return *this; 515 }*/ 516 Vector& operator+=(const Vector& p ){ 517 518 ((Base&)*this)+=p; 519 return *this; 520 } 521 Vector& operator+=(const Base& p ){ 522 523 ((Base&)*this)+=p; 524 return *this; 525 } 526 friend inline bool operator==(const Vector& p1, const Vector& p2); 527 }; 528 529 //typedef Poly PolyBase<POLY_VARIANT_RING>::create_type; 530 /*template <poly_variant v, class c> inline PolyBase<v> operator+(const PolyBase<v,c>& p1, const PolyBase<v,c>& p2){ 531 PolyImpl* res=new PolyImpl(*p1.ptr); 532 *res+=*p2.ptr; 533 return(PolyBase<v,c>(*res)); 534 }*/ 535 /*template <poly_variant v> inline PolyBase<v> operator*(const PolyBase<v>& p1, const PolyBase<v>& p2){ 536 PolyImpl* res=new PolyImpl(*p1.ptr); 537 *res *= *p2.ptr; 538 return(PolyBase<v> (*res)); 539 }*/ 540 /*template <class c> inline PolyBase<POLY_VARIANT_MODUL> operator*(const PolyBase<POLY_VARIANT_MODUL>& p1, const Number& n){ 541 PolyBase<POLY_VARIANT_MODUL> erg(p1); 542 erg*=n; 543 return erg; 544 }*/ 545 inline Poly operator*(const Poly& p, const Poly& p2){ 546 Poly erg=p; 547 erg*=p2; 548 return erg; 549 } 550 inline Vector operator*(const Number& n, const Vector& v){ 551 Vector res=v; 552 res*=n; 553 return res; 554 } 555 556 //assumes monomials commute with numbers 557 template <poly_variant variant, class create_type, class error_traits> 558 inline typename PolyBase<variant,create_type, error_traits>::create_type 559 operator* 560 (const Number& n, 561 const PolyBase<variant,create_type, class error_tratis>& p) 562 { 563 typename PolyBase<variant, create_type,error_traits>::create_type erg(p); 564 erg*=n; 565 return erg; 566 } 567 568 inline Vector operator*(const Poly& p, const Vector& v){ 569 Vector res(v); 570 res*=p; 571 return res; 572 } 573 inline Poly operator+(const Poly& p1, const Number& n){ 574 Poly f(p1); 575 f+=n; 576 return f; 577 } 578 inline bool operator==(const Poly& p1, const Poly& p2){ 579 ring r1=p1.getRing(); 580 ring r2=p2.getRing(); 581 if (r1!=r2) return false; 582 return p_EqualPolys(p1.ptr->p,p2.ptr->p,r1); 583 } 584 inline bool operator==(const Vector& p1, const Vector& p2){ 585 ring r1=p1.getRing(); 586 ring r2=p2.getRing(); 587 if (r1!=r2) return false; 588 return p_EqualPolys(p1.ptr->p,p2.ptr->p,r1); 589 } 590 template <poly_variant variant, class create_type,class error_traits> 591 inline typename PolyBase<variant,create_type,error_traits>::create_type 592 operator+ 593 (const PolyBase<variant,create_type,error_traits>& b1, 594 const PolyBase<variant,create_type,error_traits>& b2) 595 { 596 typename PolyBase<variant, create_type, error_traits>::create_type erg(b1); 597 erg+=b2; 598 return erg; 599 } 600 inline Vector unitVector(int i,ring r=currRing){ 601 poly p=p_ISet(1,r); 602 p_SetComp(p,i,r); 603 return Vector(p,r,0); 604 } 605 inline Poly operator*(const Number& n, const Poly & p){ 606 Poly res=p; 607 res*=n; 608 return res; 609 } 610 template <poly_variant variant, class create_type, class error_traits> 611 612 inline PolyBase<variant, create_type, error_traits>& 613 PolyBase<variant, create_type, error_traits>::operator*=(const Poly & p2){ 614 copy_on_write(); 615 *ptr *= *p2.ptr; 616 617 return *this; 618 } 619 #endif 620