AST Inspector: speed-up things while executing
[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_map>
24 #include <unordered_set>
25 #include <vector>
26
27 #include "core_math.h"
28 #include "tools.hxx"
29 #include "MultivariateMonomial.hxx"
30 #include "dynlib_ast.h"
31
32 namespace analysis
33 {
34 /**
35  * \struct MultivariatePolynomial
36  * \brief Represents a multivariate polynomial
37  */
38 struct EXTERN_AST MultivariatePolynomial
39 {
40     typedef std::unordered_set<MultivariateMonomial, MultivariateMonomial::Hash, MultivariateMonomial::Eq> Polynomial;
41     int64_t constant;
42     bool valid;
43     Polynomial polynomial;
44
45     /**
46      * \brief constructor
47      * \param var to init polynomial
48      */
49     MultivariatePolynomial(const uint64_t var) : constant(0), valid(true)
50     {
51         polynomial.emplace(var);
52     }
53
54     /**
55      * \brief constructor
56      * \param _constant to init polynomial
57      */
58     MultivariatePolynomial(const int64_t _constant = 0, const bool _valid = true) : constant(_constant), valid(_valid) { }
59
60     /**
61      * \brief constructor
62      * \param _size the size of the polynomial (used to reserve the unordered_set)
63      * \param _constant to init polynomial
64      */
65     MultivariatePolynomial(const unsigned int _size, const int64_t _constant) : constant(_constant), valid(true), polynomial(_size) { }
66
67     /**
68      * \brief copy constructor
69      */
70     MultivariatePolynomial(const MultivariatePolynomial & mp) : constant(mp.constant), valid(mp.valid), polynomial(mp.polynomial) { }
71
72     /**
73      * \brief Get an invalid polynomial (i.e. constant == NaN)
74      * \return invalid polynomial
75      */
76     static MultivariatePolynomial getInvalid();
77
78     /**
79      * \brief Check if it is valid
80      * \return true if it is valid
81      */
82     bool isValid() const;
83
84     /**
85      * \brief Check if it is invalid
86      * \return true if it is invalid
87      */
88     bool isInvalid() const;
89
90     /**
91      * \brief Invalidate the polynomial
92      */
93     void invalid();
94
95     /**
96      * \brief Check if a variable is contained in the polynomial
97      * \param var an id
98      * \return true if the polynomial contains the var
99      */
100     bool contains(const uint64_t var) const;
101
102     /**
103      * \brief Check if the variables of the polynomial have an id lower or equal to max
104      * \param max an id
105      * \return true if all the variables have an id leq to max
106      */
107     bool checkVariable(const uint64_t max) const;
108
109     /**
110      * \brief Check if the variables of the polynomial have an id greater or equal to min
111      * \param min an id
112      * \return true if the polynomial contains a var with an id geq than min
113      */
114     bool containsVarsGEq(const uint64_t min) const;
115
116     /**
117      * \brief Translate the variables of the polynomial which have an id greater or equal to min
118      * \param min an id
119      * \return a translated polynomial
120      */
121     MultivariatePolynomial translateVariables(const uint64_t t, const uint64_t min) const;
122
123     /**
124      * \brief Add a monomial in the polynomial
125      * \param m the monomial to add
126      * \param coeff the multiplcative coefficient to applicate to the monomial
127      * \return *this
128      */
129     MultivariatePolynomial & add(const MultivariateMonomial & m, const int64_t coeff = 1);
130
131     /**
132      * \brief Subtract a monomial in the polynomial
133      * \param m the monomial to add
134      */
135     void sub(const MultivariateMonomial & m);
136
137     /**
138      * \brief Overload of the + operator
139      */
140     MultivariatePolynomial operator+(const MultivariateMonomial & R) const;
141
142     /**
143      * \brief Overload of the += operator
144      */
145     MultivariatePolynomial & operator+=(const MultivariateMonomial & R);
146
147     /**
148      * \brief Overload of the + operator
149      */
150     MultivariatePolynomial operator+(const int64_t R) const;
151
152     /**
153      * \brief Overload of the += operator
154      */
155     MultivariatePolynomial & operator+=(const int64_t R);
156
157     /**
158      * \brief Overload of the - operator
159      */
160     MultivariatePolynomial operator-(const MultivariateMonomial & R) const;
161
162     /**
163      * \brief Overload of the -= operator
164      */
165     MultivariatePolynomial & operator-=(const MultivariateMonomial & R);
166
167     /**
168      * \brief Overload of the - operator
169      */
170     MultivariatePolynomial operator-(const int64_t R) const;
171
172     /**
173      * \brief Overload of the - operator (unary)
174      */
175     MultivariatePolynomial operator-() const;
176
177     /**
178      * \brief Overload of the -= operator
179      */
180     MultivariatePolynomial & operator-=(const int64_t R);
181
182     /**
183      * \brief Overload of the + operator
184      */
185     MultivariatePolynomial operator+(const MultivariatePolynomial & R) const;
186
187     /**
188      * \brief Overload of the += operator
189      */
190     MultivariatePolynomial & operator+=(const MultivariatePolynomial & R);
191
192     /**
193      * \brief Overload of the - operator
194      */
195     MultivariatePolynomial operator-(const MultivariatePolynomial & R) const;
196
197     /**
198      * \brief Overload of the -= operator
199      */
200     MultivariatePolynomial & operator-=(const MultivariatePolynomial & R);
201
202     /**
203      * \brief Overload of the * operator
204      */
205     MultivariatePolynomial operator*(const MultivariatePolynomial & R) const;
206
207     /**
208      * \brief Overload of the / operator
209      */
210     MultivariatePolynomial operator/(const MultivariatePolynomial & R) const;
211
212     /**
213      * \brief Overload of the /= operator
214      */
215     MultivariatePolynomial & operator/=(const MultivariatePolynomial & R);
216
217     /**
218      * \brief Overload of the *= operator
219      */
220     MultivariatePolynomial & operator*=(const MultivariatePolynomial & R);
221
222     /**
223      * \brief Overload of the * operator
224      */
225     MultivariatePolynomial operator*(const MultivariateMonomial & R) const;
226
227     /**
228      * \brief Overload of the *= operator
229      */
230     MultivariatePolynomial & operator*=(const MultivariateMonomial & R);
231
232     /**
233      * \brief Overload of the * operator
234      */
235     MultivariatePolynomial operator*(const int64_t R) const;
236
237     /**
238      * \brief Overload of the *= operator
239      */
240     MultivariatePolynomial & operator*=(const int64_t R);
241
242     /**
243      * \brief Overload of the / operator
244      */
245     MultivariatePolynomial operator/(const int64_t R) const;
246
247     /**
248      * \brief Overload of the /= operator
249      */
250     MultivariatePolynomial & operator/=(const int64_t R);
251
252     /**
253      * \brief Overload of the ^ operator (exponentiation)
254      */
255     MultivariatePolynomial operator^(unsigned int R) const;
256
257     /**
258      * \brief Overload of the ^ operator (exponentiation)
259      */
260     MultivariatePolynomial operator^(const MultivariatePolynomial & R) const;
261
262     /**
263      * \brief Evaluate a polynomial
264      * \param values an unordered_map<ULL, const MultivariatePolynomial *> containing mapping between var number and polynomial
265      *        or a std::vector<const MultivariatePolynomial *> (0->vector[0], 1->vector[1], ...)
266      * \return the result of the evaluation
267      */
268     template<typename T>
269     inline MultivariatePolynomial eval(T && values) const
270     {
271         if (isInvalid())
272         {
273             return getInvalid();
274         }
275         if (!MultivariatePolynomial::__isValid(values))
276         {
277             return getInvalid();
278         }
279
280         std::unordered_map<uint64_t, std::set<unsigned int>> expected_exps;
281         for (const auto & m : polynomial)
282         {
283             for (const auto & ve : m.monomial)
284             {
285                 if (ve.exp >= 2 && MultivariatePolynomial::__contains(values, ve.var))
286                 {
287                     expected_exps[ve.var].emplace(ve.exp);
288                 }
289             }
290         }
291
292         std::unordered_map<uint64_t, std::unordered_map<unsigned int, MultivariatePolynomial>> exps;
293         for (const auto & p : expected_exps)
294         {
295             if (p.second.size() == 1)
296             {
297                 const unsigned int expo = *p.second.begin();
298                 auto & map = exps.emplace(p.first, std::unordered_map<unsigned int, MultivariatePolynomial>()).first->second;
299                 map.emplace(expo, (*__getSafe(values, p.first)) ^ expo);
300             }
301             else
302             {
303                 unsigned int estim = 0;
304                 for (const auto i : p.second)
305                 {
306                     estim += tools::popcount(i) + tools::log2(i) - 1;
307                 }
308                 unsigned int max = *std::prev(p.second.end());
309                 auto & map = exps.emplace(p.first, std::unordered_map<unsigned int, MultivariatePolynomial>()).first->second;
310
311                 // if we have p^3, p^5, p^6 and p^9 to compute:
312                 //   i) cost of p^3 is 2 mults
313                 //   ii) cost of p^5 is 3 mults
314                 //   iii) cost of p^6 is 3 mults
315                 //   iv) cost of p^9 is 4 mults
316                 //   v) total cost is 12
317                 //   vi) 12 > 9 so we compute p^2, p^3, p^4, ..., p^9 and we retains only the exponents 3, 5, 6, 9
318                 // if we just have p^3, p^9:
319                 //   i) total cost is 6
320                 //   ii) we compute p^3 and p^9 separatly in using fast exponentiation
321                 if (estim > max)
322                 {
323                     MultivariatePolynomial mp(*__getSafe(values, p.first));
324                     auto it = p.second.begin();
325                     for (unsigned int i = 2; i <= max; ++i)
326                     {
327                         mp *= *__getSafe(values, p.first);
328                         if (i == *it)
329                         {
330                             map.emplace(i, mp);
331                             ++it;
332                         }
333                     }
334                 }
335                 else
336                 {
337                     for (const auto expo : p.second)
338                     {
339                         map.emplace(expo, (*__getSafe(values, p.first)) ^ expo);
340                     }
341                 }
342             }
343         }
344
345         MultivariatePolynomial res(constant);
346         for (const auto & m : polynomial)
347         {
348             MultivariatePolynomial r(m.coeff);
349             for (const auto & ve : m.monomial)
350             {
351                 const MultivariatePolynomial * mp = MultivariatePolynomial::__get(values, ve.var);
352                 if (mp)
353                 {
354                     if (ve.exp == 1)
355                     {
356                         r *= *mp;
357                     }
358                     else if (ve.exp > 1)
359                     {
360                         r *= exps[ve.var][ve.exp];
361                     }
362                 }
363                 else
364                 {
365                     MultivariateMonomial mm(int64_t(1));
366                     r *= mm.add(ve);
367                 }
368             }
369             res += r;
370         }
371
372         return res;
373     }
374
375     /**
376      * \brief Check divisibility by an integer
377      * \return true if all the coeffs are divisible by n
378      */
379     bool isDivisibleBy(const int64_t n) const;
380
381     /**
382      * \brief Check divisibility by a polynomial
383      * For now the polynomial must be constant.
384      * \return true if this is divisible by mp
385      */
386     bool isDivisibleBy(const MultivariatePolynomial & mp) const;
387
388     /**
389      * \brief Check positivity
390      * \return true if all the coeffs are positive and the exponents are even
391      */
392     bool isPositive() const;
393
394     /**
395      * \brief Check strict positivity
396      * \param checkConstant if true the constant is checked too
397      * \return true if all the coeffs are strict positive
398      */
399     bool isCoeffStrictPositive(const bool checkConstant = true) const;
400
401     /**
402      * \brief Check positivity
403      * \param checkConstant if true the constant is checked too
404      * \return true if all the coeffs are positive
405      */
406     bool isCoeffPositive(const bool checkConstant = true) const;
407
408     /**
409      * \brief Check strict negativity
410      * \param checkConstant if true the constant is checked too
411      * \return true if all the coeffs are strict negative
412      */
413     bool isCoeffStrictNegative(const bool checkConstant = true) const;
414
415     /**
416      * \brief Check negativity
417      * \param checkConstant if true the constant is checked too
418      * \return true if all the coeffs are negative
419      */
420     bool isCoeffNegative(const bool checkConstant = true) const;
421
422     /**
423      * \brief Helper function to print a polynomial
424      * \param vars a mapping between vars numbers and wstring representation
425      * \return a wstring representing this polynomial
426      */
427     const std::wstring print(const std::map<uint64_t, std::wstring> & vars) const;
428
429     /**
430      * \brief Overload of << operator
431      */
432     friend std::wostream & operator<<(std::wostream & out, const MultivariatePolynomial & p);
433
434     /**
435      * \return true if the two polynomials are differing only by a constant
436      */
437     bool isDiffConstant(const MultivariatePolynomial & R) const;
438
439     /**
440      * \return true if the polynomial is constant
441      */
442     bool isConstant() const;
443
444     /**
445      * \brief Check if a polynomial is constant and equal to a value
446      * \param val the constant value to check
447      * \return true if the polynomial is constant and equal to val
448      */
449     bool isConstant(const int64_t val) const;
450
451     /**
452      * \brief Get the common coeff for all the monomials
453      * \param[out] common the common value
454      * \return true if there is a common coeff
455      */
456     bool getCommonCoeff(int64_t & common) const;
457
458     /**
459      * \brief Overload of == operator
460      */
461     bool operator==(const MultivariatePolynomial & R) const;
462
463     /**
464      * \brief Overload of != operator
465      */
466     bool operator!=(const MultivariatePolynomial & R) const;
467
468     /**
469      * \brief Overload of == operator
470      */
471     bool operator==(const int64_t R) const;
472
473     /**
474      * \brief Overload of != operator
475      */
476     bool operator!=(const int64_t R) const;
477
478     /**
479      * \brief Overload of == operator
480      */
481     friend bool operator==(const int64_t L, const MultivariatePolynomial & R);
482
483     /**
484      * \brief Overload of != operator
485      */
486     friend bool operator!=(const int64_t L, const MultivariatePolynomial & R);
487
488     /**
489      * \return a hash
490      */
491     std::size_t hash() const;
492
493     /**
494      * \struct Hash
495      * \brief To be used in an unordered container
496      */
497     struct Hash
498     {
499         inline std::size_t operator()(const MultivariatePolynomial & mp) const
500         {
501             return mp.hash();
502         }
503     };
504
505     /**
506      * \struct Eq
507      * \brief To be used in an unordered container
508      */
509     struct Eq
510     {
511         inline bool operator()(const MultivariatePolynomial & L, const MultivariatePolynomial & R) const
512         {
513             return L == R;
514         }
515     };
516
517 private:
518
519     // Helper functions to use with eval
520     static bool __isValid(const std::unordered_map<uint64_t, const MultivariatePolynomial *> & values);
521     static bool __isValid(const std::vector<const MultivariatePolynomial *> & values);
522     static bool __isValid(const std::pair<uint64_t, const MultivariatePolynomial *> & values);
523     static bool __contains(const std::unordered_map<uint64_t, const MultivariatePolynomial *> & values, const uint64_t val);
524     static bool __contains(const std::vector<const MultivariatePolynomial *> & values, const uint64_t val);
525     static bool __contains(const std::pair<uint64_t, const MultivariatePolynomial *> & values, const uint64_t val);
526     static const MultivariatePolynomial * __get(const std::unordered_map<uint64_t, const MultivariatePolynomial *> & values, const uint64_t val);
527     static const MultivariatePolynomial * __getSafe(const std::unordered_map<uint64_t, const MultivariatePolynomial *> & values, const uint64_t val);
528     static const MultivariatePolynomial * __get(const std::vector<const MultivariatePolynomial *> & values, const uint64_t val);
529     static const MultivariatePolynomial * __getSafe(const std::vector<const MultivariatePolynomial *> & values, const uint64_t val);
530     static const MultivariatePolynomial * __get(const std::pair<uint64_t, const MultivariatePolynomial *> & values, const uint64_t val);
531     static const MultivariatePolynomial * __getSafe(const std::pair<uint64_t, const MultivariatePolynomial *> & values, const uint64_t val);
532 };
533
534 } // namespace analysis
535
536 #endif // __MULTIVARIATE_POLYNOMIAL_HXX__