f1402849
dmayerich
renewed commit
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
|
/*RTS Complex number class. This class is CUDA compatible,
and can therefore be used in CUDA code and on CUDA devices.
*/
#ifndef RTS_COMPLEX
#define RTS_COMPLEX
#include "cuda_callable.h"
#include <cmath>
#include <string>
#include <sstream>
#include <iostream>
namespace rts
{
template <class T>
|
7006df5f
David Mayerich
reformat of direc...
|
18
|
struct complex
|
f1402849
dmayerich
renewed commit
|
19
20
21
22
|
{
T r, i;
//default constructor
|
7006df5f
David Mayerich
reformat of direc...
|
23
|
CUDA_CALLABLE complex()
|
f1402849
dmayerich
renewed commit
|
24
25
26
27
28
29
|
{
r = 0.0;
i = 0.0;
}
//access methods
|
a47a23a9
dmayerich
added ENVI functions
|
30
|
CUDA_CALLABLE T real()
|
f1402849
dmayerich
renewed commit
|
31
32
33
34
|
{
return r;
}
|
a47a23a9
dmayerich
added ENVI functions
|
35
|
CUDA_CALLABLE T real(T r_val)
|
f1402849
dmayerich
renewed commit
|
36
37
38
39
40
|
{
r = r_val;
return r_val;
}
|
a47a23a9
dmayerich
added ENVI functions
|
41
|
CUDA_CALLABLE T imag()
|
f1402849
dmayerich
renewed commit
|
42
43
44
|
{
return i;
}
|
a47a23a9
dmayerich
added ENVI functions
|
45
|
CUDA_CALLABLE T imag(T i_val)
|
f1402849
dmayerich
renewed commit
|
46
47
48
49
50
|
{
i = i_val;
return i_val;
}
|
f1402849
dmayerich
renewed commit
|
51
|
//constructor when given real and imaginary values
|
7006df5f
David Mayerich
reformat of direc...
|
52
|
CUDA_CALLABLE complex(T r, T i)
|
f1402849
dmayerich
renewed commit
|
53
54
55
56
57
58
|
{
this->r = r;
this->i = i;
}
//return the current value multiplied by i
|
7006df5f
David Mayerich
reformat of direc...
|
59
|
CUDA_CALLABLE complex<T> imul()
|
f1402849
dmayerich
renewed commit
|
60
|
{
|
7006df5f
David Mayerich
reformat of direc...
|
61
|
complex<T> result;
|
f1402849
dmayerich
renewed commit
|
62
63
64
65
66
67
68
69
70
|
result.r = -i;
result.i = r;
return result;
}
//ARITHMETIC OPERATORS--------------------
//binary + operator (returns the result of adding two complex values)
|
7006df5f
David Mayerich
reformat of direc...
|
71
|
CUDA_CALLABLE complex<T> operator+ (const complex<T> rhs)
|
f1402849
dmayerich
renewed commit
|
72
|
{
|
7006df5f
David Mayerich
reformat of direc...
|
73
|
complex<T> result;
|
f1402849
dmayerich
renewed commit
|
74
75
76
77
78
|
result.r = r + rhs.r;
result.i = i + rhs.i;
return result;
}
|
7006df5f
David Mayerich
reformat of direc...
|
79
|
CUDA_CALLABLE complex<T> operator+ (const T rhs)
|
f1402849
dmayerich
renewed commit
|
80
|
{
|
7006df5f
David Mayerich
reformat of direc...
|
81
|
complex<T> result;
|
f1402849
dmayerich
renewed commit
|
82
83
84
85
86
87
|
result.r = r + rhs;
result.i = i;
return result;
}
//binary - operator (returns the result of adding two complex values)
|
7006df5f
David Mayerich
reformat of direc...
|
88
|
CUDA_CALLABLE complex<T> operator- (const complex<T> rhs)
|
f1402849
dmayerich
renewed commit
|
89
|
{
|
7006df5f
David Mayerich
reformat of direc...
|
90
|
complex<T> result;
|
f1402849
dmayerich
renewed commit
|
91
92
93
94
95
96
|
result.r = r - rhs.r;
result.i = i - rhs.i;
return result;
}
//binary - operator (returns the result of adding two complex values)
|
7006df5f
David Mayerich
reformat of direc...
|
97
|
CUDA_CALLABLE complex<T> operator- (const T rhs)
|
f1402849
dmayerich
renewed commit
|
98
|
{
|
7006df5f
David Mayerich
reformat of direc...
|
99
|
complex<T> result;
|
f1402849
dmayerich
renewed commit
|
100
101
102
103
104
105
|
result.r = r - rhs;
result.i = i;
return result;
}
//binary MULTIPLICATION operators (returns the result of multiplying complex values)
|
7006df5f
David Mayerich
reformat of direc...
|
106
|
CUDA_CALLABLE complex<T> operator* (const complex<T> rhs)
|
f1402849
dmayerich
renewed commit
|
107
|
{
|
7006df5f
David Mayerich
reformat of direc...
|
108
|
complex<T> result;
|
f1402849
dmayerich
renewed commit
|
109
110
111
112
|
result.r = r * rhs.r - i * rhs.i;
result.i = r * rhs.i + i * rhs.r;
return result;
}
|
7006df5f
David Mayerich
reformat of direc...
|
113
|
CUDA_CALLABLE complex<T> operator* (const T rhs)
|
f1402849
dmayerich
renewed commit
|
114
|
{
|
7006df5f
David Mayerich
reformat of direc...
|
115
|
return complex<T>(r * rhs, i * rhs);
|
f1402849
dmayerich
renewed commit
|
116
117
118
|
}
//binary DIVISION operators (returns the result of dividing complex values)
|
7006df5f
David Mayerich
reformat of direc...
|
119
|
CUDA_CALLABLE complex<T> operator/ (const complex<T> rhs)
|
f1402849
dmayerich
renewed commit
|
120
|
{
|
7006df5f
David Mayerich
reformat of direc...
|
121
|
complex<T> result;
|
f1402849
dmayerich
renewed commit
|
122
123
124
125
126
127
|
T denom = rhs.r * rhs.r + rhs.i * rhs.i;
result.r = (r * rhs.r + i * rhs.i) / denom;
result.i = (- r * rhs.i + i * rhs.r) / denom;
return result;
}
|
7006df5f
David Mayerich
reformat of direc...
|
128
|
CUDA_CALLABLE complex<T> operator/ (const T rhs)
|
f1402849
dmayerich
renewed commit
|
129
|
{
|
7006df5f
David Mayerich
reformat of direc...
|
130
|
return complex<T>(r / rhs, i / rhs);
|
f1402849
dmayerich
renewed commit
|
131
132
133
|
}
//ASSIGNMENT operators-----------------------------------
|
7006df5f
David Mayerich
reformat of direc...
|
134
|
CUDA_CALLABLE complex<T> & operator=(const complex<T> &rhs)
|
f1402849
dmayerich
renewed commit
|
135
136
137
138
139
140
141
142
143
|
{
//check for self-assignment
if(this != &rhs)
{
this->r = rhs.r;
this->i = rhs.i;
}
return *this;
}
|
7006df5f
David Mayerich
reformat of direc...
|
144
|
CUDA_CALLABLE complex<T> & operator=(const T &rhs)
|
f1402849
dmayerich
renewed commit
|
145
146
147
148
149
150
151
152
|
{
this->r = rhs;
this->i = 0;
return *this;
}
//arithmetic assignment operators
|
7006df5f
David Mayerich
reformat of direc...
|
153
|
CUDA_CALLABLE complex<T> operator+=(const complex<T> &rhs)
|
f1402849
dmayerich
renewed commit
|
154
155
156
157
|
{
*this = *this + rhs;
return *this;
}
|
7006df5f
David Mayerich
reformat of direc...
|
158
|
CUDA_CALLABLE complex<T> operator+=(const T &rhs)
|
f1402849
dmayerich
renewed commit
|
159
160
161
162
163
|
{
*this = *this + rhs;
return *this;
}
|
7006df5f
David Mayerich
reformat of direc...
|
164
|
CUDA_CALLABLE complex<T> operator*=(const complex<T> &rhs)
|
f1402849
dmayerich
renewed commit
|
165
166
167
168
|
{
*this = *this * rhs;
return *this;
}
|
7006df5f
David Mayerich
reformat of direc...
|
169
|
CUDA_CALLABLE complex<T> operator*=(const T &rhs)
|
f1402849
dmayerich
renewed commit
|
170
171
172
173
174
|
{
*this = *this * rhs;
return *this;
}
//divide and assign
|
7006df5f
David Mayerich
reformat of direc...
|
175
|
CUDA_CALLABLE complex<T> operator/=(const complex<T> &rhs)
|
f1402849
dmayerich
renewed commit
|
176
177
178
179
|
{
*this = *this / rhs;
return *this;
}
|
7006df5f
David Mayerich
reformat of direc...
|
180
|
CUDA_CALLABLE complex<T> operator/=(const T &rhs)
|
f1402849
dmayerich
renewed commit
|
181
182
183
184
185
186
187
188
189
190
191
|
{
*this = *this / rhs;
return *this;
}
//absolute value operator (returns the absolute value of the complex number)
CUDA_CALLABLE T abs()
{
return std::sqrt(r * r + i * i);
}
|
7006df5f
David Mayerich
reformat of direc...
|
192
|
CUDA_CALLABLE complex<T> log()
|
f1402849
dmayerich
renewed commit
|
193
|
{
|
7006df5f
David Mayerich
reformat of direc...
|
194
|
complex<T> result;
|
f1402849
dmayerich
renewed commit
|
195
196
197
198
199
200
201
|
result.r = std::log(std::sqrt(r * r + i * i));
result.i = std::atan2(i, r);
return result;
}
|
7006df5f
David Mayerich
reformat of direc...
|
202
|
CUDA_CALLABLE complex<T> exp()
|
f1402849
dmayerich
renewed commit
|
203
|
{
|
7006df5f
David Mayerich
reformat of direc...
|
204
|
complex<T> result;
|
f1402849
dmayerich
renewed commit
|
205
206
207
208
209
210
211
212
213
214
215
216
217
218
|
T e_r = std::exp(r);
result.r = e_r * std::cos(i);
result.i = e_r * std::sin(i);
return result;
}
/*CUDA_CALLABLE complex<T> pow(int y)
{
return pow((double)y);
}*/
|
7006df5f
David Mayerich
reformat of direc...
|
219
|
CUDA_CALLABLE complex<T> pow(T y)
|
f1402849
dmayerich
renewed commit
|
220
|
{
|
7006df5f
David Mayerich
reformat of direc...
|
221
|
complex<T> result;
|
f1402849
dmayerich
renewed commit
|
222
223
224
225
226
227
|
result = log() * y;
return result.exp();
}
|
7006df5f
David Mayerich
reformat of direc...
|
228
|
CUDA_CALLABLE complex<T> sqrt()
|
f1402849
dmayerich
renewed commit
|
229
|
{
|
7006df5f
David Mayerich
reformat of direc...
|
230
|
complex<T> result;
|
f1402849
dmayerich
renewed commit
|
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
|
//convert to polar coordinates
T a = std::sqrt(r*r + i*i);
T theta = std::atan2(i, r);
//find the square root
T a_p = std::sqrt(a);
T theta_p = theta/2.0;
//convert back to cartesian coordinates
result.r = a_p * std::cos(theta_p);
result.i = a_p * std::sin(theta_p);
return result;
}
std::string toStr()
{
std::stringstream ss;
ss<<"("<<r<<","<<i<<")";
return ss.str();
}
//COMPARISON operators
|
7006df5f
David Mayerich
reformat of direc...
|
256
|
CUDA_CALLABLE bool operator==(complex<T> rhs)
|
f1402849
dmayerich
renewed commit
|
257
258
259
260
261
262
263
264
265
266
267
268
269
|
{
if(r == rhs.r && i == rhs.i)
return true;
return false;
}
CUDA_CALLABLE bool operator==(T rhs)
{
if(r == rhs && i == (T)0.0)
return true;
return false;
}
|
f1402849
dmayerich
renewed commit
|
270
271
|
};
|
a47a23a9
dmayerich
added ENVI functions
|
272
273
|
} //end RTS namespace
|
f1402849
dmayerich
renewed commit
|
274
275
|
//addition
template<typename T>
|
7006df5f
David Mayerich
reformat of direc...
|
276
|
CUDA_CALLABLE static rts::complex<T> operator+(const double a, const rts::complex<T> b)
|
f1402849
dmayerich
renewed commit
|
277
|
{
|
7006df5f
David Mayerich
reformat of direc...
|
278
|
return rts::complex<T>(a + b.r, b.i);
|
f1402849
dmayerich
renewed commit
|
279
280
281
282
|
}
//subtraction with a real value
template<typename T>
|
7006df5f
David Mayerich
reformat of direc...
|
283
|
CUDA_CALLABLE static rts::complex<T> operator-(const double a, const rts::complex<T> b)
|
f1402849
dmayerich
renewed commit
|
284
|
{
|
7006df5f
David Mayerich
reformat of direc...
|
285
|
return rts::complex<T>(a - b.r, -b.i);
|
f1402849
dmayerich
renewed commit
|
286
287
288
289
|
}
//minus sign
template<typename T>
|
7006df5f
David Mayerich
reformat of direc...
|
290
|
CUDA_CALLABLE static rts::complex<T> operator-(const rts::complex<T> &rhs)
|
f1402849
dmayerich
renewed commit
|
291
|
{
|
7006df5f
David Mayerich
reformat of direc...
|
292
|
return rts::complex<T>(-rhs.r, -rhs.i);
|
f1402849
dmayerich
renewed commit
|
293
294
295
296
|
}
//multiply a T value by a complex value
template<typename T>
|
7006df5f
David Mayerich
reformat of direc...
|
297
|
CUDA_CALLABLE static rts::complex<T> operator*(const double a, const rts::complex<T> b)
|
f1402849
dmayerich
renewed commit
|
298
|
{
|
7006df5f
David Mayerich
reformat of direc...
|
299
|
return rts::complex<T>((T)a * b.r, (T)a * b.i);
|
f1402849
dmayerich
renewed commit
|
300
301
302
303
|
}
//divide a T value by a complex value
template<typename T>
|
7006df5f
David Mayerich
reformat of direc...
|
304
|
CUDA_CALLABLE static rts::complex<T> operator/(const double a, const rts::complex<T> b)
|
f1402849
dmayerich
renewed commit
|
305
306
|
{
//return complex<T>(a * b.r, a * b.i);
|
7006df5f
David Mayerich
reformat of direc...
|
307
|
rts::complex<T> result;
|
f1402849
dmayerich
renewed commit
|
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
|
T denom = b.r * b.r + b.i * b.i;
result.r = (a * b.r) / denom;
result.i = -(a * b.i) / denom;
return result;
}
//POW function
/*template<typename T>
CUDA_CALLABLE static complex<T> pow(complex<T> x, int y)
{
return x.pow(y);
}*/
template<typename T>
|
7006df5f
David Mayerich
reformat of direc...
|
325
|
CUDA_CALLABLE static rts::complex<T> pow(rts::complex<T> x, T y)
|
f1402849
dmayerich
renewed commit
|
326
327
328
329
330
331
|
{
return x.pow(y);
}
//log function
template<typename T>
|
7006df5f
David Mayerich
reformat of direc...
|
332
|
CUDA_CALLABLE static rts::complex<T> log(rts::complex<T> x)
|
f1402849
dmayerich
renewed commit
|
333
334
335
336
337
338
|
{
return x.log();
}
//exp function
template<typename T>
|
7006df5f
David Mayerich
reformat of direc...
|
339
|
CUDA_CALLABLE static rts::complex<T> exp(rts::complex<T> x)
|
f1402849
dmayerich
renewed commit
|
340
341
342
343
344
345
|
{
return x.exp();
}
//sqrt function
template<typename T>
|
7006df5f
David Mayerich
reformat of direc...
|
346
|
CUDA_CALLABLE static rts::complex<T> sqrt(rts::complex<T> x)
|
f1402849
dmayerich
renewed commit
|
347
348
349
350
351
352
|
{
return x.sqrt();
}
template <typename T>
|
7006df5f
David Mayerich
reformat of direc...
|
353
|
CUDA_CALLABLE static T abs(rts::complex<T> a)
|
f1402849
dmayerich
renewed commit
|
354
355
356
357
358
|
{
return a.abs();
}
template <typename T>
|
7006df5f
David Mayerich
reformat of direc...
|
359
|
CUDA_CALLABLE static T real(rts::complex<T> a)
|
f1402849
dmayerich
renewed commit
|
360
361
362
363
|
{
return a.r;
}
|
a47a23a9
dmayerich
added ENVI functions
|
364
365
366
367
368
369
|
//template <typename T>
CUDA_CALLABLE static float real(float a)
{
return a;
}
|
f1402849
dmayerich
renewed commit
|
370
|
template <typename T>
|
7006df5f
David Mayerich
reformat of direc...
|
371
|
CUDA_CALLABLE static T imag(rts::complex<T> a)
|
f1402849
dmayerich
renewed commit
|
372
373
374
375
376
377
|
{
return a.i;
}
//trigonometric functions
template<class A>
|
7006df5f
David Mayerich
reformat of direc...
|
378
|
CUDA_CALLABLE rts::complex<A> sin(const rts::complex<A> x)
|
f1402849
dmayerich
renewed commit
|
379
|
{
|
7006df5f
David Mayerich
reformat of direc...
|
380
|
rts::complex<A> result;
|
f1402849
dmayerich
renewed commit
|
381
382
383
384
385
386
387
|
result.r = std::sin(x.r) * std::cosh(x.i);
result.i = std::cos(x.r) * std::sinh(x.i);
return result;
}
template<class A>
|
7006df5f
David Mayerich
reformat of direc...
|
388
|
CUDA_CALLABLE rts::complex<A> cos(const rts::complex<A> x)
|
f1402849
dmayerich
renewed commit
|
389
|
{
|
7006df5f
David Mayerich
reformat of direc...
|
390
|
rts::complex<A> result;
|
f1402849
dmayerich
renewed commit
|
391
392
393
394
395
396
397
|
result.r = std::cos(x.r) * std::cosh(x.i);
result.i = -(std::sin(x.r) * std::sinh(x.i));
return result;
}
|
f1402849
dmayerich
renewed commit
|
398
|
template<class A>
|
7006df5f
David Mayerich
reformat of direc...
|
399
|
std::ostream& operator<<(std::ostream& os, rts::complex<A> x)
|
f1402849
dmayerich
renewed commit
|
400
401
402
403
404
|
{
os<<x.toStr();
return os;
}
|
7006df5f
David Mayerich
reformat of direc...
|
405
406
407
408
|
#if __GNUC__ > 3 && __GNUC_MINOR__ > 7
template<class T> using rtsComplex = rts::complex<T>;
#endif
|
f1402849
dmayerich
renewed commit
|
409
410
411
|
#endif
|