1 #include <symengine/complex.h>
2 #include <symengine/ntheory.h>
3
4 namespace SymEngine
5 {
6
is_re_zero() const7 bool ComplexBase::is_re_zero() const
8 {
9 return this->real_part()->is_zero();
10 }
11
Complex(rational_class real,rational_class imaginary)12 Complex::Complex(rational_class real, rational_class imaginary)
13 : real_{real}, imaginary_{imaginary}
14 {
15 SYMENGINE_ASSIGN_TYPEID()
16 SYMENGINE_ASSERT(is_canonical(this->real_, this->imaginary_))
17 }
18
is_canonical(const rational_class & real,const rational_class & imaginary) const19 bool Complex::is_canonical(const rational_class &real,
20 const rational_class &imaginary) const
21 {
22 rational_class re = real;
23 rational_class im = imaginary;
24 canonicalize(re);
25 canonicalize(im);
26 // If 'im' is 0, it should not be Complex:
27 if (get_num(im) == 0)
28 return false;
29 // if 'real' or `imaginary` are not in canonical form:
30 if (get_num(re) != get_num(real))
31 return false;
32 if (get_den(re) != get_den(real))
33 return false;
34 if (get_num(im) != get_num(imaginary))
35 return false;
36 if (get_den(im) != get_den(imaginary))
37 return false;
38 return true;
39 }
40
__hash__() const41 hash_t Complex::__hash__() const
42 {
43 // only the least significant bits that fit into "signed long int" are
44 // hashed:
45 hash_t seed = SYMENGINE_COMPLEX;
46 hash_combine<long long int>(seed, mp_get_si(get_num(this->real_)));
47 hash_combine<long long int>(seed, mp_get_si(get_den(this->real_)));
48 hash_combine<long long int>(seed, mp_get_si(get_num(this->imaginary_)));
49 hash_combine<long long int>(seed, mp_get_si(get_den(this->imaginary_)));
50 return seed;
51 }
52
__eq__(const Basic & o) const53 bool Complex::__eq__(const Basic &o) const
54 {
55 if (is_a<Complex>(o)) {
56 const Complex &s = down_cast<const Complex &>(o);
57 return ((this->real_ == s.real_)
58 and (this->imaginary_ == s.imaginary_));
59 }
60 return false;
61 }
62
compare(const Basic & o) const63 int Complex::compare(const Basic &o) const
64 {
65 SYMENGINE_ASSERT(is_a<Complex>(o))
66 const Complex &s = down_cast<const Complex &>(o);
67 if (real_ == s.real_) {
68 if (imaginary_ == s.imaginary_) {
69 return 0;
70 } else {
71 return imaginary_ < s.imaginary_ ? -1 : 1;
72 }
73 } else {
74 return real_ < s.real_ ? -1 : 1;
75 }
76 }
77
real_part() const78 RCP<const Number> Complex::real_part() const
79 {
80 return Rational::from_mpq(real_);
81 }
82
imaginary_part() const83 RCP<const Number> Complex::imaginary_part() const
84 {
85 return Rational::from_mpq(imaginary_);
86 }
87
conjugate() const88 RCP<const Basic> Complex::conjugate() const
89 {
90 return Complex::from_mpq(real_, -imaginary_);
91 }
92
from_mpq(const rational_class re,const rational_class im)93 RCP<const Number> Complex::from_mpq(const rational_class re,
94 const rational_class im)
95 {
96 // It is assumed that `re` and `im` are already in canonical form.
97 if (get_num(im) == 0) {
98 return Rational::from_mpq(re);
99 } else {
100 return make_rcp<const Complex>(re, im);
101 }
102 }
103
from_two_rats(const Rational & re,const Rational & im)104 RCP<const Number> Complex::from_two_rats(const Rational &re, const Rational &im)
105 {
106 return Complex::from_mpq(re.as_rational_class(), im.as_rational_class());
107 }
108
from_two_nums(const Number & re,const Number & im)109 RCP<const Number> Complex::from_two_nums(const Number &re, const Number &im)
110 {
111 if (is_a<Integer>(re) and is_a<Integer>(im)) {
112 rational_class re_mpq(
113 down_cast<const Integer &>(re).as_integer_class(),
114 down_cast<const Integer &>(*one).as_integer_class());
115 rational_class im_mpq(
116 down_cast<const Integer &>(im).as_integer_class(),
117 down_cast<const Integer &>(*one).as_integer_class());
118 return Complex::from_mpq(re_mpq, im_mpq);
119 } else if (is_a<Rational>(re) and is_a<Integer>(im)) {
120 rational_class re_mpq
121 = down_cast<const Rational &>(re).as_rational_class();
122 rational_class im_mpq(
123 down_cast<const Integer &>(im).as_integer_class(),
124 down_cast<const Integer &>(*one).as_integer_class());
125 return Complex::from_mpq(re_mpq, im_mpq);
126 } else if (is_a<Integer>(re) and is_a<Rational>(im)) {
127 rational_class re_mpq(
128 down_cast<const Integer &>(re).as_integer_class(),
129 down_cast<const Integer &>(*one).as_integer_class());
130 rational_class im_mpq
131 = down_cast<const Rational &>(im).as_rational_class();
132 return Complex::from_mpq(re_mpq, im_mpq);
133 } else if (is_a<Rational>(re) and is_a<Rational>(im)) {
134 rational_class re_mpq
135 = down_cast<const Rational &>(re).as_rational_class();
136 rational_class im_mpq
137 = down_cast<const Rational &>(im).as_rational_class();
138 return Complex::from_mpq(re_mpq, im_mpq);
139 } else {
140 throw SymEngineException(
141 "Invalid Format: Expected Integer or Rational");
142 }
143 }
144
pow_number(const Complex & x,unsigned long n)145 RCP<const Number> pow_number(const Complex &x, unsigned long n)
146 {
147 unsigned long mask = 1;
148 rational_class r_re(1);
149 rational_class r_im(0);
150
151 rational_class p_re = x.real_;
152 rational_class p_im = x.imaginary_;
153
154 rational_class tmp;
155
156 while (true) {
157 if (n & mask) {
158 // Multiply r by p
159 tmp = r_re * p_re - r_im * p_im;
160 r_im = r_re * p_im + r_im * p_re;
161 r_re = tmp;
162 }
163 mask = mask << 1;
164 if (not(mask > 0 and n >= mask)) {
165 break;
166 }
167 // Multiply p by p
168 tmp = p_re * p_re - p_im * p_im;
169 p_im = 2 * p_re * p_im;
170 p_re = tmp;
171 }
172 return Complex::from_mpq(r_re, r_im);
173 }
174
powcomp(const Integer & other) const175 RCP<const Number> Complex::powcomp(const Integer &other) const
176 {
177 if (this->is_re_zero()) {
178 // Imaginary Number raised to an integer power.
179 RCP<const Number> im = Rational::from_mpq(this->imaginary_);
180 long rem = mod_f(other, *integer(4))->as_int();
181 RCP<const Number> res;
182 if (rem == 0) {
183 res = one;
184 } else if (rem == 1) {
185 res = I;
186 } else if (rem == 2) {
187 res = minus_one;
188 } else {
189 res = mulnum(I, minus_one);
190 }
191 return mulnum(im->pow(other), res);
192 } else if (other.is_positive()) {
193 return pow_number(*this, other.as_int());
194 } else {
195 return one->div(*pow_number(*this, -1 * other.as_int()));
196 }
197 }
198 } // namespace SymEngine
199