Overload << operator for ckcomplex
[charm.git] / src / util / ckcomplex.h
1 #ifndef __CKCOMPLEX_H__
2 #define __CKCOMPLEX_H__
3
4 #include "pup.h"
5
6 #if USE_FFTW_DECLS
7 #include "fftw.h"
8 typedef fftw_real  RealType;
9 #else
10 typedef double     RealType;
11 #endif
12
13 #include <cmath>
14 #include <ostream>
15
16 struct ckcomplex {
17     RealType  re;
18     RealType  im;   
19
20  
21     inline ckcomplex(RealType _re=0., RealType _im=0.): re(_re), im(_im){}
22     //    inline ckcomplex(RealType r) {re=r; im=0;}
23     //    inline ckcomplex(RealType r,RealType i) {re=r; im=i;}
24     
25     inline ~ckcomplex() {}
26
27     inline RealType getMagSqr(void) const { 
28       return re*re+im*im; 
29     }
30
31     inline ckcomplex operator+(ckcomplex a) { 
32       return ckcomplex(re+a.re,im+a.im); 
33     }
34
35     // note: not a member
36     inline friend ckcomplex operator-(ckcomplex lhs, ckcomplex rhs) {
37       return ckcomplex(lhs.re - rhs.re, lhs.im - rhs.im);
38     }
39
40     inline ckcomplex conj(void) { 
41         return ckcomplex(re, -im); 
42     }
43
44     inline void operator+=(ckcomplex a) { 
45       re+=a.re; im+=a.im; 
46     }
47     
48     // note: not a member
49     inline friend ckcomplex operator*(RealType lhs, ckcomplex rhs) {
50       return ckcomplex(rhs.re*lhs, rhs.im*lhs);
51     }
52
53     inline friend ckcomplex operator/(ckcomplex lhs, RealType rhs) {
54         return ckcomplex(lhs.re/rhs, lhs.im/rhs);
55     }
56
57     inline ckcomplex operator*(RealType a) { 
58       return ckcomplex(re*a, im*a); 
59     } 
60
61     inline bool notzero() const { return( (0.0 != re) ? true : (0.0 != im)); }
62
63     inline void operator*=(ckcomplex a) {        
64         RealType treal, tim;
65         treal = re * a.re - im * a.im;
66         tim = re * a.im + im * a.re;
67         re = treal;
68         im = tim;
69     }
70
71     inline ckcomplex operator*(ckcomplex a) {
72       return ckcomplex( re * a.re - im * a.im, re * a.im + im * a.re); 
73     }
74
75
76     inline void operator -= (ckcomplex a) {
77       re -= a.re; im -= a.im;
78     }
79
80
81     inline ckcomplex multiplyByi () {
82         return ckcomplex(-im, re);
83     }
84         
85     inline void * operator new[] (size_t size){
86         void *buf = malloc(size);
87         return buf;
88     }
89     
90     inline void operator delete[] (void *buf){
91         free(buf);
92     }
93
94     inline friend std::ostream& operator<< (std::ostream &out, const ckcomplex &aNum) {
95         out<<"("<<aNum.re<<","<<aNum.im<<")";
96         return out;
97     }
98 };
99
100 typedef ckcomplex complex;
101
102 PUPbytes(ckcomplex)
103
104
105 /// Overload std::isfinite for complex numbers. @note: Doesn't seem to be part of the standard
106 inline int isfinite(complex aNum) { return ( std::isfinite(aNum.re) && std::isfinite(aNum.im) ); }
107
108 /// Like std::norm, return the square of the distance from (0,0) in the complex number plane
109 inline RealType norm(complex aNum) { return ( aNum.re*aNum.re + aNum.im*aNum.im ); }
110 /// Return the distance from (0,0) in the complex plane
111 inline RealType abs(complex aNum) { return ( sqrt(norm(aNum)) ); }
112
113 #endif
114