Analysis: fix bugs
[scilab.git] / scilab / modules / ast / src / cpp / analysis / MultivariatePolynomial.cpp
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 #include "gvn/MultivariatePolynomial.hxx"
14
15 namespace analysis
16 {
17
18 MultivariatePolynomial MultivariatePolynomial::getInvalid()
19 {
20     return MultivariatePolynomial(0, false);
21 }
22
23 bool MultivariatePolynomial::isValid() const
24 {
25     return valid;
26 }
27
28 bool MultivariatePolynomial::isInvalid() const
29 {
30     return !valid;
31 }
32
33 void MultivariatePolynomial::invalid()
34 {
35     constant = 0;
36     valid = false;
37     polynomial.clear();
38 }
39
40 bool MultivariatePolynomial::contains(const uint64_t var) const
41 {
42     for (const auto & m : polynomial)
43     {
44         if (m.contains(var))
45         {
46             return true;
47         }
48     }
49
50     return false;
51 }
52
53 bool MultivariatePolynomial::checkVariable(const uint64_t max) const
54 {
55     for (const auto & m : polynomial)
56     {
57         if (!m.checkVariable(max))
58         {
59             return false;
60         }
61     }
62     return true;
63 }
64
65 bool MultivariatePolynomial::containsVarsGEq(const uint64_t min) const
66 {
67     for (const auto & m : polynomial)
68     {
69         if (m.monomial.lower_bound(min) != m.monomial.end())
70         {
71             return true;
72         }
73     }
74
75     return false;
76 }
77
78 MultivariatePolynomial MultivariatePolynomial::translateVariables(const uint64_t t, const uint64_t min) const
79 {
80     MultivariatePolynomial mp(polynomial.size(), constant);
81     for (const auto & m : polynomial)
82     {
83         MultivariateMonomial mm(m);
84         MultivariateMonomial::Monomial::iterator i = mm.monomial.lower_bound(min);
85         if (i != mm.monomial.end())
86         {
87             // We don't modify the order in the set, so we can const_cast
88             for (MultivariateMonomial::Monomial::iterator j = std::prev(mm.monomial.end()); j != i; --j)
89             {
90                 const_cast<VarExp &>(*j).var += t;
91             }
92             const_cast<VarExp &>(*i).var += t;
93         }
94         mp.add(mm);
95     }
96
97     return mp;
98 }
99
100 MultivariatePolynomial & MultivariatePolynomial::add(const MultivariateMonomial & m, const int64_t coeff)
101 {
102     const int64_t c = m.coeff * coeff;
103     if (c)
104     {
105         Polynomial::iterator i = polynomial.find(m);
106         if (i == polynomial.end())
107         {
108             Polynomial::iterator j = polynomial.insert(m).first;
109             j->coeff = c;
110         }
111         else
112         {
113             if (i->coeff == -c)
114             {
115                 polynomial.erase(i);
116             }
117             else
118             {
119                 i->coeff += c;
120             }
121         }
122     }
123     return *this;
124 }
125
126 void MultivariatePolynomial::sub(const MultivariateMonomial & m)
127 {
128     Polynomial::iterator i = polynomial.find(m);
129     if (i == polynomial.end())
130     {
131         if (m.coeff)
132         {
133             polynomial.insert(m).first->coeff = -m.coeff;
134         }
135     }
136     else
137     {
138         if (i->coeff == m.coeff)
139         {
140             polynomial.erase(i);
141         }
142         else
143         {
144             i->coeff -= m.coeff;
145         }
146     }
147 }
148
149 MultivariatePolynomial MultivariatePolynomial::operator+(const MultivariateMonomial & R) const
150 {
151     if (isValid())
152     {
153         MultivariatePolynomial res(*this);
154         res.add(R);
155         return res;
156     }
157     return getInvalid();
158 }
159
160 MultivariatePolynomial & MultivariatePolynomial::operator+=(const MultivariateMonomial & R)
161 {
162     if (isValid())
163     {
164         add(R);
165     }
166     return *this;
167 }
168
169 MultivariatePolynomial MultivariatePolynomial::operator+(const int64_t R) const
170 {
171     if (isValid())
172     {
173         MultivariatePolynomial res(*this);
174         res.constant += R;
175         return res;
176     }
177     return *this;
178 }
179
180 MultivariatePolynomial & MultivariatePolynomial::operator+=(const int64_t R)
181 {
182     if (isValid())
183     {
184         constant += R;
185     }
186     return *this;
187 }
188
189 MultivariatePolynomial MultivariatePolynomial::operator-(const MultivariateMonomial & R) const
190 {
191     if (isValid())
192     {
193         MultivariatePolynomial res(*this);
194         res.sub(R);
195         return res;
196     }
197     return *this;
198 }
199
200 MultivariatePolynomial & MultivariatePolynomial::operator-=(const MultivariateMonomial & R)
201 {
202     if (isValid())
203     {
204         sub(R);
205     }
206     return *this;
207 }
208
209 MultivariatePolynomial MultivariatePolynomial::operator-(const int64_t R) const
210 {
211     if (isValid())
212     {
213         MultivariatePolynomial res(*this);
214         res.constant -= R;
215         return res;
216     }
217     return *this;
218 }
219
220 MultivariatePolynomial MultivariatePolynomial::operator-() const
221 {
222     if (isValid())
223     {
224         MultivariatePolynomial res(*this);
225         res.constant = -res.constant;
226         for (auto & m : res.polynomial)
227         {
228             m.coeff = -m.coeff;
229         }
230         return res;
231     }
232     return *this;
233 }
234
235 MultivariatePolynomial & MultivariatePolynomial::operator-=(const int64_t R)
236 {
237     if (isValid())
238     {
239         constant -= R;
240     }
241     return *this;
242 }
243
244 MultivariatePolynomial MultivariatePolynomial::operator+(const MultivariatePolynomial & R) const
245 {
246     if (isValid() && R.isValid())
247     {
248         MultivariatePolynomial res(*this);
249         res.constant += R.constant;
250         for (const auto & m : R.polynomial)
251         {
252             res.add(m);
253         }
254         return res;
255     }
256     return getInvalid();
257 }
258
259 MultivariatePolynomial & MultivariatePolynomial::operator+=(const MultivariatePolynomial & R)
260 {
261     if (isValid() && R.isValid())
262     {
263         constant += R.constant;
264         for (const auto & m : R.polynomial)
265         {
266             add(m);
267         }
268     }
269     else
270     {
271         invalid();
272     }
273
274     return *this;
275 }
276
277 MultivariatePolynomial MultivariatePolynomial::operator-(const MultivariatePolynomial & R) const
278 {
279     if (isValid() && R.isValid())
280     {
281         MultivariatePolynomial res(*this);
282         res.constant -= R.constant;
283         for (const auto & m : R.polynomial)
284         {
285             res.sub(m);
286         }
287         return res;
288     }
289     return getInvalid();
290 }
291
292 MultivariatePolynomial & MultivariatePolynomial::operator-=(const MultivariatePolynomial & R)
293 {
294     if (isValid() && R.isValid())
295     {
296         constant -= R.constant;
297         for (const auto & m : R.polynomial)
298         {
299             sub(m);
300         }
301     }
302     else
303     {
304         invalid();
305     }
306     return *this;
307 }
308
309 MultivariatePolynomial MultivariatePolynomial::operator*(const MultivariatePolynomial & R) const
310 {
311     if (isValid() && R.isValid())
312     {
313         if (R.isConstant())
314         {
315             return *this * R.constant;
316         }
317         else if (isConstant())
318         {
319             return R * constant;
320         }
321         else
322         {
323             MultivariatePolynomial res(static_cast<unsigned int>((polynomial.size() + 1) * (R.polynomial.size() + 1) - 1), constant * R.constant);
324             for (const auto & mR : R.polynomial)
325             {
326                 res.add(mR, constant);
327             }
328             for (const auto & mL : polynomial)
329             {
330                 res.add(mL, R.constant);
331                 for (const auto & mR : R.polynomial)
332                 {
333                     res.add(mL * mR);
334                 }
335             }
336             return res;
337         }
338     }
339     return getInvalid();
340 }
341
342 MultivariatePolynomial MultivariatePolynomial::operator/(const MultivariatePolynomial & R) const
343 {
344     if (isValid() && R.isValid())
345     {
346         if (R.isConstant())
347         {
348             if (R.constant != 1)
349             {
350                 return *this / R.constant;
351             }
352             else
353             {
354                 return *this;
355             }
356         }
357     }
358     return getInvalid();
359 }
360
361 MultivariatePolynomial & MultivariatePolynomial::operator/=(const MultivariatePolynomial & R)
362 {
363     if (isValid() && R.isValid())
364     {
365         if (R.polynomial.empty())
366         {
367             constant /= R.constant;
368             for (auto & m : polynomial)
369             {
370                 m.coeff /= R.constant;
371             }
372         }
373         else
374         {
375             MultivariatePolynomial res(*this / R);
376             polynomial = res.polynomial;
377             constant = res.constant;
378         }
379     }
380     else
381     {
382         invalid();
383     }
384     return *this;
385 }
386
387 MultivariatePolynomial & MultivariatePolynomial::operator*=(const MultivariatePolynomial & R)
388 {
389     if (isValid() && R.isValid())
390     {
391         if (R.polynomial.empty())
392         {
393             constant *= R.constant;
394             for (auto & m : polynomial)
395             {
396                 m.coeff *= R.constant;
397             }
398         }
399         else
400         {
401             MultivariatePolynomial res(*this * R);
402             polynomial = res.polynomial;
403             constant = res.constant;
404         }
405     }
406     else
407     {
408         invalid();
409     }
410     return *this;
411 }
412
413 MultivariatePolynomial MultivariatePolynomial::operator*(const MultivariateMonomial & R) const
414 {
415     if (isValid())
416     {
417         MultivariatePolynomial res(static_cast<unsigned int>(polynomial.size() + 1), int64_t(0));
418         res.add(constant * R);
419         for (const auto & mL : polynomial)
420         {
421             res.add(mL * R);
422         }
423         return res;
424     }
425     else
426     {
427         return getInvalid();
428     }
429 }
430
431 MultivariatePolynomial & MultivariatePolynomial::operator*=(const MultivariateMonomial & R)
432 {
433     if (isValid())
434     {
435         MultivariatePolynomial res = *this * R;
436         polynomial = res.polynomial;
437         constant = res.constant;
438     }
439     return *this;
440 }
441
442 MultivariatePolynomial MultivariatePolynomial::operator*(const int64_t R) const
443 {
444     if (isValid())
445     {
446         if (R)
447         {
448             if (R == 1)
449             {
450                 return *this;
451             }
452             else
453             {
454                 MultivariatePolynomial res(*this);
455                 res.constant *= R;
456                 for (auto & m : res.polynomial)
457                 {
458                     m.coeff *= R;
459                 }
460                 return res;
461             }
462         }
463         else
464         {
465             return MultivariatePolynomial(int64_t(0));
466         }
467     }
468     return getInvalid();
469 }
470
471 MultivariatePolynomial & MultivariatePolynomial::operator*=(const int64_t R)
472 {
473     if (isValid())
474     {
475         if (R == 0)
476         {
477             constant = 0;
478             polynomial.clear();
479         }
480         else if (R != 1)
481         {
482             constant *= R;
483             for (auto & m : polynomial)
484             {
485                 m.coeff *= R;
486             }
487         }
488     }
489     return *this;
490 }
491
492 MultivariatePolynomial MultivariatePolynomial::operator/(const int64_t R) const
493 {
494     if (isValid())
495     {
496         if (R != 1)
497         {
498             MultivariatePolynomial res(*this);
499             res.constant /= R;
500             for (auto & m : res.polynomial)
501             {
502                 m.coeff /= R;
503             }
504             return res;
505         }
506     }
507     return *this;
508 }
509
510 MultivariatePolynomial & MultivariatePolynomial::operator/=(const int64_t R)
511 {
512     if (isValid())
513     {
514         if (R != 1)
515         {
516             constant /= R;
517             for (auto & m : polynomial)
518             {
519                 m.coeff /= R;
520             }
521         }
522     }
523     return *this;
524 }
525
526 MultivariatePolynomial MultivariatePolynomial::operator^(unsigned int R) const
527 {
528     if (isValid())
529     {
530         if (R == 0)
531         {
532             return MultivariatePolynomial(int64_t(1));
533         }
534         else if (R == 1)
535         {
536             return *this;
537         }
538         else
539         {
540             if (constant == 0)
541             {
542                 if (polynomial.empty())
543                 {
544                     return MultivariatePolynomial(int64_t(0));
545                 }
546                 else if (polynomial.size() == 1)
547                 {
548                     const MultivariateMonomial & m = *polynomial.begin();
549                     MultivariatePolynomial res;
550                     res.polynomial.emplace(m ^ R);
551
552                     return res;
553                 }
554             }
555
556             if (polynomial.empty())
557             {
558                 return MultivariatePolynomial(tools::powui(constant, R));
559             }
560
561             MultivariatePolynomial p = *this;
562             MultivariatePolynomial y = (R & 1) ? *this : MultivariatePolynomial(int64_t(1));
563
564             while (R >>= 1)
565             {
566                 p *= p;
567                 if (R & 1)
568                 {
569                     y *= p;
570                 }
571             }
572
573             return y;
574         }
575     }
576     return getInvalid();
577 }
578
579 MultivariatePolynomial MultivariatePolynomial::operator^(const MultivariatePolynomial & R) const
580 {
581     if (isValid() && R.isValid())
582     {
583         if (R.isConstant() && R.constant == (unsigned int)R.constant)
584         {
585             return (*this) ^ ((unsigned int)R.constant);
586         }
587     }
588     return getInvalid();
589 }
590
591 bool MultivariatePolynomial::isDivisibleBy(const int64_t n) const
592 {
593     if (constant % n == 0)
594     {
595         for (const auto & m : polynomial)
596         {
597             if (m.coeff % n != 0)
598             {
599                 return false;
600             }
601         }
602         return true;
603     }
604     return false;
605 }
606
607 bool MultivariatePolynomial::isDivisibleBy(const MultivariatePolynomial & mp) const
608 {
609     if (mp.polynomial.empty())
610     {
611         return isDivisibleBy(mp.constant);
612     }
613
614     return false;
615 }
616
617 bool MultivariatePolynomial::isPositive() const
618 {
619     if (constant >= 0)
620     {
621         for (const auto & m : polynomial)
622         {
623             if (m.coeff >= 0)
624             {
625                 for (const auto & ve : m.monomial)
626                 {
627                     if (ve.exp % 2) // exp is odd
628                     {
629                         return false;
630                     }
631                 }
632             }
633             else
634             {
635                 return false;
636             }
637         }
638         return true;
639     }
640     return false;
641 }
642
643 bool MultivariatePolynomial::isCoeffStrictPositive(const bool checkConstant) const
644 {
645     if (!checkConstant || (constant > 0))
646     {
647         for (const auto & m : polynomial)
648         {
649             if (m.coeff <= 0)
650             {
651                 return false;
652             }
653         }
654         return true;
655     }
656     return false;
657 }
658
659 bool MultivariatePolynomial::isCoeffPositive(const bool checkConstant) const
660 {
661     if (!checkConstant || (constant >= 0))
662     {
663         for (const auto & m : polynomial)
664         {
665             if (m.coeff < 0)
666             {
667                 return false;
668             }
669         }
670         return true;
671     }
672     return false;
673 }
674
675 bool MultivariatePolynomial::isCoeffStrictNegative(const bool checkConstant) const
676 {
677     if (!checkConstant || (constant < 0))
678     {
679         for (const auto & m : polynomial)
680         {
681             if (m.coeff >= 0)
682             {
683                 return false;
684             }
685         }
686         return true;
687     }
688     return false;
689 }
690
691 bool MultivariatePolynomial::isCoeffNegative(const bool checkConstant) const
692 {
693     if (!checkConstant || (constant <= 0))
694     {
695         for (const auto & m : polynomial)
696         {
697             if (m.coeff > 0)
698             {
699                 return false;
700             }
701         }
702         return true;
703     }
704     return false;
705 }
706
707 const std::wstring MultivariatePolynomial::print(const std::map<uint64_t, std::wstring> & vars) const
708 {
709     std::wostringstream wos;
710     wos << constant;
711     std::set<MultivariateMonomial, MultivariateMonomial::Compare> s(polynomial.begin(), polynomial.end());
712     for (const auto & m : s)
713     {
714         wos << L" + " << m.print(vars);
715     }
716     return wos.str();
717 }
718
719 std::wostream & operator<<(std::wostream & out, const MultivariatePolynomial & p)
720 {
721     const std::map<uint64_t, std::wstring> vars;
722     out << p.constant;
723     std::set<MultivariateMonomial, MultivariateMonomial::Compare> s(p.polynomial.begin(), p.polynomial.end());
724     for (const auto & m : s)
725     {
726         out << L" + " << m.print(vars);
727     }
728     return out;
729 }
730
731 bool MultivariatePolynomial::isDiffConstant(const MultivariatePolynomial & R) const
732 {
733     return polynomial == R.polynomial;
734 }
735
736 bool MultivariatePolynomial::isConstant() const
737 {
738     return polynomial.empty();
739 }
740
741 bool MultivariatePolynomial::isConstant(const int64_t val) const
742 {
743     return isConstant() && constant == val;
744 }
745
746 bool MultivariatePolynomial::getCommonCoeff(int64_t & common) const
747 {
748     if (constant != 0)
749     {
750         return false;
751     }
752     if (polynomial.empty())
753     {
754         common = constant;
755         return true;
756     }
757
758     common = polynomial.begin()->coeff;
759     for (Polynomial::const_iterator i = std::next(polynomial.begin()), e = polynomial.end(); i != e; ++i)
760     {
761         if (i->coeff != common)
762         {
763             return false;
764         }
765     }
766     return true;
767 }
768
769 bool MultivariatePolynomial::operator==(const MultivariatePolynomial & R) const
770 {
771     return constant == R.constant && polynomial == R.polynomial;
772 }
773
774 bool MultivariatePolynomial::operator!=(const MultivariatePolynomial & R) const
775 {
776     return !(*this == R);
777 }
778
779 bool MultivariatePolynomial::operator==(const int64_t R) const
780 {
781     return polynomial.empty() && constant == R;
782 }
783
784 bool MultivariatePolynomial::operator!=(const int64_t R) const
785 {
786     return !(*this == R);
787 }
788
789 bool operator==(const int64_t L, const MultivariatePolynomial & R)
790 {
791     return R == L;
792 }
793
794 bool operator!=(const int64_t L, const MultivariatePolynomial & R)
795 {
796     return R != L;
797 }
798
799 std::size_t MultivariatePolynomial::hash() const
800 {
801     std::size_t h = std::hash<int64_t>()(constant);
802     for (const auto & m : polynomial)
803     {
804         // since the order of the monomials is not always the same
805         // we must use a commutative operation to combine the monomial's hashes
806         h += tools::hash_combine(std::hash<int64_t>()(m.coeff), MultivariateMonomial::Hash()(m));
807     }
808
809     return h;
810 }
811
812 bool MultivariatePolynomial::__isValid(const std::unordered_map<uint64_t, const MultivariatePolynomial *> & values)
813 {
814     for (const auto & p : values)
815     {
816         if (p.second->isInvalid())
817         {
818             return false;
819         }
820     }
821     return true;
822 }
823
824 bool MultivariatePolynomial::__isValid(const std::vector<const MultivariatePolynomial *> & values)
825 {
826     for (const auto & v : values)
827     {
828         if (v->isInvalid())
829         {
830             return false;
831         }
832     }
833     return true;
834 }
835
836 bool MultivariatePolynomial::__isValid(const std::pair<uint64_t, const MultivariatePolynomial *> & values)
837 {
838     return values.second->isValid();
839 }
840
841 bool MultivariatePolynomial::__contains(const std::unordered_map<uint64_t, const MultivariatePolynomial *> & values, const uint64_t val)
842 {
843     return values.find(val) != values.end();
844 }
845
846 bool MultivariatePolynomial::__contains(const std::vector<const MultivariatePolynomial *> & values, const uint64_t val)
847 {
848     return val < values.size();
849 }
850
851 bool MultivariatePolynomial::__contains(const std::pair<uint64_t, const MultivariatePolynomial *> & values, const uint64_t val)
852 {
853     return values.first == val;
854 }
855
856 const MultivariatePolynomial * MultivariatePolynomial::__get(const std::unordered_map<uint64_t, const MultivariatePolynomial *> & values, const uint64_t val)
857 {
858     const auto i = values.find(val);
859     if (i != values.end())
860     {
861         return i->second;
862     }
863     return nullptr;
864 }
865
866 const MultivariatePolynomial * MultivariatePolynomial::__getSafe(const std::unordered_map<uint64_t, const MultivariatePolynomial *> & values, const uint64_t val)
867 {
868     return values.find(val)->second;
869 }
870
871 const MultivariatePolynomial * MultivariatePolynomial::__get(const std::vector<const MultivariatePolynomial *> & values, const uint64_t val)
872 {
873     if (val < values.size())
874     {
875         return values[val];
876     }
877     return nullptr;
878 }
879
880 const MultivariatePolynomial * MultivariatePolynomial::__getSafe(const std::vector<const MultivariatePolynomial *> & values, const uint64_t val)
881 {
882     return values[val];
883 }
884
885 const MultivariatePolynomial * MultivariatePolynomial::__get(const std::pair<uint64_t, const MultivariatePolynomial *> & values, const uint64_t val)
886 {
887     if (values.first == val)
888     {
889         return values.second;
890     }
891     return nullptr;
892 }
893
894 const MultivariatePolynomial * MultivariatePolynomial::__getSafe(const std::pair<uint64_t, const MultivariatePolynomial *> & values, const uint64_t val)
895 {
896     return values.second;
897 }
898
899 } // namespace analysis