/[escript]/branches/diaplayground/cusplibrary/cusp/complex.h
ViewVC logotype

Annotation of /branches/diaplayground/cusplibrary/cusp/complex.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4955 - (hide annotations)
Tue May 20 04:33:15 2014 UTC (6 years, 2 months ago) by caltinay
File MIME type: text/plain
File size: 35442 byte(s)
added pristine copy of cusplibrary (apache license) to be used by ripley.

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

  ViewVC Help
Powered by ViewVC 1.1.26