TBTK
Need a break? Support the development by playing Polarity Puzzles
Complex.h
1 /* Copyright 2019 Kristofer Björnson
2  *
3  * Licensed under the Apache License, Version 2.0 (the "License");
4  * you may not use this file except in compliance with the License.
5  * You may obtain a copy of the License at
6  *
7  * http://www.apache.org/licenses/LICENSE-2.0
8  *
9  * Unless required by applicable law or agreed to in writing, software
10  * distributed under the License is distributed on an "AS IS" BASIS,
11  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12  * See the License for the specific language governing permissions and
13  * limitations under the License.
14  */
15 
17 
24 #ifndef COM_DAFER45_TBTK_ARBITRARY_PRECISION_COMPLEX
25 #define COM_DAFER45_TBTK_ARBITRARY_PRECISION_COMPLEX
26 
27 #include "TBTK/ArbitraryPrecision/Real.h"
28 #include "TBTK/Streams.h"
29 #include "TBTK/TBTKMacros.h"
30 
31 #include <complex>
32 #include <sstream>
33 #include <string>
34 
35 #include <gmpxx.h>
36 
37 namespace TBTK{
38 namespace ArbitraryPrecision{
39 
40 class Complex{
41 public:
47  Complex();
48 
53  Complex(unsigned int precision);
54 
61  Complex(unsigned int precision, const std::complex<double> &value);
62 
70  Complex(unsigned int precision, double real, double imag);
71 
76  Complex(const Real &real, const Real &imag);
77 
83  Complex& operator=(const std::complex<double> &rhs);
84 
92  Complex& operator=(const std::string &rhs);
93 
99  Complex& operator+=(const Complex &rhs);
100 
106  Complex operator+(const Complex &rhs) const;
107 
113  Complex& operator-=(const Complex &rhs);
114 
120  Complex operator-(const Complex &rhs) const;
121 
127  Complex& operator*=(const Complex &rhs);
128 
134  Complex operator*(const Complex &rhs) const;
135 
141  Complex& operator/=(const Complex &rhs);
142 
148  Complex operator/(const Complex &rhs) const;
149 
153  Complex operator-() const;
154 
159  friend std::ostream& operator<<(std::ostream &os, const Complex &complex);
160 
164  std::complex<double> getComplexDouble() const;
165 
169  const Real& getReal() const;
170 
174  const Real& getImag() const;
175 private:
177  Real real;
178 
180  Real imag;
181 };
182 
183 inline Complex::Complex(){
184 }
185 
186 inline Complex::Complex(
187  unsigned int precision
188 ):
189  real(precision),
190  imag(precision)
191 {
192 }
193 
194 inline Complex::Complex(
195  unsigned int precision,
196  const std::complex<double> &value
197 ) :
198  real(precision, value.real()),
199  imag(precision, value.imag())
200 {
201 }
202 
203 inline Complex::Complex(
204  unsigned int precision,
205  double real,
206  double imag
207 ) :
208  real(precision, real),
209  imag(precision, imag)
210 {
211 }
212 
213 inline Complex::Complex(
214  const Real &real,
215  const Real &imag
216 ) :
217  real(real),
218  imag(imag)
219 {
220 }
221 
222 inline Complex& Complex::operator=(const std::complex<double> &rhs){
223  real = rhs.real();
224  imag = rhs.imag();
225 
226  return *this;
227 }
228 
229 inline Complex& Complex::operator=(const std::string &rhs){
230  std::stringstream ss(rhs);
231  const int STATE_BEFORE_FIRST_TERM = 0;
232  const int STATE_READING_FIRST_TERM = 1;
233  const int STATE_AFTER_FIRST_TERM = 2;
234  const int STATE_BEFORE_SECOND_TERM = 3;
235  const int STATE_READING_SECOND_TERM = 4;
236  const int STATE_AFTER_SECOND_TERM = 5;
237  int state = STATE_BEFORE_FIRST_TERM;
238  char c;
239  std::string firstTerm;
240  std::string secondTerm;
241  bool firstTermIsReal = true;
242  bool secondTermIsReal = true;
243  std::string firstTermSign;
244  std::string secondTermSign;
245  int position = 0;
246  while(ss >> c){
247  switch(c){
248  case '0':
249  case '1':
250  case '2':
251  case '3':
252  case '4':
253  case '5':
254  case '6':
255  case '7':
256  case '8':
257  case '9':
258  case '.':
259  switch(state){
260  case STATE_BEFORE_FIRST_TERM:
261  state = STATE_READING_FIRST_TERM;
262  case STATE_READING_FIRST_TERM:
263  firstTerm += c;
264  break;
265  case STATE_BEFORE_SECOND_TERM:
266  state = STATE_READING_SECOND_TERM;
267  case STATE_READING_SECOND_TERM:
268  secondTerm += c;
269  break;
270  default:
271  TBTKExit(
272  "ArbitraryPrecision::Complex::operator=()",
273  "Unable to parse '" << rhs << "' as a"
274  << " complex number. Found unexpected"
275  << " tocken '" << c << "' at position"
276  << " '" << position << "'.",
277  ""
278  );
279  }
280  break;
281  case '+':
282  switch(state){
283  case STATE_BEFORE_FIRST_TERM:
284  if(firstTermSign.length() != 0){
285  TBTKExit(
286  "ArbitraryPrecision::Complex::operator=()",
287  "Unable to parse '" << rhs
288  << "' as a complex number."
289  << " Found unexpected tocken '"
290  << c << "' at position" << " '"
291  << position << "'.",
292  ""
293  );
294  }
295  firstTermSign = "+";
296  break;
297  case STATE_READING_FIRST_TERM:
298  case STATE_AFTER_FIRST_TERM:
299  state = STATE_BEFORE_SECOND_TERM;
300  if(secondTermSign.length() != 0){
301  TBTKExit(
302  "ArbitraryPrecision::Complex::operator=()",
303  "Unable to parse '" << rhs
304  << "' as a complex number."
305  << " Found unexpected tocken '"
306  << c << "' at position" << " '"
307  << position << "'.",
308  ""
309  );
310  }
311  secondTermSign = "+";
312  break;
313  default:
314  TBTKExit(
315  "ArbitraryPrecision::Complex::operator=()",
316  "Unable to parse '" << rhs << "' as a"
317  << " complex number. Found unexpected"
318  << " tocken '" << c << "' at position"
319  << " '" << position << "'.",
320  ""
321  );
322  }
323  break;
324  case '-':
325  switch(state){
326  case STATE_BEFORE_FIRST_TERM:
327  if(firstTermSign.length() != 0){
328  TBTKExit(
329  "ArbitraryPrecision::Complex::operator=()",
330  "Unable to parse '" << rhs
331  << "' as a complex number."
332  << " Found unexpected tocken '"
333  << c << "' at position" << " '"
334  << position << "'.",
335  ""
336  );
337  }
338  firstTermSign = "-";
339  break;
340  case STATE_READING_FIRST_TERM:
341  case STATE_AFTER_FIRST_TERM:
342  state = STATE_BEFORE_SECOND_TERM;
343  if(secondTermSign.length() != 0){
344  TBTKExit(
345  "ArbitraryPrecision::Complex::operator=()",
346  "Unable to parse '" << rhs
347  << "' as a complex number."
348  << " Found unexpected tocken '"
349  << c << "' at position" << " '"
350  << position << "'.",
351  ""
352  );
353  }
354  secondTermSign = "-";
355  break;
356  default:
357  TBTKExit(
358  "ArbitraryPrecision::Complex::operator=()",
359  "Unable to parse '" << rhs << "' as a"
360  << " complex number. Found unexpected"
361  << " tocken '" << c << "' at position"
362  << " '" << position << "'.",
363  ""
364  );
365  }
366  break;
367  case 'i':
368  switch(state){
369  case STATE_READING_FIRST_TERM:
370  state = STATE_AFTER_FIRST_TERM;
371  case STATE_AFTER_FIRST_TERM:
372  case STATE_BEFORE_FIRST_TERM:
373  if(!firstTermIsReal){
374  TBTKExit(
375  "ArbitraryPrecision::Complex::operator=()",
376  "Unable to parse '" << rhs
377  << "' as a complex number."
378  << " Found unexpected tocken '"
379  << c << "' at position '"
380  << position << "'.",
381  ""
382  );
383  }
384  firstTermIsReal = false;
385  break;
386  case STATE_READING_SECOND_TERM:
387  state = STATE_AFTER_SECOND_TERM;
388  case STATE_AFTER_SECOND_TERM:
389  case STATE_BEFORE_SECOND_TERM:
390  if(!secondTermIsReal){
391  TBTKExit(
392  "ArbitraryPrecision::Complex::operator=()",
393  "Unable to parse '" << rhs
394  << "' as a complex number."
395  << " Found unexpected tocken '"
396  << c << "' at position '"
397  << position << "'.",
398  ""
399  );
400  }
401  secondTermIsReal = false;
402  break;
403  }
404  break;
405  case ' ':
406  switch(state){
407  case STATE_BEFORE_FIRST_TERM:
408  case STATE_AFTER_FIRST_TERM:
409  case STATE_BEFORE_SECOND_TERM:
410  case STATE_AFTER_SECOND_TERM:
411  break;
412  case STATE_READING_FIRST_TERM:
413  state = STATE_AFTER_FIRST_TERM;
414  break;
415  case STATE_READING_SECOND_TERM:
416  state = STATE_AFTER_SECOND_TERM;
417  break;
418  }
419  break;
420  default:
421  TBTKExit(
422  "ArbitraryPrecision::Complex::operator=()",
423  "Unable to parse '" << rhs << "' as a complex"
424  << " number. Found unexpected tocken '" << c
425  << "' at position '" << position << "'.",
426  ""
427  );
428  }
429 
430  position++;
431  }
432 
433  if(secondTerm.length() == 0 && firstTermIsReal)
434  secondTermIsReal = false;
435 
436  TBTKAssert(
437  firstTermIsReal != secondTermIsReal,
438  "ArbitraryPrecision::Complex::operator=()",
439  "Unable to parse '" << rhs << "' as a complex number. The two"
440  << " terms cannot both be real or both be imaginary.",
441  ""
442  );
443 
444  if(firstTermSign.compare("-") == 0)
445  firstTerm = "-" + firstTerm;
446  if(secondTermSign.compare("-") == 0)
447  secondTerm = "-" + secondTerm;
448 
449  if(firstTermIsReal){
450  real = firstTerm;
451  imag = secondTerm;
452  }
453  else{
454  real = secondTerm;
455  imag = firstTerm;
456  }
457 
458  return *this;
459 }
460 
461 inline Complex& Complex::operator+=(const Complex &rhs){
462  real += rhs.real;
463  imag += rhs.imag;
464 
465  return *this;
466 }
467 
468 inline Complex Complex::operator+(const Complex &rhs) const{
469  Complex complex = *this;
470 
471  return complex += rhs;
472 }
473 
474 inline Complex& Complex::operator-=(const Complex &rhs){
475  real -= rhs.real;
476  imag -= rhs.imag;
477 
478  return *this;
479 }
480 
481 inline Complex Complex::operator-(const Complex &rhs) const{
482  Complex complex = *this;
483 
484  return complex -= rhs;
485 }
486 
487 inline Complex& Complex::operator*=(const Complex &rhs){
488  Real r = real;
489  Real i = imag;
490  real = r*rhs.real - i*rhs.imag;
491  imag = r*rhs.imag + i*rhs.real;
492 
493  return *this;
494 }
495 
496 inline Complex Complex::operator*(const Complex &rhs) const{
497  Complex complex = *this;
498 
499  return complex *= rhs;
500 }
501 
502 inline Complex& Complex::operator/=(const Complex &rhs){
503  Real denominator = rhs.real*rhs.real + rhs.imag*rhs.imag;
504 
505  Real r = real;
506  Real i = imag;
507  real = (r*rhs.real + i*rhs.imag)/denominator;
508  imag = (i*rhs.real - r*rhs.imag)/denominator;
509 
510  return *this;
511 }
512 
513 inline Complex Complex::operator/(const Complex &rhs) const{
514  Complex complex = *this;
515 
516  return complex /= rhs;
517 }
518 
519 inline Complex Complex::operator-() const{
520  return Complex(-real, -imag);
521 }
522 
523 inline std::ostream& operator<<(std::ostream &os, const Complex &complex){
524  if(complex.imag.getDouble() >= 0)
525  os << complex.real << " + i" << complex.imag;
526  else
527  os << complex.real << " - i" << -complex.imag;
528 
529  return os;
530 }
531 
532 inline std::complex<double> Complex::getComplexDouble() const{
533  return std::complex<double>(real.getDouble(), imag.getDouble());
534 }
535 
536 inline const Real& Complex::getReal() const{
537  return real;
538 }
539 
540 inline const Real& Complex::getImag() const{
541  return imag;
542 }
543 
544 //Complex pow(Complex base, long unsigned int power);
545 inline Complex pow(Complex base, long unsigned int power){
546  Complex result;
547  if(power%2)
548  result = base;
549  else{
550  mp_bitcnt_t realPrecision = base.getReal().getPrecision();
551  mp_bitcnt_t imagPrecision = base.getReal().getPrecision();
552  if(realPrecision > imagPrecision)
553  result = Complex(realPrecision, 1);
554  else
555  result = Complex(realPrecision, 1);
556  }
557 
558  while(power >>= 1){
559  //Update to *= when this has been implemented.
560  base = base*base;
561  if(power%2){
562  //Update to *= when this has been implemented.
563  result = result*base;
564  }
565  }
566 
567  return result;
568 }
569 
570 inline Complex pow(const Complex &base, long int power){
571  if(power < 0){
572  Complex result = pow(base, (long unsigned int)(-power));
573 
574  return Complex(result.getReal().getPrecision(), 1)/result;
575  }
576  else{
577  return pow(base, (long unsigned int)power);
578  }
579 }
580 
581 inline Complex pow(const Complex &base, int power){
582  return pow(base, (long int)power);
583 }
584 
585 }; //End of namespace ArbitraryPrecision
586 }; //End of namesapce TBTK
587 
588 #endif
589 
TBTK::Math::real
Array< double > real(const Array< std::complex< double >> &array)
Definition: ArrayAlgorithms.h:802
TBTK::Complex::operator+
friend Complex operator+(const Complex &lhs, const double &rhs)
Definition: Complex.h:114
TBTK::Complex::operator/=
Complex & operator/=(const Complex &rhs)
Definition: Complex.h:231
TBTK::Complex::operator*=
Complex & operator*=(const Complex &rhs)
Definition: Complex.h:194
Streams.h
Streams for TBTK output.
TBTK::Complex::real
constexpr friend Real real(const Complex &complex)
Definition: Complex.h:305
TBTK::Complex::operator*
friend Complex operator*(const Complex &lhs, const double &rhs)
Definition: Complex.h:206
TBTK::Complex::operator-
Complex operator-() const
Definition: Complex.h:185
TBTK::Complex::operator=
Complex & operator=(const std::complex< double > &rhs)
Definition: Complex.h:76
TBTK::Complex::operator+=
Complex & operator+=(const Complex &rhs)
Definition: Complex.h:87
TBTK::Complex::imag
constexpr friend Real imag(const Complex &complex)
Definition: Complex.h:314
TBTK::Math::pow
Array< DataType > pow(const Array< DataType > &array, double exponent)
Definition: ArrayAlgorithms.h:747
TBTK::Complex::operator/
friend Complex operator/(const Complex &lhs, const double &rhs)
Definition: Complex.h:242
TBTK::Math::imag
Array< double > imag(const Array< std::complex< double >> &array)
Definition: ArrayAlgorithms.h:812
TBTK::Complex::Complex
Complex()
Definition: Complex.h:41
TBTKMacros.h
Precompiler macros.
TBTK::operator<<
std::ostream & operator<<(std::ostream &stream, const HoppingAmplitude &hoppingAmplitude)
Definition: HoppingAmplitude.h:315
TBTK::operator*
const Vector2d operator*(double lhs, const Vector2d &rhs)
Definition: Vector2d.h:129
TBTK::Complex::operator-=
Complex & operator-=(const Complex &rhs)
Definition: Complex.h:139