DGtal  0.9.4beta
MPolynomial.h
1 
17 #pragma once
18 
36 #if defined(MPolynomial_RECURSES)
37 #error Recursive header files inclusion detected in MPolynomial.h
38 #else // defined(MPolynomial_RECURSES)
39 
40 #define MPolynomial_RECURSES
41 
42 #if !defined MPolynomial_h
43 
44 #define MPolynomial_h
45 
47 // Inclusions
48 #include <iostream>
49 #include <vector>
50 #include "DGtal/base/Common.h"
52 
53 namespace DGtal
54 {
55 
56  template < int n, typename TRing,
57  typename TAlloc = std::allocator<TRing> >
58  class MPolynomial;
59  template < int N, int n, typename TRing,
60  typename TAlloc = std::allocator<TRing> >
62  template < int n, typename TRing,
63  typename TAlloc = std::allocator<TRing>,
64  typename TX = TRing >
66 
67  template < int n, typename TRing, typename TOwner,
68  typename TAlloc, typename TX >
70 
74  template < typename TRing, typename TAlloc >
79 
88  template < typename TRing, typename TOwner,
89  typename TAlloc, typename TX >
90  class MPolynomialEvaluatorImpl< 1, TRing, TOwner, TAlloc, TX >
91  {
92  public:
93  typedef TRing Ring;
94  typedef TOwner Owner;
95  typedef TAlloc Alloc;
96  typedef TX X;
97 
98  template <int nn, class TT, class AA, class SS>
99  friend class MPolynomialEvaluator;
100 
101  template <int nn, class TT, class HLHL, class AA, class SS>
103 
104 private:
105  const Owner & myOwner;
106  const X & myEvalPoint;
107 
108  inline MPolynomialEvaluatorImpl( const Owner & owner, const X & evalpoint)
109  : myOwner(owner), myEvalPoint(evalpoint)
110  {}
111 
115  class EvalFun
116  {
118 
119  public:
120  inline
122  owner )
123  : myOwner(owner)
124  {}
125 
131  inline
133  {
134  X res = (X) 0;
135  X xx = (X) 1;
136  for ( int i = 0; i < (int) p.myValue.size(); ++i )
137  {
138  res += (X)(Ring) p.myValue[i] * xx;
139  xx = xx * myOwner.myEvalPoint;
140  }
141  return res;
142  }
143  };
144 
145 public:
149  inline operator X() const
150  {
151  X res = (X) 0;
152  myOwner.evaluate( res, EvalFun( *this ) );
153  return res;
154  }
155 
159  inline X operator() () const
160  {
161  return (X)(*this);
162  }
163 
164  };
165 
180  template < int n, typename TRing, typename TOwner,
181  typename TAlloc, typename TX >
182  class MPolynomialEvaluatorImpl
183  {
184  public:
185  typedef TRing Ring;
186  typedef TOwner Owner;
187  typedef TAlloc Alloc;
188  typedef TX X;
191 
198  typedef MPolynomial< n - 1, X,
199  typename Alloc::template rebind<X>::other > MPolyNM1;
200 
201  template<int nn, class TT, class AA, class SS>
202  friend class MPolynomialEvaluator;
203 
204  template<int nn, class TT, class HLHL, class AA, class SS>
206 
207  private:
208  const Owner & myOwner; // The "owner"
209  const X & myEvalPoint; // The evaluation point on this level
210 
211  inline
212  MPolynomialEvaluatorImpl(const Owner & owner, const X & evalpoint)
213  : myOwner(owner), myEvalPoint(evalpoint)
214  {}
215 
221  template < typename XX, typename Fun >
222  class EvalFun
223  {
225  const Fun & myEvalFun;
226 
227  public:
228  inline
230  const Fun & evalfun)
231  : myOwner(owner), myEvalFun(evalfun)
232  {}
233 
240  inline XX operator() ( const MPolynomial<n, Ring, Alloc> & p ) const
241  {
242  XX res = (XX) 0;
243  X xx = (X) 1;
244  for ( int i = 0; i < (int) p.myValue.size(); ++i )
245  {
246  res += myEvalFun(p.myValue[i]) * (XX) xx;
247  xx = xx * myOwner.myEvalPoint;
248  }
249  return res;
250  }
251  };
252 
259  template < typename XX, typename Fun >
260  inline
261  void evaluate( XX & res, const Fun & evalfun ) const
262  {
263  // We have to pass evaluation on to our owner, but give a new
264  // functor which now evaluates polynomials of type poly<n, T>.
265  myOwner.evaluate( res, EvalFun< XX, Fun >( *this, evalfun ) );
266  }
267 
274  class EvalFun2
275  {
277 
278  public:
279  inline
281  : myOwner(owner)
282  {
283  }
284 
291  inline MPolyNM1 operator() (const MPolyN & p) const
292  {
293  MPolyNM1 res;
294  X xx = (X) 1;
295  for ( int i = 0; i < (int) p.myValue.size(); ++i )
296  {
297  res += ( (MPolyNM1) p.myValue[i] ) * xx;
298  xx = xx * myOwner.myEvalPoint;
299  }
300  return res;
301  }
302  };
303 
304  public:
308  inline
309  operator MPolyNM1() const
310  // operator poly<n - 1, S, typename Alloc::template rebind<S>::other>() const
311  {
312  MPolyNM1 res; // missing: determine allocator object
313  // We need to pass evaluation on to our owner
314  myOwner.evaluate( res, EvalFun2( *this ) );
315  return res;
316  }
317 
323  template < typename XX >
324  inline
326  < n - 1, Ring,
328  Alloc, XX > operator() ( const XX & x ) const
329  {
332  Alloc, XX >( *this, x );
333  }
334  };
335 
344  template < typename TRing, typename TAlloc, typename TX >
345  class MPolynomialEvaluator< 1, TRing, TAlloc, TX>
346  {
347  friend class MPolynomial< 1, TRing, TAlloc >;
348  public:
349  typedef TRing Ring;
350  typedef TAlloc Alloc;
351  typedef TX X;
353 
354  private:
355  const MPoly1 & myPoly;
356  const X & myEvalPoint;
357 
358  inline
359  MPolynomialEvaluator( const MPoly1 & poly, const X & evalpoint)
360  : myPoly( poly ), myEvalPoint( evalpoint )
361  {}
362 
363  public:
367  inline
368  operator X() const
369  {
370  X res = (X) 0;
371  X xx = (X) 1;
372  for ( int i = 0; i < (int) myPoly.myValue.size(); ++i )
373  {
374  res += ( (X) (Ring) myPoly.myValue[i] ) * xx;
375  xx = xx * myEvalPoint;
376  }
377  return res;
378  }
379 
383  inline
384  X operator() () const
385  {
386  return (X) (*this);
387  }
388  };
389 
390 
399  template < int n, typename TRing, typename TAlloc, typename TX >
400  class MPolynomialEvaluator
401  {
402  friend class MPolynomial< n, TRing, TAlloc>;
403 
404  template<int nn, class TT, class HLHL, class AA, class SS>
406 
407  public:
408  typedef TRing Ring;
409  typedef TAlloc Alloc;
410  typedef TX X;
412  typedef MPolynomial< n - 1, X,
413  typename Alloc::template rebind<X>::other > MPolyNM1;
414  private:
415  const MPolyN & myPoly;
416  const X & myEvalPoint;
417 
418  inline
419  MPolynomialEvaluator( const MPolyN & poly, const X & evalpoint)
420  : myPoly( poly ), myEvalPoint( evalpoint )
421  {}
422 
429  template < typename XX, typename Fun >
430  inline
431  void evaluate( XX & res, const Fun & evalfun ) const
432  {
433  X xx = (X) 1;
434  for ( int i = 0; i < (int) myPoly.myValue.size(); ++i )
435  {
436  res += evalfun( myPoly.myValue[i] ) * xx;
437  xx = xx * myEvalPoint;
438  }
439  }
440 
441 public:
445  inline
446  operator MPolyNM1() const
447  {
448  MPolyNM1 res( myPoly.getAllocator() );
449  X xx = (X) 1;
450  for ( int i = 0; i < (int) myPoly.myValue.size(); ++i )
451  {
452  res += MPolyNM1( myPoly.myValue[i] ) * xx;
453  xx = xx * myEvalPoint;
454  }
455  return res;
456  }
457 
461  template <typename XX>
462  inline
465  operator() (const XX & x) const
466  {
468  < n - 1, Ring, MPolynomialEvaluator< n, Ring, Alloc, X >, Alloc, XX >
469  ( *this, x );
470  }
471  };
472 
494  template < typename TRing, typename TAlloc >
495  class MPolynomial< 0, TRing, TAlloc >
496  {
497  public:
498  typedef TRing Ring;
499  typedef TAlloc Alloc;
500 
501  private:
502  Alloc myAllocator;
503  Ring myValue;
504 
505  public:
513  inline
514  MPolynomial( const Ring & v = 0, const Alloc & allocator = Alloc() )
515  : myAllocator(allocator), myValue(v)
516  {}
517 
524  inline
525  MPolynomial( const Alloc & allocator )
526  : myAllocator(allocator), myValue( Ring() )
527  {}
528 
532  inline bool isZero() const
533  {
534  return myValue == 0;
535  }
536 
541  inline operator const Ring & () const
542  {
543  return myValue;
544  }
545 
551  inline MPolynomial & operator = (const Ring & v)
552  {
553  myValue = v;
554  return *this;
555  }
556 
561  inline Ring operator() () const
562  {
563  return myValue;
564  }
565 
571  inline MPolynomial operator * (const Ring & v) const
572  {
573  return MPolynomial(myValue * v);
574  }
575 
581  inline MPolynomial operator / (const Ring & v) const
582  {
583  return MPolynomial(myValue / v);
584  }
585 
591  inline MPolynomial operator + (const Ring & v) const
592  {
593  return MPolynomial(myValue + v);
594  }
595 
601  inline MPolynomial operator - (const Ring & v) const
602  {
603  return MPolynomial(myValue - v);
604  }
605 
610  inline MPolynomial operator - () const
611  {
612  return MPolynomial(-myValue);
613  }
614 
620  inline MPolynomial & operator *= (const Ring & v)
621  {
622  myValue *= v;
623  return *this;
624  }
625 
631  inline MPolynomial & operator /= (const Ring & v)
632  {
633  myValue /= v;
634  return *this;
635  }
636 
642  inline MPolynomial & operator += (const Ring & v)
643  {
644  myValue += v;
645  return *this;
646  }
647 
653  inline MPolynomial & operator -= (const Ring & v)
654  {
655  myValue -= v;
656  return *this;
657  }
658 
664  inline bool operator == (const Ring & v) const
665  {
666  return myValue == v;
667  }
668 
674  inline bool operator != (const Ring & v) const
675  {
676  return myValue != v;
677  }
678 
683  void selfDisplay( std::ostream & s, int /*N = 0*/ ) const
684  {
685  s << myValue;
686  }
687 
692  inline void swap( MPolynomial & p )
693  {
694  std::swap(myValue, p.myValue);
695  }
696 
700  Alloc getAllocator() const
701  {
702  return myAllocator;
703  }
704  };
705 
706 
719  template < typename T, typename TAlloc = std::allocator<T>,
720  bool usePointers = false>
721  class IVector
722  {
723  public:
724  typedef TAlloc Alloc;
725  typedef typename std::vector<T, Alloc>::size_type Size;
726  private:
727  std::vector<T, Alloc> myVec;
728 
729  public:
730  IVector( const Alloc & allocator = Alloc() )
731  : myVec(allocator)
732  {}
733 
734  IVector( Size aSize, const Alloc & allocator = Alloc() )
735  : myVec( aSize, T(), allocator )
736  {}
737 
738  IVector( Size aSize, const T & entry, const Alloc & allocator = Alloc() )
739  : myVec(aSize, entry, allocator)
740  {}
741 
742  Size size() const
743  {
744  return myVec.size();
745  }
746 
747  void resize( Size aSize, const T & entry = T() )
748  {
749  myVec.resize(aSize, entry);
750  }
751 
752  const T & operator[] ( Size i ) const
753  {
754  return myVec[i];
755  }
756 
757  T & operator[] ( Size i )
758  {
759  return myVec[i];
760  }
761 
762  const T & back() const
763  {
764  return myVec.back();
765  }
766 
767  T & back()
768  {
769  return myVec.back();
770  }
771 
772  void swap(IVector & v)
773  {
774  myVec.swap( v.myVec );
775  }
776 
777  Alloc get_allocator() const
778  {
779  return myVec.get_allocator();
780  }
781 
782  Alloc getAllocator() const
783  {
784  return myVec.get_allocator();
785  }
786  };
787 
792  template < typename T, typename TAlloc>
793  class IVector< T, TAlloc, true >
794  {
795  public:
796  typedef TAlloc Alloc;
797  typedef typename std::vector<typename Alloc::pointer, typename Alloc::template rebind<typename Alloc::pointer>::other>::size_type Size;
798 
799  private:
800  Alloc myAllocator;
801  std::vector<typename Alloc::pointer, typename Alloc::template rebind<typename Alloc::pointer>::other> myVec;
802 
803  void create( Size begin, Size end, typename Alloc::const_reference entry)
804  {
805  for (Size i = begin; i < end; ++i)
806  {
807  myVec[i] = myAllocator.allocate(sizeof(T));
808  myAllocator.construct(myVec[i], entry);
809  }
810  }
811 
812  void free(Size begin, Size end)
813  {
814  for (Size i = begin; i < end; ++i)
815  {
816  myAllocator.destroy(myVec[i]);
817  myAllocator.deallocate(myVec[i], sizeof(T));
818  }
819  }
820 
821  template<class A>
822  void copy_from(const std::vector<typename Alloc::pointer, A> & source)
823  {
824  for (Size i = 0; i < myVec.size(); ++i)
825  {
826  myVec[i] = myAllocator.allocate(sizeof(T));
827  myAllocator.construct(myVec[i], *source[i]);
828  }
829  }
830 
831  public:
832  IVector(const Alloc & allocator = Alloc())
833  : myAllocator(allocator), myVec(allocator)
834  {}
835 
836  IVector(Size aSize, const Alloc & allocator = Alloc())
837  : myAllocator(allocator), myVec(aSize, 0, allocator)
838  {
839  create(0, aSize, T());
840  }
841 
842  IVector(Size aSize, const T & entry, const Alloc & allocator = Alloc())
843  : myAllocator(allocator), myVec(aSize, 0, allocator)
844  {
845  create(0, aSize, entry);
846  }
847 
848  IVector(const IVector & v)
849  : myVec(v.size())
850  {
851  copy_from(v.myVec);
852  }
853 
855  {
856  free(0, (Size)myVec.size());
857  }
858 
859  IVector & operator = (const IVector & v)
860  {
861  if (&v != this)
862  {
863  free(0, (Size)myVec.size());
864  myVec.resize(v.size());
865  copy_from(v.myVec);
866  }
867  return *this;
868  }
869 
870  Size size() const
871  {
872  return (Size)myVec.size();
873  }
874 
875  void resize(Size aSize, const T & entry = T())
876  {
877  Size oldsize = (Size)myVec.size();
878  if (oldsize > aSize)
879  free(aSize, oldsize);
880  myVec.resize(aSize);
881  if (oldsize < aSize)
882  create(oldsize, aSize, entry);
883  }
884 
885  const T & operator[] (Size i) const
886  {
887  return *myVec[i];
888  }
889 
890  T & operator[] (Size i)
891  {
892  return *myVec[i];
893  }
894 
895  const T & back() const
896  {
897  return *myVec.back();
898  }
899 
900  T & back()
901  {
902  return *myVec.back();
903  }
904 
905  void swap(IVector & v)
906  {
907  myVec.swap(v.myVec);
908  }
909 
910  Alloc get_allocator() const
911  {
912  return myVec.get_allocator();
913  }
914 
915  Alloc getAllocator() const
916  {
917  return myVec.get_allocator();
918  }
919  };
920 
921 
923  // template class MPolynomial
953  template < int n, typename TRing, class TAlloc >
954  class MPolynomial
955  {
956  // ----------------------- friends ---------------------------------------
957 
958  template<int NN, int nn, typename TT, typename AA>
960 
961  friend void euclidDiv<TRing, TAlloc>
963  const MPolynomial<1, TRing, TAlloc> &,
964  MPolynomial<1, TRing, TAlloc> &,
965  MPolynomial<1, TRing, TAlloc> & );
966 
967  template<int nn, typename TT, typename AA, typename SS>
968  friend class MPolynomialEvaluator;
969 
970  template<int nn, typename TT, typename HLHL, typename AA, typename SS>
972 
973  public:
974  typedef TRing Ring;
975  typedef TAlloc Alloc;
976  typedef MPolynomial< n - 1, Ring, Alloc > MPolyNM1;
983  typedef IVector< MPolyNM1,
984  typename Alloc::template rebind<MPolyNM1 >::other, (n>1) >
986  typedef typename Storage::Size Size;
987 
988  // ----------------------- Datas ----------------------------
989 private:
991  static MPolyNM1 myZeroPolynomial;
999 
1003  inline
1004  MPolynomial( bool, Size s, const Alloc & /*allocator*/)
1005  : myValue(s)
1006  {}
1007 
1008 public:
1015  inline void normalize()
1016  {
1017  Size dp1 = myValue.size();
1018  while ( dp1 )
1019  {
1020  if (myValue[dp1 - 1].isZero())
1021  --dp1;
1022  else
1023  break;
1024  }
1025  if (dp1 < myValue.size())
1026  myValue.resize(dp1);
1027  }
1028 
1029  // ----------------------- Standard services ----------------------------
1030  public:
1031 
1035  inline MPolynomial( const Alloc & allocator = Alloc() )
1036  : myValue( allocator )
1037  {}
1038 
1042  inline MPolynomial( const Ring & v, const Alloc & allocator = Alloc() )
1043  : myValue( 1, MPolyNM1( v ), allocator )
1044  {}
1045 
1050  inline MPolynomial( const MPolyNM1 & pdm1,
1051  const Alloc & /*allocator = Alloc()*/ )
1052  : myValue( 1, pdm1 )
1053  {}
1054 
1059  template < typename Ring2, typename Alloc2 >
1061  const Alloc & allocator = Alloc() )
1062  : myValue( p.degree() + 1, MPolyNM1(), allocator )
1063  {
1064  for ( Size i = 0; i < myValue.size(); ++i )
1065  myValue[i] = p[i];
1066  normalize();
1067  }
1068 
1073  template < typename Ring2, typename Alloc2 >
1074  inline
1075  MPolynomial & operator=
1077  {
1078  myValue.resize(p.degree() + 1);
1079  for ( Size i = 0; i < myValue.size(); ++i )
1080  myValue[i] = p[i];
1081  normalize();
1082  return *this;
1083  }
1084 
1089  inline void swap( MPolynomial & p )
1090  {
1091  myValue.swap(p.myValue);
1092  }
1093 
1097  Alloc getAllocator() const
1098  {
1099  return myValue.getAllocator();
1100  }
1101 
1102  // ----------------------- MPolynomial services ----------------------------
1103  public:
1104 
1109  inline int degree() const
1110  {
1111  return (int)(myValue.size() - 1);
1112  }
1113 
1119  inline const MPolyNM1 & leading() const
1120  {
1121  return myValue.size() ? myValue.back() : myZeroPolynomial;
1122  }
1123 
1127  inline bool isZero() const
1128  {
1129  return myValue.size() == 0;
1130  }
1131 
1137  inline MPolyNM1 & operator[] ( Size i )
1138  {
1139  if (i >= myValue.size())
1140  myValue.resize(i + 1);
1141  return myValue[i];
1142  }
1143 
1147  inline const MPolyNM1 & operator[] (Size i) const
1148  {
1149  return i < myValue.size() ? myValue[i] : myZeroPolynomial;
1150  }
1151 
1158  inline
1160  {
1161  return MPolynomialEvaluator<n, Ring, Alloc, Ring>( *this, x );
1162  }
1163 
1171  template <typename Ring2>
1172  inline
1174  {
1176  }
1177 
1178  // ----------------------- Operators --------------------------------------
1179  public:
1180 
1186  inline MPolynomial operator * (const Ring & v) const
1187  {
1188  MPolynomial r(*this);
1189  for ( Size i = 0; i < myValue.size(); ++i )
1190  r[i] *= v;
1191  return r;
1192  }
1193 
1199  inline MPolynomial operator / (const Ring & v) const
1200  {
1201  MPolynomial r(*this);
1202  for ( Size i = 0; i < myValue.size(); ++i )
1203  r[i] /= v;
1204  return r;
1205  }
1206 
1213  {
1214  return *this = *this * p;
1215  }
1216 
1222  inline MPolynomial & operator *= (const Ring & v)
1223  {
1224  MPolynomial r( *this );
1225  for ( Size i = 0; i < myValue.size(); ++i )
1226  myValue[i] *= v;
1227  return *this;
1228  }
1229 
1235  inline MPolynomial & operator /= (const Ring & v)
1236  {
1237  for ( Size i = 0; i < myValue.size(); ++i )
1238  myValue[i] /= v;
1239  return *this;
1240  }
1241 
1248  friend
1249  inline MPolynomial operator * (const Ring & v, const MPolynomial & p)
1250  {
1251  MPolynomial r(p);
1252  for ( Size i = 0; i < p.myValue.size(); ++i )
1253  r[i] *= v;
1254  return r;
1255  }
1256 
1260  inline MPolynomial operator - () const
1261  {
1262  MPolynomial r( true, myValue.size(), getAllocator() );
1263  for ( Size i = 0; i < myValue.size(); ++i )
1264  r[i] = -myValue[i];
1265  return r;
1266  }
1267 
1273  inline MPolynomial operator + (const MPolynomial & q) const
1274  {
1275  MPolynomial r(*this);
1276  if (r.myValue.size() < q.myValue.size())
1277  r.myValue.resize(q.myValue.size());
1278  for ( Size i = 0; i < q.myValue.size(); ++i )
1279  r[i] += q[i];
1280  r.normalize();
1281  return r;
1282  }
1283 
1289  inline MPolynomial operator - (const MPolynomial & q) const
1290  {
1291  MPolynomial r(*this);
1292  if (r.myValue.size() < q.myValue.size())
1293  r.myValue.resize(q.myValue.size());
1294  for ( Size i = 0; i < q.myValue.size(); ++i )
1295  r[i] -= q[i];
1296  r.normalize();
1297  return r;
1298  }
1299 
1306  {
1307  if (myValue.size() < q.myValue.size())
1308  myValue.resize(q.myValue.size());
1309  for ( Size i = 0; i < q.myValue.size(); ++i )
1310  myValue[i] += q[i];
1311  normalize();
1312  return *this;
1313  }
1314 
1321  {
1322  if (myValue.size() < q.myValue.size())
1323  myValue.resize(q.myValue.size());
1324  for ( Size i = 0; i < q.myValue.size(); ++i )
1325  myValue[i] -= q[i];
1326  normalize();
1327  return *this;
1328  }
1329 
1336  inline MPolynomial operator + (const MPolyNM1 & q) const
1337  {
1338  MPolynomial r(*this);
1339  if (r.myValue.size() < 1)
1340  r.myValue.resize(1);
1341  r[0] += q;
1342  r.normalize();
1343  return r;
1344  }
1345 
1352  inline MPolynomial operator - (const MPolyNM1 & q) const
1353  {
1354  MPolynomial r(*this);
1355  if (r.myValue.size() < 1)
1356  r.myValue.resize(1);
1357  r[0] -= q;
1358  r.normalize();
1359  return r;
1360  }
1361 
1367  inline MPolynomial operator + (const Ring & v) const
1368  {
1369  MPolynomial r(*this);
1370  if (r.myValue.size() < 1)
1371  r.myValue.resize(1);
1372  r[0] += v;
1373  r.normalize();
1374  return r;
1375  }
1376 
1382  inline MPolynomial operator - (const Ring & v) const
1383  {
1384  MPolynomial r(*this);
1385  if (r.myValue.size() < 1)
1386  r.myValue.resize(1);
1387  r[0] -= v;
1388  r.normalize();
1389  return r;
1390  }
1391 
1398  inline MPolynomial & operator += (const MPolyNM1 & q)
1399  {
1400  if (myValue.size() < 1)
1401  myValue.resize(1);
1402  myValue[0] += q;
1403  normalize();
1404  return *this;
1405  }
1406 
1413  inline MPolynomial & operator -= (const MPolyNM1 & q)
1414  {
1415  if (myValue.size() < 1)
1416  myValue.resize(1);
1417  myValue[0] -= q;
1418  normalize();
1419  return *this;
1420  }
1421 
1427  inline MPolynomial & operator += (const Ring & v)
1428  {
1429  if (myValue.size() < 1)
1430  myValue.resize(1);
1431  myValue[0] += v;
1432  normalize();
1433  return *this;
1434  }
1435 
1441  inline MPolynomial & operator -= (const Ring & v)
1442  {
1443  if (myValue.size() < 1)
1444  myValue.resize(1);
1445  myValue[0] -= v;
1446  normalize();
1447  return *this;
1448  }
1449 
1457  inline MPolynomial operator * (const MPolynomial & p) const
1458  {
1459  int d = p.degree() + degree();
1460  if (d < -1)
1461  d = -1;
1462  MPolynomial r( true, d + 1, getAllocator() );
1463  if (!isZero() && !p.isZero())
1464  for ( Size i = 0; i < r.myValue.size(); ++i )
1465  for ( Size j = 0; j < myValue.size(); ++j )
1466  if (i < j + p.myValue.size())
1467  r[i] += myValue[j] * p[i - j];
1468  r.normalize();
1469  return r;
1470  }
1471 
1472  // ----------------------- Comparison operators ---------------------------
1473  public:
1474 
1480  inline bool operator == (const MPolynomial & q) const
1481  {
1482  if (myValue.size() != q.myValue.size())
1483  return false;
1484  for (Size i = 0; i < myValue.size(); ++i)
1485  if (myValue[i] != q[i])
1486  return false;
1487  return true;
1488  }
1489 
1495  inline bool operator != (const MPolynomial & q) const
1496  {
1497  return !(*this == q);
1498  }
1499 
1505  inline bool operator == (const Ring & v) const
1506  {
1507  if ((v == 0) && (myValue.size() == 0))
1508  return true;
1509  if (myValue.size() != 1)
1510  return false;
1511  return myValue[0] == v;
1512  }
1513 
1519  inline bool operator != (const Ring & v) const
1520  {
1521  return !(*this == v);
1522  }
1523 
1524 
1525 
1526  // ----------------------- Interface --------------------------------------
1527  public:
1528 
1536  void selfDisplay ( std::ostream & s, int N = n ) const
1537  {
1538  if (isZero())
1539  s << (Ring) 0;
1540  else
1541  {
1542  Size nonzero = 0;
1543  for (Size i = 0; i < myValue.size(); ++i)
1544  if (!myValue[i].isZero())
1545  ++nonzero;
1546  if (nonzero > 1) s << "(";
1547  bool first = true;
1548  for (Size i = 0; i < myValue.size(); ++i)
1549  if (!myValue[i].isZero())
1550  {
1551  if (first) first = false;
1552  else s << " + ";
1553  myValue[i].selfDisplay(s, N);
1554  if (i > 0)
1555  {
1556  s << " ";
1557  s << "X_" << N - n;
1558  if (i > 1) s << "^" << i;
1559  }
1560  }
1561  if (nonzero > 1)
1562  s << ")";
1563  }
1564  }
1565 
1570  bool isValid() const;
1571 
1572 
1573  // ------------------------- Datas --------------------------------------
1574  private:
1575  // ------------------------- Hidden services ----------------------------
1576  protected:
1577 
1578 
1579  }; // end of class MPolynomial
1580 
1581 
1588  template <int N, typename TRing, class TAlloc>
1589  std::ostream&
1590  operator<< ( std::ostream & out,
1591  const MPolynomial<N, TRing, TAlloc> & object );
1592 
1593 
1594  // ------------------------- monomial services ----------------------------
1595 
1602  template <int n, typename Ring, typename Alloc>
1604  {
1605  public:
1607 
1614  get( unsigned int k, unsigned int e )
1615  {
1617  if ( k == 0 )
1618  p[e] = Xe_kComputer<n-1,Ring,Alloc>().get( k-1, e );
1619  else
1620  p[0] = Xe_kComputer<n-1,Ring,Alloc>().get( k-1, e );
1621  p.normalize();
1622  //std::cerr << "Xe_k(" << k << "," << e << ")=" << p << std::endl;
1623  return p;
1624  }
1625 
1626  };
1627 
1628  template <typename Ring, typename Alloc>
1629  class Xe_kComputer<0,Ring,Alloc>
1630  {
1631  public:
1633 
1635  get( unsigned int /*k*/, unsigned int /*e*/ )
1636  {
1638  return p;
1639  }
1640  };
1641 
1651  template <int n, typename Ring, typename Alloc>
1652  inline
1653  MPolynomial<n, Ring, Alloc>
1654  Xe_k( unsigned int k, unsigned int e )
1655  {
1656  return Xe_kComputer<n, Ring, Alloc>().get( k, e );
1657  }
1658 
1667  template <int n, typename Ring>
1668  inline
1669  MPolynomial<n, Ring, std::allocator<Ring> >
1670  Xe_k( unsigned int k, unsigned int e )
1671  {
1672  return Xe_kComputer<n, Ring, std::allocator<Ring> >().get( k, e );
1673  }
1674 
1675 
1683  template <typename Ring, typename Alloc>
1684  inline
1685  MPolynomial<1, Ring, Alloc>
1686  mmonomial(unsigned int e)
1687  {
1689  p[e] = 1;
1690  return p;
1691  }
1692 
1701  template <typename Ring, typename Alloc>
1702  inline
1703  MPolynomial<2, Ring, Alloc>
1704  mmonomial(unsigned int e, unsigned int f)
1705  {
1707  p[e][f] = 1;
1708  return p;
1709  }
1710 
1720  template <typename Ring, typename Alloc>
1721  inline MPolynomial<3, Ring, Alloc>
1722  mmonomial(unsigned int e, unsigned int f, unsigned int g)
1723  {
1725  p[e][f][g] = 1;
1726  return p;
1727  }
1728 
1739  template <typename Ring, typename Alloc>
1740  inline
1741  MPolynomial<4, Ring, Alloc>
1742  mmonomial(unsigned int e, unsigned int f, unsigned int g, unsigned int h)
1743  {
1745  p[e][f][g][h] = 1;
1746  return p;
1747  }
1748 
1755  template <typename Ring>
1756  inline
1757  MPolynomial<1, Ring, std::allocator<Ring> >
1758  mmonomial(unsigned int e)
1759  {
1761  p[e] = 1;
1762  return p;
1763  }
1764 
1772  template <typename Ring>
1773  inline
1774  MPolynomial<2, Ring, std::allocator<Ring> >
1775  mmonomial(unsigned int e, unsigned int f)
1776  {
1778  p[e][f] = 1;
1779  return p;
1780  }
1781 
1790  template <typename Ring>
1791  inline
1792  MPolynomial<3, Ring, std::allocator<Ring> >
1793  mmonomial(unsigned int e, unsigned int f, unsigned int g)
1794  {
1796  p[e][f][g] = 1;
1797  return p;
1798  }
1799 
1809  template <typename Ring>
1810  inline
1811  MPolynomial<4, Ring, std::allocator<Ring> >
1812  mmonomial
1813  (unsigned int e, unsigned int f, unsigned int g, unsigned int h)
1814  {
1816  p[e][f][g][h] = 1;
1817  return p;
1818  }
1819 
1820 
1841  template <int n, typename Ring, class Alloc>
1842  class MPolynomialDerivativeComputer<0, n, Ring, Alloc>
1843  {
1844  public:
1847 
1855  static inline
1856  void computeDerivative( const MPolyN & src, MPolyN & dest)
1857  {
1858  dest.myValue.resize(src.degree() >= 0 ? src.degree() : 0);
1859  for ( int i = 1; i <= src.degree(); ++i )
1860  dest[i - 1] = src[i] * (Ring)i;
1861  dest.normalize();
1862  }
1863  };
1864 
1886  template <int N, int n, typename Ring, typename Alloc>
1887  class MPolynomialDerivativeComputer
1888  {
1889  public:
1892 
1900  static inline
1901  void computeDerivative(const MPolyN & src, MPolyN & dest)
1902  {
1903  dest.myValue.resize(src.degree() + 1);
1904  for ( int i = 0; i <= src.degree(); ++i )
1906  ::computeDerivative( src[i], dest[i] );
1907  dest.normalize();
1908  }
1909  };
1910 
1914  template<typename Ring, class Alloc>
1915  class MPolynomialDerivativeComputer<0, 0, Ring, Alloc>
1916  {
1917  public:
1918  class ERROR_N_must_be_strictly_less_than_n_in_derivative_template;
1920 
1921  static inline
1922  void computeDerivative(const MPoly0 &, MPoly0 &)
1923  {
1924  ERROR_N_must_be_strictly_less_than_n_in_derivative_template();
1925  }
1926  };
1927 
1931  template<int N, typename Ring, class Alloc>
1932  class MPolynomialDerivativeComputer<N, 0, Ring, Alloc>
1933  {
1934  public:
1935  class ERROR_N_must_be_strictly_less_than_n_in_derivative_template;
1937 
1938  static inline
1939  void computeDerivative(const MPoly0 &, MPoly0 &)
1940  {
1941  ERROR_N_must_be_strictly_less_than_n_in_derivative_template();
1942  }
1943  };
1944 
1963  template <int N, int n, typename Ring, typename Alloc>
1964  inline
1965  MPolynomial<n, Ring, Alloc>
1967  {
1970  return res;
1971  }
1972 
1986  template<int N, int n, typename Ring>
1987  inline
1988  MPolynomial<n, Ring, std::allocator<Ring> >
1989  derivative
1990  (const MPolynomial<n,Ring,std::allocator<Ring> > & p)
1991  {
1992  MPolynomial<n, Ring, std::allocator<Ring> > res( p.getAllocator() );
1994  ::computeDerivative( p, res );
1995  return res;
1996  }
1997 
2001  template<typename Ring, typename Alloc>
2002  void
2004  const MPolynomial<1, Ring, Alloc> & g,
2007  {
2008  if (f.degree() < g.degree())
2009  {
2010  // Ignore the trivial case
2011  q = MPolynomial<1, Ring, Alloc>(f.getAllocator());
2012  r = f;
2013  return;
2014  }
2015  q = MPolynomial<1, Ring, Alloc>( true, f.degree() - g.degree() + 1,
2016  f.getAllocator() );
2017  r = f;
2018  for (int i = q.degree(); i >= 0; --i)
2019  {
2020  q[i] = r[i + g.degree()] / g.leading();
2021  for (int j = g.degree(); j >= 0; --j)
2022  r[i + j] -= q[i] * g[j];
2023  }
2024  r.normalize();
2025  // Note that the degree of q is already correct.
2026  }
2027 
2031  template <typename Ring>
2032  void
2033  euclidDiv( const MPolynomial<1, Ring, std::allocator<Ring> > & f,
2034  const MPolynomial<1, Ring, std::allocator<Ring> > & g,
2035  MPolynomial<1, Ring, std::allocator<Ring> > & q,
2036  MPolynomial<1, Ring, std::allocator<Ring> > & r )
2037  {
2038  euclidDiv<Ring, std::allocator<Ring> >(f, g, q, r);
2039  }
2040 
2045  template<typename Ring, typename Alloc>
2046  MPolynomial<1, Ring, Alloc>
2048  const MPolynomial<1, Ring, Alloc> & g )
2049  {
2050  if (f.isZero())
2051  {
2052  if (g.isZero()) return f; // both are zero
2053  else return g / g.leading(); // make g monic
2054  }
2056  d1(f / f.leading()),
2057  d2(g / g.leading()),
2058  q(f.getAllocator()),
2059  r(f.getAllocator());
2060  while (!d2.isZero())
2061  {
2062  euclidDiv(d1, d2, q, r);
2063  d1.swap(d2);
2064  d2 = r;
2065  d2 /= r.leading(); // make r monic
2066  }
2067  return d1;
2068  }
2069 
2074  template<typename Ring>
2075  MPolynomial<1, Ring, std::allocator<Ring> >
2076  gcd( const MPolynomial<1, Ring, std::allocator<Ring> > & f,
2077  const MPolynomial<1, Ring, std::allocator<Ring> > & g )
2078  {
2079  return gcd<Ring, std::allocator<Ring> >(f, g);
2080  }
2081 
2082 } // namespace DGtal
2083 
2084 
2086 // Includes inline functions.
2087 #include "DGtal/math/MPolynomial.ih"
2088 
2089 // //
2091 
2092 #endif // !defined MPolynomial_h
2093 
2094 #undef MPolynomial_RECURSES
2095 #endif // else defined(MPolynomial_RECURSES)
IVector(Size aSize, const Alloc &allocator=Alloc())
Definition: MPolynomial.h:836
IVector< MPolyNM1, typename Alloc::template rebind< MPolyNM1 >::other,(n >1) > Storage
Definition: MPolynomial.h:985
XX operator()(const MPolynomial< n, Ring, Alloc > &p) const
Definition: MPolynomial.h:240
MPolynomial< n, Ring, Alloc > MPolyN
Type for the multivariate polynomial.
Definition: MPolynomial.h:190
const MPoly1 & myPoly
The polynomial in question.
Definition: MPolynomial.h:355
bool isZero() const
Definition: MPolynomial.h:1127
Alloc getAllocator() const
Definition: MPolynomial.h:1097
std::vector< T, Alloc >::size_type Size
Definition: MPolynomial.h:725
const MPolyN & myPoly
The polynomial in question.
Definition: MPolynomial.h:415
void resize(Size aSize, const T &entry=T())
Definition: MPolynomial.h:875
void evaluate(XX &res, const Fun &evalfun) const
Definition: MPolynomial.h:431
MPolynomial operator-() const
Definition: MPolynomial.h:1260
static MPolyNM1 myZeroPolynomial
The zero polynomial of n-1 variables for a n-multivariate polynomial.
Definition: MPolynomial.h:991
MPolynomial< n, Ring, Alloc > Xe_k(unsigned int k, unsigned int e)
Definition: MPolynomial.h:1654
const X & myEvalPoint
the evaluation point
Definition: MPolynomial.h:416
EvalFun2(const MPolynomialEvaluatorImpl< n, Ring, Owner, Alloc, X > &owner)
Definition: MPolynomial.h:280
MPolynomial operator/(const Ring &v) const
Definition: MPolynomial.h:1199
const MPolyNM1 & leading() const
Definition: MPolynomial.h:1119
MPolynomial< 1, Ring, Alloc > mmonomial(unsigned int e)
Definition: MPolynomial.h:1686
void resize(Size aSize, const T &entry=T())
Definition: MPolynomial.h:747
MPolynomial(const Ring &v=0, const Alloc &allocator=Alloc())
Definition: MPolynomial.h:514
MPolynomialEvaluatorImpl< n-1, Ring, MPolynomialEvaluator< n, Ring, Alloc, X >, Alloc, XX > operator()(const XX &x) const
Definition: MPolynomial.h:465
MPolynomial< n-1, X, typename Alloc::template rebind< X >::other > MPolyNM1
Definition: MPolynomial.h:413
Aim: Represents a multivariate polynomial, i.e. an element of , where K is some ring or field...
Definition: MPolynomial.h:58
Alloc get_allocator() const
Definition: MPolynomial.h:777
std::vector< typename Alloc::pointer, typename Alloc::template rebind< typename Alloc::pointer >::other >::size_type Size
Definition: MPolynomial.h:797
EvalFun(const MPolynomialEvaluatorImpl< n, Ring, Owner, Alloc, X > &owner, const Fun &evalfun)
Definition: MPolynomial.h:229
bool isValid() const
MPolynomialEvaluatorImpl(const Owner &owner, const X &evalpoint)
Definition: MPolynomial.h:108
const X & myEvalPoint
The evaluation point on this level.
Definition: MPolynomial.h:106
void selfDisplay(std::ostream &s, int) const
Definition: MPolynomial.h:683
void swap(MPolynomial &p)
Definition: MPolynomial.h:1089
MPolynomial & operator*=(const MPolynomial &p)
Definition: MPolynomial.h:1212
IVector(Size aSize, const Alloc &allocator=Alloc())
Definition: MPolynomial.h:734
MPolynomialEvaluator< n, Ring, Alloc, Ring > operator()(const Ring &x) const
Definition: MPolynomial.h:1159
MPolynomialEvaluatorImpl< n-1, Ring, MPolynomialEvaluatorImpl< n, Ring, Owner, Alloc, X >, Alloc, XX > operator()(const XX &x) const
Definition: MPolynomial.h:328
MPolynomial operator+(const MPolynomial &q) const
Definition: MPolynomial.h:1273
IVector(Size aSize, const T &entry, const Alloc &allocator=Alloc())
Definition: MPolynomial.h:842
void free(Size begin, Size end)
Definition: MPolynomial.h:812
static void computeDerivative(const MPolyN &src, MPolyN &dest)
Definition: MPolynomial.h:1856
MPolynomial(const Alloc &allocator=Alloc())
Definition: MPolynomial.h:1035
MPolynomial(const Ring &v, const Alloc &allocator=Alloc())
Definition: MPolynomial.h:1042
const T & back() const
Definition: MPolynomial.h:762
bool operator==(const MPolynomial &q) const
Definition: MPolynomial.h:1480
MPolynomial< n, Ring, Alloc > MPolyN
Type for polynomial with n variable in the ring Ring.
Definition: MPolynomial.h:1891
MPolynomialEvaluatorImpl(const Owner &owner, const X &evalpoint)
Definition: MPolynomial.h:212
MPolynomial< n-1, X, typename Alloc::template rebind< X >::other > MPolyNM1
Definition: MPolynomial.h:199
void euclidDiv(const MPolynomial< 1, TRing, TAlloc > &f, const MPolynomial< 1, TRing, TAlloc > &g, MPolynomial< 1, TRing, TAlloc > &q, MPolynomial< 1, TRing, TAlloc > &r)
std::ostream & operator<<(std::ostream &out, const ClosedIntegerHalfPlane< TSpace > &object)
const MPolynomialEvaluatorImpl< 1, Ring, Owner, Alloc, X > & myOwner
Definition: MPolynomial.h:117
MPolynomial & operator-=(const MPolynomial &q)
Definition: MPolynomial.h:1320
MPolynomial & operator+=(const MPolynomial &q)
Definition: MPolynomial.h:1305
MPolynomial & operator=(const MPolynomial< n, Ring2, Alloc2 > &p)
Definition: MPolynomial.h:1076
void selfDisplay(std::ostream &s, int N=n) const
Definition: MPolynomial.h:1536
MPolynomial(const MPolynomial< n, Ring2, Alloc2 > &p, const Alloc &allocator=Alloc())
Definition: MPolynomial.h:1060
const X & myEvalPoint
The evaluation point.
Definition: MPolynomial.h:356
MPolynomial< n, Ring, Alloc > MPolyN
Definition: MPolynomial.h:411
static void computeDerivative(const MPolyN &src, MPolyN &dest)
Definition: MPolynomial.h:1901
static void computeDerivative(const MPoly0 &, MPoly0 &)
Definition: MPolynomial.h:1939
IVector(Size aSize, const T &entry, const Alloc &allocator=Alloc())
Definition: MPolynomial.h:738
MPolynomial & operator/=(const Ring &v)
Definition: MPolynomial.h:1235
MPolyNM1 & operator[](Size i)
Definition: MPolynomial.h:1137
MPolynomial< 1, Ring, Alloc > gcd(const MPolynomial< 1, Ring, Alloc > &f, const MPolynomial< 1, Ring, Alloc > &g)
Definition: MPolynomial.h:2047
DGtal is the top-level namespace which contains all DGtal functions and types.
EvalFun(const MPolynomialEvaluatorImpl< 1, Ring, Owner, Alloc, X > &owner)
Definition: MPolynomial.h:121
Size size() const
Definition: MPolynomial.h:742
MPolynomialEvaluator(const MPoly1 &poly, const X &evalpoint)
Definition: MPolynomial.h:359
MPolynomial(const MPolyNM1 &pdm1, const Alloc &)
Definition: MPolynomial.h:1050
void evaluate(XX &res, const Fun &evalfun) const
Definition: MPolynomial.h:261
const MPolynomialEvaluatorImpl< n, Ring, Owner, Alloc, X > & myOwner
Definition: MPolynomial.h:276
const MPolynomialEvaluatorImpl< n, Ring, Owner, Alloc, X > & myOwner
Definition: MPolynomial.h:224
MPolynomial(bool, Size s, const Alloc &)
Definition: MPolynomial.h:1004
MPolynomial< n, Ring, Alloc > MPolyN
Type for polynomial with n variable in the ring Ring.
Definition: MPolynomial.h:1846
Alloc getAllocator() const
Definition: MPolynomial.h:782
void create(Size begin, Size end, typename Alloc::const_reference entry)
Definition: MPolynomial.h:803
int degree() const
Definition: MPolynomial.h:1109
void swap(IVector &v)
Definition: MPolynomial.h:772
MPolynomial< n, Ring, Alloc > derivative(const MPolynomial< n, Ring, Alloc > &p)
Definition: MPolynomial.h:1966
const T & operator[](Size i) const
Definition: MPolynomial.h:752
bool operator!=(const MPolynomial &q) const
Definition: MPolynomial.h:1495
Storage::Size Size
Definition: MPolynomial.h:986
MPolynomial< n-1, Ring, Alloc > MPolyNM1
Definition: MPolynomial.h:976
MPolynomial operator*(const Ring &v) const
Definition: MPolynomial.h:1186
MPolynomial(const Alloc &allocator)
Definition: MPolynomial.h:525
void copy_from(const std::vector< typename Alloc::pointer, A > &source)
Definition: MPolynomial.h:822
IVector(const Alloc &allocator=Alloc())
Definition: MPolynomial.h:832
MPolynomialEvaluator(const MPolyN &poly, const X &evalpoint)
Definition: MPolynomial.h:419
std::vector< T, Alloc > myVec
Definition: MPolynomial.h:727
static void computeDerivative(const MPoly0 &, MPoly0 &)
Definition: MPolynomial.h:1922
std::vector< typename Alloc::pointer, typename Alloc::template rebind< typename Alloc::pointer >::other > myVec
Definition: MPolynomial.h:801
MPolyNM1 operator()(const MPolyN &p) const
Definition: MPolynomial.h:291
IVector(const Alloc &allocator=Alloc())
Definition: MPolynomial.h:730