#ifndef POLYCPP_HEADER #define POLYCPP_HEADER #include #include #include #include #include #include #include #include #include #include using std::exception; #define BOOST_DISABLE_THREADS class DifferentDomainException: public exception{ }; class ExceptionBasedErrorHandler{ public: static const bool handleErrors=true; static void handleDifferentRing(ring r, ring s){ PrintS("throwing"); throw DifferentDomainException(); } }; //PolyImpl is a 08/15 poly wrapper //Poly wraps around PolyImpl with reference counting using boost class TrivialErrorHandler{ public: static const bool handleErrors=false; static void handleDifferentRing(ring r, ring s){ } }; typedef TrivialErrorHandler MyErrorHandler; class PolyImpl{ template friend class PolyBase; //friend class PolyBase; //friend class PolyBase; //friend class PolyBase; friend class Poly; friend class Vector; //friend class Number; protected: poly getInternalReference() const{ return p; } public: ring getRing() const{ return r.get(); } friend inline bool operator==(const Poly& p1, const Poly& p2); friend inline bool operator==(const Vector& p1, const Vector& p2); friend PolyImpl operator+(const PolyImpl& p1, const PolyImpl& n2); friend PolyImpl operator-(const PolyImpl& p1, const PolyImpl& n2); friend PolyImpl operator/(const PolyImpl& p1, const PolyImpl& n2); friend PolyImpl operator*(const PolyImpl& p1, const PolyImpl& n2); friend bool operator==(const PolyImpl& p1, const PolyImpl& n2); friend PolyImpl operator+(const PolyImpl& p1, int n2); friend PolyImpl operator-(const PolyImpl& p1, int n2); friend PolyImpl operator/(const PolyImpl& p1, int n2); friend PolyImpl operator*(const PolyImpl& p1, int n2); friend bool operator==(const PolyImpl& p1, int n2); Number leadCoef(){ return Number(p->coef,r.get()); } PolyImpl& operator=(const PolyImpl& p2){ //durch Reihenfolge Selbstzuweisungen berücksichtigt if (this==&p2) return *this; poly pc=p_Copy(p2.p,p2.r.get()); if(r!=NULL) p_Delete(&p,r.get()); r=p2.r; p=pc; return *this; } PolyImpl operator-(){ PolyImpl t(*this); t.p=p_Copy(p,r.get()); t.p=p_Neg(t.p,r.get()); return t; } PolyImpl& operator+=(const PolyImpl & p2){ if (r!=p2.r){ Werror("not the same ring"); return *this; } if (this==&p2){ number two=n_Init(2,r.get()); p_Mult_nn(p,two,r.get()); return *this; } p=p_Add_q(p,p_Copy(p2.p,p2.r.get()),r.get()); return *this; } PolyImpl& operator*=(const PolyImpl & p2){ if (r!=p2.r){ Werror("not the same ring"); return *this; } if (this==&p2){ poly pc=p_Copy(p,r.get()); p=p_Mult_q(p,p2.p,r.get()); return *this; } p=p_Mult_q(p,p_Copy(p2.p,p2.r.get()),r.get()); return *this; } PolyImpl& operator*=(const Number & n){ if (r!=n.r){ Werror("not the same ring"); return *this; } p=p_Mult_nn(p,n.n,r.get()); return *this; } PolyImpl& operator-=(const PolyImpl & p2){ if (r!=p2.r){ Werror("not the same ring"); return *this; } if (this==&p2){ p_Delete(&p,r.get()); p=NULL; return *this; } poly pc=p_Copy(p2.p,p2.r.get()); pc=p_Neg(pc,r.get()); p=p_Add_q(p,pc,r.get()); return *this; } PolyImpl& operator=(int n){ p_Delete(&p,r.get()); p=p_ISet(n,r.get()); return *this; } PolyImpl(){ r=currRing; p=NULL; } PolyImpl(const PolyImpl & p){ r=p.r; this->p=p_Copy(p.p,r.get()); } PolyImpl(poly p, intrusive_ptr r){ this->p=p_Copy(p,r.get()); this->r=r; } PolyImpl(poly p, intrusive_ptr r,int){ this->p=p; this->r=r; } PolyImpl(int n, intrusive_ptr r){ this->p=p_ISet(n,r.get()); this->r=r; } PolyImpl(const Number & n){ r=n.r.get(); this->p=p_NSet(n_Copy(n.n,r.get()),r.get()); } explicit PolyImpl(int n){ r=currRing; this->p=p_ISet(n,r.get()); } void print(){ p_Write(p,r.get(),r.get()); } virtual ~PolyImpl(){ if (r!=NULL) p_Delete(&p,r.get()); } protected: poly p; intrusive_ptr r; }; inline PolyImpl operator+(const PolyImpl &p1, const PolyImpl& p2){ PolyImpl erg(p1); erg+=p2; return erg; } inline PolyImpl operator*(const PolyImpl &p1, const PolyImpl& p2){ PolyImpl erg(p1); erg*=p2; return erg; } inline PolyImpl operator-(const PolyImpl &p1, const PolyImpl& p2){ PolyImpl erg(p1); erg-=p2; return erg; } inline PolyImpl operator+(const PolyImpl &p1, int p2){ PolyImpl erg(p1); erg+=PolyImpl(p2,p1.r.get()); return erg; } inline PolyImpl operator*(const PolyImpl &p1, int p2){ PolyImpl erg(p1); erg*=PolyImpl(p2,p1.r.get()); return erg; } inline PolyImpl operator-(const PolyImpl &p1, int p2){ PolyImpl erg(p1); erg-=PolyImpl(p2,p1.r); return erg; } inline PolyImpl operator+(int p1, const PolyImpl& p2){ PolyImpl erg(p2); return erg+=PolyImpl(p1,p2.getRing()); } inline PolyImpl operator*(int p1, const PolyImpl& p2){ PolyImpl erg(p2); return erg*=PolyImpl(p1,p2.getRing()); } using namespace boost; template class ConstTermReference{ private: ring r; poly t; public: operator T() const { return T(p_Head(t,r),r); } ConstTermReference(poly p, ring r){ this->t=p; this->r=r; } bool isConstant() const{ return p_LmIsConstant(t,r); } }; template class PolyInputIterator: public std::iterator,ConstTermReference > { private: poly t; ring r; public: bool operator==(const PolyInputIterator& t2){ return t2.t==t; } bool operator!=(const PolyInputIterator& t2){ return t2.t!=t; } PolyInputIterator(poly p, ring r){ t=p; this->r=r; } PolyInputIterator(const PolyInputIterator& it){ t=it.t; r=it.r; } PolyInputIterator& operator++(){ t=t->next; return *this; } PolyInputIterator operator++(int){ PolyInputIterator it(*this); ++(*this); return it; } const ConstTermReference operator*(){ return ConstTermReference (t,r); } shared_ptr operator->(){ return shared_ptr(new T(p_Head(t,r),r,0)); } }; template class PolyBase{ private: typedef PolyBase ThisType; public: poly as_poly() const{ return p_Copy(ptr->p,ptr->getRing()); } template void checkIsSameRing(T& p){ if (error_handle_traits::handleErrors){ if (p.getRing()!=this->getRing()){ error_handle_traits::handleDifferentRing( this->getRing(), p.getRing() ); } } } typedef create_type_input create_type; typedef PolyInputIterator iterator; Intvec leadExp(){ int nvars=rVar(ptr->r); Intvec res(nvars); for(int i=0;ip,i+1,ptr->getRing()); } return res; } void copy_on_write(){ if (!ptr.unique()){ ptr.reset(new PolyImpl(*ptr)); } } void print() const { ptr->print(); } //* ressource managed by Singular char* c_string() const{ return p_String(ptr->p,ptr->getRing(),ptr->getRing()); } PolyBase(ring r=currRing):ptr(new PolyImpl((poly) NULL,r)){ } PolyBase(const char* c, ring r=currRing):ptr(new PolyImpl((poly)NULL,r)){ //p_Read takes no const so do char* cp=(char*) omalloc((strlen(c)+1)*sizeof(char)); strcpy(cp,c); p_Read(cp,ptr->p,r); omfree(cp); } PolyBase(const PolyBase&p):ptr(p.ptr){ } PolyBase& operator+=(const PolyBase& p2){ checkIsSameRing(p2); copy_on_write(); *ptr += *p2.ptr; return *this; } PolyBase& operator*=(const Poly & p2); PolyBase& operator*=(Number n){ copy_on_write(); *ptr *=n; return *this; } /* void print(){ StringSetS(""); write(); { char* s = StringEndS(); PrintS(s); omFree(s); } }*/ virtual ~PolyBase(){} PolyBase(poly p, ring r):ptr(new PolyImpl(p_Copy(p,r),r)){ } PolyBase(poly p, ring r,int):ptr(new PolyImpl(p,r,0)){ } /*Poly(Poly& p){ ptr=p.ptr; }*/ PolyInputIterator begin(){ return PolyInputIterator(ptr->p,ptr->getRing()); } PolyInputIterator end(){ return PolyInputIterator(NULL, ptr->getRing()); } ring getRing() const{ return ptr->getRing(); } int lmTotalDegree() const{ return pTotaldegree(ptr->p); } Number leadCoef(){ return ptr->leadCoef(); } create_type operator-(){ create_type erg(*this); erg*=Number(-1,ptr->getRing()); return erg; } protected: PolyBase(PolyImpl& impl):ptr(&impl){ } poly getInternalReference(){ return ptr->getInternalReference(); } protected: shared_ptr ptr; }; class Poly: public PolyBase{ private: typedef PolyBase Base; friend class Vector; friend class PolyBase; public: Poly(ring r=currRing):Base ((poly)NULL,r,0){ } Poly(int n, ring r=currRing):Base(*(new PolyImpl(n,r))){ } Poly(const char* c, ring r=currRing):Base(c,r){ } Poly(const Base& p):Base(p){ } Poly(const Number& n):Base(*(new PolyImpl(n))){ } Poly(poly p, ring r):Base(p,r){ } Poly(poly p, ring r, int):Base(p,r,0){ } Poly(const std::vector& v, ring r=currRing):Base(*(new PolyImpl((poly) NULL,r))){ unsigned int i; int s=v.size(); poly p=p_ISet(1,r); for(i=0;i&) (*this))+=p2; return *this; }*/ Poly& operator+=(const Poly& p ){ ((Base&)*this)+=p; return *this; } Poly& operator+=(const Base& p ){ ((Base&)*this)+=p; return *this; } friend inline bool operator==(const Poly& p1, const Poly& p2); }; class Vector: public PolyBase{ private: typedef PolyBase Base; public: Vector(ring r=currRing):Base ((poly)NULL,r,0){ } Vector(int n, ring r=currRing):Base(*(new PolyImpl(n,r))){ } Vector(const char* c, ring r=currRing):Base(c,r){ } Vector(const Base& p):Base(p){ } Vector(poly p, ring r):Base(p,r){ } Vector(poly p, ring r, int):Base(p,r,0){ } Vector(std::vector v, ring r=currRing):Base(*(new PolyImpl((poly) NULL,r))){ unsigned int i; int s=v.size(); poly p=p_ISet(1,r); for(i=0;i&) (*this))+=p2; return *this; }*/ Vector& operator+=(const Vector& p ){ ((Base&)*this)+=p; return *this; } Vector& operator+=(const Base& p ){ ((Base&)*this)+=p; return *this; } friend inline bool operator==(const Vector& p1, const Vector& p2); }; //typedef Poly PolyBase::create_type; /*template inline PolyBase operator+(const PolyBase& p1, const PolyBase& p2){ PolyImpl* res=new PolyImpl(*p1.ptr); *res+=*p2.ptr; return(PolyBase(*res)); }*/ /*template inline PolyBase operator*(const PolyBase& p1, const PolyBase& p2){ PolyImpl* res=new PolyImpl(*p1.ptr); *res *= *p2.ptr; return(PolyBase (*res)); }*/ /*template inline PolyBase operator*(const PolyBase& p1, const Number& n){ PolyBase erg(p1); erg*=n; return erg; }*/ inline Poly operator*(const Poly& p, const Poly& p2){ Poly erg=p; erg*=p2; return erg; } inline Vector operator*(const Number& n, const Vector& v){ Vector res=v; res*=n; return res; } //assumes monomials commute with numbers template inline typename PolyBase::create_type operator* (const Number& n, const PolyBase& p) { typename PolyBase::create_type erg(p); erg*=n; return erg; } inline Vector operator*(const Poly& p, const Vector& v){ Vector res(v); res*=p; return res; } inline Poly operator+(const Poly& p1, const Number& n){ Poly f(p1); f+=n; return f; } inline bool operator==(const Poly& p1, const Poly& p2){ ring r1=p1.getRing(); ring r2=p2.getRing(); if (r1!=r2) return false; return p_EqualPolys(p1.ptr->p,p2.ptr->p,r1); } inline bool operator==(const Vector& p1, const Vector& p2){ ring r1=p1.getRing(); ring r2=p2.getRing(); if (r1!=r2) return false; return p_EqualPolys(p1.ptr->p,p2.ptr->p,r1); } template inline typename PolyBase::create_type operator+ (const PolyBase& b1, const PolyBase& b2) { typename PolyBase::create_type erg(b1); erg+=b2; return erg; } inline Vector unitVector(int i,ring r=currRing){ poly p=p_ISet(1,r); p_SetComp(p,i,r); return Vector(p,r,0); } inline Poly operator*(const Number& n, const Poly & p){ Poly res=p; res*=n; return res; } template inline PolyBase& PolyBase::operator*=(const Poly & p2){ copy_on_write(); *ptr *= *p2.ptr; return *this; } #endif