source: git/kernel/old.Poly.h @ f7f084

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