DGtal  1.5.beta
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)
40 #define MPolynomial_RECURSES
41 
42 #if !defined MPolynomial_h
44 #define MPolynomial_h
45 
46 
47 //David Coeurjolly: too many deprecation warnings
48 //There are pragma pop at the end of this file.
49 #pragma GCC diagnostic push
50 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
51 #pragma clang diagnostic push
52 #pragma clang diagnostic ignored "-Wdeprecated-declarations"
53 
54 
55 
57 // Inclusions
58 #include <iostream>
59 #include <vector>
60 #include "DGtal/base/Common.h"
62 
63 namespace DGtal
64 {
65 
66  template < int n, typename TRing,
67  typename TAlloc = std::allocator<TRing> >
68  class MPolynomial;
69  template < int N, int n, typename TRing,
70  typename TAlloc = std::allocator<TRing> >
71  class MPolynomialDerivativeComputer;
72  template < int n, typename TRing,
73  typename TAlloc = std::allocator<TRing>,
74  typename TX = TRing >
75  class MPolynomialEvaluator;
76 
77  template < int n, typename TRing, typename TOwner,
78  typename TAlloc, typename TX >
79  class MPolynomialEvaluatorImpl;
80 
84  template < typename TRing, typename TAlloc >
89 
98  template < typename TRing, typename TOwner,
99  typename TAlloc, typename TX >
100  class MPolynomialEvaluatorImpl< 1, TRing, TOwner, TAlloc, TX >
101  {
102  public:
103  typedef TRing Ring;
104  typedef TOwner Owner;
105  typedef TAlloc Alloc;
106  typedef TX X;
107 
108  template <int nn, class TT, class AA, class SS>
109  friend class MPolynomialEvaluator;
110 
111  template <int nn, class TT, class HLHL, class AA, class SS>
113 
114 private:
115  const Owner & myOwner;
116  const X & myEvalPoint;
117 
118  inline MPolynomialEvaluatorImpl( const Owner & owner, const X & evalpoint)
119  : myOwner(owner), myEvalPoint(evalpoint)
120  {}
121 
125  class EvalFun
126  {
128 
129  public:
130  inline
132  owner )
133  : myOwner(owner)
134  {}
135 
141  inline
143  {
144  X res = (X) 0;
145  X xx = (X) 1;
146  for ( int i = 0; i < (int) p.myValue.size(); ++i )
147  {
148  res += (X)(Ring) p.myValue[i] * xx;
149  xx = xx * myOwner.myEvalPoint;
150  }
151  return res;
152  }
153  };
154 
155 public:
159  inline operator X() const
160  {
161  X res = (X) 0;
162  myOwner.evaluate( res, EvalFun( *this ) );
163  return res;
164  }
165 
169  inline X operator() () const
170  {
171  return (X)(*this);
172  }
173 
174  };
175 
190  template < int n, typename TRing, typename TOwner,
191  typename TAlloc, typename TX >
193  {
194  public:
195  typedef TRing Ring;
196  typedef TOwner Owner;
197  typedef TAlloc Alloc;
198  typedef TX X;
201 
208  typedef MPolynomial< n - 1, X,
209  typename Alloc::template rebind<X>::other > MPolyNM1;
210 
211  template<int nn, class TT, class AA, class SS>
212  friend class MPolynomialEvaluator;
213 
214  template<int nn, class TT, class HLHL, class AA, class SS>
216 
217  private:
218  const Owner & myOwner; // The "owner"
219  const X & myEvalPoint; // The evaluation point on this level
220 
221  inline
222  MPolynomialEvaluatorImpl(const Owner & owner, const X & evalpoint)
223  : myOwner(owner), myEvalPoint(evalpoint)
224  {}
225 
231  template < typename XX, typename Fun >
232  class EvalFun
233  {
235  const Fun & myEvalFun;
236 
237  public:
238  inline
240  const Fun & evalfun)
241  : myOwner(owner), myEvalFun(evalfun)
242  {}
243 
250  inline XX operator() ( const MPolynomial<n, Ring, Alloc> & p ) const
251  {
252  XX res = (XX) 0;
253  X xx = (X) 1;
254  for ( int i = 0; i < (int) p.myValue.size(); ++i )
255  {
256  res += myEvalFun(p.myValue[i]) * (XX) xx;
257  xx = xx * myOwner.myEvalPoint;
258  }
259  return res;
260  }
261  };
262 
269  template < typename XX, typename Fun >
270  inline
271  void evaluate( XX & res, const Fun & evalfun ) const
272  {
273  // We have to pass evaluation on to our owner, but give a new
274  // functor which now evaluates polynomials of type poly<n, T>.
275  myOwner.evaluate( res, EvalFun< XX, Fun >( *this, evalfun ) );
276  }
277 
284  class EvalFun2
285  {
287 
288  public:
289  inline
291  : myOwner(owner)
292  {
293  }
294 
301  inline MPolyNM1 operator() (const MPolyN & p) const
302  {
303  MPolyNM1 res;
304  X xx = (X) 1;
305  for ( int i = 0; i < (int) p.myValue.size(); ++i )
306  {
307  res += ( (MPolyNM1) p.myValue[i] ) * xx;
308  xx = xx * myOwner.myEvalPoint;
309  }
310  return res;
311  }
312  };
313 
314  public:
318  inline
319  operator MPolyNM1() const
320  // operator poly<n - 1, S, typename Alloc::template rebind<S>::other>() const
321  {
322  MPolyNM1 res; // missing: determine allocator object
323  // We need to pass evaluation on to our owner
324  myOwner.evaluate( res, EvalFun2( *this ) );
325  return res;
326  }
327 
333  template < typename XX >
334  inline
336  < n - 1, Ring,
338  Alloc, XX > operator() ( const XX & x ) const
339  {
342  Alloc, XX >( *this, x );
343  }
344  };
345 
354  template < typename TRing, typename TAlloc, typename TX >
355  class MPolynomialEvaluator< 1, TRing, TAlloc, TX>
356  {
357  friend class MPolynomial< 1, TRing, TAlloc >;
358  public:
359  typedef TRing Ring;
360  typedef TAlloc Alloc;
361  typedef TX X;
363 
364  private:
365  const MPoly1 & myPoly;
366  const X & myEvalPoint;
367 
368  inline
369  MPolynomialEvaluator( const MPoly1 & poly, const X & evalpoint)
370  : myPoly( poly ), myEvalPoint( evalpoint )
371  {}
372 
373  public:
377  inline
378  operator X() const
379  {
380  X res = (X) 0;
381  X xx = (X) 1;
382  for ( int i = 0; i < (int) myPoly.myValue.size(); ++i )
383  {
384  res += ( (X) (Ring) myPoly.myValue[i] ) * xx;
385  xx = xx * myEvalPoint;
386  }
387  return res;
388  }
389 
393  inline
394  X operator() () const
395  {
396  return (X) (*this);
397  }
398  };
399 
400 
409  template < int n, typename TRing, typename TAlloc, typename TX >
411  {
412  friend class MPolynomial< n, TRing, TAlloc>;
413 
414  template<int nn, class TT, class HLHL, class AA, class SS>
416 
417  public:
418  typedef TRing Ring;
419  typedef TAlloc Alloc;
420  typedef TX X;
422  typedef MPolynomial< n - 1, X,
423  typename Alloc::template rebind<X>::other > MPolyNM1;
424  private:
425  const MPolyN & myPoly;
426  const X & myEvalPoint;
427 
428  inline
429  MPolynomialEvaluator( const MPolyN & poly, const X & evalpoint)
430  : myPoly( poly ), myEvalPoint( evalpoint )
431  {}
432 
439  template < typename XX, typename Fun >
440  inline
441  void evaluate( XX & res, const Fun & evalfun ) const
442  {
443  X xx = (X) 1;
444  for ( int i = 0; i < (int) myPoly.myValue.size(); ++i )
445  {
446  res += evalfun( myPoly.myValue[i] ) * xx;
447  xx = xx * myEvalPoint;
448  }
449  }
450 
451 public:
455  inline
456  operator MPolyNM1() const
457  {
458  MPolyNM1 res( myPoly.getAllocator() );
459  X xx = (X) 1;
460  for ( int i = 0; i < (int) myPoly.myValue.size(); ++i )
461  {
462  res += MPolyNM1( myPoly.myValue[i] ) * xx;
463  xx = xx * myEvalPoint;
464  }
465  return res;
466  }
467 
471  template <typename XX>
472  inline
475  operator() (const XX & x) const
476  {
479  ( *this, x );
480  }
481  };
482 
504  template < typename TRing, typename TAlloc >
505  class MPolynomial< 0, TRing, TAlloc >
506  {
507  public:
508  typedef TRing Ring;
509  typedef TAlloc Alloc;
510 
511  private:
514 
515  public:
523  inline
524  MPolynomial( const Ring & v = 0, const Alloc & allocator = Alloc() )
525  : myAllocator(allocator), myValue(v)
526  {}
527 
534  inline
535  MPolynomial( const Alloc & allocator )
536  : myAllocator(allocator), myValue( Ring() )
537  {}
538 
542  inline bool isZero() const
543  {
544  return myValue == 0;
545  }
546 
551  inline operator const Ring & () const
552  {
553  return myValue;
554  }
555 
561  inline MPolynomial & operator = (const Ring & v)
562  {
563  myValue = v;
564  return *this;
565  }
566 
571  inline Ring operator() () const
572  {
573  return myValue;
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 
611  inline MPolynomial operator - (const Ring & v) const
612  {
613  return MPolynomial(myValue - v);
614  }
615 
620  inline MPolynomial operator - () const
621  {
622  return MPolynomial(-myValue);
623  }
624 
630  inline MPolynomial & operator *= (const Ring & v)
631  {
632  myValue *= v;
633  return *this;
634  }
635 
641  inline MPolynomial & operator /= (const Ring & v)
642  {
643  myValue /= v;
644  return *this;
645  }
646 
652  inline MPolynomial & operator += (const Ring & v)
653  {
654  myValue += v;
655  return *this;
656  }
657 
663  inline MPolynomial & operator -= (const Ring & v)
664  {
665  myValue -= v;
666  return *this;
667  }
668 
674  inline bool operator == (const Ring & v) const
675  {
676  return myValue == v;
677  }
678 
684  inline bool operator != (const Ring & v) const
685  {
686  return myValue != v;
687  }
688 
693  void selfDisplay( std::ostream & s, int /*N = 0*/ ) const
694  {
695  s << myValue;
696  }
697 
702  inline void swap( MPolynomial & p )
703  {
704  std::swap(myValue, p.myValue);
705  }
706 
711  {
712  return myAllocator;
713  }
714  };
715 
716 
729  template < typename T, typename TAlloc = std::allocator<T>,
730  bool usePointers = false>
731  class IVector
732  {
733  public:
734  typedef TAlloc Alloc;
735  typedef typename std::vector<T, Alloc>::size_type Size;
736  private:
737  std::vector<T, Alloc> myVec;
738 
739  public:
740  IVector( const Alloc & allocator = Alloc() )
741  : myVec(allocator)
742  {}
743 
744  IVector( Size aSize, const Alloc & allocator = Alloc() )
745  : myVec( aSize, T(), allocator )
746  {}
747 
748  IVector( Size aSize, const T & entry, const Alloc & allocator = Alloc() )
749  : myVec(aSize, entry, allocator)
750  {}
751 
752  Size size() const
753  {
754  return myVec.size();
755  }
756 
757  void resize( Size aSize, const T & entry = T() )
758  {
759  myVec.resize(aSize, entry);
760  }
761 
762  const T & operator[] ( Size i ) const
763  {
764  return myVec[i];
765  }
766 
767  T & operator[] ( Size i )
768  {
769  return myVec[i];
770  }
771 
772  const T & back() const
773  {
774  return myVec.back();
775  }
776 
777  T & back()
778  {
779  return myVec.back();
780  }
781 
782  void swap(IVector & v)
783  {
784  myVec.swap( v.myVec );
785  }
786 
788  {
789  return myVec.get_allocator();
790  }
791 
793  {
794  return myVec.get_allocator();
795  }
796  };
797 
802  template < typename T, typename TAlloc>
803  class IVector< T, TAlloc, true >
804  {
805  public:
806  typedef TAlloc Alloc;
807  typedef typename std::vector<typename Alloc::pointer, typename Alloc::template rebind<typename Alloc::pointer>::other>::size_type Size;
808 
809  private:
811  std::vector<typename Alloc::pointer, typename Alloc::template rebind<typename Alloc::pointer>::other> myVec;
812 
813  void create( Size begin, Size end, typename Alloc::const_reference entry)
814  {
815  for (Size i = begin; i < end; ++i)
816  {
817  myVec[i] = myAllocator.allocate(sizeof(T));
818  myAllocator.construct(myVec[i], entry);
819  }
820  }
821 
822  void free(Size begin, Size end)
823  {
824  for (Size i = begin; i < end; ++i)
825  {
826  myAllocator.destroy(myVec[i]);
827  myAllocator.deallocate(myVec[i], sizeof(T));
828  }
829  }
830 
831  template<class A>
832  void copy_from(const std::vector<typename Alloc::pointer, A> & source)
833  {
834  for (Size i = 0; i < myVec.size(); ++i)
835  {
836  myVec[i] = myAllocator.allocate(sizeof(T));
837  myAllocator.construct(myVec[i], *source[i]);
838  }
839  }
840 
841  public:
842  IVector(const Alloc & allocator = Alloc())
843  : myAllocator(allocator), myVec(allocator)
844  {}
845 
846  IVector(Size aSize, const Alloc & allocator = Alloc())
847  : myAllocator(allocator), myVec(aSize, 0, allocator)
848  {
849  create(0, aSize, T());
850  }
851 
852  IVector(Size aSize, const T & entry, const Alloc & allocator = Alloc())
853  : myAllocator(allocator), myVec(aSize, 0, allocator)
854  {
855  create(0, aSize, entry);
856  }
857 
858  IVector(const IVector & v)
859  : myVec(v.size())
860  {
861  copy_from(v.myVec);
862  }
863 
865  {
866  free(0, (Size)myVec.size());
867  }
868 
869  IVector & operator = (const IVector & v)
870  {
871  if (&v != this)
872  {
873  free(0, (Size)myVec.size());
874  myVec.resize(v.size());
875  copy_from(v.myVec);
876  }
877  return *this;
878  }
879 
880  Size size() const
881  {
882  return (Size)myVec.size();
883  }
884 
885  void resize(Size aSize, const T & entry = T())
886  {
887  Size oldsize = (Size)myVec.size();
888  if (oldsize > aSize)
889  free(aSize, oldsize);
890  myVec.resize(aSize);
891  if (oldsize < aSize)
892  create(oldsize, aSize, entry);
893  }
894 
895  const T & operator[] (Size i) const
896  {
897  return *myVec[i];
898  }
899 
901  {
902  return *myVec[i];
903  }
904 
905  const T & back() const
906  {
907  return *myVec.back();
908  }
909 
910  T & back()
911  {
912  return *myVec.back();
913  }
914 
915  void swap(IVector & v)
916  {
917  myVec.swap(v.myVec);
918  }
919 
921  {
922  return myVec.get_allocator();
923  }
924 
926  {
927  return myVec.get_allocator();
928  }
929  };
930 
931 
933  // template class MPolynomial
963  template < int n, typename TRing, class TAlloc >
965  {
966  // ----------------------- friends ---------------------------------------
967 
968  template<int NN, int nn, typename TT, typename AA>
970 
971  friend void euclidDiv<TRing, TAlloc>
976 
977  template<int nn, typename TT, typename AA, typename SS>
978  friend class MPolynomialEvaluator;
979 
980  template<int nn, typename TT, typename HLHL, typename AA, typename SS>
982 
983  public:
984  typedef TRing Ring;
985  typedef TAlloc Alloc;
986  typedef MPolynomial< n - 1, Ring, Alloc > MPolyNM1;
993  typedef IVector< MPolyNM1,
994  typename Alloc::template rebind<MPolyNM1 >::other, (n>1) >
996  typedef typename Storage::Size Size;
997 
998  // ----------------------- Datas ----------------------------
999 private:
1009 
1013  inline
1014  MPolynomial( bool, Size s, const Alloc & /*allocator*/)
1015  : myValue(s)
1016  {}
1017 
1018 public:
1025  inline void normalize()
1026  {
1027  Size dp1 = myValue.size();
1028  while ( dp1 )
1029  {
1030  if (myValue[dp1 - 1].isZero())
1031  --dp1;
1032  else
1033  break;
1034  }
1035  if (dp1 < myValue.size())
1036  myValue.resize(dp1);
1037  }
1038 
1039  // ----------------------- Standard services ----------------------------
1040  public:
1041 
1045  inline MPolynomial( const Alloc & allocator = Alloc() )
1046  : myValue( allocator )
1047  {}
1048 
1052  inline MPolynomial( const Ring & v, const Alloc & allocator = Alloc() )
1053  : myValue( 1, MPolyNM1( v ), allocator )
1054  {}
1055 
1060  inline MPolynomial( const MPolyNM1 & pdm1,
1061  const Alloc & /*allocator = Alloc()*/ )
1062  : myValue( 1, pdm1 )
1063  {}
1064 
1069  template < typename Ring2, typename Alloc2 >
1071  const Alloc & allocator = Alloc() )
1072  : myValue( p.degree() + 1, MPolyNM1(), allocator )
1073  {
1074  for ( Size i = 0; i < myValue.size(); ++i )
1075  myValue[i] = p[i];
1076  normalize();
1077  }
1078 
1083  template < typename Ring2, typename Alloc2 >
1084  inline
1085  MPolynomial & operator=
1086  (const MPolynomial< n, Ring2, Alloc2> & p )
1087  {
1088  myValue.resize(p.degree() + 1);
1089  for ( Size i = 0; i < myValue.size(); ++i )
1090  myValue[i] = p[i];
1091  normalize();
1092  return *this;
1093  }
1094 
1099  inline void swap( MPolynomial & p )
1100  {
1101  myValue.swap(p.myValue);
1102  }
1103 
1108  {
1109  return myValue.getAllocator();
1110  }
1111 
1112  // ----------------------- MPolynomial services ----------------------------
1113  public:
1114 
1119  inline int degree() const
1120  {
1121  return (int)(myValue.size() - 1);
1122  }
1123 
1129  inline const MPolyNM1 & leading() const
1130  {
1131  return myValue.size() ? myValue.back() : myZeroPolynomial;
1132  }
1133 
1137  inline bool isZero() const
1138  {
1139  return myValue.size() == 0;
1140  }
1141 
1147  inline MPolyNM1 & operator[] ( Size i )
1148  {
1149  if (i >= myValue.size())
1150  myValue.resize(i + 1);
1151  return myValue[i];
1152  }
1153 
1157  inline const MPolyNM1 & operator[] (Size i) const
1158  {
1159  return i < myValue.size() ? myValue[i] : myZeroPolynomial;
1160  }
1161 
1168  inline
1170  {
1171  return MPolynomialEvaluator<n, Ring, Alloc, Ring>( *this, x );
1172  }
1173 
1181  template <typename Ring2>
1182  inline
1184  {
1186  }
1187 
1188  // ----------------------- Operators --------------------------------------
1189  public:
1190 
1196  inline MPolynomial operator * (const Ring & v) const
1197  {
1198  MPolynomial r(*this);
1199  for ( Size i = 0; i < myValue.size(); ++i )
1200  r[i] *= v;
1201  return r;
1202  }
1203 
1209  inline MPolynomial operator / (const Ring & v) const
1210  {
1211  MPolynomial r(*this);
1212  for ( Size i = 0; i < myValue.size(); ++i )
1213  r[i] /= v;
1214  return r;
1215  }
1216 
1223  {
1224  return *this = *this * p;
1225  }
1226 
1232  inline MPolynomial & operator *= (const Ring & v)
1233  {
1234  MPolynomial r( *this );
1235  for ( Size i = 0; i < myValue.size(); ++i )
1236  myValue[i] *= v;
1237  return *this;
1238  }
1239 
1245  inline MPolynomial & operator /= (const Ring & v)
1246  {
1247  for ( Size i = 0; i < myValue.size(); ++i )
1248  myValue[i] /= v;
1249  return *this;
1250  }
1251 
1258  friend
1259  inline MPolynomial operator * (const Ring & v, const MPolynomial & p)
1260  {
1261  MPolynomial r(p);
1262  for ( Size i = 0; i < p.myValue.size(); ++i )
1263  r[i] *= v;
1264  return r;
1265  }
1266 
1270  inline MPolynomial operator - () const
1271  {
1272  MPolynomial r( true, myValue.size(), getAllocator() );
1273  for ( Size i = 0; i < myValue.size(); ++i )
1274  r[i] = -myValue[i];
1275  return r;
1276  }
1277 
1283  inline MPolynomial operator + (const MPolynomial & q) const
1284  {
1285  MPolynomial r(*this);
1286  if (r.myValue.size() < q.myValue.size())
1287  r.myValue.resize(q.myValue.size());
1288  for ( Size i = 0; i < q.myValue.size(); ++i )
1289  r[i] += q[i];
1290  r.normalize();
1291  return r;
1292  }
1293 
1299  inline MPolynomial operator - (const MPolynomial & q) const
1300  {
1301  MPolynomial r(*this);
1302  if (r.myValue.size() < q.myValue.size())
1303  r.myValue.resize(q.myValue.size());
1304  for ( Size i = 0; i < q.myValue.size(); ++i )
1305  r[i] -= q[i];
1306  r.normalize();
1307  return r;
1308  }
1309 
1316  {
1317  if (myValue.size() < q.myValue.size())
1318  myValue.resize(q.myValue.size());
1319  for ( Size i = 0; i < q.myValue.size(); ++i )
1320  myValue[i] += q[i];
1321  normalize();
1322  return *this;
1323  }
1324 
1331  {
1332  if (myValue.size() < q.myValue.size())
1333  myValue.resize(q.myValue.size());
1334  for ( Size i = 0; i < q.myValue.size(); ++i )
1335  myValue[i] -= q[i];
1336  normalize();
1337  return *this;
1338  }
1339 
1346  inline MPolynomial operator + (const MPolyNM1 & q) const
1347  {
1348  MPolynomial r(*this);
1349  if (r.myValue.size() < 1)
1350  r.myValue.resize(1);
1351  r[0] += q;
1352  r.normalize();
1353  return r;
1354  }
1355 
1362  inline MPolynomial operator - (const MPolyNM1 & q) const
1363  {
1364  MPolynomial r(*this);
1365  if (r.myValue.size() < 1)
1366  r.myValue.resize(1);
1367  r[0] -= q;
1368  r.normalize();
1369  return r;
1370  }
1371 
1377  inline MPolynomial operator + (const Ring & v) const
1378  {
1379  MPolynomial r(*this);
1380  if (r.myValue.size() < 1)
1381  r.myValue.resize(1);
1382  r[0] += v;
1383  r.normalize();
1384  return r;
1385  }
1386 
1392  inline MPolynomial operator - (const Ring & v) const
1393  {
1394  MPolynomial r(*this);
1395  if (r.myValue.size() < 1)
1396  r.myValue.resize(1);
1397  r[0] -= v;
1398  r.normalize();
1399  return r;
1400  }
1401 
1408  inline MPolynomial & operator += (const MPolyNM1 & q)
1409  {
1410  if (myValue.size() < 1)
1411  myValue.resize(1);
1412  myValue[0] += q;
1413  normalize();
1414  return *this;
1415  }
1416 
1423  inline MPolynomial & operator -= (const MPolyNM1 & q)
1424  {
1425  if (myValue.size() < 1)
1426  myValue.resize(1);
1427  myValue[0] -= q;
1428  normalize();
1429  return *this;
1430  }
1431 
1437  inline MPolynomial & operator += (const Ring & v)
1438  {
1439  if (myValue.size() < 1)
1440  myValue.resize(1);
1441  myValue[0] += v;
1442  normalize();
1443  return *this;
1444  }
1445 
1451  inline MPolynomial & operator -= (const Ring & v)
1452  {
1453  if (myValue.size() < 1)
1454  myValue.resize(1);
1455  myValue[0] -= v;
1456  normalize();
1457  return *this;
1458  }
1459 
1467  inline MPolynomial operator * (const MPolynomial & p) const
1468  {
1469  int d = p.degree() + degree();
1470  if (d < -1)
1471  d = -1;
1472  MPolynomial r( true, d + 1, getAllocator() );
1473  if (!isZero() && !p.isZero())
1474  for ( Size i = 0; i < r.myValue.size(); ++i )
1475  for ( Size j = 0; j < myValue.size(); ++j )
1476  if (i < j + p.myValue.size())
1477  r[i] += myValue[j] * p[i - j];
1478  r.normalize();
1479  return r;
1480  }
1481 
1482  // ----------------------- Comparison operators ---------------------------
1483  public:
1484 
1490  inline bool operator == (const MPolynomial & q) const
1491  {
1492  if (myValue.size() != q.myValue.size())
1493  return false;
1494  for (Size i = 0; i < myValue.size(); ++i)
1495  if (myValue[i] != q[i])
1496  return false;
1497  return true;
1498  }
1499 
1505  inline bool operator != (const MPolynomial & q) const
1506  {
1507  return !(*this == q);
1508  }
1509 
1515  inline bool operator == (const Ring & v) const
1516  {
1517  if ((v == 0) && (myValue.size() == 0))
1518  return true;
1519  if (myValue.size() != 1)
1520  return false;
1521  return myValue[0] == v;
1522  }
1523 
1529  inline bool operator != (const Ring & v) const
1530  {
1531  return !(*this == v);
1532  }
1533 
1534 
1535 
1536  // ----------------------- Interface --------------------------------------
1537  public:
1538 
1546  void selfDisplay ( std::ostream & s, int N = n ) const
1547  {
1548  if (isZero())
1549  s << (Ring) 0;
1550  else
1551  {
1552  Size nonzero = 0;
1553  for (Size i = 0; i < myValue.size(); ++i)
1554  if (!myValue[i].isZero())
1555  ++nonzero;
1556  if (nonzero > 1) s << "(";
1557  bool first = true;
1558  for (Size i = 0; i < myValue.size(); ++i)
1559  if (!myValue[i].isZero())
1560  {
1561  if (first) first = false;
1562  else s << " + ";
1563  myValue[i].selfDisplay(s, N);
1564  if (i > 0)
1565  {
1566  s << " ";
1567  s << "X_" << N - n;
1568  if (i > 1) s << "^" << i;
1569  }
1570  }
1571  if (nonzero > 1)
1572  s << ")";
1573  }
1574  }
1575 
1580  bool isValid() const;
1581 
1582 
1583  // ------------------------- Datas --------------------------------------
1584  private:
1585  // ------------------------- Hidden services ----------------------------
1586  protected:
1587 
1588 
1589  }; // end of class MPolynomial
1590 
1591 
1598  template <int N, typename TRing, class TAlloc>
1599  std::ostream&
1600  operator<< ( std::ostream & out,
1601  const MPolynomial<N, TRing, TAlloc> & object );
1602 
1603 
1604  // ------------------------- monomial services ----------------------------
1605 
1612  template <int n, typename Ring, typename Alloc>
1614  {
1615  public:
1617 
1624  get( unsigned int k, unsigned int e )
1625  {
1627  if ( k == 0 )
1628  p[e] = Xe_kComputer<n-1,Ring,Alloc>().get( k-1, e );
1629  else
1630  p[0] = Xe_kComputer<n-1,Ring,Alloc>().get( k-1, e );
1631  p.normalize();
1632  //std::cerr << "Xe_k(" << k << "," << e << ")=" << p << std::endl;
1633  return p;
1634  }
1635 
1636  };
1637 
1638  template <typename Ring, typename Alloc>
1639  class Xe_kComputer<0,Ring,Alloc>
1640  {
1641  public:
1643 
1645  get( unsigned int /*k*/, unsigned int /*e*/ )
1646  {
1648  return p;
1649  }
1650  };
1651 
1661  template <int n, typename Ring, typename Alloc>
1662  inline
1663  MPolynomial<n, Ring, Alloc>
1664  Xe_k( unsigned int k, unsigned int e )
1665  {
1666  return Xe_kComputer<n, Ring, Alloc>().get( k, e );
1667  }
1668 
1677  template <int n, typename Ring>
1678  inline
1679  MPolynomial<n, Ring, std::allocator<Ring> >
1680  Xe_k( unsigned int k, unsigned int e )
1681  {
1682  return Xe_kComputer<n, Ring, std::allocator<Ring> >().get( k, e );
1683  }
1684 
1685 
1693  template <typename Ring, typename Alloc>
1694  inline
1695  MPolynomial<1, Ring, Alloc>
1696  mmonomial(unsigned int e)
1697  {
1699  p[e] = 1;
1700  return p;
1701  }
1702 
1711  template <typename Ring, typename Alloc>
1712  inline
1713  MPolynomial<2, Ring, Alloc>
1714  mmonomial(unsigned int e, unsigned int f)
1715  {
1717  p[e][f] = 1;
1718  return p;
1719  }
1720 
1730  template <typename Ring, typename Alloc>
1731  inline MPolynomial<3, Ring, Alloc>
1732  mmonomial(unsigned int e, unsigned int f, unsigned int g)
1733  {
1735  p[e][f][g] = 1;
1736  return p;
1737  }
1738 
1749  template <typename Ring, typename Alloc>
1750  inline
1751  MPolynomial<4, Ring, Alloc>
1752  mmonomial(unsigned int e, unsigned int f, unsigned int g, unsigned int h)
1753  {
1755  p[e][f][g][h] = 1;
1756  return p;
1757  }
1758 
1765  template <typename Ring>
1766  inline
1767  MPolynomial<1, Ring, std::allocator<Ring> >
1768  mmonomial(unsigned int e)
1769  {
1771  p[e] = 1;
1772  return p;
1773  }
1774 
1782  template <typename Ring>
1783  inline
1784  MPolynomial<2, Ring, std::allocator<Ring> >
1785  mmonomial(unsigned int e, unsigned int f)
1786  {
1788  p[e][f] = 1;
1789  return p;
1790  }
1791 
1800  template <typename Ring>
1801  inline
1802  MPolynomial<3, Ring, std::allocator<Ring> >
1803  mmonomial(unsigned int e, unsigned int f, unsigned int g)
1804  {
1806  p[e][f][g] = 1;
1807  return p;
1808  }
1809 
1819  template <typename Ring>
1820  inline
1821  MPolynomial<4, Ring, std::allocator<Ring> >
1823  (unsigned int e, unsigned int f, unsigned int g, unsigned int h)
1824  {
1826  p[e][f][g][h] = 1;
1827  return p;
1828  }
1829 
1830 
1851  template <int n, typename Ring, class Alloc>
1852  class MPolynomialDerivativeComputer<0, n, Ring, Alloc>
1853  {
1854  public:
1857 
1865  static inline
1866  void computeDerivative( const MPolyN & src, MPolyN & dest)
1867  {
1868  dest.myValue.resize(src.degree() >= 0 ? src.degree() : 0);
1869  for ( int i = 1; i <= src.degree(); ++i )
1870  dest[i - 1] = src[i] * (Ring)i;
1871  dest.normalize();
1872  }
1873  };
1874 
1896  template <int N, int n, typename Ring, typename Alloc>
1898  {
1899  public:
1902 
1910  static inline
1911  void computeDerivative(const MPolyN & src, MPolyN & dest)
1912  {
1913  dest.myValue.resize(src.degree() + 1);
1914  for ( int i = 0; i <= src.degree(); ++i )
1916  ::computeDerivative( src[i], dest[i] );
1917  dest.normalize();
1918  }
1919  };
1920 
1924  template<typename Ring, class Alloc>
1925  class MPolynomialDerivativeComputer<0, 0, Ring, Alloc>
1926  {
1927  public:
1928  class ERROR_N_must_be_strictly_less_than_n_in_derivative_template;
1930 
1931  static inline
1933  {
1934  ERROR_N_must_be_strictly_less_than_n_in_derivative_template();
1935  }
1936  };
1937 
1941  template<int N, typename Ring, class Alloc>
1942  class MPolynomialDerivativeComputer<N, 0, Ring, Alloc>
1943  {
1944  public:
1945  class ERROR_N_must_be_strictly_less_than_n_in_derivative_template;
1947 
1948  static inline
1950  {
1951  ERROR_N_must_be_strictly_less_than_n_in_derivative_template();
1952  }
1953  };
1954 
1973  template <int N, int n, typename Ring, typename Alloc>
1974  inline
1975  MPolynomial<n, Ring, Alloc>
1977  {
1978  MPolynomial<n, Ring, Alloc> res( p.getAllocator() );
1980  return res;
1981  }
1982 
1996  template<int N, int n, typename Ring>
1997  inline
1998  MPolynomial<n, Ring, std::allocator<Ring> >
2000  (const MPolynomial<n,Ring,std::allocator<Ring> > & p)
2001  {
2002  MPolynomial<n, Ring, std::allocator<Ring> > res( p.getAllocator() );
2004  ::computeDerivative( p, res );
2005  return res;
2006  }
2007 
2011  template<typename Ring, typename Alloc>
2012  void
2014  const MPolynomial<1, Ring, Alloc> & g,
2017  {
2018  if (f.degree() < g.degree())
2019  {
2020  // Ignore the trivial case
2021  q = MPolynomial<1, Ring, Alloc>(f.getAllocator());
2022  r = f;
2023  return;
2024  }
2025  q = MPolynomial<1, Ring, Alloc>( true, f.degree() - g.degree() + 1,
2026  f.getAllocator() );
2027  r = f;
2028  for (int i = q.degree(); i >= 0; --i)
2029  {
2030  q[i] = r[i + g.degree()] / g.leading();
2031  for (int j = g.degree(); j >= 0; --j)
2032  r[i + j] -= q[i] * g[j];
2033  }
2034  r.normalize();
2035  // Note that the degree of q is already correct.
2036  }
2037 
2041  template <typename Ring>
2042  void
2043  euclidDiv( const MPolynomial<1, Ring, std::allocator<Ring> > & f,
2044  const MPolynomial<1, Ring, std::allocator<Ring> > & g,
2045  MPolynomial<1, Ring, std::allocator<Ring> > & q,
2046  MPolynomial<1, Ring, std::allocator<Ring> > & r )
2047  {
2048  euclidDiv<Ring, std::allocator<Ring> >(f, g, q, r);
2049  }
2050 
2055  template<typename Ring, typename Alloc>
2056  MPolynomial<1, Ring, Alloc>
2058  const MPolynomial<1, Ring, Alloc> & g )
2059  {
2060  if (f.isZero())
2061  {
2062  if (g.isZero()) return f; // both are zero
2063  else return g / g.leading(); // make g monic
2064  }
2066  d1(f / f.leading()),
2067  d2(g / g.leading()),
2068  q(f.getAllocator()),
2069  r(f.getAllocator());
2070  while (!d2.isZero())
2071  {
2072  euclidDiv(d1, d2, q, r);
2073  d1.swap(d2);
2074  d2 = r;
2075  d2 /= r.leading(); // make r monic
2076  }
2077  return d1;
2078  }
2079 
2084  template<typename Ring>
2085  MPolynomial<1, Ring, std::allocator<Ring> >
2086  gcd( const MPolynomial<1, Ring, std::allocator<Ring> > & f,
2087  const MPolynomial<1, Ring, std::allocator<Ring> > & g )
2088  {
2089  return gcd<Ring, std::allocator<Ring> >(f, g);
2090  }
2091 
2092 } // namespace DGtal
2093 
2094 
2096 // Includes inline functions.
2097 #include "DGtal/math/MPolynomial.ih"
2098 
2099 #pragma GCC diagnostic pop
2100 #pragma clang diagnostic pop
2101 
2102 
2103 // //
2105 
2106 #endif // !defined MPolynomial_h
2107 
2108 #undef MPolynomial_RECURSES
2109 #endif // else defined(MPolynomial_RECURSES)
IVector(Size aSize, const Alloc &allocator=Alloc())
Definition: MPolynomial.h:846
IVector(Size aSize, const T &entry, const Alloc &allocator=Alloc())
Definition: MPolynomial.h:852
void copy_from(const std::vector< typename Alloc::pointer, A > &source)
Definition: MPolynomial.h:832
IVector(const Alloc &allocator=Alloc())
Definition: MPolynomial.h:842
void create(Size begin, Size end, typename Alloc::const_reference entry)
Definition: MPolynomial.h:813
void resize(Size aSize, const T &entry=T())
Definition: MPolynomial.h:885
std::vector< typename Alloc::pointer, typename Alloc::template rebind< typename Alloc::pointer >::other > myVec
Definition: MPolynomial.h:811
void free(Size begin, Size end)
Definition: MPolynomial.h:822
std::vector< typename Alloc::pointer, typename Alloc::template rebind< typename Alloc::pointer >::other >::size_type Size
Definition: MPolynomial.h:807
Size size() const
Definition: MPolynomial.h:752
std::vector< T, Alloc >::size_type Size
Definition: MPolynomial.h:735
IVector(const Alloc &allocator=Alloc())
Definition: MPolynomial.h:740
IVector(Size aSize, const Alloc &allocator=Alloc())
Definition: MPolynomial.h:744
void swap(IVector &v)
Definition: MPolynomial.h:782
Alloc get_allocator() const
Definition: MPolynomial.h:787
std::vector< T, Alloc > myVec
Definition: MPolynomial.h:737
const T & operator[](Size i) const
Definition: MPolynomial.h:762
IVector(Size aSize, const T &entry, const Alloc &allocator=Alloc())
Definition: MPolynomial.h:748
Alloc getAllocator() const
Definition: MPolynomial.h:792
const T & back() const
Definition: MPolynomial.h:772
void resize(Size aSize, const T &entry=T())
Definition: MPolynomial.h:757
static void computeDerivative(const MPoly0 &, MPoly0 &)
Definition: MPolynomial.h:1932
static void computeDerivative(const MPolyN &src, MPolyN &dest)
Definition: MPolynomial.h:1866
MPolynomial< n, Ring, Alloc > MPolyN
Type for polynomial with n variable in the ring Ring.
Definition: MPolynomial.h:1856
static void computeDerivative(const MPoly0 &, MPoly0 &)
Definition: MPolynomial.h:1949
MPolynomial< n, Ring, Alloc > MPolyN
Type for polynomial with n variable in the ring Ring.
Definition: MPolynomial.h:1901
static void computeDerivative(const MPolyN &src, MPolyN &dest)
Definition: MPolynomial.h:1911
const MPolynomialEvaluatorImpl< n, Ring, Owner, Alloc, X > & myOwner
Definition: MPolynomial.h:286
MPolyNM1 operator()(const MPolyN &p) const
Definition: MPolynomial.h:301
EvalFun2(const MPolynomialEvaluatorImpl< n, Ring, Owner, Alloc, X > &owner)
Definition: MPolynomial.h:290
const MPolynomialEvaluatorImpl< n, Ring, Owner, Alloc, X > & myOwner
Definition: MPolynomial.h:234
XX operator()(const MPolynomial< n, Ring, Alloc > &p) const
Definition: MPolynomial.h:250
EvalFun(const MPolynomialEvaluatorImpl< n, Ring, Owner, Alloc, X > &owner, const Fun &evalfun)
Definition: MPolynomial.h:239
const MPolynomialEvaluatorImpl< 1, Ring, Owner, Alloc, X > & myOwner
Definition: MPolynomial.h:127
EvalFun(const MPolynomialEvaluatorImpl< 1, Ring, Owner, Alloc, X > &owner)
Definition: MPolynomial.h:131
MPolynomialEvaluatorImpl(const Owner &owner, const X &evalpoint)
Definition: MPolynomial.h:118
const X & myEvalPoint
The evaluation point on this level.
Definition: MPolynomial.h:116
void evaluate(XX &res, const Fun &evalfun) const
Definition: MPolynomial.h:271
MPolynomial< n - 1, X, typename Alloc::template rebind< X >::other > MPolyNM1
Definition: MPolynomial.h:209
MPolynomialEvaluatorImpl< n - 1, Ring, MPolynomialEvaluatorImpl< n, Ring, Owner, Alloc, X >, Alloc, XX > operator()(const XX &x) const
Definition: MPolynomial.h:338
MPolynomial< n, Ring, Alloc > MPolyN
Type for the multivariate polynomial.
Definition: MPolynomial.h:200
MPolynomialEvaluatorImpl(const Owner &owner, const X &evalpoint)
Definition: MPolynomial.h:222
const X & myEvalPoint
The evaluation point.
Definition: MPolynomial.h:366
const MPoly1 & myPoly
The polynomial in question.
Definition: MPolynomial.h:365
MPolynomialEvaluator(const MPoly1 &poly, const X &evalpoint)
Definition: MPolynomial.h:369
MPolynomial< n, Ring, Alloc > MPolyN
Definition: MPolynomial.h:421
const X & myEvalPoint
the evaluation point
Definition: MPolynomial.h:426
MPolynomial< n - 1, X, typename Alloc::template rebind< X >::other > MPolyNM1
Definition: MPolynomial.h:423
MPolynomialEvaluator(const MPolyN &poly, const X &evalpoint)
Definition: MPolynomial.h:429
MPolynomialEvaluatorImpl< n - 1, Ring, MPolynomialEvaluator< n, Ring, Alloc, X >, Alloc, XX > operator()(const XX &x) const
Definition: MPolynomial.h:475
void evaluate(XX &res, const Fun &evalfun) const
Definition: MPolynomial.h:441
const MPolyN & myPoly
The polynomial in question.
Definition: MPolynomial.h:425
MPolynomial(const Alloc &allocator)
Definition: MPolynomial.h:535
MPolynomial(const Ring &v=0, const Alloc &allocator=Alloc())
Definition: MPolynomial.h:524
void selfDisplay(std::ostream &s, int) const
Definition: MPolynomial.h:693
Aim: Represents a multivariate polynomial, i.e. an element of , where K is some ring or field.
Definition: MPolynomial.h:965
MPolynomial(const Ring &v, const Alloc &allocator=Alloc())
Definition: MPolynomial.h:1052
bool operator==(const MPolynomial &q) const
Definition: MPolynomial.h:1490
Storage::Size Size
Definition: MPolynomial.h:996
static MPolyNM1 myZeroPolynomial
The zero polynomial of n-1 variables for a n-multivariate polynomial.
Definition: MPolynomial.h:1001
MPolyNM1 & operator[](Size i)
Definition: MPolynomial.h:1147
void selfDisplay(std::ostream &s, int N=n) const
Definition: MPolynomial.h:1546
MPolynomial(const Alloc &allocator=Alloc())
Definition: MPolynomial.h:1045
MPolynomial & operator+=(const MPolynomial &q)
Definition: MPolynomial.h:1315
MPolynomial & operator/=(const Ring &v)
Definition: MPolynomial.h:1245
MPolynomial operator/(const Ring &v) const
Definition: MPolynomial.h:1209
MPolynomial< n - 1, Ring, Alloc > MPolyNM1
Definition: MPolynomial.h:986
MPolynomial operator-() const
Definition: MPolynomial.h:1270
bool isValid() const
const MPolyNM1 & leading() const
Definition: MPolynomial.h:1129
MPolynomial & operator*=(const MPolynomial &p)
Definition: MPolynomial.h:1222
MPolynomial(const MPolyNM1 &pdm1, const Alloc &)
Definition: MPolynomial.h:1060
Alloc getAllocator() const
Definition: MPolynomial.h:1107
MPolynomialEvaluator< n, Ring, Alloc, Ring > operator()(const Ring &x) const
Definition: MPolynomial.h:1169
MPolynomial operator*(const Ring &v) const
Definition: MPolynomial.h:1196
bool operator!=(const MPolynomial &q) const
Definition: MPolynomial.h:1505
MPolynomial & operator-=(const MPolynomial &q)
Definition: MPolynomial.h:1330
void swap(MPolynomial &p)
Definition: MPolynomial.h:1099
int degree() const
Definition: MPolynomial.h:1119
bool isZero() const
Definition: MPolynomial.h:1137
MPolynomial & operator=(const MPolynomial< n, Ring2, Alloc2 > &p)
Definition: MPolynomial.h:1086
MPolynomial(bool, Size s, const Alloc &)
Definition: MPolynomial.h:1014
MPolynomial(const MPolynomial< n, Ring2, Alloc2 > &p, const Alloc &allocator=Alloc())
Definition: MPolynomial.h:1070
MPolynomial operator+(const MPolynomial &q) const
Definition: MPolynomial.h:1283
MPolynomial< 0, Ring, Alloc > get(unsigned int, unsigned int)
Definition: MPolynomial.h:1645
MPolynomial< n, Ring, Alloc > get(unsigned int k, unsigned int e)
Definition: MPolynomial.h:1624
DGtal is the top-level namespace which contains all DGtal functions and types.
void euclidDiv(const MPolynomial< 1, TRing, TAlloc > &f, const MPolynomial< 1, TRing, TAlloc > &g, MPolynomial< 1, TRing, TAlloc > &q, MPolynomial< 1, TRing, TAlloc > &r)
MPolynomial< 1, Ring, Alloc > gcd(const MPolynomial< 1, Ring, Alloc > &f, const MPolynomial< 1, Ring, Alloc > &g)
Definition: MPolynomial.h:2057
std::ostream & operator<<(std::ostream &out, const ATu0v1< TKSpace, TLinearAlgebra > &object)
MPolynomial< 1, Ring, Alloc > mmonomial(unsigned int e)
Definition: MPolynomial.h:1696
MPolynomial< n, Ring, Alloc > Xe_k(unsigned int k, unsigned int e)
Definition: MPolynomial.h:1664
MPolynomial< n, Ring, Alloc > derivative(const MPolynomial< n, Ring, Alloc > &p)
Definition: MPolynomial.h:1976