Loading [MathJax]/extensions/TeX/AMSsymbols.js
DGtal 2.0.0
MPolynomial.h
1
16
17#pragma once
18
35
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// std::allocator_traits (as a replacement of removes Alloc::rebind)
61#include <memory>
62#include "DGtal/base/Common.h"
64
65namespace DGtal
66{
67
68 template < int n, typename TRing,
69 typename TAlloc = std::allocator<TRing> >
70 class MPolynomial;
71 template < int N, int n, typename TRing,
72 typename TAlloc = std::allocator<TRing> >
74 template < int n, typename TRing,
75 typename TAlloc = std::allocator<TRing>,
76 typename TX = TRing >
78
79 template < int n, typename TRing, typename TOwner,
80 typename TAlloc, typename TX >
82
86 template < typename TRing, typename TAlloc >
91
100 template < typename TRing, typename TOwner,
101 typename TAlloc, typename TX >
102 class MPolynomialEvaluatorImpl< 1, TRing, TOwner, TAlloc, TX >
103 {
104 public:
105 typedef TRing Ring;
106 typedef TOwner Owner;
107 typedef TAlloc Alloc;
108 typedef TX X;
109
110 template <int nn, class TT, class AA, class SS>
112
113 template <int nn, class TT, class HLHL, class AA, class SS>
115
116private:
117 const Owner & myOwner;
118 const X & myEvalPoint;
119
120 inline MPolynomialEvaluatorImpl( const Owner & owner, const X & evalpoint)
121 : myOwner(owner), myEvalPoint(evalpoint)
122 {}
123
128 {
130
131 public:
132 inline
137
143 inline
145 {
146 X res = (X) 0;
147 X xx = (X) 1;
148 for ( int i = 0; i < (int) p.myValue.size(); ++i )
149 {
150 res += (X)(Ring) p.myValue[i] * xx;
151 xx = xx * myOwner.myEvalPoint;
152 }
153 return res;
154 }
155 };
156
157public:
161 inline operator X() const
162 {
163 X res = (X) 0;
164 myOwner.evaluate( res, EvalFun( *this ) );
165 return res;
166 }
167
171 inline X operator() () const
172 {
173 return (X)(*this);
174 }
175
176 };
177
192 template < int n, typename TRing, typename TOwner,
193 typename TAlloc, typename TX >
195 {
196 public:
197 typedef TRing Ring;
198 typedef TOwner Owner;
199 typedef TAlloc Alloc;
200 typedef TX X;
203
210 typedef MPolynomial< n - 1, X,
211 typename std::allocator_traits<Alloc>::template rebind_alloc<X> > MPolyNM1;
212
213 template<int nn, class TT, class AA, class SS>
215
216 template<int nn, class TT, class HLHL, class AA, class SS>
218
219 private:
220 const Owner & myOwner; // The "owner"
221 const X & myEvalPoint; // The evaluation point on this level
222
223 inline
224 MPolynomialEvaluatorImpl(const Owner & owner, const X & evalpoint)
225 : myOwner(owner), myEvalPoint(evalpoint)
226 {}
227
233 template < typename XX, typename Fun >
235 {
237 const Fun & myEvalFun;
238
239 public:
240 inline
242 const Fun & evalfun)
243 : myOwner(owner), myEvalFun(evalfun)
244 {}
245
252 inline XX operator() ( const MPolynomial<n, Ring, Alloc> & p ) const
253 {
254 XX res = (XX) 0;
255 X xx = (X) 1;
256 for ( int i = 0; i < (int) p.myValue.size(); ++i )
257 {
258 res += myEvalFun(p.myValue[i]) * (XX) xx;
259 xx = xx * myOwner.myEvalPoint;
260 }
261 return res;
262 }
263 };
264
271 template < typename XX, typename Fun >
272 inline
273 void evaluate( XX & res, const Fun & evalfun ) const
274 {
275 // We have to pass evaluation on to our owner, but give a new
276 // functor which now evaluates polynomials of type poly<n, T>.
277 myOwner.evaluate( res, EvalFun< XX, Fun >( *this, evalfun ) );
278 }
279
287 {
289
290 public:
291 inline
296
303 inline MPolyNM1 operator() (const MPolyN & p) const
304 {
305 MPolyNM1 res;
306 X xx = (X) 1;
307 for ( int i = 0; i < (int) p.myValue.size(); ++i )
308 {
309 res += ( (MPolyNM1) p.myValue[i] ) * xx;
310 xx = xx * myOwner.myEvalPoint;
311 }
312 return res;
313 }
314 };
315
316 public:
320 inline
321 operator MPolyNM1() const
322 // operator poly<n - 1, S, typename Alloc::template rebind<S>::other>() const
323 {
324 MPolyNM1 res; // missing: determine allocator object
325 // We need to pass evaluation on to our owner
326 myOwner.evaluate( res, EvalFun2( *this ) );
327 return res;
328 }
329
335 template < typename XX >
336 inline
338 < n - 1, Ring,
340 Alloc, XX > operator() ( const XX & x ) const
341 {
344 Alloc, XX >( *this, x );
345 }
346 };
347
356 template < typename TRing, typename TAlloc, typename TX >
357 class MPolynomialEvaluator< 1, TRing, TAlloc, TX>
358 {
359 friend class MPolynomial< 1, TRing, TAlloc >;
360 public:
361 typedef TRing Ring;
362 typedef TAlloc Alloc;
363 typedef TX X;
365
366 private:
367 const MPoly1 & myPoly;
368 const X & myEvalPoint;
369
370 inline
371 MPolynomialEvaluator( const MPoly1 & poly, const X & evalpoint)
372 : myPoly( poly ), myEvalPoint( evalpoint )
373 {}
374
375 public:
379 inline
380 operator X() const
381 {
382 X res = (X) 0;
383 X xx = (X) 1;
384 for ( int i = 0; i < (int) myPoly.myValue.size(); ++i )
385 {
386 res += ( (X) (Ring) myPoly.myValue[i] ) * xx;
387 xx = xx * myEvalPoint;
388 }
389 return res;
390 }
391
395 inline
397 {
398 return (X) (*this);
399 }
400 };
401
402
411 template < int n, typename TRing, typename TAlloc, typename TX >
413 {
414 friend class MPolynomial< n, TRing, TAlloc>;
415
416 template<int nn, class TT, class HLHL, class AA, class SS>
418
419 public:
420 typedef TRing Ring;
421 typedef TAlloc Alloc;
422 typedef TX X;
424 typedef MPolynomial< n - 1, X,
425 typename std::allocator_traits<Alloc>::template rebind_alloc<X> > MPolyNM1;
426 private:
427 const MPolyN & myPoly;
428 const X & myEvalPoint;
429
430 inline
431 MPolynomialEvaluator( const MPolyN & poly, const X & evalpoint)
432 : myPoly( poly ), myEvalPoint( evalpoint )
433 {}
434
441 template < typename XX, typename Fun >
442 inline
443 void evaluate( XX & res, const Fun & evalfun ) const
444 {
445 X xx = (X) 1;
446 for ( int i = 0; i < (int) myPoly.myValue.size(); ++i )
447 {
448 res += evalfun( myPoly.myValue[i] ) * xx;
449 xx = xx * myEvalPoint;
450 }
451 }
452
453public:
457 inline
458 operator MPolyNM1() const
459 {
460 MPolyNM1 res( myPoly.getAllocator() );
461 X xx = (X) 1;
462 for ( int i = 0; i < (int) myPoly.myValue.size(); ++i )
463 {
464 res += MPolyNM1( myPoly.myValue[i] ) * xx;
465 xx = xx * myEvalPoint;
466 }
467 return res;
468 }
469
473 template <typename XX>
474 inline
477 operator() (const XX & x) const
478 {
481 ( *this, x );
482 }
483 };
484
506 template < typename TRing, typename TAlloc >
507 class MPolynomial< 0, TRing, TAlloc >
508 {
509 public:
510 typedef TRing Ring;
511 typedef TAlloc Alloc;
512
513 private:
516
517 public:
525 inline
526 MPolynomial( const Ring & v = 0, const Alloc & allocator = Alloc() )
527 : myAllocator(allocator), myValue(v)
528 {}
529
536 inline
537 MPolynomial( const Alloc & allocator )
538 : myAllocator(allocator), myValue( Ring() )
539 {}
540
544 inline bool isZero() const
545 {
546 return myValue == 0;
547 }
548
553 inline operator const Ring & () const
554 {
555 return myValue;
556 }
557
563 inline MPolynomial & operator = (const Ring & v)
564 {
565 myValue = v;
566 return *this;
567 }
568
573 inline Ring operator() () const
574 {
575 return myValue;
576 }
577
583 inline MPolynomial operator * (const Ring & v) const
584 {
585 return MPolynomial(myValue * v);
586 }
587
593 inline MPolynomial operator / (const Ring & v) const
594 {
595 return MPolynomial(myValue / v);
596 }
597
603 inline MPolynomial operator + (const Ring & v) const
604 {
605 return MPolynomial(myValue + v);
606 }
607
613 inline MPolynomial operator - (const Ring & v) const
614 {
615 return MPolynomial(myValue - v);
616 }
617
623 {
624 return MPolynomial(-myValue);
625 }
626
632 inline MPolynomial & operator *= (const Ring & v)
633 {
634 myValue *= v;
635 return *this;
636 }
637
643 inline MPolynomial & operator /= (const Ring & v)
644 {
645 myValue /= v;
646 return *this;
647 }
648
654 inline MPolynomial & operator += (const Ring & v)
655 {
656 myValue += v;
657 return *this;
658 }
659
665 inline MPolynomial & operator -= (const Ring & v)
666 {
667 myValue -= v;
668 return *this;
669 }
670
676 inline bool operator == (const Ring & v) const
677 {
678 return myValue == v;
679 }
680
686 inline bool operator != (const Ring & v) const
687 {
688 return myValue != v;
689 }
690
695 void selfDisplay( std::ostream & s, int /*N = 0*/ ) const
696 {
697 s << myValue;
698 }
699
704 inline void swap( MPolynomial & p )
705 {
706 std::swap(myValue, p.myValue);
707 }
708
713 {
714 return myAllocator;
715 }
716 };
717
718
731 template < typename T, typename TAlloc = std::allocator<T>,
732 bool usePointers = false>
734 {
735 public:
736 typedef TAlloc Alloc;
737 typedef typename std::vector<T, Alloc>::size_type Size;
738 private:
739 std::vector<T, Alloc> myVec;
740
741 public:
742 IVector( const Alloc & allocator = Alloc() )
743 : myVec(allocator)
744 {}
745
746 IVector( Size aSize, const Alloc & allocator = Alloc() )
747 : myVec( aSize, T(), allocator )
748 {}
749
750 IVector( Size aSize, const T & entry, const Alloc & allocator = Alloc() )
751 : myVec(aSize, entry, allocator)
752 {}
753
754 Size size() const
755 {
756 return myVec.size();
757 }
758
759 void resize( Size aSize, const T & entry = T() )
760 {
761 myVec.resize(aSize, entry);
762 }
763
764 const T & operator[] ( Size i ) const
765 {
766 return myVec[i];
767 }
768
770 {
771 return myVec[i];
772 }
773
774 const T & back() const
775 {
776 return myVec.back();
777 }
778
779 T & back()
780 {
781 return myVec.back();
782 }
783
784 void swap(IVector & v)
785 {
786 myVec.swap( v.myVec );
787 }
788
790 {
791 return myVec.get_allocator();
792 }
793
795 {
796 return myVec.get_allocator();
797 }
798 };
799
804 template < typename T, typename TAlloc>
805 class IVector< T, TAlloc, true >
806 {
807 public:
808 typedef TAlloc Alloc;
809 typedef typename std::allocator_traits<Alloc>::pointer TPointer;
810 typedef typename std::vector<TPointer, typename std::allocator_traits<Alloc>::template rebind_alloc<TPointer> >::size_type Size;
811
812 private:
814 std::vector<TPointer, typename std::allocator_traits<Alloc>::template rebind_alloc<TPointer> > myVec;
815
816 // Previous : typename Alloc::const_reference
817 // const value_type& should be equivalent according to : https://stackoverflow.com/questions/16360068/why-is-allocatorreference-being-phased-out
818 void create( Size begin, Size end, const typename Alloc::value_type& entry)
819 {
820 for (Size i = begin; i < end; ++i)
821 {
822 myVec[i] = myAllocator.allocate(sizeof(T));
823 std::allocator_traits<Alloc>::construct(myAllocator, myVec[i], entry);
824 }
825 }
826
827 void free(Size begin, Size end)
828 {
829 for (Size i = begin; i < end; ++i)
830 {
831 std::allocator_traits<Alloc>::destroy(myAllocator, myVec[i]);
832 myAllocator.deallocate(myVec[i], sizeof(T));
833 }
834 }
835
836 template<class A>
837 void copy_from(const std::vector<TPointer, A> & source)
838 {
839 for (Size i = 0; i < myVec.size(); ++i)
840 {
841 myVec[i] = myAllocator.allocate(sizeof(T));
842 std::allocator_traits<Alloc>::construct(myAllocator, myVec[i], *source[i]);
843 }
844 }
845
846 public:
847 IVector(const Alloc & allocator = Alloc())
848 : myAllocator(allocator), myVec(allocator)
849 {}
850
851 IVector(Size aSize, const Alloc & allocator = Alloc())
852 : myAllocator(allocator), myVec(aSize, 0, allocator)
853 {
854 create(0, aSize, T());
855 }
856
857 IVector(Size aSize, const T & entry, const Alloc & allocator = Alloc())
858 : myAllocator(allocator), myVec(aSize, 0, allocator)
859 {
860 create(0, aSize, entry);
861 }
862
863 IVector(const IVector & v)
864 : myVec(v.size())
865 {
866 copy_from(v.myVec);
867 }
868
870 {
871 free(0, (Size)myVec.size());
872 }
873
874 IVector & operator = (const IVector & v)
875 {
876 if (&v != this)
877 {
878 free(0, (Size)myVec.size());
879 myVec.resize(v.size());
880 copy_from(v.myVec);
881 }
882 return *this;
883 }
884
885 Size size() const
886 {
887 return (Size)myVec.size();
888 }
889
890 void resize(Size aSize, const T & entry = T())
891 {
892 Size oldsize = (Size)myVec.size();
893 if (oldsize > aSize)
894 free(aSize, oldsize);
895 myVec.resize(aSize);
896 if (oldsize < aSize)
897 create(oldsize, aSize, entry);
898 }
899
900 const T & operator[] (Size i) const
901 {
902 return *myVec[i];
903 }
904
906 {
907 return *myVec[i];
908 }
909
910 const T & back() const
911 {
912 return *myVec.back();
913 }
914
915 T & back()
916 {
917 return *myVec.back();
918 }
919
920 void swap(IVector & v)
921 {
922 myVec.swap(v.myVec);
923 }
924
926 {
927 return myVec.get_allocator();
928 }
929
931 {
932 return myVec.get_allocator();
933 }
934 };
935
936
938 // template class MPolynomial
968 template < int n, typename TRing, class TAlloc >
970 {
971 // ----------------------- friends ---------------------------------------
972
973 template<int NN, int nn, typename TT, typename AA>
975
981
982 template<int nn, typename TT, typename AA, typename SS>
984
985 template<int nn, typename TT, typename HLHL, typename AA, typename SS>
987
988 public:
989 typedef TRing Ring;
990 typedef TAlloc Alloc;
991 typedef MPolynomial< n - 1, Ring, Alloc > MPolyNM1;
998 typedef IVector< MPolyNM1,
999 typename std::allocator_traits<Alloc>::template rebind_alloc<MPolyNM1>, (n > 1) >
1001 typedef typename Storage::Size Size;
1002
1003 // ----------------------- Datas ----------------------------
1004private:
1014
1018 inline
1019 MPolynomial( bool, Size s, const Alloc & /*allocator*/)
1020 : myValue(s)
1021 {}
1022
1023public:
1030 inline void normalize()
1031 {
1032 Size dp1 = myValue.size();
1033 while ( dp1 )
1034 {
1035 if (myValue[dp1 - 1].isZero())
1036 --dp1;
1037 else
1038 break;
1039 }
1040 if (dp1 < myValue.size())
1041 myValue.resize(dp1);
1042 }
1043
1044 // ----------------------- Standard services ----------------------------
1045 public:
1046
1050 inline MPolynomial( const Alloc & allocator = Alloc() )
1051 : myValue( allocator )
1052 {}
1053
1057 inline MPolynomial( const Ring & v, const Alloc & allocator = Alloc() )
1058 : myValue( 1, MPolyNM1( v ), allocator )
1059 {}
1060
1065 inline MPolynomial( const MPolyNM1 & pdm1,
1066 const Alloc & /*allocator = Alloc()*/ )
1067 : myValue( 1, pdm1 )
1068 {}
1069
1074 template < typename Ring2, typename Alloc2 >
1076 const Alloc & allocator = Alloc() )
1077 : myValue( p.degree() + 1, MPolyNM1(), allocator )
1078 {
1079 for ( Size i = 0; i < myValue.size(); ++i )
1080 myValue[i] = p[i];
1081 normalize();
1082 }
1083
1088 template < typename Ring2, typename Alloc2 >
1089 inline
1090 MPolynomial & operator=
1092 {
1093 myValue.resize(p.degree() + 1);
1094 for ( Size i = 0; i < myValue.size(); ++i )
1095 myValue[i] = p[i];
1096 normalize();
1097 return *this;
1098 }
1099
1104 inline void swap( MPolynomial & p )
1105 {
1106 myValue.swap(p.myValue);
1107 }
1108
1113 {
1114 return myValue.getAllocator();
1115 }
1116
1117 // ----------------------- MPolynomial services ----------------------------
1118 public:
1119
1124 inline int degree() const
1125 {
1126 return (int)(myValue.size() - 1);
1127 }
1128
1134 inline const MPolyNM1 & leading() const
1135 {
1136 return myValue.size() ? myValue.back() : myZeroPolynomial;
1137 }
1138
1142 inline bool isZero() const
1143 {
1144 return myValue.size() == 0;
1145 }
1146
1153 {
1154 if (i >= myValue.size())
1155 myValue.resize(i + 1);
1156 return myValue[i];
1157 }
1158
1162 inline const MPolyNM1 & operator[] (Size i) const
1163 {
1164 return i < myValue.size() ? myValue[i] : myZeroPolynomial;
1165 }
1166
1173 inline
1178
1186 template <typename Ring2>
1187 inline
1192
1193 // ----------------------- Operators --------------------------------------
1194 public:
1195
1201 inline MPolynomial operator * (const Ring & v) const
1202 {
1203 MPolynomial r(*this);
1204 for ( Size i = 0; i < myValue.size(); ++i )
1205 r[i] *= v;
1206 return r;
1207 }
1208
1214 inline MPolynomial operator / (const Ring & v) const
1215 {
1216 MPolynomial r(*this);
1217 for ( Size i = 0; i < myValue.size(); ++i )
1218 r[i] /= v;
1219 return r;
1220 }
1221
1228 {
1229 return *this = *this * p;
1230 }
1231
1237 inline MPolynomial & operator *= (const Ring & v)
1238 {
1239 MPolynomial r( *this );
1240 for ( Size i = 0; i < myValue.size(); ++i )
1241 myValue[i] *= v;
1242 return *this;
1243 }
1244
1250 inline MPolynomial & operator /= (const Ring & v)
1251 {
1252 for ( Size i = 0; i < myValue.size(); ++i )
1253 myValue[i] /= v;
1254 return *this;
1255 }
1256
1263 friend
1264 inline MPolynomial operator * (const Ring & v, const MPolynomial & p)
1265 {
1266 MPolynomial r(p);
1267 for ( Size i = 0; i < p.myValue.size(); ++i )
1268 r[i] *= v;
1269 return r;
1270 }
1271
1276 {
1277 MPolynomial r( true, myValue.size(), getAllocator() );
1278 for ( Size i = 0; i < myValue.size(); ++i )
1279 r[i] = -myValue[i];
1280 return r;
1281 }
1282
1288 inline MPolynomial operator + (const MPolynomial & q) const
1289 {
1290 MPolynomial r(*this);
1291 if (r.myValue.size() < q.myValue.size())
1292 r.myValue.resize(q.myValue.size());
1293 for ( Size i = 0; i < q.myValue.size(); ++i )
1294 r[i] += q[i];
1295 r.normalize();
1296 return r;
1297 }
1298
1304 inline MPolynomial operator - (const MPolynomial & q) const
1305 {
1306 MPolynomial r(*this);
1307 if (r.myValue.size() < q.myValue.size())
1308 r.myValue.resize(q.myValue.size());
1309 for ( Size i = 0; i < q.myValue.size(); ++i )
1310 r[i] -= q[i];
1311 r.normalize();
1312 return r;
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 {
1337 if (myValue.size() < q.myValue.size())
1338 myValue.resize(q.myValue.size());
1339 for ( Size i = 0; i < q.myValue.size(); ++i )
1340 myValue[i] -= q[i];
1341 normalize();
1342 return *this;
1343 }
1344
1351 inline MPolynomial operator + (const MPolyNM1 & q) const
1352 {
1353 MPolynomial r(*this);
1354 if (r.myValue.size() < 1)
1355 r.myValue.resize(1);
1356 r[0] += q;
1357 r.normalize();
1358 return r;
1359 }
1360
1367 inline MPolynomial operator - (const MPolyNM1 & q) const
1368 {
1369 MPolynomial r(*this);
1370 if (r.myValue.size() < 1)
1371 r.myValue.resize(1);
1372 r[0] -= q;
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
1397 inline MPolynomial operator - (const Ring & v) const
1398 {
1399 MPolynomial r(*this);
1400 if (r.myValue.size() < 1)
1401 r.myValue.resize(1);
1402 r[0] -= v;
1403 r.normalize();
1404 return r;
1405 }
1406
1414 {
1415 if (myValue.size() < 1)
1416 myValue.resize(1);
1417 myValue[0] += q;
1418 normalize();
1419 return *this;
1420 }
1421
1429 {
1430 if (myValue.size() < 1)
1431 myValue.resize(1);
1432 myValue[0] -= q;
1433 normalize();
1434 return *this;
1435 }
1436
1442 inline MPolynomial & operator += (const Ring & v)
1443 {
1444 if (myValue.size() < 1)
1445 myValue.resize(1);
1446 myValue[0] += v;
1447 normalize();
1448 return *this;
1449 }
1450
1456 inline MPolynomial & operator -= (const Ring & v)
1457 {
1458 if (myValue.size() < 1)
1459 myValue.resize(1);
1460 myValue[0] -= v;
1461 normalize();
1462 return *this;
1463 }
1464
1472 inline MPolynomial operator * (const MPolynomial & p) const
1473 {
1474 int d = p.degree() + degree();
1475 if (d < -1)
1476 d = -1;
1477 MPolynomial r( true, d + 1, getAllocator() );
1478 if (!isZero() && !p.isZero())
1479 for ( Size i = 0; i < r.myValue.size(); ++i )
1480 for ( Size j = 0; j < myValue.size(); ++j )
1481 if (i < j + p.myValue.size())
1482 r[i] += myValue[j] * p[i - j];
1483 r.normalize();
1484 return r;
1485 }
1486
1487 // ----------------------- Comparison operators ---------------------------
1488 public:
1489
1495 inline bool operator == (const MPolynomial & q) const
1496 {
1497 if (myValue.size() != q.myValue.size())
1498 return false;
1499 for (Size i = 0; i < myValue.size(); ++i)
1500 if (myValue[i] != q[i])
1501 return false;
1502 return true;
1503 }
1504
1510 inline bool operator != (const MPolynomial & q) const
1511 {
1512 return !(*this == q);
1513 }
1514
1520 inline bool operator == (const Ring & v) const
1521 {
1522 if ((v == 0) && (myValue.size() == 0))
1523 return true;
1524 if (myValue.size() != 1)
1525 return false;
1526 return myValue[0] == v;
1527 }
1528
1534 inline bool operator != (const Ring & v) const
1535 {
1536 return !(*this == v);
1537 }
1538
1539
1540
1541 // ----------------------- Interface --------------------------------------
1542 public:
1543
1551 void selfDisplay ( std::ostream & s, int N = n ) const
1552 {
1553 if (isZero())
1554 s << (Ring) 0;
1555 else
1556 {
1557 Size nonzero = 0;
1558 for (Size i = 0; i < myValue.size(); ++i)
1559 if (!myValue[i].isZero())
1560 ++nonzero;
1561 if (nonzero > 1) s << "(";
1562 bool first = true;
1563 for (Size i = 0; i < myValue.size(); ++i)
1564 if (!myValue[i].isZero())
1565 {
1566 if (first) first = false;
1567 else s << " + ";
1568 myValue[i].selfDisplay(s, N);
1569 if (i > 0)
1570 {
1571 s << " ";
1572 s << "X_" << N - n;
1573 if (i > 1) s << "^" << i;
1574 }
1575 }
1576 if (nonzero > 1)
1577 s << ")";
1578 }
1579 }
1580
1585 bool isValid() const;
1586
1587
1588 // ------------------------- Datas --------------------------------------
1589 private:
1590 // ------------------------- Hidden services ----------------------------
1591 protected:
1592
1593
1594 }; // end of class MPolynomial
1595
1596
1603 template <int N, typename TRing, class TAlloc>
1604 std::ostream&
1605 operator<< ( std::ostream & out,
1606 const MPolynomial<N, TRing, TAlloc> & object );
1607
1608
1609 // ------------------------- monomial services ----------------------------
1610
1617 template <int n, typename Ring, typename Alloc>
1619 {
1620 public:
1622
1629 get( unsigned int k, unsigned int e )
1630 {
1632 if ( k == 0 )
1633 p[e] = Xe_kComputer<n-1,Ring,Alloc>().get( k-1, e );
1634 else
1635 p[0] = Xe_kComputer<n-1,Ring,Alloc>().get( k-1, e );
1636 p.normalize();
1637 //std::cerr << "Xe_k(" << k << "," << e << ")=" << p << std::endl;
1638 return p;
1639 }
1640
1641 };
1642
1643 template <typename Ring, typename Alloc>
1644 class Xe_kComputer<0,Ring,Alloc>
1645 {
1646 public:
1648
1650 get( unsigned int /*k*/, unsigned int /*e*/ )
1651 {
1653 return p;
1654 }
1655 };
1656
1666 template <int n, typename Ring, typename Alloc>
1667 inline
1668 MPolynomial<n, Ring, Alloc>
1669 Xe_k( unsigned int k, unsigned int e )
1670 {
1671 return Xe_kComputer<n, Ring, Alloc>().get( k, e );
1672 }
1673
1682 template <int n, typename Ring>
1683 inline
1684 MPolynomial<n, Ring, std::allocator<Ring> >
1685 Xe_k( unsigned int k, unsigned int e )
1686 {
1687 return Xe_kComputer<n, Ring, std::allocator<Ring> >().get( k, e );
1688 }
1689
1690
1698 template <typename Ring, typename Alloc>
1699 inline
1700 MPolynomial<1, Ring, Alloc>
1701 mmonomial(unsigned int e)
1702 {
1704 p[e] = 1;
1705 return p;
1706 }
1707
1716 template <typename Ring, typename Alloc>
1717 inline
1718 MPolynomial<2, Ring, Alloc>
1719 mmonomial(unsigned int e, unsigned int f)
1720 {
1722 p[e][f] = 1;
1723 return p;
1724 }
1725
1735 template <typename Ring, typename Alloc>
1736 inline MPolynomial<3, Ring, Alloc>
1737 mmonomial(unsigned int e, unsigned int f, unsigned int g)
1738 {
1740 p[e][f][g] = 1;
1741 return p;
1742 }
1743
1754 template <typename Ring, typename Alloc>
1755 inline
1756 MPolynomial<4, Ring, Alloc>
1757 mmonomial(unsigned int e, unsigned int f, unsigned int g, unsigned int h)
1758 {
1760 p[e][f][g][h] = 1;
1761 return p;
1762 }
1763
1770 template <typename Ring>
1771 inline
1772 MPolynomial<1, Ring, std::allocator<Ring> >
1773 mmonomial(unsigned int e)
1774 {
1776 p[e] = 1;
1777 return p;
1778 }
1779
1787 template <typename Ring>
1788 inline
1789 MPolynomial<2, Ring, std::allocator<Ring> >
1790 mmonomial(unsigned int e, unsigned int f)
1791 {
1793 p[e][f] = 1;
1794 return p;
1795 }
1796
1805 template <typename Ring>
1806 inline
1807 MPolynomial<3, Ring, std::allocator<Ring> >
1808 mmonomial(unsigned int e, unsigned int f, unsigned int g)
1809 {
1811 p[e][f][g] = 1;
1812 return p;
1813 }
1814
1824 template <typename Ring>
1825 inline
1826 MPolynomial<4, Ring, std::allocator<Ring> >
1828 (unsigned int e, unsigned int f, unsigned int g, unsigned int h)
1829 {
1831 p[e][f][g][h] = 1;
1832 return p;
1833 }
1834
1835
1856 template <int n, typename Ring, class Alloc>
1858 {
1859 public:
1862
1870 static inline
1871 void computeDerivative( const MPolyN & src, MPolyN & dest)
1872 {
1873 dest.myValue.resize(src.degree() >= 0 ? src.degree() : 0);
1874 for ( int i = 1; i <= src.degree(); ++i )
1875 dest[i - 1] = src[i] * (Ring)i;
1876 dest.normalize();
1877 }
1878 };
1879
1901 template <int N, int n, typename Ring, typename Alloc>
1903 {
1904 public:
1907
1915 static inline
1916 void computeDerivative(const MPolyN & src, MPolyN & dest)
1917 {
1918 dest.myValue.resize(src.degree() + 1);
1919 for ( int i = 0; i <= src.degree(); ++i )
1921 ::computeDerivative( src[i], dest[i] );
1922 dest.normalize();
1923 }
1924 };
1925
1929 template<typename Ring, class Alloc>
1931 {
1932 public:
1933 class ERROR_N_must_be_strictly_less_than_n_in_derivative_template;
1935
1936 static inline
1938 {
1939 ERROR_N_must_be_strictly_less_than_n_in_derivative_template();
1940 }
1941 };
1942
1946 template<int N, typename Ring, class Alloc>
1948 {
1949 public:
1950 class ERROR_N_must_be_strictly_less_than_n_in_derivative_template;
1952
1953 static inline
1955 {
1956 ERROR_N_must_be_strictly_less_than_n_in_derivative_template();
1957 }
1958 };
1959
1978 template <int N, int n, typename Ring, typename Alloc>
1979 inline
1980 MPolynomial<n, Ring, Alloc>
1987
2001 template<int N, int n, typename Ring>
2002 inline
2003 MPolynomial<n, Ring, std::allocator<Ring> >
2005 (const MPolynomial<n,Ring,std::allocator<Ring> > & p)
2006 {
2007 MPolynomial<n, Ring, std::allocator<Ring> > res( p.getAllocator() );
2009 ::computeDerivative( p, res );
2010 return res;
2011 }
2012
2016 template<typename Ring, typename Alloc>
2017 void
2022 {
2023 if (f.degree() < g.degree())
2024 {
2025 // Ignore the trivial case
2027 r = f;
2028 return;
2029 }
2030 q = MPolynomial<1, Ring, Alloc>( true, f.degree() - g.degree() + 1,
2031 f.getAllocator() );
2032 r = f;
2033 for (int i = q.degree(); i >= 0; --i)
2034 {
2035 q[i] = r[i + g.degree()] / g.leading();
2036 for (int j = g.degree(); j >= 0; --j)
2037 r[i + j] -= q[i] * g[j];
2038 }
2039 r.normalize();
2040 // Note that the degree of q is already correct.
2041 }
2042
2046 template <typename Ring>
2047 void
2048 euclidDiv( const MPolynomial<1, Ring, std::allocator<Ring> > & f,
2049 const MPolynomial<1, Ring, std::allocator<Ring> > & g,
2050 MPolynomial<1, Ring, std::allocator<Ring> > & q,
2051 MPolynomial<1, Ring, std::allocator<Ring> > & r )
2052 {
2054 }
2055
2060 template<typename Ring, typename Alloc>
2061 MPolynomial<1, Ring, Alloc>
2063 const MPolynomial<1, Ring, Alloc> & g )
2064 {
2065 if (f.isZero())
2066 {
2067 if (g.isZero()) return f; // both are zero
2068 else return g / g.leading(); // make g monic
2069 }
2071 d1(f / f.leading()),
2072 d2(g / g.leading()),
2073 q(f.getAllocator()),
2074 r(f.getAllocator());
2075 while (!d2.isZero())
2076 {
2077 euclidDiv(d1, d2, q, r);
2078 d1.swap(d2);
2079 d2 = r;
2080 d2 /= r.leading(); // make r monic
2081 }
2082 return d1;
2083 }
2084
2089 template<typename Ring>
2090 MPolynomial<1, Ring, std::allocator<Ring> >
2091 gcd( const MPolynomial<1, Ring, std::allocator<Ring> > & f,
2092 const MPolynomial<1, Ring, std::allocator<Ring> > & g )
2093 {
2094 return gcd<Ring, std::allocator<Ring> >(f, g);
2095 }
2096
2097} // namespace DGtal
2098
2099
2101// Includes inline functions.
2102#include "DGtal/math/MPolynomial.ih"
2103
2104#pragma GCC diagnostic pop
2105#pragma clang diagnostic pop
2106
2107
2108// //
2110
2111#endif // !defined MPolynomial_h
2112
2113#undef MPolynomial_RECURSES
2114#endif // else defined(MPolynomial_RECURSES)
IVector(Size aSize, const Alloc &allocator=Alloc())
IVector(Size aSize, const T &entry, const Alloc &allocator=Alloc())
IVector(const Alloc &allocator=Alloc())
void copy_from(const std::vector< TPointer, A > &source)
void resize(Size aSize, const T &entry=T())
std::vector< TPointer, typenamestd::allocator_traits< Alloc >::templaterebind_alloc< TPointer > >::size_type Size
void create(Size begin, Size end, const typename Alloc::value_type &entry)
std::allocator_traits< Alloc >::pointer TPointer
void free(Size begin, Size end)
std::vector< TPointer, typename std::allocator_traits< Alloc >::template rebind_alloc< TPointer > > myVec
Size size() const
std::vector< T, Alloc >::size_type Size
IVector(const Alloc &allocator=Alloc())
IVector(Size aSize, const Alloc &allocator=Alloc())
const T & back() const
void swap(IVector &v)
Alloc get_allocator() const
const T & operator[](Size i) const
IVector(Size aSize, const T &entry, const Alloc &allocator=Alloc())
Alloc getAllocator() const
void resize(Size aSize, const T &entry=T())
static void computeDerivative(const MPoly0 &, MPoly0 &)
static void computeDerivative(const MPolyN &src, MPolyN &dest)
MPolynomial< n, Ring, Alloc > MPolyN
Type for polynomial with n variable in the ring Ring.
static void computeDerivative(const MPoly0 &, MPoly0 &)
MPolynomial< n, Ring, Alloc > MPolyN
Type for polynomial with n variable in the ring Ring.
static void computeDerivative(const MPolyN &src, MPolyN &dest)
const MPolynomialEvaluatorImpl< n, Ring, Owner, Alloc, X > & myOwner
MPolyNM1 operator()(const MPolyN &p) const
EvalFun2(const MPolynomialEvaluatorImpl< n, Ring, Owner, Alloc, X > &owner)
const MPolynomialEvaluatorImpl< n, Ring, Owner, Alloc, X > & myOwner
XX operator()(const MPolynomial< n, Ring, Alloc > &p) const
EvalFun(const MPolynomialEvaluatorImpl< n, Ring, Owner, Alloc, X > &owner, const Fun &evalfun)
const MPolynomialEvaluatorImpl< 1, Ring, Owner, Alloc, X > & myOwner
EvalFun(const MPolynomialEvaluatorImpl< 1, Ring, Owner, Alloc, X > &owner)
MPolynomialEvaluatorImpl(const Owner &owner, const X &evalpoint)
const X & myEvalPoint
The evaluation point on this level.
void evaluate(XX &res, const Fun &evalfun) const
MPolynomialEvaluatorImpl< n - 1, Ring, MPolynomialEvaluatorImpl< n, Ring, Owner, Alloc, X >, Alloc, XX > operator()(const XX &x) const
MPolynomial< n, Ring, Alloc > MPolyN
Type for the multivariate polynomial.
MPolynomial< n - 1, X, typename std::allocator_traits< Alloc >::template rebind_alloc< X > > MPolyNM1
MPolynomialEvaluatorImpl(const Owner &owner, const X &evalpoint)
const MPoly1 & myPoly
The polynomial in question.
MPolynomialEvaluator(const MPoly1 &poly, const X &evalpoint)
friend class MPolynomialEvaluatorImpl
const X & myEvalPoint
the evaluation point
MPolynomial< n - 1, X, typename std::allocator_traits< Alloc >::template rebind_alloc< X > > MPolyNM1
MPolynomialEvaluator(const MPolyN &poly, const X &evalpoint)
MPolynomialEvaluatorImpl< n - 1, Ring, MPolynomialEvaluator< n, Ring, Alloc, X >, Alloc, XX > operator()(const XX &x) const
void evaluate(XX &res, const Fun &evalfun) const
const MPolyN & myPoly
The polynomial in question.
MPolynomial(const Alloc &allocator)
MPolynomial(const Ring &v=0, const Alloc &allocator=Alloc())
void selfDisplay(std::ostream &s, int) const
Aim: Represents a multivariate polynomial, i.e. an element of , where K is some ring or field.
MPolynomial(const Ring &v, const Alloc &allocator=Alloc())
bool operator==(const MPolynomial &q) const
MPolynomialEvaluator< n, Ring, Alloc, Ring > operator()(const Ring &x) const
friend void euclidDiv(const MPolynomial< 1, TRing, TAlloc > &, const MPolynomial< 1, TRing, TAlloc > &, MPolynomial< 1, TRing, TAlloc > &, MPolynomial< 1, TRing, TAlloc > &)
void selfDisplay(std::ostream &s, int N=n) const
MPolynomial(const Alloc &allocator=Alloc())
MPolynomial & operator-=(const MPolynomial &q)
MPolynomial operator/(const Ring &v) const
MPolynomial & operator/=(const Ring &v)
MPolynomial< n - 1, Ring, Alloc > MPolyNM1
MPolynomial operator-() const
MPolynomial & operator+=(const MPolynomial &q)
bool isValid() const
MPolynomial(const MPolyNM1 &pdm1, const Alloc &)
IVector< MPolyNM1, typename std::allocator_traits< Alloc >::template rebind_alloc< MPolyNM1 >,(n > 1) > Storage
MPolyNM1 & operator[](Size i)
MPolynomial & operator*=(const MPolynomial &p)
Alloc getAllocator() const
MPolynomial operator*(const Ring &v) const
bool operator!=(const MPolynomial &q) const
void swap(MPolynomial &p)
const MPolyNM1 & leading() const
bool isZero() const
MPolynomial(bool, Size s, const Alloc &)
MPolynomial(const MPolynomial< n, Ring2, Alloc2 > &p, const Alloc &allocator=Alloc())
MPolynomial & operator=(const MPolynomial< n, Ring2, Alloc2 > &p)
MPolynomial operator+(const MPolynomial &q) const
MPolynomial< 0, Ring, Alloc > get(unsigned int, unsigned int)
MPolynomial< n, Ring, Alloc > get(unsigned int k, unsigned int e)
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 > mmonomial(unsigned int e)
MPolynomial< n, Ring, Alloc > Xe_k(unsigned int k, unsigned int e)
std::ostream & operator<<(std::ostream &out, const ClosedIntegerHalfPlane< TSpace > &object)
MPolynomial< 1, Ring, Alloc > gcd(const MPolynomial< 1, Ring, Alloc > &f, const MPolynomial< 1, Ring, Alloc > &g)
MPolynomial< n, Ring, Alloc > derivative(const MPolynomial< n, Ring, Alloc > &p)
RealPointT::Coordinate Ring