57729e5b
David Mayerich
changed directory...
|
1
2
3
4
5
6
7
|
/*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
|
7731cf24
David Mayerich
fixed precision p...
|
8
|
#include "../cuda/callable.h"
|
57729e5b
David Mayerich
changed directory...
|
9
10
11
12
13
|
#include <cmath>
#include <string>
#include <sstream>
#include <iostream>
|
8a86bd56
David Mayerich
changed rts names...
|
14
|
namespace stim
|
57729e5b
David Mayerich
changed directory...
|
15
16
17
18
19
20
21
22
23
24
|
{
template <class T>
struct complex
{
T r, i;
//default constructor
CUDA_CALLABLE complex()
{
|
7731cf24
David Mayerich
fixed precision p...
|
25
|
r = 0;
|
559d0fcb
David Mayerich
wave interactions...
|
26
27
28
29
30
31
32
33
|
i = 0;
}
//constructor when given real and imaginary values
CUDA_CALLABLE complex(T r, T i = 0)
{
this->r = r;
this->i = i;
|
57729e5b
David Mayerich
changed directory...
|
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
|
}
//access methods
CUDA_CALLABLE T real()
{
return r;
}
CUDA_CALLABLE T real(T r_val)
{
r = r_val;
return r_val;
}
CUDA_CALLABLE T imag()
{
return i;
}
CUDA_CALLABLE T imag(T i_val)
{
i = i_val;
return i_val;
}
|
559d0fcb
David Mayerich
wave interactions...
|
58
|
|
57729e5b
David Mayerich
changed directory...
|
59
60
61
62
63
64
65
66
67
68
69
|
//return the current value multiplied by i
CUDA_CALLABLE complex<T> imul()
{
complex<T> result;
result.r = -i;
result.i = r;
return result;
}
|
ecfd14df
David Mayerich
implemented D fie...
|
70
71
72
73
74
75
76
|
//returns the complex signum (-1, 0, 1)
CUDA_CALLABLE int sgn(){
if(r > 0) return 1;
else if(r < 0) return -1;
else return (0 < i - i < 0);
}
|
57729e5b
David Mayerich
changed directory...
|
77
78
79
|
//ARITHMETIC OPERATORS--------------------
//binary + operator (returns the result of adding two complex values)
|
ecfd14df
David Mayerich
implemented D fie...
|
80
|
CUDA_CALLABLE complex<T> operator+ (const complex<T> rhs) const
|
57729e5b
David Mayerich
changed directory...
|
81
82
83
84
85
86
87
|
{
complex<T> result;
result.r = r + rhs.r;
result.i = i + rhs.i;
return result;
}
|
ecfd14df
David Mayerich
implemented D fie...
|
88
|
CUDA_CALLABLE complex<T> operator+ (const T rhs) const
|
57729e5b
David Mayerich
changed directory...
|
89
90
91
92
93
94
95
96
|
{
complex<T> result;
result.r = r + rhs;
result.i = i;
return result;
}
//binary - operator (returns the result of adding two complex values)
|
ecfd14df
David Mayerich
implemented D fie...
|
97
|
CUDA_CALLABLE complex<T> operator- (const complex<T> rhs) const
|
57729e5b
David Mayerich
changed directory...
|
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
|
{
complex<T> result;
result.r = r - rhs.r;
result.i = i - rhs.i;
return result;
}
//binary - operator (returns the result of adding two complex values)
CUDA_CALLABLE complex<T> operator- (const T rhs)
{
complex<T> result;
result.r = r - rhs;
result.i = i;
return result;
}
//binary MULTIPLICATION operators (returns the result of multiplying complex values)
|
ecfd14df
David Mayerich
implemented D fie...
|
115
|
CUDA_CALLABLE complex<T> operator* (const complex<T> rhs) const
|
57729e5b
David Mayerich
changed directory...
|
116
117
118
119
120
121
122
123
124
125
126
127
|
{
complex<T> result;
result.r = r * rhs.r - i * rhs.i;
result.i = r * rhs.i + i * rhs.r;
return result;
}
CUDA_CALLABLE complex<T> operator* (const T rhs)
{
return complex<T>(r * rhs, i * rhs);
}
//binary DIVISION operators (returns the result of dividing complex values)
|
ecfd14df
David Mayerich
implemented D fie...
|
128
|
CUDA_CALLABLE complex<T> operator/ (const complex<T> rhs) const
|
57729e5b
David Mayerich
changed directory...
|
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
|
{
complex<T> result;
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;
}
CUDA_CALLABLE complex<T> operator/ (const T rhs)
{
return complex<T>(r / rhs, i / rhs);
}
//ASSIGNMENT operators-----------------------------------
CUDA_CALLABLE complex<T> & operator=(const complex<T> &rhs)
{
//check for self-assignment
if(this != &rhs)
{
this->r = rhs.r;
this->i = rhs.i;
}
return *this;
}
CUDA_CALLABLE complex<T> & operator=(const T &rhs)
{
|
3b012a80
David Mayerich
added code for co...
|
155
156
|
this->r = rhs;
this->i = 0;
|
57729e5b
David Mayerich
changed directory...
|
157
158
159
160
161
162
163
164
|
return *this;
}
//arithmetic assignment operators
CUDA_CALLABLE complex<T> operator+=(const complex<T> &rhs)
{
*this = *this + rhs;
|
3b012a80
David Mayerich
added code for co...
|
165
|
return *this;
|
57729e5b
David Mayerich
changed directory...
|
166
167
168
169
|
}
CUDA_CALLABLE complex<T> operator+=(const T &rhs)
{
*this = *this + rhs;
|
3b012a80
David Mayerich
added code for co...
|
170
|
return *this;
|
57729e5b
David Mayerich
changed directory...
|
171
172
|
}
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
173
174
175
176
177
178
179
180
181
182
183
|
CUDA_CALLABLE complex<T> operator-=(const complex<T> &rhs)
{
*this = *this - rhs;
return *this;
}
CUDA_CALLABLE complex<T> operator-=(const T &rhs)
{
*this = *this - rhs;
return *this;
}
|
57729e5b
David Mayerich
changed directory...
|
184
185
186
|
CUDA_CALLABLE complex<T> operator*=(const complex<T> &rhs)
{
*this = *this * rhs;
|
3b012a80
David Mayerich
added code for co...
|
187
|
return *this;
|
57729e5b
David Mayerich
changed directory...
|
188
189
190
191
|
}
CUDA_CALLABLE complex<T> operator*=(const T &rhs)
{
*this = *this * rhs;
|
3b012a80
David Mayerich
added code for co...
|
192
|
return *this;
|
57729e5b
David Mayerich
changed directory...
|
193
194
195
196
197
|
}
//divide and assign
CUDA_CALLABLE complex<T> operator/=(const complex<T> &rhs)
{
*this = *this / rhs;
|
3b012a80
David Mayerich
added code for co...
|
198
|
return *this;
|
57729e5b
David Mayerich
changed directory...
|
199
200
201
202
|
}
CUDA_CALLABLE complex<T> operator/=(const T &rhs)
{
*this = *this / rhs;
|
3b012a80
David Mayerich
added code for co...
|
203
|
return *this;
|
57729e5b
David Mayerich
changed directory...
|
204
205
206
207
208
209
210
211
212
213
|
}
//absolute value operator (returns the absolute value of the complex number)
CUDA_CALLABLE T abs()
{
return std::sqrt(r * r + i * i);
}
CUDA_CALLABLE complex<T> log()
{
|
3b012a80
David Mayerich
added code for co...
|
214
215
216
|
complex<T> result;
result.r = (T)std::log(std::sqrt(r * r + i * i));
result.i = (T)std::atan2(i, r);
|
57729e5b
David Mayerich
changed directory...
|
217
218
|
|
3b012a80
David Mayerich
added code for co...
|
219
|
return result;
|
57729e5b
David Mayerich
changed directory...
|
220
221
222
223
|
}
CUDA_CALLABLE complex<T> exp()
{
|
3b012a80
David Mayerich
added code for co...
|
224
|
complex<T> result;
|
57729e5b
David Mayerich
changed directory...
|
225
|
|
3b012a80
David Mayerich
added code for co...
|
226
227
228
|
T e_r = std::exp(r);
result.r = e_r * (T)std::cos(i);
result.i = e_r * (T)std::sin(i);
|
57729e5b
David Mayerich
changed directory...
|
229
|
|
3b012a80
David Mayerich
added code for co...
|
230
|
return result;
|
57729e5b
David Mayerich
changed directory...
|
231
232
233
234
235
236
237
238
239
240
|
}
/*CUDA_CALLABLE complex<T> pow(int y)
{
return pow((double)y);
}*/
CUDA_CALLABLE complex<T> pow(T y)
{
|
3b012a80
David Mayerich
added code for co...
|
241
|
complex<T> result;
|
57729e5b
David Mayerich
changed directory...
|
242
|
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
243
|
result = log() * y;
|
57729e5b
David Mayerich
changed directory...
|
244
|
|
3b012a80
David Mayerich
added code for co...
|
245
|
return result.exp();
|
57729e5b
David Mayerich
changed directory...
|
246
247
248
249
250
251
252
253
254
255
256
257
|
}
CUDA_CALLABLE complex<T> sqrt()
{
complex<T> result;
//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);
|
7731cf24
David Mayerich
fixed precision p...
|
258
|
T theta_p = theta/2.0f;
|
57729e5b
David Mayerich
changed directory...
|
259
260
261
262
263
264
265
266
|
//convert back to cartesian coordinates
result.r = a_p * std::cos(theta_p);
result.i = a_p * std::sin(theta_p);
return result;
}
|
ecfd14df
David Mayerich
implemented D fie...
|
267
|
std::string str()
|
57729e5b
David Mayerich
changed directory...
|
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
|
{
std::stringstream ss;
ss<<"("<<r<<","<<i<<")";
return ss.str();
}
//COMPARISON operators
CUDA_CALLABLE bool operator==(complex<T> rhs)
{
if(r == rhs.r && i == rhs.i)
return true;
return false;
}
CUDA_CALLABLE bool operator==(T rhs)
{
|
7731cf24
David Mayerich
fixed precision p...
|
285
|
if(r == rhs && i == 0)
|
57729e5b
David Mayerich
changed directory...
|
286
287
288
289
|
return true;
return false;
}
|
ecfd14df
David Mayerich
implemented D fie...
|
290
291
292
293
294
295
296
|
CUDA_CALLABLE bool operator!=(T rhs)
{
if(r != rhs || i != 0)
return true;
return false;
}
|
0ef519a4
David Mayerich
optimized materia...
|
297
298
299
300
301
302
303
304
305
306
307
308
309
|
CUDA_CALLABLE bool operator<(complex<T> rhs){
return abs() < rhs.abs();
}
CUDA_CALLABLE bool operator<=(complex<T> rhs){
return abs() <= rhs.abs();
}
CUDA_CALLABLE bool operator>(complex<T> rhs){
return abs() > rhs.abs();
}
CUDA_CALLABLE bool operator >=(complex<T> rhs){
return abs() >= rhs.abs();
}
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
|
//CASTING operators
template < typename otherT >
operator complex<otherT>()
{
complex<otherT> result((otherT)r, (otherT)i);
return result;
}
template< typename otherT >
complex( const complex<otherT> &rhs)
{
r = (T)rhs.r;
i = (T)rhs.i;
}
template< typename otherT >
complex& operator=(const complex<otherT> &rhs)
{
r = (T)rhs.r;
i = (T)rhs.i;
return *this;
}
|
57729e5b
David Mayerich
changed directory...
|
331
332
333
334
335
336
|
};
} //end RTS namespace
//addition
template<typename T>
|
ef70330e
David Mayerich
changed rts to stim
|
337
|
CUDA_CALLABLE static stim::complex<T> operator+(const double a, const stim::complex<T> b)
|
57729e5b
David Mayerich
changed directory...
|
338
|
{
|
ef70330e
David Mayerich
changed rts to stim
|
339
|
return stim::complex<T>((T)a + b.r, b.i);
|
57729e5b
David Mayerich
changed directory...
|
340
341
342
343
|
}
//subtraction with a real value
template<typename T>
|
ef70330e
David Mayerich
changed rts to stim
|
344
|
CUDA_CALLABLE static stim::complex<T> operator-(const double a, const stim::complex<T> b)
|
57729e5b
David Mayerich
changed directory...
|
345
|
{
|
ef70330e
David Mayerich
changed rts to stim
|
346
|
return stim::complex<T>((T)a - b.r, -b.i);
|
57729e5b
David Mayerich
changed directory...
|
347
348
349
350
|
}
//minus sign
template<typename T>
|
ef70330e
David Mayerich
changed rts to stim
|
351
|
CUDA_CALLABLE static stim::complex<T> operator-(const stim::complex<T> &rhs)
|
57729e5b
David Mayerich
changed directory...
|
352
|
{
|
ef70330e
David Mayerich
changed rts to stim
|
353
|
return stim::complex<T>(-rhs.r, -rhs.i);
|
57729e5b
David Mayerich
changed directory...
|
354
355
356
357
|
}
//multiply a T value by a complex value
template<typename T>
|
ef70330e
David Mayerich
changed rts to stim
|
358
|
CUDA_CALLABLE static stim::complex<T> operator*(const double a, const stim::complex<T> b)
|
57729e5b
David Mayerich
changed directory...
|
359
|
{
|
ef70330e
David Mayerich
changed rts to stim
|
360
|
return stim::complex<T>((T)a * b.r, (T)a * b.i);
|
57729e5b
David Mayerich
changed directory...
|
361
362
363
364
|
}
//divide a T value by a complex value
template<typename T>
|
ef70330e
David Mayerich
changed rts to stim
|
365
|
CUDA_CALLABLE static stim::complex<T> operator/(const double a, const stim::complex<T> b)
|
57729e5b
David Mayerich
changed directory...
|
366
|
{
|
ef70330e
David Mayerich
changed rts to stim
|
367
|
stim::complex<T> result;
|
57729e5b
David Mayerich
changed directory...
|
368
369
370
|
T denom = b.r * b.r + b.i * b.i;
|
3b012a80
David Mayerich
added code for co...
|
371
372
|
result.r = ((T)a * b.r) / denom;
result.i = -((T)a * b.i) / denom;
|
57729e5b
David Mayerich
changed directory...
|
373
374
375
376
|
return result;
}
|
57729e5b
David Mayerich
changed directory...
|
377
378
|
template<typename T>
|
ef70330e
David Mayerich
changed rts to stim
|
379
|
CUDA_CALLABLE static stim::complex<T> pow(stim::complex<T> x, T y)
|
57729e5b
David Mayerich
changed directory...
|
380
381
382
|
{
return x.pow(y);
}
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
383
|
template<typename T>
|
ef70330e
David Mayerich
changed rts to stim
|
384
|
CUDA_CALLABLE static stim::complex<T> pow(stim::complex<T> x, int y)
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
385
386
387
|
{
return x.pow(y);
}
|
57729e5b
David Mayerich
changed directory...
|
388
389
390
|
//log function
template<typename T>
|
ef70330e
David Mayerich
changed rts to stim
|
391
|
CUDA_CALLABLE static stim::complex<T> log(stim::complex<T> x)
|
57729e5b
David Mayerich
changed directory...
|
392
393
394
395
396
397
|
{
return x.log();
}
//exp function
template<typename T>
|
ef70330e
David Mayerich
changed rts to stim
|
398
|
CUDA_CALLABLE static stim::complex<T> exp(stim::complex<T> x)
|
57729e5b
David Mayerich
changed directory...
|
399
400
401
402
403
404
|
{
return x.exp();
}
//sqrt function
template<typename T>
|
ef70330e
David Mayerich
changed rts to stim
|
405
|
CUDA_CALLABLE static stim::complex<T> sqrt(stim::complex<T> x)
|
57729e5b
David Mayerich
changed directory...
|
406
407
408
409
410
411
|
{
return x.sqrt();
}
template <typename T>
|
ef70330e
David Mayerich
changed rts to stim
|
412
|
CUDA_CALLABLE static T abs(stim::complex<T> a)
|
57729e5b
David Mayerich
changed directory...
|
413
414
415
416
417
|
{
return a.abs();
}
template <typename T>
|
ef70330e
David Mayerich
changed rts to stim
|
418
|
CUDA_CALLABLE static T real(stim::complex<T> a)
|
57729e5b
David Mayerich
changed directory...
|
419
420
421
422
423
424
425
426
427
428
429
|
{
return a.r;
}
//template <typename T>
CUDA_CALLABLE static float real(float a)
{
return a;
}
template <typename T>
|
ef70330e
David Mayerich
changed rts to stim
|
430
|
CUDA_CALLABLE static T imag(stim::complex<T> a)
|
57729e5b
David Mayerich
changed directory...
|
431
432
433
434
435
|
{
return a.i;
}
//trigonometric functions
|
3b012a80
David Mayerich
added code for co...
|
436
|
//template<class A>
|
ef70330e
David Mayerich
changed rts to stim
|
437
|
/*CUDA_CALLABLE static stim::complex<float> sinf(const stim::complex<float> x)
|
3b012a80
David Mayerich
added code for co...
|
438
|
{
|
ef70330e
David Mayerich
changed rts to stim
|
439
|
stim::complex<float> result;
|
3b012a80
David Mayerich
added code for co...
|
440
441
442
443
444
445
|
result.r = sinf(x.r) * coshf(x.i);
result.i = cosf(x.r) * sinhf(x.i);
return result;
}*/
|
57729e5b
David Mayerich
changed directory...
|
446
|
template<class A>
|
ef70330e
David Mayerich
changed rts to stim
|
447
|
CUDA_CALLABLE stim::complex<A> sin(const stim::complex<A> x)
|
57729e5b
David Mayerich
changed directory...
|
448
|
{
|
ef70330e
David Mayerich
changed rts to stim
|
449
|
stim::complex<A> result;
|
3b012a80
David Mayerich
added code for co...
|
450
451
|
result.r = (A)std::sin(x.r) * (A)std::cosh(x.i);
result.i = (A)std::cos(x.r) * (A)std::sinh(x.i);
|
57729e5b
David Mayerich
changed directory...
|
452
453
454
455
|
return result;
}
|
3b012a80
David Mayerich
added code for co...
|
456
457
|
//floating point template
//template<class A>
|
ef70330e
David Mayerich
changed rts to stim
|
458
|
/*CUDA_CALLABLE static stim::complex<float> cosf(const stim::complex<float> x)
|
3b012a80
David Mayerich
added code for co...
|
459
|
{
|
ef70330e
David Mayerich
changed rts to stim
|
460
|
stim::complex<float> result;
|
3b012a80
David Mayerich
added code for co...
|
461
462
463
464
465
466
|
result.r = cosf(x.r) * coshf(x.i);
result.i = -(sinf(x.r) * sinhf(x.i));
return result;
}*/
|
57729e5b
David Mayerich
changed directory...
|
467
|
template<class A>
|
ef70330e
David Mayerich
changed rts to stim
|
468
|
CUDA_CALLABLE stim::complex<A> cos(const stim::complex<A> x)
|
57729e5b
David Mayerich
changed directory...
|
469
|
{
|
ef70330e
David Mayerich
changed rts to stim
|
470
|
stim::complex<A> result;
|
3b012a80
David Mayerich
added code for co...
|
471
472
|
result.r = (A)std::cos(x.r) * (A)std::cosh(x.i);
result.i = -((A)std::sin(x.r) * (A)std::sinh(x.i));
|
57729e5b
David Mayerich
changed directory...
|
473
474
475
476
477
478
|
return result;
}
template<class A>
|
ef70330e
David Mayerich
changed rts to stim
|
479
|
std::ostream& operator<<(std::ostream& os, stim::complex<A> x)
|
57729e5b
David Mayerich
changed directory...
|
480
|
{
|
ecfd14df
David Mayerich
implemented D fie...
|
481
|
os<<x.str();
|
57729e5b
David Mayerich
changed directory...
|
482
483
484
|
return os;
}
|
0ef519a4
David Mayerich
optimized materia...
|
485
|
template<class A>
|
ef70330e
David Mayerich
changed rts to stim
|
486
|
std::istream& operator>>(std::istream& is, stim::complex<A>& x)
|
0ef519a4
David Mayerich
optimized materia...
|
487
488
489
490
491
492
493
494
495
496
497
498
|
{
A r, i;
r = i = 0; //initialize the real and imaginary parts to zero
is>>r; //parse
is>>i;
x.real(r); //assign the parsed values to x
x.imag(i);
return is; //return the stream
}
|
57729e5b
David Mayerich
changed directory...
|
499
|
//#if __GNUC__ > 3 && __GNUC_MINOR__ > 7
|
ef70330e
David Mayerich
changed rts to stim
|
500
|
//template<class T> using rtsComplex = stim::complex<T>;
|
57729e5b
David Mayerich
changed directory...
|
501
502
503
504
505
|
//#endif
#endif
|