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