source: git/kernel/Poly.h @ a41f3aa

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