mul.cpp
1 #include <symengine/add.h>
2 #include <symengine/pow.h>
3 #include <symengine/complex.h>
4 #include <symengine/symengine_exception.h>
5 #include <symengine/test_visitors.h>
6 
7 namespace SymEngine
8 {
9 
10 Mul::Mul(const RCP<const Number> &coef, map_basic_basic &&dict)
11  : coef_{coef}, dict_{std::move(dict)}
12 {
13  SYMENGINE_ASSIGN_TYPEID()
14  SYMENGINE_ASSERT(is_canonical(coef, dict_))
15 }
16 
17 bool Mul::is_canonical(const RCP<const Number> &coef,
18  const map_basic_basic &dict) const
19 {
20  if (coef == null)
21  return false;
22  // e.g. 0*x*y
23  if (coef->is_zero())
24  return false;
25  if (dict.size() == 0)
26  return false;
27  if (dict.size() == 1) {
28  // e.g. 1*x, 1*x**2
29  if (coef->is_one())
30  return false;
31  }
32  // Check that each term in 'dict' is in canonical form
33  for (const auto &p : dict) {
34  if (p.first == null)
35  return false;
36  if (p.second == null)
37  return false;
38  // e.g. 2**3, (2/3)**4
39  // However for Complex no simplification is done
40  if ((is_a<Integer>(*p.first) or is_a<Rational>(*p.first))
41  and is_a<Integer>(*p.second))
42  return false;
43  // e.g. 0**x
44  if (is_a<Integer>(*p.first)
45  and down_cast<const Integer &>(*p.first).is_zero())
46  return false;
47  // e.g. 1**x
48  if (is_a<Integer>(*p.first)
49  and down_cast<const Integer &>(*p.first).is_one())
50  return false;
51  // e.g. x**0
52  if (is_a_Number(*p.second)
53  and down_cast<const Number &>(*p.second).is_zero())
54  return false;
55  // e.g. (x*y)**2 (={xy:2}), which should be represented as x**2*y**2
56  // (={x:2, y:2})
57  if (is_a<Mul>(*p.first)) {
58  if (is_a<Integer>(*p.second))
59  return false;
60  if (is_a_Number(*p.second)
61  and neq(*down_cast<const Mul &>(*p.first).coef_, *one)
62  and neq(*down_cast<const Mul &>(*p.first).coef_, *minus_one))
63  return false;
64  }
65  // e.g. x**2**y (={x**2:y}), which should be represented as x**(2y)
66  // (={x:2y})
67  if (is_a<Pow>(*p.first) && is_a<Integer>(*p.second))
68  return false;
69  // e.g. 0.5^2.0 should be represented as 0.25
70  if (is_a_Number(*p.first)
71  and not down_cast<const Number &>(*p.first).is_exact()
72  and is_a_Number(*p.second)
73  and not down_cast<const Number &>(*p.second).is_exact())
74  return false;
75  }
76  return true;
77 }
78 
79 hash_t Mul::__hash__() const
80 {
81  hash_t seed = SYMENGINE_MUL;
82  hash_combine<Basic>(seed, *coef_);
83  for (const auto &p : dict_) {
84  hash_combine<Basic>(seed, *(p.first));
85  hash_combine<Basic>(seed, *(p.second));
86  }
87  return seed;
88 }
89 
90 bool Mul::__eq__(const Basic &o) const
91 {
92  if (is_a<Mul>(o) and eq(*coef_, *(down_cast<const Mul &>(o).coef_))
93  and unified_eq(dict_, down_cast<const Mul &>(o).dict_))
94  return true;
95 
96  return false;
97 }
98 
99 int Mul::compare(const Basic &o) const
100 {
101  SYMENGINE_ASSERT(is_a<Mul>(o))
102  const Mul &s = down_cast<const Mul &>(o);
103  // # of elements
104  if (dict_.size() != s.dict_.size())
105  return (dict_.size() < s.dict_.size()) ? -1 : 1;
106 
107  // coef
108  int cmp = coef_->__cmp__(*s.coef_);
109  if (cmp != 0)
110  return cmp;
111 
112  // Compare dictionaries:
113  return unified_compare(dict_, s.dict_);
114 }
115 
116 RCP<const SymEngine::Basic> Mul::from_dict(const RCP<const Number> &coef,
117  map_basic_basic &&d)
118 {
119  if (coef->is_zero())
120  return coef;
121  if (d.size() == 0) {
122  return coef;
123  } else if (d.size() == 1) {
124  auto p = d.begin();
125  if (is_a<Integer>(*(p->second))) {
126  if (coef->is_one()) {
127  if ((down_cast<const Integer &>(*(p->second))).is_one()) {
128  // For x**1 we simply return "x":
129  return p->first;
130  }
131  } else {
132  // For coef*x or coef*x**3 we simply return Mul:
133  return make_rcp<const Mul>(coef, std::move(d));
134  }
135  }
136  if (coef->is_one()) {
137  // Create a Pow() here:
138  if (eq(*p->second, *one)) {
139  return p->first;
140  }
141  return make_rcp<const Pow>(p->first, p->second);
142  } else {
143  return make_rcp<const Mul>(coef, std::move(d));
144  }
145  } else {
146  return make_rcp<const Mul>(coef, std::move(d));
147  }
148 }
149 
150 // Mul (t**exp) to the dict "d"
151 void Mul::dict_add_term(map_basic_basic &d, const RCP<const Basic> &exp,
152  const RCP<const Basic> &t)
153 {
154  auto it = d.find(t);
155  if (it == d.end()) {
156  insert(d, t, exp);
157  } else {
158  // Very common case, needs to be fast:
159  if (is_a_Number(*it->second) and is_a_Number(*exp)) {
160  RCP<const Number> tmp = rcp_static_cast<const Number>(it->second);
161  iaddnum(outArg(tmp), rcp_static_cast<const Number>(exp));
162  if (tmp->is_zero()) {
163  d.erase(it);
164  } else {
165  it->second = tmp;
166  }
167  } else {
168  // General case:
169  it->second = add(it->second, exp);
170  if (is_number_and_zero(*it->second)) {
171  d.erase(it);
172  }
173  }
174  }
175 }
176 
177 // Mul (t**exp) to the dict "d"
178 void Mul::dict_add_term_new(const Ptr<RCP<const Number>> &coef,
179  map_basic_basic &d, const RCP<const Basic> &exp,
180  const RCP<const Basic> &t)
181 {
182  auto it = d.find(t);
183  if (it == d.end()) {
184  // Don't check for `exp = 0` here
185  // `pow` for Complex is not expanded by default
186  if (is_a<Integer>(*t) or is_a<Rational>(*t)) {
187  if (is_a<Integer>(*exp)) {
188  imulnum(outArg(*coef),
189  pownum(rcp_static_cast<const Number>(t),
190  rcp_static_cast<const Number>(exp)));
191  } else if (is_a<Rational>(*exp)) {
192  RCP<const Basic> res;
193  if (is_a<Integer>(*t)) {
194  res = down_cast<const Rational &>(*exp).rpowrat(
195  down_cast<const Integer &>(*t));
196  } else {
197  res = down_cast<const Rational &>(*t).powrat(
198  down_cast<const Rational &>(*exp));
199  }
200  if (is_a_Number(*res)) {
201  imulnum(outArg(*coef), rcp_static_cast<const Number>(res));
202  } else if (is_a<Mul>(*res)) {
203  RCP<const Mul> m = rcp_static_cast<const Mul>(res);
204  imulnum(outArg(*coef), m->coef_);
205  for (auto &p : m->dict_) {
206  Mul::dict_add_term_new(coef, d, p.second, p.first);
207  }
208  } else {
209  insert(d, t, exp);
210  }
211  } else {
212  insert(d, t, exp);
213  }
214  } else if (is_a<Integer>(*exp) and is_a<Complex>(*t)) {
215  if (down_cast<const Integer &>(*exp).is_one()) {
216  imulnum(outArg(*coef), rcp_static_cast<const Number>(t));
217  } else if (down_cast<const Integer &>(*exp).is_minus_one()) {
218  idivnum(outArg(*coef), rcp_static_cast<const Number>(t));
219  } else {
220  insert(d, t, exp);
221  }
222  } else {
223  insert(d, t, exp);
224  }
225  } else {
226  // Very common case, needs to be fast:
227  if (is_a_Number(*exp) and is_a_Number(*it->second)) {
228  RCP<const Number> tmp = rcp_static_cast<const Number>(it->second);
229  iaddnum(outArg(tmp), rcp_static_cast<const Number>(exp));
230  it->second = tmp;
231  } else
232  it->second = add(it->second, exp);
233 
234  if (is_a<Integer>(*it->second)) {
235  // `pow` for Complex is not expanded by default
236  if (is_a<Integer>(*t) or is_a<Rational>(*t)) {
237  if (not down_cast<const Integer &>(*(it->second)).is_zero()) {
238  imulnum(outArg(*coef),
239  pownum(rcp_static_cast<const Number>(t),
240  rcp_static_cast<const Number>(it->second)));
241  }
242  d.erase(it);
243  return;
244  } else if (down_cast<const Integer &>(*(it->second)).is_zero()) {
245  d.erase(it);
246  return;
247  } else if (is_a<Complex>(*t)) {
248  if (down_cast<const Integer &>(*(it->second)).is_one()) {
249  imulnum(outArg(*coef), rcp_static_cast<const Number>(t));
250  d.erase(it);
251  } else if (down_cast<const Integer &>(*(it->second))
252  .is_minus_one()) {
253  idivnum(outArg(*coef), rcp_static_cast<const Number>(t));
254  d.erase(it);
255  }
256  return;
257  }
258  } else if (is_a<Rational>(*it->second)) {
259  if (is_a<Integer>(*t) or is_a<Rational>(*t)) {
260  RCP<const Basic> res;
261  if (is_a<Integer>(*t)) {
262  res = down_cast<const Rational &>(*it->second)
263  .rpowrat(down_cast<const Integer &>(*t));
264  } else {
265  res = down_cast<const Rational &>(*t).powrat(
266  down_cast<const Rational &>(*it->second));
267  }
268  if (is_a_Number(*res)) {
269  d.erase(it);
270  imulnum(outArg(*coef), rcp_static_cast<const Number>(res));
271  return;
272  } else if (is_a<Mul>(*res)) {
273  d.erase(it);
274  RCP<const Mul> m = rcp_static_cast<const Mul>(res);
275  imulnum(outArg(*coef), m->coef_);
276  for (auto &p : m->dict_) {
277  Mul::dict_add_term_new(coef, d, p.second, p.first);
278  }
279  return;
280  }
281  }
282  }
283  if (is_a_Number(*it->second)) {
284  if (down_cast<const Number &>(*it->second).is_zero()) {
285  // In 1*x**0.0, result should be 1.0
286  imulnum(
287  outArg(*coef),
288  pownum(rcp_static_cast<const Number>(it->second), zero));
289  d.erase(it);
290  } else if (is_a<Mul>(*it->first)) {
291  RCP<const Mul> m = rcp_static_cast<const Mul>(it->first);
292  if (is_a<Integer>(*it->second)
293  or (neq(*m->coef_, *one) and neq(*m->coef_, *minus_one))) {
294  RCP<const Number> exp_
295  = rcp_static_cast<const Number>(it->second);
296  d.erase(it);
297  m->power_num(outArg(*coef), d, exp_);
298  }
299  }
300  }
301  }
302 }
303 
304 void Mul::as_two_terms(const Ptr<RCP<const Basic>> &a,
305  const Ptr<RCP<const Basic>> &b) const
306 {
307  // Example: if this=3*x**2*y**2*z**2, then a=x**2 and b=3*y**2*z**2
308  auto p = dict_.begin();
309  *a = pow(p->first, p->second);
311  d.erase(p->first);
312  *b = Mul::from_dict(coef_, std::move(d));
313 }
314 
315 void Mul::as_base_exp(const RCP<const Basic> &self,
316  const Ptr<RCP<const Basic>> &exp,
317  const Ptr<RCP<const Basic>> &base)
318 {
319  if (is_a_Number(*self)) {
320  // Always ensure it is of form |num| > |den|
321  // in case of Integers den = 1
322  if (is_a<Rational>(*self)) {
323  RCP<const Rational> self_new
324  = rcp_static_cast<const Rational>(self);
325  if (mp_abs(get_num(self_new->as_rational_class()))
326  < mp_abs(get_den(self_new->as_rational_class()))) {
327  *exp = minus_one;
328  *base = self_new->rdiv(*rcp_static_cast<const Number>(one));
329  } else {
330  *exp = one;
331  *base = self;
332  }
333  } else {
334  *exp = one;
335  *base = self;
336  }
337  } else if (is_a<Pow>(*self)) {
338  *exp = down_cast<const Pow &>(*self).get_exp();
339  *base = down_cast<const Pow &>(*self).get_base();
340  } else {
341  SYMENGINE_ASSERT(not is_a<Mul>(*self));
342  *exp = one;
343  *base = self;
344  }
345 }
346 
347 RCP<const Basic> mul(const RCP<const Basic> &a, const RCP<const Basic> &b)
348 {
350  RCP<const Number> coef = one;
351  if (is_a<Mul>(*a) and is_a<Mul>(*b)) {
352  RCP<const Mul> A = rcp_static_cast<const Mul>(a);
353  RCP<const Mul> B = rcp_static_cast<const Mul>(b);
354  // This is important optimization, as coef=1 if Mul is inside an Add.
355  // To further speed this up, the upper level code could tell us that we
356  // are inside an Add, then we don't even have can simply skip the
357  // following two lines.
358  if (not(A->get_coef()->is_one()) or not(B->get_coef()->is_one()))
359  coef = mulnum(A->get_coef(), B->get_coef());
360  d = A->get_dict();
361  for (const auto &p : B->get_dict())
362  Mul::dict_add_term_new(outArg(coef), d, p.second, p.first);
363  } else if (is_a<Mul>(*a)) {
364  RCP<const Basic> exp;
365  RCP<const Basic> t;
366  coef = (down_cast<const Mul &>(*a)).get_coef();
367  d = (down_cast<const Mul &>(*a)).get_dict();
368  if (is_a_Number(*b)) {
369  imulnum(outArg(coef), rcp_static_cast<const Number>(b));
370  } else {
371  Mul::as_base_exp(b, outArg(exp), outArg(t));
372  Mul::dict_add_term_new(outArg(coef), d, exp, t);
373  }
374  } else if (is_a<Mul>(*b)) {
375  RCP<const Basic> exp;
376  RCP<const Basic> t;
377  coef = (down_cast<const Mul &>(*b)).get_coef();
378  d = (down_cast<const Mul &>(*b)).get_dict();
379  if (is_a_Number(*a)) {
380  imulnum(outArg(coef), rcp_static_cast<const Number>(a));
381  } else {
382  Mul::as_base_exp(a, outArg(exp), outArg(t));
383  Mul::dict_add_term_new(outArg(coef), d, exp, t);
384  }
385  } else {
386  RCP<const Basic> exp;
387  RCP<const Basic> t;
388  if (is_a_Number(*a)) {
389  imulnum(outArg(coef), rcp_static_cast<const Number>(a));
390  } else {
391  Mul::as_base_exp(a, outArg(exp), outArg(t));
392  Mul::dict_add_term_new(outArg(coef), d, exp, t);
393  }
394  if (is_a_Number(*b)) {
395  imulnum(outArg(coef), rcp_static_cast<const Number>(b));
396  } else {
397  Mul::as_base_exp(b, outArg(exp), outArg(t));
398  Mul::dict_add_term_new(outArg(coef), d, exp, t);
399  }
400  }
401  return Mul::from_dict(coef, std::move(d));
402 }
403 
404 RCP<const Basic> mul(const vec_basic &a)
405 {
407  RCP<const Number> coef = one;
408  for (const auto &i : a) {
409  if (is_a<Mul>(*i)) {
410  RCP<const Mul> A = rcp_static_cast<const Mul>(i);
411  imulnum(outArg(coef), A->get_coef());
412  for (const auto &p : A->get_dict())
413  Mul::dict_add_term_new(outArg(coef), d, p.second, p.first);
414  } else if (is_a_Number(*i)) {
415  imulnum(outArg(coef), rcp_static_cast<const Number>(i));
416  } else {
417  RCP<const Basic> exp;
418  RCP<const Basic> t;
419  Mul::as_base_exp(i, outArg(exp), outArg(t));
420  Mul::dict_add_term_new(outArg(coef), d, exp, t);
421  }
422  }
423  return Mul::from_dict(coef, std::move(d));
424 }
425 
426 RCP<const Basic> div(const RCP<const Basic> &a, const RCP<const Basic> &b)
427 {
428  if (is_number_and_zero(*b)) {
429  if (is_number_and_zero(*a)) {
430  return Nan;
431  } else {
432  return ComplexInf;
433  }
434  }
435  return mul(a, pow(b, minus_one));
436 }
437 
438 RCP<const Basic> neg(const RCP<const Basic> &a)
439 {
440  return mul(minus_one, a);
441 }
442 
443 void Mul::power_num(const Ptr<RCP<const Number>> &coef, map_basic_basic &d,
444  const RCP<const Number> &exp) const
445 {
446  if (exp->is_zero()) {
447  // (x*y)**(0.0) should return 1.0
448  imulnum(coef, pownum(rcp_static_cast<const Number>(exp), zero));
449  return;
450  }
451  RCP<const Basic> new_coef;
452  RCP<const Basic> new_exp;
453  if (is_a<Integer>(*exp)) {
454  // For eg. (3*y*(x**(1/2))**2 should be expanded to 9*x*y**2
455  new_coef = pow(coef_, exp);
456  for (const auto &p : dict_) {
457  new_exp = mul(p.second, exp);
458  if (is_a<Integer>(*new_exp) and is_a<Mul>(*p.first)) {
459  down_cast<const Mul &>(*p.first).power_num(
460  coef, d, rcp_static_cast<const Number>(new_exp));
461  } else {
462  // No need for additional dict checks here.
463  // The dict should be of standard form before this is
464  // called.
465  Mul::dict_add_term_new(coef, d, new_exp, p.first);
466  }
467  }
468  } else {
469  if (coef_->is_negative()) {
470  // (-3*x*y)**(1/2) -> 3**(1/2)*(-x*y)**(1/2)
471  new_coef = pow(coef_->mul(*minus_one), exp);
472  map_basic_basic d1 = dict_;
473  Mul::dict_add_term_new(coef, d, exp,
474  Mul::from_dict(minus_one, std::move(d1)));
475  } else if (coef_->is_positive() and not coef_->is_one()) {
476  // (3*x*y)**(1/2) -> 3**(1/2)*(x*y)**(1/2)
477  new_coef = pow(coef_, exp);
478  map_basic_basic d1 = dict_;
479  Mul::dict_add_term_new(coef, d, exp,
480  Mul::from_dict(one, std::move(d1)));
481  } else {
482  // ((1+2*I)*x*y)**(1/2) is kept as it is
483  new_coef = one;
484  Mul::dict_add_term_new(coef, d, exp, this->rcp_from_this());
485  }
486  }
487  if (is_a_Number(*new_coef)) {
488  imulnum(coef, rcp_static_cast<const Number>(new_coef));
489  } else if (is_a<Mul>(*new_coef)) {
490  RCP<const Mul> tmp = rcp_static_cast<const Mul>(new_coef);
491  imulnum(coef, tmp->coef_);
492  for (const auto &q : tmp->dict_) {
493  Mul::dict_add_term_new(coef, d, q.second, q.first);
494  }
495  } else {
496  RCP<const Basic> _exp, t;
497  Mul::as_base_exp(new_coef, outArg(_exp), outArg(t));
498  Mul::dict_add_term_new(coef, d, _exp, t);
499  }
500 }
501 
503 {
504  vec_basic args;
505  if (not coef_->is_one()) {
506  args.reserve(dict_.size() + 1);
507  args.push_back(coef_);
508  } else {
509  args.reserve(dict_.size());
510  }
511  for (const auto &p : dict_) {
512  if (eq(*p.second, *one)) {
513  args.push_back(p.first);
514  } else {
515  args.push_back(make_rcp<const Pow>(p.first, p.second));
516  }
517  }
518  return args;
519 }
520 
521 } // namespace SymEngine
Classes and functions relating to the binary operation of addition.
T begin(T... args)
The lowest unit of symbolic representation.
Definition: basic.h:95
RCP< Basic > rcp_from_this()
Get RCP<T> pointer to self (it will cast the pointer to T)
bool is_canonical(const RCP< const Number > &coef, const map_basic_basic &dict) const
Definition: mul.cpp:17
Mul(const RCP< const Number > &coef, map_basic_basic &&dict)
Definition: mul.cpp:10
static void as_base_exp(const RCP< const Basic > &self, const Ptr< RCP< const Basic >> &exp, const Ptr< RCP< const Basic >> &base)
Convert to a base and exponent form.
Definition: mul.cpp:315
static void dict_add_term(map_basic_basic &d, const RCP< const Basic > &exp, const RCP< const Basic > &t)
Add terms to dict.
Definition: mul.cpp:151
virtual vec_basic get_args() const
Returns the list of arguments.
Definition: mul.cpp:502
virtual bool __eq__(const Basic &o) const
Definition: mul.cpp:90
void as_two_terms(const Ptr< RCP< const Basic >> &a, const Ptr< RCP< const Basic >> &b) const
Rewrite as 2 terms.
Definition: mul.cpp:304
static RCP< const Basic > from_dict(const RCP< const Number > &coef, map_basic_basic &&d)
Create a Mul from a dict.
Definition: mul.cpp:116
virtual hash_t __hash__() const
Definition: mul.cpp:79
void power_num(const Ptr< RCP< const Number >> &coef, map_basic_basic &d, const RCP< const Number > &exp) const
Power all terms with the exponent exp
Definition: mul.cpp:443
map_basic_basic dict_
The coefficient (e.g. 2 in 2*x*y)
Definition: mul.h:78
virtual int compare(const Basic &o) const
Definition: mul.cpp:99
T end(T... args)
T erase(T... args)
T find(T... args)
T move(T... args)
Main namespace for SymEngine package.
Definition: add.cpp:19
bool is_a_Number(const Basic &b)
Definition: number.h:130
RCP< const Basic > div(const RCP< const Basic > &a, const RCP< const Basic > &b)
Division.
Definition: mul.cpp:426
bool is_number_and_zero(const Basic &b)
Definition: number.h:139
RCP< const Number > mulnum(const RCP< const Number > &self, const RCP< const Number > &other)
Multiply self and other
Definition: number.h:93
RCP< const Number > pownum(const RCP< const Number > &self, const RCP< const Number > &other)
Raise self to power other
Definition: number.h:105
bool eq(const Basic &a, const Basic &b)
Checks equality for a and b
Definition: basic-inl.h:21
RCP< const Basic > exp(const RCP< const Basic > &x)
Returns the natural exponential function E**x = pow(E, x)
Definition: pow.cpp:270
RCP< const Basic > mul(const RCP< const Basic > &a, const RCP< const Basic > &b)
Multiplication.
Definition: mul.cpp:347
void insert(T1 &m, const T2 &first, const T3 &second)
Definition: dict.h:83
RCP< const Basic > add(const RCP< const Basic > &a, const RCP< const Basic > &b)
Adds two objects (safely).
Definition: add.cpp:425
int unified_compare(const T &a, const T &b)
Definition: dict.h:205
bool neq(const Basic &a, const Basic &b)
Checks inequality for a and b
Definition: basic-inl.h:29
RCP< const Basic > neg(const RCP< const Basic > &a)
Negation.
Definition: mul.cpp:438
STL namespace.
T push_back(T... args)
T reserve(T... args)
T size(T... args)