complex.cpp
1 #include <symengine/complex.h>
2 #include <symengine/ntheory.h>
3 
4 namespace SymEngine
5 {
6 
7 bool ComplexBase::is_re_zero() const
8 {
9  return this->real_part()->is_zero();
10 }
11 
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 
19 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 
41 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 
53 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 
63 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 
78 RCP<const Number> Complex::real_part() const
79 {
80  return Rational::from_mpq(real_);
81 }
82 
83 RCP<const Number> Complex::imaginary_part() const
84 {
85  return Rational::from_mpq(imaginary_);
86 }
87 
88 RCP<const Basic> Complex::conjugate() const
89 {
90  return Complex::from_mpq(real_, -imaginary_);
91 }
92 
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 
104 RCP<const Number> Complex::from_two_rats(const Rational &re, const Rational &im)
105 {
107 }
108 
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 
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 
175 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
The lowest unit of symbolic representation.
Definition: basic.h:95
Complex Class.
Definition: complex.h:33
Complex(rational_class real, rational_class imaginary)
Constructor of Complex class.
Definition: complex.cpp:12
static RCP< const Number > from_mpq(const rational_class re, const rational_class im)
Definition: complex.cpp:93
virtual RCP< const Number > real_part() const
Get the real part of the complex number.
Definition: complex.cpp:78
virtual hash_t __hash__() const
Definition: complex.cpp:41
rational_class real_
Definition: complex.h:38
RCP< const Number > powcomp(const Integer &other) const
Definition: complex.cpp:175
virtual RCP< const Number > imaginary_part() const
Get the imaginary part of the complex number.
Definition: complex.cpp:83
virtual RCP< const Basic > conjugate() const
Get the conjugate of the complex number.
Definition: complex.cpp:88
virtual bool __eq__(const Basic &o) const
Definition: complex.cpp:53
static RCP< const Number > from_two_rats(const Rational &re, const Rational &im)
Definition: complex.cpp:104
bool is_canonical(const rational_class &real, const rational_class &imaginary) const
Definition: complex.cpp:19
virtual int compare(const Basic &o) const
Definition: complex.cpp:63
static RCP< const Number > from_two_nums(const Number &re, const Number &im)
Definition: complex.cpp:109
Integer Class.
Definition: integer.h:19
signed long int as_int() const
Convert to int, raise an exception if it does not fit.
Definition: integer.cpp:34
virtual bool is_positive() const
Definition: integer.h:68
Rational Class.
Definition: rational.h:16
static RCP< const Number > from_mpq(const rational_class &i)
Definition: rational.cpp:23
const rational_class & as_rational_class() const
Convert to rational_class.
Definition: rational.h:50
Main namespace for SymEngine package.
Definition: add.cpp:19
std::enable_if< std::is_integral< T >::value, RCP< const Integer > >::type integer(T i)
Definition: integer.h:200
RCP< const Number > mulnum(const RCP< const Number > &self, const RCP< const Number > &other)
Multiply self and other
Definition: number.h:93
RCP< const Integer > mod_f(const Integer &n, const Integer &d)
modulo round toward -inf
Definition: ntheory.cpp:84