d288424f3cf170e4ffa23bbd829dcfb9c2ca9753
[scilab.git] / scilab / modules / ast / includes / analysis / gvn / MultivariatePolynomial.hxx
1 /*
2  *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3  *  Copyright (C) 2015 - Scilab Enterprises - Calixte DENIZET
4  *
5  *  This file must be used under the terms of the CeCILL.
6  *  This source file is licensed as described in the file COPYING, which
7  *  you should have received as part of this distribution.  The terms
8  *  are also available at
9  *  http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
10  *
11  */
12
13 #ifndef __MULTIVARIATE_POLYNOMIAL_HXX__
14 #define __MULTIVARIATE_POLYNOMIAL_HXX__
15
16 #include <cmath>
17 #include <functional>
18 #include <iostream>
19 #include <map>
20 #include <set>
21 #include <sstream>
22 #include <string>
23 #include <unordered_set>
24 #include <vector>
25
26 #include "core_math.h"
27 #include "tools.hxx"
28 #include "MultivariateMonomial.hxx"
29
30 namespace analysis
31 {
32 /**
33  * \struct MultivariatePolynomial
34  * \brief Represents a multivariate polynomial
35  */
36 struct MultivariatePolynomial
37 {
38     typedef std::unordered_set<MultivariateMonomial, MultivariateMonomial::Hash, MultivariateMonomial::Eq> Polynomial;
39     double constant;
40     Polynomial polynomial;
41
42     /**
43      * \brief constructor
44      * \param var to init polynomial
45      */
46     MultivariatePolynomial(const unsigned long long var) : constant(0)
47     {
48         polynomial.emplace(var);
49     }
50
51     /**
52      * \brief constructor
53      * \param _constant to init polynomial
54      */
55     MultivariatePolynomial(const double _constant = 0) : constant(_constant) { }
56
57     /**
58      * \brief constructor
59      * \param _size the size of the polynomial (used to reserve the unordered_set)
60      * \param _constant to init polynomial
61      */
62     MultivariatePolynomial(const unsigned int _size, const double _constant) : constant(_constant), polynomial(_size) { }
63
64     /**
65      * \brief copy constructor
66      */
67     MultivariatePolynomial(const MultivariatePolynomial & mp) : constant(mp.constant), polynomial(mp.polynomial) { }
68
69     /**
70      * \brief Get an invalid polynomial (i.e. constant == NaN)
71      * \return invalid polynomial
72      */
73     inline static MultivariatePolynomial getInvalid()
74     {
75         return MultivariatePolynomial(tools::NaN());
76     }
77
78     /**
79      * \brief Check if it is valid
80      * \return true if it is valid
81      */
82     inline bool isValid() const
83     {
84         return !tools::isNaN(constant);
85     }
86
87     /**
88      * \brief Check if it is invalid
89      * \return true if it is invalid
90      */
91     inline bool isInvalid() const
92     {
93         return tools::isNaN(constant);
94     }
95
96     /**
97      * \brief Invalidate the polynomial
98      */
99     inline void invalid()
100     {
101         constant = tools::NaN();
102         polynomial.clear();
103     }
104
105     /**
106      * \brief Check if the variables of the polynomial have an id lower or equal to max
107      * \param max an id
108      * \return true if all the variables have an id leq to max
109      */
110     inline bool checkVariable(const unsigned long long max) const
111     {
112         for (const auto & m : polynomial)
113         {
114             if (!m.checkVariable(max))
115             {
116                 return false;
117             }
118         }
119         return true;
120     }
121
122     /**
123      * \brief Add a monomial in the polynomial
124      * \param m the monomial to add
125      * \param coeff the multiplcative coefficient to applicate to the monomial
126      * \return *this
127      */
128     inline MultivariatePolynomial & add(const MultivariateMonomial & m, const double coeff = 1)
129     {
130         const double c = m.coeff * coeff;
131         if (c)
132         {
133             Polynomial::iterator i = polynomial.find(m);
134             if (i == polynomial.end())
135             {
136                 Polynomial::iterator j = polynomial.insert(m).first;
137                 j->coeff = c;
138             }
139             else
140             {
141                 if (i->coeff == -c)
142                 {
143                     polynomial.erase(i);
144                 }
145                 else
146                 {
147                     i->coeff += c;
148                 }
149             }
150         }
151         return *this;
152     }
153
154     /**
155      * \brief Subtract a monomial in the polynomial
156      * \param m the monomial to add
157      */
158     inline void sub(const MultivariateMonomial & m)
159     {
160         Polynomial::iterator i = polynomial.find(m);
161         if (i == polynomial.end())
162         {
163             if (m.coeff)
164             {
165                 polynomial.insert(m).first->coeff = -m.coeff;
166             }
167         }
168         else
169         {
170             if (i->coeff == m.coeff)
171             {
172                 polynomial.erase(i);
173             }
174             else
175             {
176                 i->coeff -= m.coeff;
177             }
178         }
179     }
180
181     /**
182      * \brief Overload of the + operator
183      */
184     inline MultivariatePolynomial operator+(const MultivariateMonomial & R) const
185     {
186         if (isValid())
187         {
188             MultivariatePolynomial res(*this);
189             res.add(R);
190             return res;
191         }
192         return getInvalid();
193     }
194
195     /**
196      * \brief Overload of the += operator
197      */
198     inline MultivariatePolynomial & operator+=(const MultivariateMonomial & R)
199     {
200         if (isValid())
201         {
202             add(R);
203         }
204         return *this;
205     }
206
207     /**
208      * \brief Overload of the + operator
209      */
210     inline MultivariatePolynomial operator+(const double R) const
211     {
212         if (isValid())
213         {
214             MultivariatePolynomial res(*this);
215             res.constant += R;
216             return res;
217         }
218         return *this;
219     }
220
221     /**
222      * \brief Overload of the += operator
223      */
224     inline MultivariatePolynomial & operator+=(const double R)
225     {
226         if (isValid())
227         {
228             constant += R;
229         }
230         return *this;
231     }
232
233     /**
234      * \brief Overload of the - operator
235      */
236     inline MultivariatePolynomial operator-(const MultivariateMonomial & R) const
237     {
238         if (isValid())
239         {
240             MultivariatePolynomial res(*this);
241             res.sub(R);
242             return res;
243         }
244         return *this;
245     }
246
247     /**
248      * \brief Overload of the -= operator
249      */
250     inline MultivariatePolynomial & operator-=(const MultivariateMonomial & R)
251     {
252         if (isValid())
253         {
254             sub(R);
255         }
256         return *this;
257     }
258
259     /**
260      * \brief Overload of the - operator
261      */
262     inline MultivariatePolynomial operator-(const double R) const
263     {
264         if (isValid())
265         {
266             MultivariatePolynomial res(*this);
267             res.constant -= R;
268             return res;
269         }
270         return *this;
271     }
272
273     /**
274      * \brief Overload of the - operator (unary)
275      */
276     inline MultivariatePolynomial operator-() const
277     {
278         if (isValid())
279         {
280             MultivariatePolynomial res(*this);
281             res.constant = -res.constant;
282             for (auto & m : res.polynomial)
283             {
284                 m.coeff = -m.coeff;
285             }
286             return res;
287         }
288         return *this;
289     }
290
291     /**
292      * \brief Overload of the -= operator
293      */
294     inline MultivariatePolynomial & operator-=(const double R)
295     {
296         if (isValid())
297         {
298             constant -= R;
299         }
300         return *this;
301     }
302
303     /**
304      * \brief Overload of the + operator
305      */
306     inline MultivariatePolynomial operator+(const MultivariatePolynomial & R) const
307     {
308         if (isValid() && R.isValid())
309         {
310             MultivariatePolynomial res(*this);
311             res.constant += R.constant;
312             for (const auto & m : R.polynomial)
313             {
314                 res.add(m);
315             }
316             return res;
317         }
318         return getInvalid();
319     }
320
321     /**
322      * \brief Overload of the += operator
323      */
324     inline MultivariatePolynomial & operator+=(const MultivariatePolynomial & R)
325     {
326         if (isValid() && R.isValid())
327         {
328             constant += R.constant;
329             for (const auto & m : R.polynomial)
330             {
331                 add(m);
332             }
333         }
334         else
335         {
336             invalid();
337         }
338
339         return *this;
340     }
341
342     /**
343      * \brief Overload of the - operator
344      */
345     inline MultivariatePolynomial operator-(const MultivariatePolynomial & R) const
346     {
347         if (isValid() && R.isValid())
348         {
349             MultivariatePolynomial res(*this);
350             res.constant -= R.constant;
351             for (const auto & m : R.polynomial)
352             {
353                 res.sub(m);
354             }
355             return res;
356         }
357         return getInvalid();
358     }
359
360     /**
361      * \brief Overload of the -= operator
362      */
363     inline MultivariatePolynomial & operator-=(const MultivariatePolynomial & R)
364     {
365         if (isValid() && R.isValid())
366         {
367             constant -= R.constant;
368             for (const auto & m : R.polynomial)
369             {
370                 sub(m);
371             }
372         }
373         else
374         {
375             invalid();
376         }
377         return *this;
378     }
379
380     /**
381      * \brief Overload of the * operator
382      */
383     inline MultivariatePolynomial operator*(const MultivariatePolynomial & R) const
384     {
385         if (isValid() && R.isValid())
386         {
387             if (R.isConstant())
388             {
389                 return *this * R.constant;
390             }
391             else if (isConstant())
392             {
393                 return R * constant;
394             }
395             else
396             {
397                 MultivariatePolynomial res((polynomial.size() + 1) * (R.polynomial.size() + 1) - 1, constant * R.constant);
398                 for (const auto & mR : R.polynomial)
399                 {
400                     res.add(mR, constant);
401                 }
402                 for (const auto & mL : polynomial)
403                 {
404                     res.add(mL, R.constant);
405                     for (const auto & mR : R.polynomial)
406                     {
407                         res.add(mL * mR);
408                     }
409                 }
410                 return res;
411             }
412         }
413         return getInvalid();
414     }
415
416     /**
417      * \brief Overload of the / operator
418      */
419     inline MultivariatePolynomial operator/(const MultivariatePolynomial & R) const
420     {
421         if (isValid() && R.isValid())
422         {
423             if (R.isConstant())
424             {
425                 if (R.constant != 1)
426                 {
427                     return *this / R.constant;
428                 }
429                 else
430                 {
431                     return *this;
432                 }
433             }
434         }
435         return getInvalid();
436     }
437
438     /**
439      * \brief Overload of the /= operator
440      */
441     inline MultivariatePolynomial & operator/=(const MultivariatePolynomial & R)
442     {
443         if (isValid() && R.isValid())
444         {
445             if (R.polynomial.empty())
446             {
447                 constant /= R.constant;
448                 for (auto & m : polynomial)
449                 {
450                     m.coeff /= R.constant;
451                 }
452             }
453             else
454             {
455                 MultivariatePolynomial res(*this / R);
456                 polynomial = res.polynomial;
457                 constant = res.constant;
458             }
459         }
460         else
461         {
462             invalid();
463         }
464         return *this;
465     }
466
467     /**
468      * \brief Overload of the *= operator
469      */
470     inline MultivariatePolynomial & operator*=(const MultivariatePolynomial & R)
471     {
472         if (isValid() && R.isValid())
473         {
474             if (R.polynomial.empty())
475             {
476                 constant *= R.constant;
477                 for (auto & m : polynomial)
478                 {
479                     m.coeff *= R.constant;
480                 }
481             }
482             else
483             {
484                 MultivariatePolynomial res(*this * R);
485                 polynomial = res.polynomial;
486                 constant = res.constant;
487             }
488         }
489         else
490         {
491             invalid();
492         }
493         return *this;
494     }
495
496     /**
497      * \brief Overload of the * operator
498      */
499     inline MultivariatePolynomial operator*(const MultivariateMonomial & R) const
500     {
501         if (isValid())
502         {
503             MultivariatePolynomial res(polynomial.size() + 1, 0.);
504             res.add(constant * R);
505             for (const auto & mL : polynomial)
506             {
507                 res.add(mL * R);
508             }
509             return res;
510         }
511         else
512         {
513             return getInvalid();
514         }
515     }
516
517     /**
518      * \brief Overload of the *= operator
519      */
520     inline MultivariatePolynomial & operator*=(const MultivariateMonomial & R)
521     {
522         if (isValid())
523         {
524             MultivariatePolynomial res = *this * R;
525             polynomial = res.polynomial;
526             constant = res.constant;
527         }
528         return *this;
529     }
530
531     /**
532      * \brief Overload of the * operator
533      */
534     inline MultivariatePolynomial operator*(const double R) const
535     {
536         if (isValid())
537         {
538             if (R)
539             {
540                 if (R == 1)
541                 {
542                     return *this;
543                 }
544                 else
545                 {
546                     MultivariatePolynomial res(*this);
547                     res.constant *= R;
548                     for (auto & m : res.polynomial)
549                     {
550                         m.coeff *= R;
551                     }
552                     return res;
553                 }
554             }
555             else
556             {
557                 return MultivariatePolynomial(0.);
558             }
559         }
560         return getInvalid();
561     }
562
563     /**
564      * \brief Overload of the *= operator
565      */
566     inline MultivariatePolynomial & operator*=(const double R)
567     {
568         if (isValid())
569         {
570             if (R == 0)
571             {
572                 constant = 0.;
573                 polynomial.clear();
574             }
575             else if (R != 1)
576             {
577                 constant *= R;
578                 for (auto & m : polynomial)
579                 {
580                     m.coeff *= R;
581                 }
582             }
583         }
584         return *this;
585     }
586
587     /**
588      * \brief Overload of the / operator
589      */
590     inline MultivariatePolynomial operator/(const double R) const
591     {
592         if (isValid())
593         {
594             if (R != 1)
595             {
596                 MultivariatePolynomial res(*this);
597                 res.constant /= R;
598                 for (auto & m : res.polynomial)
599                 {
600                     m.coeff /= R;
601                 }
602                 return res;
603             }
604         }
605         return *this;
606     }
607
608     /**
609      * \brief Overload of the /= operator
610      */
611     inline MultivariatePolynomial & operator/=(const double R)
612     {
613         if (isValid())
614         {
615             if (R != 1)
616             {
617                 constant /= R;
618                 for (auto & m : polynomial)
619                 {
620                     m.coeff /= R;
621                 }
622             }
623         }
624         return *this;
625     }
626
627     /**
628      * \brief Overload of the ^ operator (exponentiation)
629      */
630     inline MultivariatePolynomial operator^(unsigned int R) const
631     {
632         if (isValid())
633         {
634             if (R == 0)
635             {
636                 return MultivariatePolynomial(1.);
637             }
638             else if (R == 1)
639             {
640                 return *this;
641             }
642             else
643             {
644                 if (constant == 0)
645                 {
646                     if (polynomial.empty())
647                     {
648                         return MultivariatePolynomial(0.);
649                     }
650                     else if (polynomial.size() == 1)
651                     {
652                         const MultivariateMonomial & m = *polynomial.begin();
653                         MultivariatePolynomial res;
654                         res.polynomial.emplace(m ^ R);
655
656                         return res;
657                     }
658                 }
659
660                 if (polynomial.empty())
661                 {
662                     return MultivariatePolynomial(std::pow(constant, R));
663                 }
664
665                 MultivariatePolynomial p = *this;
666                 MultivariatePolynomial y = (R & 1) ? *this : MultivariatePolynomial(1.);
667
668                 while (R >>= 1)
669                 {
670                     p *= p;
671                     if (R & 1)
672                     {
673                         y *= p;
674                     }
675                 }
676
677                 return y;
678             }
679         }
680         return getInvalid();
681     }
682
683     /**
684      * \brief Overload of the ^ operator (exponentiation)
685      */
686     inline MultivariatePolynomial operator^(const MultivariatePolynomial & R) const
687     {
688         if (isValid() && R.isValid())
689         {
690             if (R.isConstant() && R.constant == (unsigned int)R.constant)
691             {
692                 return (*this) ^ ((unsigned int)R.constant);
693             }
694         }
695         return getInvalid();
696     }
697
698     /**
699      * \brief Evaluate a polynomial
700      * \param values an unordered_map<ULL, const MultivariatePolynomial *> containing mapping between var number and polynomial
701      *        or a std::vector<const MultivariatePolynomial *> (0->vector[0], 1->vector[1], ...)
702      * \return the result of the evaluation
703      */
704     template<typename T>
705     inline MultivariatePolynomial eval(T && values) const
706     {
707         if (isInvalid())
708         {
709             return getInvalid();
710         }
711         if (!MultivariatePolynomial::__isValid(values))
712         {
713             return getInvalid();
714         }
715
716         std::unordered_map<unsigned long long, std::set<unsigned int>> expected_exps;
717         for (const auto & m : polynomial)
718         {
719             for (const auto & ve : m.monomial)
720             {
721                 if (ve.exp >= 2 && MultivariatePolynomial::__contains(values, ve.var))
722                 {
723                     expected_exps[ve.var].emplace(ve.exp);
724                 }
725             }
726         }
727
728         std::unordered_map<unsigned long long, std::unordered_map<unsigned int, MultivariatePolynomial>> exps;
729         for (const auto & p : expected_exps)
730         {
731             if (p.second.size() == 1)
732             {
733                 const unsigned int expo = *p.second.begin();
734                 auto & map = exps.emplace(p.first, std::unordered_map<unsigned int, MultivariatePolynomial>()).first->second;
735                 map.emplace(expo, (*values[p.first]) ^ expo);
736             }
737             else
738             {
739                 unsigned int estim = 0;
740                 for (const auto i : p.second)
741                 {
742                     estim += tools::popcount(i) + tools::log2(i) - 1;
743                 }
744                 unsigned int max = *std::prev(p.second.end());
745                 auto & map = exps.emplace(p.first, std::unordered_map<unsigned int, MultivariatePolynomial>()).first->second;
746
747                 // if we have p^3, p^5, p^6 and p^9 to compute:
748                 //   i) cost of p^3 is 2 mults
749                 //   ii) cost of p^5 is 3 mults
750                 //   iii) cost of p^6 is 3 mults
751                 //   iv) cost of p^9 is 4 mults
752                 //   v) total cost is 12
753                 //   vi) 12 > 9 so we compute p^2, p^3, p^4, ..., p^9 and we retains only the exponents 3, 5, 6, 9
754                 // if we just have p^3, p^9:
755                 //   i) total cost is 6
756                 //   ii) we compute p^3 and p^9 separatly in using fast exponentiation
757                 if (estim > max)
758                 {
759                     MultivariatePolynomial mp(*values[p.first]);
760                     auto it = p.second.begin();
761                     for (unsigned int i = 2; i <= max; ++i)
762                     {
763                         mp *= *values[p.first];
764                         if (i == *it)
765                         {
766                             map.emplace(i, mp);
767                             ++it;
768                         }
769                     }
770                 }
771                 else
772                 {
773                     for (const auto expo : p.second)
774                     {
775                         map.emplace(expo, (*values[p.first]) ^ expo);
776                     }
777                 }
778             }
779         }
780
781         MultivariatePolynomial res(constant);
782         for (const auto & m : polynomial)
783         {
784             MultivariatePolynomial r(m.coeff);
785             for (const auto & ve : m.monomial)
786             {
787                 const MultivariatePolynomial * mp = MultivariatePolynomial::__get(values, ve.var);
788                 if (mp)
789                 {
790                     if (ve.exp == 1)
791                     {
792                         r *= *mp;
793                     }
794                     else if (ve.exp > 1)
795                     {
796                         r *= exps[ve.var][ve.exp];
797                     }
798                 }
799                 else
800                 {
801                     MultivariateMonomial mm(1.);
802                     r *= mm.add(ve);
803                 }
804             }
805             res += r;
806         }
807
808         return res;
809     }
810
811     /**
812      * \brief Check positivity
813      * \return true if all the coeffs are positive and the exponents are even
814      */
815     inline bool isPositive() const
816     {
817         if (constant >= 0)
818         {
819             for (const auto & m : polynomial)
820             {
821                 if (m.coeff >= 0)
822                 {
823                     for (const auto & ve : m.monomial)
824                     {
825                         if (ve.exp % 2) // exp is odd
826                         {
827                             return false;
828                         }
829                     }
830                 }
831                 else
832                 {
833                     return false;
834                 }
835             }
836             return true;
837         }
838         return false;
839     }
840
841     /**
842      * \brief Check strict positivity
843      * \param checkConstant if true the constant is checked too
844      * \return true if all the coeffs are strict positive
845      */
846     inline bool isCoeffStrictPositive(const bool checkConstant = true) const
847     {
848         if (!checkConstant || (constant > 0))
849         {
850             for (const auto & m : polynomial)
851             {
852                 if (m.coeff <= 0)
853                 {
854                     return false;
855                 }
856             }
857             return true;
858         }
859         return false;
860     }
861
862     /**
863      * \brief Check positivity
864      * \param checkConstant if true the constant is checked too
865      * \return true if all the coeffs are positive
866      */
867     inline bool isCoeffPositive(const bool checkConstant = true) const
868     {
869         if (!checkConstant || (constant >= 0))
870         {
871             for (const auto & m : polynomial)
872             {
873                 if (m.coeff < 0)
874                 {
875                     return false;
876                 }
877             }
878             return true;
879         }
880         return false;
881     }
882
883     /**
884      * \brief Check strict negativity
885      * \param checkConstant if true the constant is checked too
886      * \return true if all the coeffs are strict negative
887      */
888     inline bool isCoeffStrictNegative(const bool checkConstant = true) const
889     {
890         if (!checkConstant || (constant < 0))
891         {
892             for (const auto & m : polynomial)
893             {
894                 if (m.coeff >= 0)
895                 {
896                     return false;
897                 }
898             }
899             return true;
900         }
901         return false;
902     }
903
904     /**
905      * \brief Check negativity
906      * \param checkConstant if true the constant is checked too
907      * \return true if all the coeffs are negative
908      */
909     inline bool isCoeffNegative(const bool checkConstant = true) const
910     {
911         if (!checkConstant || (constant <= 0))
912         {
913             for (const auto & m : polynomial)
914             {
915                 if (m.coeff > 0)
916                 {
917                     return false;
918                 }
919             }
920             return true;
921         }
922         return false;
923     }
924
925     /**
926      * \brief Helper function to print a polynomial
927      * \param vars a mapping between vars numbers and wstring representation
928      * \return a wstring representing this polynomial
929      */
930     inline const std::wstring print(const std::map<unsigned long long, std::wstring> & vars) const
931     {
932         std::wostringstream wos;
933         wos << constant;
934         std::set<MultivariateMonomial, MultivariateMonomial::Compare> s(polynomial.begin(), polynomial.end());
935         for (const auto & m : s)
936         {
937             wos << L" + " << m.print(vars);
938         }
939         return wos.str();
940     }
941
942     /**
943      * \brief Overload of << operator
944      */
945     friend inline std::wostream & operator<<(std::wostream & out, const MultivariatePolynomial & p)
946     {
947         const std::map<unsigned long long, std::wstring> vars;
948         out << p.constant;
949         std::set<MultivariateMonomial, MultivariateMonomial::Compare> s(p.polynomial.begin(), p.polynomial.end());
950         for (const auto & m : s)
951         {
952             out << L" + " << m.print(vars);
953         }
954         return out;
955     }
956
957     /**
958      * \return true if the two polynomials are differing only by a constant
959      */
960     inline bool isDiffConstant(const MultivariatePolynomial & R) const
961     {
962         return polynomial == R.polynomial;
963     }
964
965     /**
966      * \return true if the polynomial is constant
967      */
968     inline bool isConstant() const
969     {
970         return polynomial.empty();
971     }
972
973     /**
974      * \brief Check if a polynomial is constant and equal to a value
975      * \param val the constant value to check
976      * \return true if the polynomial is constant and equal to val
977      */
978     inline bool isConstant(const double val) const
979     {
980         return isConstant() && constant == val;
981     }
982
983     /**
984      * \brief Get the common coeff for all the monomials
985      * \param[out] common the common value
986      * \return true if there is a common coeff
987      */
988     inline bool getCommonCoeff(double & common) const
989     {
990         if (constant != 0)
991         {
992             return false;
993         }
994         if (polynomial.empty())
995         {
996             common = constant;
997             return true;
998         }
999
1000         common = polynomial.begin()->coeff;
1001         for (Polynomial::const_iterator i = std::next(polynomial.begin()), e = polynomial.end(); i != e; ++i)
1002         {
1003             if (i->coeff != common)
1004             {
1005                 return false;
1006             }
1007         }
1008         return true;
1009     }
1010
1011     /**
1012      * \brief Overload of == operator
1013      */
1014     inline bool operator==(const MultivariatePolynomial & R) const
1015     {
1016         return constant == R.constant && polynomial == R.polynomial;
1017     }
1018
1019     /**
1020      * \brief Overload of != operator
1021      */
1022     inline bool operator!=(const MultivariatePolynomial & R) const
1023     {
1024         return !(*this == R);
1025     }
1026
1027     /**
1028      * \brief Overload of == operator
1029      */
1030     inline bool operator==(const double R) const
1031     {
1032         return polynomial.empty() && constant == R;
1033     }
1034
1035     /**
1036      * \brief Overload of != operator
1037      */
1038     inline bool operator!=(const double R) const
1039     {
1040         return !(*this == R);
1041     }
1042
1043     /**
1044      * \brief Overload of == operator
1045      */
1046     friend inline bool operator==(const double L, const MultivariatePolynomial & R)
1047     {
1048         return R == L;
1049     }
1050
1051     /**
1052      * \brief Overload of != operator
1053      */
1054     friend inline bool operator!=(const double L, const MultivariatePolynomial & R)
1055     {
1056         return R != L;
1057     }
1058
1059     /**
1060      * \return a hash
1061      */
1062     inline std::size_t hash() const
1063     {
1064         std::size_t h = std::hash<double>()(constant);
1065         for (const auto & m : polynomial)
1066         {
1067             // since the order of the monomials is not always the same
1068             // we must use a commutative operation to combine the monomial's hashes
1069             h += tools::hash_combine(std::hash<double>()(m.coeff), MultivariateMonomial::Hash()(m));
1070         }
1071
1072         return h;
1073     }
1074
1075     /**
1076      * \struct Hash
1077      * \brief To be used in an unordered container
1078      */
1079     struct Hash
1080     {
1081         inline std::size_t operator()(const MultivariatePolynomial & mp) const
1082         {
1083             return mp.hash();
1084         }
1085     };
1086
1087     /**
1088      * \struct Eq
1089      * \brief To be used in an unordered container
1090      */
1091     struct Eq
1092     {
1093         inline bool operator()(const MultivariatePolynomial & L, const MultivariatePolynomial & R) const
1094         {
1095             return L == R;
1096         }
1097     };
1098
1099 private:
1100
1101     // Helper function to use with eval
1102     inline static bool __isValid(const std::unordered_map<unsigned long long, const MultivariatePolynomial *> & values)
1103     {
1104         for (const auto & p : values)
1105         {
1106             if (p.second->isInvalid())
1107             {
1108                 return false;
1109             }
1110         }
1111         return true;
1112     }
1113
1114     // Helper function to use with eval
1115     inline static bool __isValid(const std::vector<const MultivariatePolynomial *> & values)
1116     {
1117         for (const auto & v : values)
1118         {
1119             if (v->isInvalid())
1120             {
1121                 return false;
1122             }
1123         }
1124         return true;
1125     }
1126
1127     // Helper function to use with eval
1128     inline static bool __contains(const std::unordered_map<unsigned long long, const MultivariatePolynomial *> & values, const unsigned long long val)
1129     {
1130         return values.find(val) != values.end();
1131     }
1132
1133     // Helper function to use with eval
1134     inline static bool __contains(const std::vector<const MultivariatePolynomial *> & values, const unsigned long long val)
1135     {
1136         return val < values.size();
1137     }
1138
1139     // Helper function to use with eval
1140     inline static const MultivariatePolynomial * __get(const std::unordered_map<unsigned long long, const MultivariatePolynomial *> & values, const unsigned long long val)
1141     {
1142         const auto i = values.find(val);
1143         if (i != values.end())
1144         {
1145             return i->second;
1146         }
1147         return nullptr;
1148     }
1149
1150     // Helper function to use with eval
1151     inline static const MultivariatePolynomial * __get(const std::vector<const MultivariatePolynomial *> & values, const unsigned long long val)
1152     {
1153         return values[val];
1154     }
1155 };
1156
1157 } // namespace analysis
1158
1159 #endif // __MULTIVARIATE_POLYNOMIAL_HXX__