Loading...
Searching...
No Matches
as_real_imag.cpp
1#include <symengine/visitor.h>
2#include <symengine/basic.h>
3
4namespace SymEngine
5{
6
7void 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
37class RealImagVisitor : public BaseVisitor<RealImagVisitor>
38{
39private:
40 Ptr<RCP<const Basic>> real_, imag_;
41
42public:
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
282void 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
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
virtual vec_basic get_args() const
Returns the arguments of the Add.
Definition: add.cpp:397
The lowest unit of symbolic representation.
Definition: basic.h:97
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:507
RCP< const Basic > get_arg() const
Definition: functions.h:36
RCP< const Basic > get_base() const
Definition: pow.h:37
RCP< const Basic > get_exp() const
Definition: pow.h:42
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:431
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:1614
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:2212
RCP< const Basic > mul(const RCP< const Basic > &a, const RCP< const Basic > &b)
Multiplication.
Definition: mul.cpp:352
RCP< const Basic > cos(const RCP< const Basic > &arg)
Canonicalize Cos:
Definition: functions.cpp:942
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:2127
RCP< const Basic > neg(const RCP< const Basic > &a)
Negation.
Definition: mul.cpp:443
RCP< const Basic > sin(const RCP< const Basic > &arg)
Canonicalize Sin:
Definition: functions.cpp:874
T tie(T... args)