as_real_imag.cpp
1 #include <symengine/visitor.h>
2 #include <symengine/basic.h>
3 
4 namespace SymEngine
5 {
6 
7 void pow_number(const RCP<const Basic> &in_re, const RCP<const Basic> &in_im,
8  unsigned long n, Ptr<RCP<const Basic>> &out_re,
9  Ptr<RCP<const Basic>> &out_im)
10 {
11  unsigned long mask = 1;
12  RCP<const Basic> tmp;
13  RCP<const Basic> p_re = in_re;
14  RCP<const Basic> p_im = in_im;
15  *out_re = one;
16  *out_im = zero;
17 
18  while (true) {
19  if (n & mask) {
20  // Multiply r by p
21  tmp = sub(mul(*out_re, p_re), mul(*out_im, p_im));
22  *out_im = add(mul(*out_re, p_im), mul(*out_im, p_re));
23  *out_re = tmp;
24  }
25  mask = mask << 1;
26  if (mask > 0 and n >= mask) {
27  // Multiply p by p
28  tmp = sub(mul(p_re, p_re), mul(p_im, p_im));
29  p_im = mul(two, mul(p_re, p_im));
30  p_re = tmp;
31  } else {
32  break;
33  }
34  }
35 }
36 
37 class RealImagVisitor : public BaseVisitor<RealImagVisitor>
38 {
39 private:
40  Ptr<RCP<const Basic>> real_, imag_;
41 
42 public:
43  RealImagVisitor(const Ptr<RCP<const Basic>> &real,
44  const Ptr<RCP<const Basic>> &imag)
45  : real_{real}, imag_{imag}
46  {
47  }
48 
49  void apply(const Basic &b)
50  {
51  b.accept(*this);
52  }
53 
54  void bvisit(const Mul &x)
55  {
56  RCP<const Basic> rest = one;
57  RCP<const Basic> fre_ = one, fim_ = zero;
58 
59  for (const auto &arg : x.get_args()) {
60  apply(*arg);
61  std::tie(fre_, fim_)
62  = std::make_tuple(sub(mul(fre_, *real_), mul(fim_, *imag_)),
63  add(mul(fre_, *imag_), mul(fim_, *real_)));
64  }
65  *real_ = fre_;
66  *imag_ = fim_;
67  }
68 
69  void bvisit(const Add &x)
70  {
71  RCP<const Basic> t;
72  umap_basic_num dr, dim;
73  RCP<const Number> coefr = zero, coefim = zero, coef;
74 
75  for (const auto &arg : x.get_args()) {
76  apply(*arg);
77  if (is_a_Number(**real_)) {
78  iaddnum(outArg(coefr), rcp_static_cast<const Number>(*real_));
79  } else {
80  Add::as_coef_term(*real_, outArg(coef), outArg(t));
81  Add::dict_add_term(dr, coef, t);
82  }
83  if (is_a_Number(**imag_)) {
84  iaddnum(outArg(coefim), rcp_static_cast<const Number>(*imag_));
85  } else {
86  Add::as_coef_term(*imag_, outArg(coef), outArg(t));
87  Add::dict_add_term(dim, coef, t);
88  }
89  }
90 
91  *real_ = Add::from_dict(coefr, std::move(dr));
92  *imag_ = Add::from_dict(coefim, std::move(dim));
93  }
94 
95  void bvisit(const Pow &x)
96  {
97  RCP<const Basic> exp_;
98  exp_ = x.get_exp();
99  apply(*x.get_base());
100 
101  if (eq(**imag_, *zero)) {
102  *real_ = x.rcp_from_this();
103  *imag_ = zero;
104  return;
105  }
106  if (is_a<Integer>(*exp_)) {
107  if (static_cast<const Integer &>(*exp_).is_negative()) {
108  auto magn = add(mul(*real_, *real_), mul(*imag_, *imag_));
109  *imag_ = neg(*imag_);
110  RCP<const Integer> expx = rcp_static_cast<const Integer>(exp_);
111  expx = static_cast<const Integer &>(*exp_).neg();
112  unsigned long n = numeric_cast<unsigned long>(
113  mp_get_ui(expx->as_integer_class()));
114  RCP<const Basic> real1 = *real_, imag1 = *imag_;
115  pow_number(real1, imag1, n, real_, imag_);
116  magn = pow(magn, expx);
117  *real_ = div(*real_, magn);
118  *imag_ = div(*imag_, magn);
119  } else {
120  RCP<const Integer> expx = rcp_static_cast<const Integer>(exp_);
121  unsigned long n = numeric_cast<unsigned long>(
122  mp_get_ui(expx->as_integer_class()));
123  RCP<const Basic> real1 = *real_, imag1 = *imag_;
124  pow_number(real1, imag1, n, real_, imag_);
125  }
126  } else if (is_a<Rational>(*exp_)) {
127  auto magn = sqrt(add(mul(*real_, *real_), mul(*imag_, *imag_)));
128  auto ang = atan2(*imag_, *real_);
129  magn = pow(magn, exp_);
130  ang = mul(ang, exp_);
131  *real_ = mul(magn, cos(ang));
132  *imag_ = mul(magn, sin(ang));
133  } else {
134  throw SymEngineException("Not Implemented");
135  }
136  }
137 
138  void bvisit(const ComplexBase &x)
139  {
140  *real_ = x.real_part();
141  *imag_ = x.imaginary_part();
142  }
143 
144  void bvisit(const Infty &x)
145  {
146  if (eq(x, *ComplexInf)) {
147  *real_ = Nan;
148  *imag_ = Nan;
149  } else {
150  *real_ = x.rcp_from_this();
151  *imag_ = zero;
152  }
153  }
154 
155  void bvisit(const Sin &x)
156  {
157  apply(*x.get_arg());
158  std::tie(*real_, *imag_) = std::make_tuple(
159  mul(sin(*real_), cosh(*imag_)), mul(sinh(*imag_), cos(*real_)));
160  }
161 
162  void bvisit(const Cos &x)
163  {
164  apply(*x.get_arg());
165  std::tie(*real_, *imag_)
166  = std::make_tuple(mul(cos(*real_), cosh(*imag_)),
167  neg(mul(sinh(*imag_), sin(*real_))));
168  }
169 
170  void bvisit(const Tan &x)
171  {
172  apply(*x.get_arg());
173  if (eq(**imag_, *zero)) {
174  *real_ = x.rcp_from_this();
175  *imag_ = zero;
176  return;
177  }
178  auto twice_real_ = mul(two, *real_), twice_imag_ = mul(two, *imag_);
179  auto den = add(cos(twice_real_), cosh(twice_imag_));
180  *real_ = div(sin(twice_real_), den);
181  *imag_ = div(sinh(twice_imag_), den);
182  }
183 
184  void bvisit(const Csc &x)
185  {
186  apply(*div(one, sin(x.get_arg())));
187  }
188 
189  void bvisit(const Sec &x)
190  {
191  apply(*div(one, cos(x.get_arg())));
192  }
193 
194  void bvisit(const Cot &x)
195  {
196  apply(*x.get_arg());
197  if (eq(**imag_, *zero)) {
198  *real_ = x.rcp_from_this();
199  return;
200  }
201  auto twice_real_ = mul(two, *real_), twice_imag_ = mul(two, *imag_);
202  auto den = sub(cos(twice_real_), cosh(twice_imag_));
203  *real_ = neg(div(sin(twice_real_), den));
204  *imag_ = neg(div(sinh(twice_imag_), den));
205  }
206 
207  void bvisit(const Sinh &x)
208  {
209  apply(*x.get_arg());
210  std::tie(*real_, *imag_) = std::make_tuple(
211  mul(sinh(*real_), cos(*imag_)), mul(sin(*imag_), cosh(*real_)));
212  }
213 
214  void bvisit(const Cosh &x)
215  {
216  apply(*x.get_arg());
217  std::tie(*real_, *imag_) = std::make_tuple(
218  mul(cosh(*real_), cos(*imag_)), mul(sin(*imag_), sinh(*real_)));
219  }
220 
221  void bvisit(const Tanh &x)
222  {
223  apply(*x.get_arg());
224  if (eq(**imag_, *zero)) {
225  *real_ = x.rcp_from_this();
226  return;
227  }
228  auto sinh_re = sinh(*real_), cos_im = cos(*imag_);
229  auto den = add(pow(sinh_re, two), pow(cos_im, two));
230  *real_ = div(mul(sinh_re, cosh(*real_)), den);
231  *imag_ = div(mul(sin(*imag_), cos_im), den);
232  }
233 
234  void bvisit(const Csch &x)
235  {
236  apply(*div(one, sinh(x.get_arg())));
237  }
238 
239  void bvisit(const Sech &x)
240  {
241  apply(*div(one, cosh(x.get_arg())));
242  }
243 
244  void bvisit(const Coth &x)
245  {
246  apply(*x.get_arg());
247  if (eq(**imag_, *zero)) {
248  *real_ = x.rcp_from_this();
249  return;
250  }
251  auto sinh_re = sinh(*real_), sin_im = sin(*imag_);
252  auto den = add(pow(sinh_re, two), pow(sin_im, two));
253  *real_ = div(mul(sinh_re, cosh(*real_)), den);
254  *imag_ = neg(div(mul(sin_im, cos(*imag_)), den));
255  }
256 
257  void bvisit(const Abs &x)
258  {
259  *real_ = x.rcp_from_this();
260  *imag_ = zero;
261  }
262 
263  void bvisit(const Function &x)
264  {
265  throw SymEngineException(
266  "Not Implemented classes for real and imaginary parts");
267  }
268 
269  void bvisit(const Symbol &x)
270  {
271  throw SymEngineException(
272  "Not Implemented classes for real and imaginary parts");
273  }
274 
275  void bvisit(const Basic &x)
276  {
277  *real_ = x.rcp_from_this();
278  *imag_ = zero;
279  }
280 };
281 
282 void as_real_imag(const RCP<const Basic> &x, const Ptr<RCP<const Basic>> &real,
283  const Ptr<RCP<const Basic>> &imag)
284 {
285  RealImagVisitor v(real, imag);
286  v.apply(*x);
287 }
288 
289 } // namespace SymEngine
The base class for SymEngine.
The base class for representing addition in symbolic expressions.
Definition: add.h:27
static RCP< const Basic > from_dict(const RCP< const Number > &coef, umap_basic_num &&d)
Create an appropriate instance from dictionary quickly.
Definition: add.cpp:140
static void dict_add_term(umap_basic_num &d, const RCP< const Number > &coef, const RCP< const Basic > &t)
Adds a new term to the expression.
Definition: add.cpp:237
virtual vec_basic get_args() const
Returns the arguments of the Add.
Definition: add.cpp:397
static void as_coef_term(const RCP< const Basic > &self, const Ptr< RCP< const Number >> &coef, const Ptr< RCP< const Basic >> &term)
Converts a Basic self into the form of coefficient * term.
Definition: add.cpp:306
The lowest unit of symbolic representation.
Definition: basic.h:95
ComplexBase Class for deriving all complex classes.
Definition: complex.h:16
RCP< T > rcp_from_this()
Get RCP<T> pointer to self (it will cast the pointer to T)
Integer Class.
Definition: integer.h:19
RCP< const Integer > neg() const
Definition: integer.h:119
virtual bool is_negative() const
Definition: integer.h:73
virtual vec_basic get_args() const
Returns the list of arguments.
Definition: mul.cpp:502
RCP< const Basic > get_arg() const
Definition: functions.h:36
RCP< const Basic > get_exp() const
Definition: pow.h:42
RCP< const Basic > get_base() const
Definition: pow.h:37
T make_tuple(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 eq(const Basic &a, const Basic &b)
Checks equality for a and b
Definition: basic-inl.h:21
RCP< const Basic > atan2(const RCP< const Basic > &num, const RCP< const Basic > &den)
Canonicalize ATan2:
Definition: functions.cpp:1571
RCP< const Basic > sub(const RCP< const Basic > &a, const RCP< const Basic > &b)
Substracts b from a.
Definition: add.cpp:495
RCP< const Basic > cosh(const RCP< const Basic > &arg)
Canonicalize Cosh:
Definition: functions.cpp:2169
RCP< const Basic > mul(const RCP< const Basic > &a, const RCP< const Basic > &b)
Multiplication.
Definition: mul.cpp:347
RCP< const Basic > cos(const RCP< const Basic > &arg)
Canonicalize Cos:
Definition: functions.cpp:899
RCP< const Basic > add(const RCP< const Basic > &a, const RCP< const Basic > &b)
Adds two objects (safely).
Definition: add.cpp:425
RCP< const Basic > sinh(const RCP< const Basic > &arg)
Canonicalize Sinh:
Definition: functions.cpp:2084
RCP< const Basic > neg(const RCP< const Basic > &a)
Negation.
Definition: mul.cpp:438
RCP< const Basic > sin(const RCP< const Basic > &arg)
Canonicalize Sin:
Definition: functions.cpp:831
T tie(T... args)