/[escript]/trunk/cusplibrary/cusp/complex.h
ViewVC logotype

Contents of /trunk/cusplibrary/cusp/complex.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5588 - (show annotations)
Tue Apr 21 05:07:16 2015 UTC (5 years, 2 months ago) by caltinay
File MIME type: text/plain
File size: 35591 byte(s)
Update to copyright headers in cusp library.

1 /*
2 * Copyright 2008-2009 NVIDIA Corporation
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16
17 /*
18 * Copyright (c) 2010, The Regents of the University of California,
19 * through Lawrence Berkeley National Laboratory (subject to receipt of
20 * any required approvals from U.S. Dept. of Energy) All rights reserved.
21 *
22 * Redistribution and use in source and binary forms, with or
23 * without modification, are permitted provided that the
24 * following conditions are met:
25 *
26 * * Redistributions of source code must retain the above
27 * copyright notice, this list of conditions and the following
28 * disclaimer.
29 *
30 * * Redistributions in binary form must reproduce the
31 * above copyright notice, this list of conditions and the
32 * following disclaimer in the documentation and/or other
33 * materials provided with the distribution.
34 *
35 * * Neither the name of the University of California,
36 * Berkeley, nor the names of its contributors may be used to
37 * endorse or promote products derived from this software
38 * without specific prior written permission.
39 *
40 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
41 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
42 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
43 * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
44 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
45 * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
46 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
47 * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
48 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
49 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
50 * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
51 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
52 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
53 * SUCH DAMAGE.
54 *
55 */
56
57 /*
58 * Modifications to this file:
59 * Copyright (c) 2014-2015, The University of Queensland
60 * Licensed under the Apache License, Version 2.0.
61 *
62 */
63
64 /*! \file complex.h
65 * \brief Complex numbers
66 */
67
68 #pragma once
69
70 #include <cusp/detail/config.h>
71
72 #if (defined THRUST_DEVICE_BACKEND && THRUST_DEVICE_BACKEND == THRUST_DEVICE_BACKEND_CUDA) || (defined THRUST_DEVICE_SYSTEM && THRUST_DEVICE_SYSTEM == THRUST_DEVICE_SYSTEM_CUDA)
73
74 #ifdef _WIN32
75 #define _USE_MATH_DEFINES 1 // make sure M_PI is defined
76 #endif
77
78 #include <math.h>
79 #include <complex>
80 #include <cuComplex.h>
81 #include <sstream>
82 #include <cusp/cmath.h>
83
84 namespace cusp
85 {
86
87 template <typename ValueType> struct complex;
88 template <> struct complex<float>;
89 template <> struct complex<double>;
90
91
92 /// Returns the magnitude of z.
93 template<typename ValueType>
94 __host__ __device__
95 ValueType abs(const complex<ValueType>& z);
96 /// Returns the phase angle of z.
97 template<typename ValueType>
98 __host__ __device__
99 ValueType arg(const complex<ValueType>& z);
100 /// Returns the magnitude of z squared.
101 template<typename ValueType>
102 __host__ __device__
103 ValueType norm(const complex<ValueType>& z);
104
105 /// Returns the complex conjugate of z.
106 template<typename ValueType>
107 __host__ __device__
108 complex<ValueType> conj(const complex<ValueType>& z);
109 /// Returns the complex with magnitude m and angle theta in radians.
110 template<typename ValueType>
111 __host__ __device__
112 complex<ValueType> polar(const ValueType& m, const ValueType& theta = 0);
113
114 // Arithmetic operators:
115 // Multiplication
116 template <typename ValueType>
117 __host__ __device__
118 inline complex<ValueType> operator*(const complex<ValueType>& lhs, const complex<ValueType>& rhs);
119 template <typename ValueType>
120 __host__ __device__
121 inline complex<ValueType> operator*(const complex<ValueType>& lhs, const ValueType & rhs);
122 template <typename ValueType>
123 __host__ __device__
124 inline complex<ValueType> operator*(const ValueType& lhs, const complex<ValueType>& rhs);
125 // Division
126 template <typename ValueType>
127 __host__ __device__
128 inline complex<ValueType> operator/(const complex<ValueType>& lhs, const complex<ValueType>& rhs);
129 template <>
130 __host__ __device__
131 inline complex<float> operator/(const complex<float>& lhs, const complex<float>& rhs);
132 template <>
133 __host__ __device__
134 inline complex<double> operator/(const complex<double>& lhs, const complex<double>& rhs);
135
136 // Addition
137 template <typename ValueType>
138 __host__ __device__
139 inline complex<ValueType> operator+(const complex<ValueType>& lhs, const complex<ValueType>& rhs);
140 template <typename ValueType>
141 __host__ __device__
142 inline complex<ValueType> operator+(const complex<ValueType>& lhs, const ValueType & rhs);
143 template <typename ValueType>
144 __host__ __device__
145 inline complex<ValueType> operator+(const ValueType& lhs, const complex<ValueType>& rhs);
146 // Subtraction
147 template <typename ValueType>
148 __host__ __device__
149 inline complex<ValueType> operator-(const complex<ValueType>& lhs, const complex<ValueType>& rhs);
150 template <typename ValueType>
151 __host__ __device__
152 inline complex<ValueType> operator-(const complex<ValueType>& lhs, const ValueType & rhs);
153 template <typename ValueType>
154 __host__ __device__
155 inline complex<ValueType> operator-(const ValueType& lhs, const complex<ValueType>& rhs);
156
157 // Unary plus and minus
158 template <typename ValueType>
159 __host__ __device__
160 inline complex<ValueType> operator+(const complex<ValueType>& rhs);
161 template <typename ValueType>
162 __host__ __device__
163 inline complex<ValueType> operator-(const complex<ValueType>& rhs);
164
165 // Transcendentals:
166 // Returns the complex cosine of z.
167 template <typename ValueType>
168 __host__ __device__
169 complex<ValueType> cos(const complex<ValueType>& z);
170 // Returns the complex hyperbolic cosine of z.
171 template <typename ValueType>
172 __host__ __device__
173 complex<ValueType> cosh(const complex<ValueType>& z);
174 // Returns the complex base e exponential of z.
175 template <typename ValueType>
176 __host__ __device__
177 complex<ValueType> exp(const complex<ValueType>& z);
178 // Returns the complex natural logarithm of z.
179 template <typename ValueType>
180 __host__ __device__
181 complex<ValueType> log(const complex<ValueType>& z);
182 // Returns the complex base 10 logarithm of z.
183 template <typename ValueType>
184 __host__ __device__
185 complex<ValueType> log10(const complex<ValueType>& z);
186 // Returns z to the n'th power.
187 template <typename ValueType>
188 __host__ __device__
189 complex<ValueType> pow(const complex<ValueType>& z, const int& n);
190 // Returns z to the x'th power.
191 template <typename ValueType>
192 __host__ __device__
193 complex<ValueType> pow(const complex<ValueType>&z, const ValueType&x);
194 // Returns z to the z2'th power.
195 template <typename ValueType>
196 __host__ __device__
197 complex<ValueType> pow(const complex<ValueType>&z,
198 const complex<ValueType>&z2);
199 // Returns x to the z'th power.
200 template <typename ValueType>
201 __host__ __device__
202 complex<ValueType> pow(const ValueType& x, const complex<ValueType>& z);
203 // Returns the complex sine of z.
204 template <typename ValueType>
205 __host__ __device__
206 complex<ValueType> sin(const complex<ValueType>&z);
207 // Returns the complex hyperbolic sine of z.
208 template <typename ValueType>
209 __host__ __device__
210 complex<ValueType> sinh(const complex<ValueType>&z);
211 // Returns the complex square root of z.
212 template <typename ValueType>
213 __host__ __device__
214 complex<ValueType> sqrt(const complex<ValueType>&z);
215 // Returns the complex tangent of z.
216 template <typename ValueType>
217 __host__ __device__
218 complex<ValueType> tan(const complex<ValueType>&z);
219 // Returns the complex hyperbolic tangent of z.
220 template <typename ValueType>
221 __host__ __device__
222 complex<ValueType> tanh(const complex<ValueType>&z);
223
224
225 // Inverse Trigonometric:
226 // Returns the complex arc cosine of z.
227 template <typename ValueType>
228 __host__ __device__
229 complex<ValueType> acos(const complex<ValueType>& z);
230 // Returns the complex arc sine of z.
231 template <typename ValueType>
232 __host__ __device__
233 complex<ValueType> asin(const complex<ValueType>& z);
234 // Returns the complex arc tangent of z.
235 template <typename ValueType>
236 __host__ __device__
237 complex<ValueType> atan(const complex<ValueType>& z);
238 // Returns the complex hyperbolic arc cosine of z.
239 template <typename ValueType>
240 __host__ __device__
241 complex<ValueType> acosh(const complex<ValueType>& z);
242 // Returns the complex hyperbolic arc sine of z.
243 template <typename ValueType>
244 __host__ __device__
245 complex<ValueType> asinh(const complex<ValueType>& z);
246 // Returns the complex hyperbolic arc tangent of z.
247 template <typename ValueType>
248 __host__ __device__
249 complex<ValueType> atanh(const complex<ValueType>& z);
250
251
252
253 // Stream operators:
254 template<typename ValueType,class charT, class traits>
255 std::basic_ostream<charT, traits>& operator<<(std::basic_ostream<charT, traits>& os, const complex<ValueType>& z);
256 template<typename ValueType, typename charT, class traits>
257 std::basic_istream<charT, traits>&
258 operator>>(std::basic_istream<charT, traits>& is, complex<ValueType>& z);
259
260
261 // Stream operators
262 template<typename ValueType,class charT, class traits>
263 std::basic_ostream<charT, traits>& operator<<(std::basic_ostream<charT, traits>& os, const complex<ValueType>& z)
264 {
265 os << '(' << z.real() << ',' << z.imag() << ')';
266 return os;
267 }
268
269 template<typename ValueType, typename charT, class traits>
270 std::basic_istream<charT, traits>&
271 operator>>(std::basic_istream<charT, traits>& is, complex<ValueType>& z)
272 {
273 ValueType re, im;
274
275 charT ch;
276 is >> ch;
277
278 if(ch == '(')
279 {
280 is >> re >> ch;
281 if (ch == ',')
282 {
283 is >> im >> ch;
284 if (ch == ')')
285 {
286 z = complex<ValueType>(re, im);
287 }
288 else
289 {
290 is.setstate(std::ios_base::failbit);
291 }
292 }
293 else if (ch == ')')
294 {
295 z = re;
296 }
297 else
298 {
299 is.setstate(std::ios_base::failbit);
300 }
301 }
302 else
303 {
304 is.putback(ch);
305 is >> re;
306 z = re;
307 }
308 return is;
309 }
310
311 template <typename T>
312 struct norm_type {
313 typedef T type;
314 };
315 template <typename T>
316 struct norm_type< complex<T> > {
317 typedef T type;
318 };
319
320 template <typename ValueType>
321 struct complex
322 {
323 public:
324 typedef ValueType value_type;
325
326 // Constructors
327 __host__ __device__
328 inline complex<ValueType>(const ValueType & re = ValueType(), const ValueType& im = ValueType())
329 {
330 real(re);
331 imag(im);
332 }
333
334 template <class X>
335 __host__ __device__
336 inline complex<ValueType>(const complex<X> & z)
337 {
338 real(z.real());
339 imag(z.imag());
340 }
341
342 template <class X>
343 __host__ __device__
344 inline complex<ValueType>(const std::complex<X> & z)
345 {
346 real(z.real());
347 imag(z.imag());
348 }
349
350 template <typename T>
351 __host__ __device__
352 inline complex<ValueType>& operator=(const complex<T> z)
353 {
354 real(z.real());
355 imag(z.imag());
356 return *this;
357 }
358
359 __host__ __device__
360 inline complex<ValueType>& operator+=(const complex<ValueType> z)
361 {
362 real(real()+z.real());
363 imag(imag()+z.imag());
364 return *this;
365 }
366
367 __host__ __device__
368 inline complex<ValueType>& operator-=(const complex<ValueType> z)
369 {
370 real(real()-z.real());
371 imag(imag()-z.imag());
372 return *this;
373 }
374
375 __host__ __device__
376 inline complex<ValueType>& operator*=(const complex<ValueType> z)
377 {
378 *this = *this * z;
379 return *this;
380 }
381
382 __host__ __device__
383 inline complex<ValueType>& operator/=(const complex<ValueType> z)
384 {
385 *this = *this / z;
386 return *this;
387 }
388
389 __host__ __device__ inline ValueType real() const volatile;
390 __host__ __device__ inline ValueType imag() const volatile;
391 __host__ __device__ inline ValueType real() const;
392 __host__ __device__ inline ValueType imag() const;
393 __host__ __device__ inline void real(ValueType) volatile;
394 __host__ __device__ inline void imag(ValueType) volatile;
395 __host__ __device__ inline void real(ValueType);
396 __host__ __device__ inline void imag(ValueType);
397 };
398
399 // TODO make cuFloatComplex and cuDoubleComplex protected
400 // TODO see if returning references is a perf hazard
401
402 template<>
403 struct complex <float> : public cuFloatComplex
404 {
405 public:
406 typedef float value_type;
407 __host__ __device__
408 inline complex<float>(){};
409 __host__ __device__
410 inline complex<float>(const float & re, const float& im = float())
411 {
412 real(re);
413 imag(im);
414 }
415
416 // For some reason having the following constructor
417 // explicitly makes things faster with at least g++
418 __host__ __device__
419 complex<float>(const complex<float> & z)
420 : cuFloatComplex(z){}
421
422 __host__ __device__
423 complex<float>(cuFloatComplex z)
424 : cuFloatComplex(z){}
425
426 template <class X>
427 inline complex<float>(const std::complex<X> & z)
428 {
429 real(z.real());
430 imag(z.imag());
431 }
432
433 // Member operators
434 template <typename T>
435 __host__ __device__
436 inline volatile complex<float>& operator=(const complex<T> z) volatile
437 {
438 real(z.real());
439 imag(z.imag());
440 return *this;
441 }
442
443 template <typename T>
444 __host__ __device__
445 inline complex<float>& operator=(const complex<T> z)
446 {
447 real(z.real());
448 imag(z.imag());
449 return *this;
450 }
451
452 __host__ __device__
453 inline complex<float>& operator+=(const complex<float> z)
454 {
455 real(real()+z.real());
456 imag(imag()+z.imag());
457 return *this;
458 }
459
460 __host__ __device__
461 inline complex<float>& operator-=(const complex<float> z)
462 {
463 real(real()-z.real());
464 imag(imag()-z.imag());
465 return *this;
466 }
467
468 __host__ __device__
469 inline complex<float>& operator*=(const complex<float> z)
470 {
471 *this = *this * z;
472 return *this;
473 }
474
475 __host__ __device__
476 inline complex<float>& operator/=(const complex<float> z)
477 {
478 *this = *this / z;
479 return *this;
480 }
481
482 // Let the compiler synthesize the copy and assignment operators.
483 __host__ __device__ inline complex<float>(const volatile complex<float> & z)
484 {
485 real(z.real());
486 imag(z.imag());
487 }
488
489 __host__ __device__ inline float real() const volatile{ return x; }
490 __host__ __device__ inline float imag() const volatile{ return y; }
491 __host__ __device__ inline float real() const{ return x; }
492 __host__ __device__ inline float imag() const{ return y; }
493 __host__ __device__ inline void real(float re)volatile{ x = re; }
494 __host__ __device__ inline void imag(float im)volatile{ y = im; }
495 __host__ __device__ inline void real(float re){ x = re; }
496 __host__ __device__ inline void imag(float im){ y = im; }
497
498 // cast operators
499 inline operator std::complex<float>() const { return std::complex<float>(real(),imag()); }
500 // inline operator float() const { return real(); }
501 };
502
503 template<>
504 struct complex <double> : public cuDoubleComplex
505 {
506 public:
507 typedef double value_type;
508 __host__ __device__
509 inline complex<double>(){};
510 __host__ __device__
511 inline complex<double>(const double & re, const double& im = double())
512 {
513 real(re);
514 imag(im);
515 }
516
517 // For some reason having the following constructor
518 // explicitly makes things faster with at least g++
519 __host__ __device__
520 inline complex<double>(const complex<double> & z)
521 : cuDoubleComplex(z) {}
522
523 __host__ __device__
524 inline complex<double>(cuDoubleComplex z)
525 : cuDoubleComplex(z) {}
526
527 template <class X>
528 inline complex<double>(const std::complex<X> & z)
529 {
530 real(z.real());
531 imag(z.imag());
532 }
533
534 // Member operators
535 template <typename T>
536 __host__ __device__
537 inline volatile complex<double>& operator=(const complex<T> z) volatile
538 {
539 real(z.real());
540 imag(z.imag());
541 return *this;
542 }
543
544 template <typename T>
545 __host__ __device__
546 inline complex<double>& operator=(const complex<T> z)
547 {
548 real(z.real());
549 imag(z.imag());
550 return *this;
551 }
552
553 __host__ __device__
554 inline complex<double>& operator+=(const complex<double> z)
555 {
556 real(real()+z.real());
557 imag(imag()+z.imag());
558 return *this;
559 }
560
561 __host__ __device__
562 inline complex<double>& operator-=(const complex<double> z)
563 {
564 real(real()-z.real());
565 imag(imag()-z.imag());
566 return *this;
567 }
568
569 __host__ __device__
570 inline complex<double>& operator*=(const complex<double> z)
571 {
572 *this = *this * z;
573 return *this;
574 }
575
576 __host__ __device__
577 inline complex<double>& operator/=(const complex<double> z)
578 {
579 *this = *this / z;
580 return *this;
581 }
582
583 __host__ __device__ inline complex<double>(const volatile complex<double> & z)
584 {
585 real(z.real());
586 imag(z.imag());
587 }
588
589 // Let the compiler synthesize the copy and assignment operators.
590 __host__ __device__ inline double real() const volatile { return x; }
591 __host__ __device__ inline double imag() const volatile { return y; }
592 __host__ __device__ inline double real() const { return x; }
593 __host__ __device__ inline double imag() const { return y; }
594 __host__ __device__ inline void real(double re)volatile{ x = re; }
595 __host__ __device__ inline void imag(double im)volatile{ y = im; }
596 __host__ __device__ inline void real(double re){ x = re; }
597 __host__ __device__ inline void imag(double im){ y = im; }
598
599 // cast operators
600 inline operator std::complex<double>() const { return std::complex<double>(real(),imag()); }
601 // inline operator double() { return real(); }
602 };
603
604
605
606 // Binary arithmetic operations
607 // At the moment I'm implementing the basic functions, and the
608 // corresponding cuComplex calls are commented.
609
610 template<typename ValueType>
611 __host__ __device__
612 inline complex<ValueType> operator+(const complex<ValueType>& lhs,
613 const complex<ValueType>& rhs){
614 return complex<ValueType>(lhs.real()+rhs.real(),lhs.imag()+rhs.imag());
615 // return cuCaddf(lhs,rhs);
616 }
617
618 template<typename ValueType>
619 __host__ __device__
620 inline complex<ValueType> operator+(const volatile complex<ValueType>& lhs,
621 const volatile complex<ValueType>& rhs){
622 return complex<ValueType>(lhs.real()+rhs.real(),lhs.imag()+rhs.imag());
623 // return cuCaddf(lhs,rhs);
624 }
625
626 template <typename ValueType>
627 __host__ __device__
628 inline complex<ValueType> operator+(const complex<ValueType>& lhs, const ValueType & rhs){
629 return complex<ValueType>(lhs.real()+rhs,lhs.imag());
630 // return cuCaddf(lhs,complex<ValueType>(rhs));
631 }
632 template <typename ValueType>
633 __host__ __device__
634 inline complex<ValueType> operator+(const ValueType& lhs, const complex<ValueType>& rhs){
635 return complex<ValueType>(rhs.real()+lhs,rhs.imag());
636 // return cuCaddf(complex<float>(lhs),rhs);
637 }
638
639 template <typename ValueType>
640 __host__ __device__
641 inline complex<ValueType> operator-(const complex<ValueType>& lhs, const complex<ValueType>& rhs){
642 return complex<ValueType>(lhs.real()-rhs.real(),lhs.imag()-rhs.imag());
643 // return cuCsubf(lhs,rhs);
644 }
645 template <typename ValueType>
646 __host__ __device__
647 inline complex<ValueType> operator-(const complex<ValueType>& lhs, const ValueType & rhs){
648 return complex<ValueType>(lhs.real()-rhs,lhs.imag());
649 // return cuCsubf(lhs,complex<float>(rhs));
650 }
651 template <typename ValueType>
652 __host__ __device__
653 inline complex<ValueType> operator-(const ValueType& lhs, const complex<ValueType>& rhs){
654 return complex<ValueType>(lhs-rhs.real(),-rhs.imag());
655 // return cuCsubf(complex<float>(lhs),rhs);
656 }
657
658 template <typename ValueType>
659 __host__ __device__
660 inline complex<ValueType> operator*(const complex<ValueType>& lhs,
661 const complex<ValueType>& rhs){
662 return complex<ValueType>(lhs.real()*rhs.real()-lhs.imag()*rhs.imag(),
663 lhs.real()*rhs.imag()+lhs.imag()*rhs.real());
664 // return cuCmulf(lhs,rhs);
665 }
666
667 template <typename ValueType>
668 __host__ __device__
669 inline complex<ValueType> operator*(const complex<ValueType>& lhs, const ValueType & rhs){
670 return complex<ValueType>(lhs.real()*rhs,lhs.imag()*rhs);
671 // return cuCmulf(lhs,complex<float>(rhs));
672 }
673
674 template <typename ValueType>
675 __host__ __device__
676 inline complex<ValueType> operator*(const ValueType& lhs, const complex<ValueType>& rhs){
677 return complex<ValueType>(rhs.real()*lhs,rhs.imag()*lhs);
678 // return cuCmulf(complex<float>(lhs),rhs);
679 }
680
681
682 template <typename ValueType>
683 __host__ __device__
684 inline complex<ValueType> operator/(const complex<ValueType>& lhs, const complex<ValueType>& rhs){
685 const ValueType cross_norm = lhs.real() * rhs.real() + lhs.imag() * rhs.imag();
686 const ValueType rhs_norm = norm(rhs);
687 return complex<ValueType>(cross_norm/rhs_norm,
688 (lhs.imag() * rhs.real() - lhs.real() * rhs.imag()) / rhs_norm);
689 }
690
691 template <>
692 __host__ __device__
693 inline complex<float> operator/(const complex<float>& lhs, const complex<float>& rhs){
694 return cuCdivf(lhs,rhs);
695 }
696
697 template <>
698 __host__ __device__
699 inline complex<double> operator/(const complex<double>& lhs, const complex<double>& rhs){
700 return cuCdiv(lhs,rhs);
701 }
702
703 template <typename ValueType>
704 __host__ __device__
705 inline complex<ValueType> operator/(const complex<ValueType>& lhs, const ValueType & rhs){
706 return complex<ValueType>(lhs.real()/rhs,lhs.imag()/rhs);
707 // return cuCdivf(lhs,complex<float>(rhs));
708 }
709
710 template <typename ValueType>
711 __host__ __device__
712 inline complex<ValueType> operator/(const ValueType& lhs, const complex<ValueType>& rhs){
713 const ValueType cross_norm = lhs * rhs.real();
714 const ValueType rhs_norm = norm(rhs);
715 return complex<ValueType>(cross_norm/rhs_norm,(-lhs.real() * rhs.imag()) / rhs_norm);
716 }
717
718 template <>
719 __host__ __device__
720 inline complex<float> operator/(const float& lhs, const complex<float>& rhs){
721 return cuCdivf(complex<float>(lhs),rhs);
722 }
723 template <>
724 __host__ __device__
725 inline complex<double> operator/(const double& lhs, const complex<double>& rhs){
726 return cuCdiv(complex<double>(lhs),rhs);
727 }
728
729
730 // Unary arithmetic operations
731 template <typename ValueType>
732 __host__ __device__
733 inline complex<ValueType> operator+(const complex<ValueType>& rhs){
734 return rhs;
735 }
736 template <typename ValueType>
737 __host__ __device__
738 inline complex<ValueType> operator-(const complex<ValueType>& rhs){
739 return rhs*-ValueType(1);
740 }
741
742 // Equality operators
743 template <typename ValueType>
744 __host__ __device__
745 inline bool operator==(const complex<ValueType>& lhs, const complex<ValueType>& rhs){
746 if(lhs.real() == rhs.real() && lhs.imag() == rhs.imag()){
747 return true;
748 }
749 return false;
750 }
751 template <typename ValueType>
752 __host__ __device__
753 inline bool operator==(const ValueType & lhs, const complex<ValueType>& rhs){
754 if(lhs == rhs.real() && rhs.imag() == 0){
755 return true;
756 }
757 return false;
758 }
759 template <typename ValueType>
760 __host__ __device__
761 inline bool operator==(const complex<ValueType> & lhs, const ValueType& rhs){
762 if(lhs.real() == rhs && lhs.imag() == 0){
763 return true;
764 }
765 return false;
766 }
767
768 template <typename ValueType>
769 __host__ __device__
770 inline bool operator!=(const complex<ValueType>& lhs, const complex<ValueType>& rhs){
771 return !(lhs == rhs);
772 }
773
774 template <typename ValueType>
775 __host__ __device__
776 inline bool operator!=(const ValueType & lhs, const complex<ValueType>& rhs){
777 return !(lhs == rhs);
778 }
779
780 template <typename ValueType>
781 __host__ __device__
782 inline bool operator!=(const complex<ValueType> & lhs, const ValueType& rhs){
783 return !(lhs == rhs);
784 }
785
786
787 template <typename ValueType>
788 __host__ __device__
789 inline complex<ValueType> conj(const complex<ValueType>& z){
790 return complex<ValueType>(z.real(),-z.imag());
791 }
792
793 template <typename ValueType>
794 __host__ __device__
795 inline ValueType abs(const complex<ValueType>& z){
796 return ::hypot(z.real(),z.imag());
797 }
798 template <>
799 __host__ __device__
800 inline float abs(const complex<float>& z){
801 return ::hypotf(z.real(),z.imag());
802 }
803 template<>
804 __host__ __device__
805 inline double abs(const complex<double>& z){
806 return ::hypot(z.real(),z.imag());
807 }
808
809 template <typename ValueType>
810 __host__ __device__
811 inline ValueType arg(const complex<ValueType>& z){
812 return atan2(z.imag(),z.real());
813 }
814 template<>
815 __host__ __device__
816 inline float arg(const complex<float>& z){
817 return atan2f(z.imag(),z.real());
818 }
819 template<>
820 __host__ __device__
821 inline double arg(const complex<double>& z){
822 return atan2(z.imag(),z.real());
823 }
824
825 template <typename ValueType>
826 __host__ __device__
827 inline ValueType norm(const complex<ValueType>& z){
828 return z.real()*z.real() + z.imag()*z.imag();
829 }
830
831 template <typename ValueType>
832 __host__ __device__
833 inline complex<ValueType> polar(const ValueType & m, const ValueType & theta){
834 return complex<ValueType>(m * ::cos(theta),m * ::sin(theta));
835 }
836
837 template <>
838 __host__ __device__
839 inline complex<float> polar(const float & magnitude, const float & angle){
840 return complex<float>(magnitude * ::cosf(angle),magnitude * ::sinf(angle));
841 }
842
843 template <>
844 __host__ __device__
845 inline complex<double> polar(const double & magnitude, const double & angle){
846 return complex<double>(magnitude * ::cos(angle),magnitude * ::sin(angle));
847 }
848
849 // Transcendental functions implementation
850 template <typename ValueType>
851 __host__ __device__
852 inline complex<ValueType> cos(const complex<ValueType>& z){
853 const ValueType re = z.real();
854 const ValueType im = z.imag();
855 return complex<ValueType>(::cos(re) * ::cosh(im), -::sin(re) * ::sinh(im));
856 }
857
858 template <>
859 __host__ __device__
860 inline complex<float> cos(const complex<float>& z){
861 const float re = z.real();
862 const float im = z.imag();
863 return complex<float>(cosf(re) * coshf(im), -sinf(re) * sinhf(im));
864 }
865
866 template <typename ValueType>
867 __host__ __device__
868 inline complex<ValueType> cosh(const complex<ValueType>& z){
869 const ValueType re = z.real();
870 const ValueType im = z.imag();
871 return complex<ValueType>(::cosh(re) * ::cos(im), ::sinh(re) * ::sin(im));
872 }
873
874 template <>
875 __host__ __device__
876 inline complex<float> cosh(const complex<float>& z){
877 const float re = z.real();
878 const float im = z.imag();
879 return complex<float>(::coshf(re) * ::cosf(im), ::sinhf(re) * ::sinf(im));
880 }
881
882
883 template <typename ValueType>
884 __host__ __device__
885 inline complex<ValueType> exp(const complex<ValueType>& z){
886 return polar(::exp(z.real()),z.imag());
887 }
888
889 template <>
890 __host__ __device__
891 inline complex<float> exp(const complex<float>& z){
892 return polar(::expf(z.real()),z.imag());
893 }
894
895 template <typename ValueType>
896 __host__ __device__
897 inline complex<ValueType> log(const complex<ValueType>& z){
898 return complex<ValueType>(::log(abs(z)),arg(z));
899 }
900
901 template <>
902 __host__ __device__
903 inline complex<float> log(const complex<float>& z){
904 return complex<float>(::logf(abs(z)),arg(z));
905 }
906
907
908 template <typename ValueType>
909 __host__ __device__
910 inline complex<ValueType> log10(const complex<ValueType>& z){
911 // Using the explicit literal prevents compile time warnings in
912 // devices that don't support doubles
913 return log(z)/ValueType(2.30258509299404568402);
914 // return log(z)/ValueType(::log(10.0));
915 }
916
917 template <typename ValueType>
918 __host__ __device__
919 inline complex<ValueType> pow(const complex<ValueType>& z, const ValueType & exponent){
920 return exp(log(z)*exponent);
921 }
922
923 template <typename ValueType>
924 __host__ __device__
925 inline complex<ValueType> pow(const complex<ValueType>& z, const complex<ValueType> & exponent){
926 return exp(log(z)*exponent);
927 }
928
929 template <typename ValueType>
930 __host__ __device__
931 inline complex<ValueType> pow(const ValueType & x, const complex<ValueType> & exponent){
932 return exp(::log(x)*exponent);
933 }
934
935 template <>
936 __host__ __device__
937 inline complex<float> pow(const float & x, const complex<float> & exponent){
938 return exp(::logf(x)*exponent);
939 }
940
941 template <typename ValueType>
942 __host__ __device__
943 inline complex<ValueType> pow(const complex<ValueType>& z,const int & exponent){
944 return exp(log(z)*ValueType(exponent));
945 }
946
947 template <typename ValueType>
948 __host__ __device__
949 inline complex<ValueType> sin(const complex<ValueType>& z){
950 const ValueType re = z.real();
951 const ValueType im = z.imag();
952 return complex<ValueType>(::sin(re) * ::cosh(im), ::cos(re) * ::sinh(im));
953 }
954
955 template <>
956 __host__ __device__
957 inline complex<float> sin(const complex<float>& z){
958 const float re = z.real();
959 const float im = z.imag();
960 return complex<float>(::sinf(re) * ::coshf(im), ::cosf(re) * ::sinhf(im));
961 }
962
963 template <typename ValueType>
964 __host__ __device__
965 inline complex<ValueType> sinh(const complex<ValueType>& z){
966 const ValueType re = z.real();
967 const ValueType im = z.imag();
968 return complex<ValueType>(::sinh(re) * ::cos(im), ::cosh(re) * ::sin(im));
969 }
970
971 template <>
972 __host__ __device__
973 inline complex<float> sinh(const complex<float>& z){
974 const float re = z.real();
975 const float im = z.imag();
976 return complex<float>(::sinhf(re) * ::cosf(im), ::coshf(re) * ::sinf(im));
977 }
978
979 template <typename ValueType>
980 __host__ __device__
981 inline complex<ValueType> sqrt(const complex<ValueType>& z){
982 return polar(::sqrt(abs(z)),arg(z)/ValueType(2));
983 }
984
985 template <typename ValueType>
986 __host__ __device__
987 inline complex<float> sqrt(const complex<float>& z){
988 return polar(::sqrtf(abs(z)),arg(z)/float(2));
989 }
990
991 template <typename ValueType>
992 __host__ __device__
993 inline complex<ValueType> tan(const complex<ValueType>& z){
994 return sin(z)/cos(z);
995 }
996
997 template <typename ValueType>
998 __host__ __device__
999 inline complex<ValueType> tanh(const complex<ValueType>& z){
1000 // This implementation seems better than the simple sin/cos
1001 return (exp(ValueType(2)*z)-ValueType(1))/(exp(ValueType(2)*z)+ValueType(1));
1002 // return sinh(z)/cosh(z);
1003 }
1004
1005 // Inverse trigonometric functions implementation
1006
1007 template <typename ValueType>
1008 __host__ __device__
1009 inline complex<ValueType> acos(const complex<ValueType>& z){
1010 const complex<ValueType> ret = asin(z);
1011 return complex<ValueType>(ValueType(M_PI/2.0) - ret.real(),-ret.imag());
1012 }
1013
1014 template <typename ValueType>
1015 __host__ __device__
1016 inline complex<ValueType> asin(const complex<ValueType>& z){
1017 const complex<ValueType> i(0,1);
1018 return -i*asinh(i*z);
1019 }
1020
1021 template <typename ValueType>
1022 __host__ __device__
1023 inline complex<ValueType> atan(const complex<ValueType>& z){
1024 const complex<ValueType> i(0,1);
1025 return -i*atanh(i*z);
1026 }
1027
1028 template <typename ValueType>
1029 __host__ __device__
1030 inline complex<ValueType> acosh(const complex<ValueType>& z){
1031 cusp::complex<ValueType> ret((z.real() - z.imag()) * (z.real() + z.imag()) - ValueType(1.0),
1032 ValueType(2.0) * z.real() * z.imag());
1033 ret = sqrt(ret);
1034 if (z.real() < ValueType(0.0)){
1035 ret = -ret;
1036 }
1037 ret += z;
1038 ret = log(ret);
1039 if (ret.real() < ValueType(0.0)){
1040 ret = -ret;
1041 }
1042 return ret;
1043
1044 /*
1045 cusp::complex<ValueType> ret = log(sqrt(z*z-ValueType(1))+z);
1046 if(ret.real() < 0){
1047 ret.real(-ret.real());
1048 }
1049 return ret;
1050 */
1051 }
1052
1053 template <typename ValueType>
1054 __host__ __device__
1055 inline complex<ValueType> asinh(const complex<ValueType>& z){
1056 return log(sqrt(z*z+ValueType(1))+z);
1057 }
1058
1059 template <typename ValueType>
1060 __host__ __device__
1061 inline complex<ValueType> atanh(const complex<ValueType>& z){
1062 ValueType imag2 = z.imag() * z.imag();
1063 ValueType n = ValueType(1.0) + z.real();
1064 n = imag2 + n * n;
1065
1066 ValueType d = ValueType(1.0) - z.real();
1067 d = imag2 + d * d;
1068 complex<ValueType> ret(ValueType(0.25) * (::log(n) - ::log(d)),0);
1069
1070 d = ValueType(1.0) - z.real() * z.real() - imag2;
1071
1072 ret.imag(ValueType(0.5) * ::atan2(ValueType(2.0) * z.imag(), d));
1073 return ret;
1074 //return (log(ValueType(1)+z)-log(ValueType(1)-z))/ValueType(2);
1075 }
1076
1077 template <typename ValueType>
1078 __host__ __device__
1079 inline complex<float> atanh(const complex<float>& z){
1080 float imag2 = z.imag() * z.imag();
1081 float n = float(1.0) + z.real();
1082 n = imag2 + n * n;
1083
1084 float d = float(1.0) - z.real();
1085 d = imag2 + d * d;
1086 complex<float> ret(float(0.25) * (::logf(n) - ::logf(d)),0);
1087
1088 d = float(1.0) - z.real() * z.real() - imag2;
1089
1090 ret.imag(float(0.5) * ::atan2f(float(2.0) * z.imag(), d));
1091 return ret;
1092 //return (log(ValueType(1)+z)-log(ValueType(1)-z))/ValueType(2);
1093
1094 }
1095
1096 } // end namespace cusp
1097
1098 #else
1099 #include <complex>
1100
1101 namespace cusp
1102 {
1103 using std::complex;
1104 using std::conj;
1105 using std::abs;
1106 using std::arg;
1107 using std::norm;
1108 using std::polar;
1109 using std::cos;
1110 using std::cosh;
1111 using std::exp;
1112 using std::log;
1113 using std::log10;
1114 using std::pow;
1115 using std::sin;
1116 using std::sinh;
1117 using std::sqrt;
1118 using std::tan;
1119 using std::tanh;
1120
1121 using std::acos;
1122 using std::asin;
1123 using std::atan;
1124 #if __cplusplus >= 201103L
1125 using std::acosh;
1126 using std::asinh;
1127 using std::atanh;
1128 #else
1129 template <typename ValueType>
1130 inline complex<ValueType> acosh(const complex<ValueType>& z){
1131 cusp::complex<ValueType> ret((z.real() - z.imag()) * (z.real() + z.imag()) - ValueType(1.0),
1132 ValueType(2.0) * z.real() * z.imag());
1133 ret = sqrt(ret);
1134 if (z.real() < ValueType(0.0)){
1135 ret = -ret;
1136 }
1137 ret += z;
1138 ret = log(ret);
1139 if (ret.real() < ValueType(0.0)){
1140 ret = -ret;
1141 }
1142 return ret;
1143 }
1144 template <typename ValueType>
1145 inline complex<ValueType> asinh(const complex<ValueType>& z){
1146 return log(sqrt(z*z+ValueType(1))+z);
1147 }
1148
1149 template <typename ValueType>
1150 inline complex<ValueType> atanh(const complex<ValueType>& z){
1151 ValueType imag2 = z.imag() * z.imag();
1152 ValueType n = ValueType(1.0) + z.real();
1153 n = imag2 + n * n;
1154
1155 ValueType d = ValueType(1.0) - z.real();
1156 d = imag2 + d * d;
1157 complex<ValueType> ret(ValueType(0.25) * (::log(n) - ::log(d)),0);
1158
1159 d = ValueType(1.0) - z.real() * z.real() - imag2;
1160
1161 ret.imag(ValueType(0.5) * ::atan2(ValueType(2.0) * z.imag(), d));
1162 return ret;
1163 }
1164 #endif
1165
1166
1167 template <typename T>
1168 struct norm_type {
1169 typedef T type;
1170 };
1171
1172 template <typename T>
1173 struct norm_type< complex<T> > {
1174 typedef T type;
1175 };
1176 }
1177 #endif

  ViewVC Help
Powered by ViewVC 1.1.26