pow.cpp
1 #include <symengine/pow.h>
2 #include <symengine/add.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 Pow::Pow(const RCP<const Basic> &base, const RCP<const Basic> &exp)
11  : base_{base}, exp_{exp}
12 {
13  SYMENGINE_ASSIGN_TYPEID()
14  SYMENGINE_ASSERT(is_canonical(*base, *exp))
15 }
16 
17 bool Pow::is_canonical(const Basic &base, const Basic &exp) const
18 {
19  // e.g. 0**x
20  if (is_a<Integer>(base) and down_cast<const Integer &>(base).is_zero()) {
21  if (is_a_Number(exp)) {
22  return false;
23  } else {
24  return true;
25  }
26  }
27  // e.g. 1**x
28  if (is_a<Integer>(base) and down_cast<const Integer &>(base).is_one())
29  return false;
30  // e.g. x**0.0
32  return false;
33  // e.g. x**1
34  if (is_a<Integer>(exp) and down_cast<const Integer &>(exp).is_one())
35  return false;
36  // e.g. 2**3, (2/3)**4
37  if ((is_a<Integer>(base) or is_a<Rational>(base)) and is_a<Integer>(exp))
38  return false;
39  // e.g. (x*y)**2, should rather be x**2*y**2
40  if (is_a<Mul>(base) and is_a<Integer>(exp))
41  return false;
42  // e.g. (x**y)**2, should rather be x**(2*y)
43  if (is_a<Pow>(base) and is_a<Integer>(exp))
44  return false;
45  // If exp is a rational, it should be between 0 and 1, i.e. we don't
46  // allow things like 2**(-1/2) or 2**(3/2)
47  if ((is_a<Rational>(base) or is_a<Integer>(base)) and is_a<Rational>(exp)
48  and (down_cast<const Rational &>(exp).as_rational_class() < 0
49  or down_cast<const Rational &>(exp).as_rational_class() > 1))
50  return false;
51  // Purely Imaginary complex numbers with integral powers are expanded
52  // e.g (2I)**3
53  if (is_a<Complex>(base) and down_cast<const Complex &>(base).is_re_zero()
54  and is_a<Integer>(exp))
55  return false;
56  // e.g. 0.5^2.0 should be represented as 0.25
57  if (is_a_Number(base) and not down_cast<const Number &>(base).is_exact()
58  and is_a_Number(exp) and not down_cast<const Number &>(exp).is_exact())
59  return false;
60  return true;
61 }
62 
63 hash_t Pow::__hash__() const
64 {
65  hash_t seed = SYMENGINE_POW;
66  hash_combine<Basic>(seed, *base_);
67  hash_combine<Basic>(seed, *exp_);
68  return seed;
69 }
70 
71 bool Pow::__eq__(const Basic &o) const
72 {
73  if (is_a<Pow>(o) and eq(*base_, *(down_cast<const Pow &>(o).base_))
74  and eq(*exp_, *(down_cast<const Pow &>(o).exp_)))
75  return true;
76 
77  return false;
78 }
79 
80 int Pow::compare(const Basic &o) const
81 {
82  SYMENGINE_ASSERT(is_a<Pow>(o))
83  const Pow &s = down_cast<const Pow &>(o);
84  int base_cmp = base_->__cmp__(*s.base_);
85  if (base_cmp == 0)
86  return exp_->__cmp__(*s.exp_);
87  else
88  return base_cmp;
89 }
90 
91 RCP<const Basic> pow(const RCP<const Basic> &a, const RCP<const Basic> &b)
92 {
93  if (is_number_and_zero(*b)) {
94  // addnum is used for converting to the type of `b`.
95  return addnum(one, rcp_static_cast<const Number>(b));
96  }
97  if (eq(*b, *one))
98  return a;
99 
100  if (eq(*a, *zero)) {
101  if (is_a_Number(*b)
102  and rcp_static_cast<const Number>(b)->is_positive()) {
103  return zero;
104  } else if (is_a_Number(*b)
105  and rcp_static_cast<const Number>(b)->is_negative()) {
106  return ComplexInf;
107  } else {
108  return make_rcp<const Pow>(a, b);
109  }
110  }
111 
112  if (eq(*a, *one) and not is_a_Number(*b))
113  return one;
114  if (eq(*a, *minus_one)) {
115  if (is_a<Integer>(*b)) {
116  return is_a<Integer>(*div(b, integer(2))) ? one : minus_one;
117  } else if (is_a<Rational>(*b) and eq(*b, *rational(1, 2))) {
118  return I;
119  }
120  }
121 
122  if (is_a_Number(*b)) {
123  if (is_a_Number(*a)) {
124  if (is_a<Integer>(*b)) {
125  return down_cast<const Number &>(*a).pow(
126  *rcp_static_cast<const Number>(b));
127  } else if (is_a<Rational>(*b)) {
128  if (is_a<Rational>(*a)) {
129  return down_cast<const Rational &>(*a).powrat(
130  down_cast<const Rational &>(*b));
131  } else if (is_a<Integer>(*a)) {
132  return down_cast<const Rational &>(*b).rpowrat(
133  down_cast<const Integer &>(*a));
134  } else if (is_a<Complex>(*a)) {
135  return make_rcp<const Pow>(a, b);
136  } else {
137  return down_cast<const Number &>(*a).pow(
138  *rcp_static_cast<const Number>(b));
139  }
140  } else if (is_a<Complex>(*b)
141  and down_cast<const Number &>(*a).is_exact()) {
142  return make_rcp<const Pow>(a, b);
143  } else {
144  return down_cast<const Number &>(*a).pow(
145  *rcp_static_cast<const Number>(b));
146  }
147  } else if (eq(*a, *E)) {
148  RCP<const Number> p = rcp_static_cast<const Number>(b);
149  if (not p->is_exact()) {
150  // Evaluate E**0.2, but not E**2
151  return p->get_eval().exp(*p);
152  }
153  } else if (is_a<Mul>(*a)) {
154  // Expand (x*y)**b = x**b*y**b
155  map_basic_basic d;
156  RCP<const Number> coef = one;
157  down_cast<const Mul &>(*a).power_num(
158  outArg(coef), d, rcp_static_cast<const Number>(b));
159  return Mul::from_dict(coef, std::move(d));
160  }
161  }
162  if (is_a<Pow>(*a) and is_a<Integer>(*b)) {
163  // Convert (x**y)**b = x**(b*y), where 'b' is an integer. This holds for
164  // any complex 'x', 'y' and integer 'b'.
165  RCP<const Pow> A = rcp_static_cast<const Pow>(a);
166  return pow(A->get_base(), mul(A->get_exp(), b));
167  }
168  if (is_a<Pow>(*a)
169  and eq(*down_cast<const Pow &>(*a).get_exp(), *minus_one)) {
170  // Convert (x**-1)**b = x**(-b)
171  RCP<const Pow> A = rcp_static_cast<const Pow>(a);
172  return pow(A->get_base(), neg(b));
173  }
174  return make_rcp<const Pow>(a, b);
175 }
176 
177 // This function can overflow, but it is fast.
178 // TODO: figure out condition for (m, n) when it overflows and raise an
179 // exception.
180 void multinomial_coefficients(unsigned m, unsigned n, map_vec_uint &r)
181 {
182  vec_uint t;
183  unsigned j, tj, start, k;
184  unsigned long long int v;
185  if (m < 2)
186  throw SymEngineException("multinomial_coefficients: m >= 2 must hold.");
187  t.assign(m, 0);
188  t[0] = n;
189  r[t] = 1;
190  if (n == 0)
191  return;
192  j = 0;
193  while (j < m - 1) {
194  tj = t[j];
195  if (j) {
196  t[j] = 0;
197  t[0] = tj;
198  }
199  if (tj > 1) {
200  t[j + 1] += 1;
201  j = 0;
202  start = 1;
203  v = 0;
204  } else {
205  j += 1;
206  start = j + 1;
207  v = r[t];
208  t[j] += 1;
209  }
210  for (k = start; k < m; k++) {
211  if (t[k]) {
212  t[k] -= 1;
213  v += r[t];
214  t[k] += 1;
215  }
216  }
217  t[0] -= 1;
218  r[t] = (v * tj) / (n - t[0]);
219  }
220 }
221 
222 // Slower, but returns exact (possibly large) integers (as mpz)
223 void multinomial_coefficients_mpz(unsigned m, unsigned n, map_vec_mpz &r)
224 {
225  vec_uint t;
226  unsigned j, tj, start, k;
227  integer_class v;
228  if (m < 2)
229  throw SymEngineException("multinomial_coefficients: m >= 2 must hold.");
230  t.assign(m, 0);
231  t[0] = n;
232  r[t] = 1;
233  if (n == 0)
234  return;
235  j = 0;
236  while (j < m - 1) {
237  tj = t[j];
238  if (j) {
239  t[j] = 0;
240  t[0] = tj;
241  }
242  if (tj > 1) {
243  t[j + 1] += 1;
244  j = 0;
245  start = 1;
246  v = 0;
247  } else {
248  j += 1;
249  start = j + 1;
250  v = r[t];
251  t[j] += 1;
252  }
253  for (k = start; k < m; k++) {
254  if (t[k]) {
255  t[k] -= 1;
256  v += r[t];
257  t[k] += 1;
258  }
259  }
260  t[0] -= 1;
261  r[t] = (v * tj) / (n - t[0]);
262  }
263 }
264 
266 {
267  return {base_, exp_};
268 }
269 
270 RCP<const Basic> exp(const RCP<const Basic> &x)
271 {
272  return pow(E, x);
273 }
274 
275 } // namespace SymEngine
Classes and functions relating to the binary operation of addition.
The lowest unit of symbolic representation.
Definition: basic.h:95
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 vec_basic get_args() const
Returns the list of arguments.
Definition: pow.cpp:265
Pow(const RCP< const Basic > &base, const RCP< const Basic > &exp)
Pow Constructor.
Definition: pow.cpp:10
virtual hash_t __hash__() const
Definition: pow.cpp:63
virtual bool __eq__(const Basic &o) const
Definition: pow.cpp:71
bool is_canonical(const Basic &base, const Basic &exp) const
Definition: pow.cpp:17
virtual int compare(const Basic &o) const
Definition: pow.cpp:80
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
std::enable_if< std::is_integral< T >::value, RCP< const Integer > >::type integer(T i)
Definition: integer.h:200
bool is_number_and_zero(const Basic &b)
Definition: number.h:139
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
tribool is_zero(const Basic &b, const Assumptions *assumptions=nullptr)
Check if a number is zero.
RCP< const Number > rational(long n, long d)
convenience creator from two longs
Definition: rational.h:328
RCP< const Number > addnum(const RCP< const Number > &self, const RCP< const Number > &other)
Add self and other
Definition: number.h:81
RCP< const Basic > neg(const RCP< const Basic > &a)
Negation.
Definition: mul.cpp:438