summaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/include/cctk_Complex.h125
-rw-r--r--src/include/cctk_Types.h18
-rw-r--r--src/main/Complex.c441
3 files changed, 467 insertions, 117 deletions
diff --git a/src/include/cctk_Complex.h b/src/include/cctk_Complex.h
index 041884f3..977896d3 100644
--- a/src/include/cctk_Complex.h
+++ b/src/include/cctk_Complex.h
@@ -12,27 +12,92 @@
#define _CCTK_COMPLEX_H_
#ifdef __cplusplus
-extern "C"
-{
+# define EXTERN_C_BEGIN extern "C" {
+# define EXTERN_C_END }
+#else
+# define EXTERN_C_BEGIN
+# define EXTERN_C_END
+#endif
+
+
+#ifndef DEFINE_CCTK_COMPLEX_EXTERN_FUNCTIONS
+# define PREFIX static inline
+#else
+# define PREFIX
+#endif
+
+
+#ifdef __cplusplus
+ /* declare C++ overloaded operators */
+# include <complex>
+# include <ostream>
+# define DECLARE_CMPLX_CXX_OPERATORS(CCTK_Cmplx, cctk_real, cctk_complex) \
+PREFIX cctk_complex operator+ (cctk_complex const & a); \
+PREFIX cctk_complex operator- (cctk_complex const & a); \
+PREFIX cctk_complex conj (cctk_complex const & a); \
+PREFIX cctk_real abs (cctk_complex const & a); \
+PREFIX cctk_real arg (cctk_complex const & a); \
+PREFIX cctk_real norm (cctk_complex const & a); \
+PREFIX cctk_complex operator+ (cctk_complex const & a, cctk_complex const & b); \
+PREFIX cctk_complex operator+ (cctk_real const & a, cctk_complex const & b); \
+PREFIX cctk_complex operator+ (cctk_complex const & a, cctk_real const & b); \
+PREFIX cctk_complex operator- (cctk_complex const & a, cctk_complex const & b); \
+PREFIX cctk_complex operator- (cctk_real const & a, cctk_complex const & b); \
+PREFIX cctk_complex operator- (cctk_complex const & a, cctk_real const & b); \
+PREFIX cctk_complex operator* (cctk_complex const & a, cctk_complex const & b); \
+PREFIX cctk_complex operator* (cctk_real const & a, cctk_complex const & b); \
+PREFIX cctk_complex operator* (cctk_complex const & a, cctk_real const & b); \
+PREFIX cctk_complex operator/ (cctk_complex const & a, cctk_complex const & b); \
+PREFIX cctk_complex operator/ (cctk_real const & a, cctk_complex const & b); \
+PREFIX cctk_complex operator/ (cctk_complex const & a, cctk_real const & b); \
+PREFIX cctk_complex operator+= (cctk_complex & a, cctk_complex const & b); \
+PREFIX cctk_complex operator+= (cctk_complex & a, cctk_real const & b); \
+PREFIX cctk_complex operator-= (cctk_complex & a, cctk_complex const & b); \
+PREFIX cctk_complex operator-= (cctk_complex & a, cctk_real const & b); \
+PREFIX cctk_complex operator*= (cctk_complex & a, cctk_complex const & b); \
+PREFIX cctk_complex operator*= (cctk_complex & a, cctk_real const & b); \
+PREFIX cctk_complex operator/= (cctk_complex & a, cctk_complex const & b); \
+PREFIX cctk_complex operator/= (cctk_complex & a, cctk_real const & b); \
+PREFIX cctk_complex pow (cctk_complex const & a, int i); \
+PREFIX cctk_complex pow (cctk_complex const & a, cctk_real const & r); \
+PREFIX cctk_complex pow (cctk_complex const & a, cctk_complex const & b); \
+PREFIX cctk_complex sin (cctk_complex const & a); \
+PREFIX cctk_complex cos (cctk_complex const & a); \
+PREFIX cctk_complex exp (cctk_complex const & a); \
+PREFIX cctk_complex log (cctk_complex const & a); \
+PREFIX cctk_complex sqrt (cctk_complex const & a); \
+PREFIX std::ostream & operator << (std::ostream & os, cctk_complex const & a);
+#else
+ /* declare no C++ overloaded operators */
+# define DECLARE_CMPLX_CXX_OPERATORS(CCTK_Cmplx, cctk_real, cctk_complex)
#endif
/* macro to declare a set of complex functions for a given precision */
-#define DECLARE_CMPLX_FUNCTIONS(CCTK_Cmplx, cctk_real, cctk_complex) \
-cctk_complex CCTK_Cmplx (cctk_real Re, cctk_real Im); \
-cctk_real CCTK_Cmplx##Real (cctk_complex complex_number); \
-cctk_real CCTK_Cmplx##Imag (cctk_complex complex_number); \
-cctk_complex CCTK_Cmplx##Conjg (cctk_complex complex_number); \
-cctk_real CCTK_Cmplx##Abs (cctk_complex complex_number); \
-cctk_complex CCTK_Cmplx##Add (cctk_complex a, cctk_complex b); \
-cctk_complex CCTK_Cmplx##Sub (cctk_complex a, cctk_complex b); \
-cctk_complex CCTK_Cmplx##Mul (cctk_complex a, cctk_complex b); \
-cctk_complex CCTK_Cmplx##Div (cctk_complex a, cctk_complex b); \
-cctk_complex CCTK_Cmplx##Sin (cctk_complex complex_number); \
-cctk_complex CCTK_Cmplx##Cos (cctk_complex complex_number); \
-cctk_complex CCTK_Cmplx##Exp (cctk_complex complex_number); \
-cctk_complex CCTK_Cmplx##Sqrt (cctk_complex complex_number); \
-cctk_complex CCTK_Cmplx##Pow (cctk_complex complex_number, cctk_real w);
+#define DECLARE_CMPLX_FUNCTIONS(CCTK_Cmplx, cctk_real, cctk_complex) \
+EXTERN_C_BEGIN \
+PREFIX cctk_complex CCTK_Cmplx (cctk_real Re, cctk_real Im); \
+PREFIX cctk_real CCTK_Cmplx##Real (cctk_complex a); \
+PREFIX cctk_real CCTK_Cmplx##Imag (cctk_complex a); \
+PREFIX cctk_complex CCTK_Cmplx##Neg (cctk_complex a); \
+PREFIX cctk_complex CCTK_Cmplx##Conjg (cctk_complex a); \
+PREFIX cctk_real CCTK_Cmplx##Abs (cctk_complex a); \
+PREFIX cctk_real CCTK_Cmplx##Arg (cctk_complex a); \
+PREFIX cctk_real CCTK_Cmplx##Norm (cctk_complex a); \
+PREFIX cctk_complex CCTK_Cmplx##Add (cctk_complex a, cctk_complex b); \
+PREFIX cctk_complex CCTK_Cmplx##Sub (cctk_complex a, cctk_complex b); \
+PREFIX cctk_complex CCTK_Cmplx##Mul (cctk_complex a, cctk_complex b); \
+PREFIX cctk_complex CCTK_Cmplx##Div (cctk_complex a, cctk_complex b); \
+PREFIX cctk_complex CCTK_Cmplx##CPow (cctk_complex a, cctk_complex b); \
+PREFIX cctk_complex CCTK_Cmplx##Sin (cctk_complex a); \
+PREFIX cctk_complex CCTK_Cmplx##Cos (cctk_complex a); \
+PREFIX cctk_complex CCTK_Cmplx##Exp (cctk_complex a); \
+PREFIX cctk_complex CCTK_Cmplx##Log (cctk_complex a); \
+PREFIX cctk_complex CCTK_Cmplx##Sqrt (cctk_complex a); \
+PREFIX cctk_complex CCTK_Cmplx##Pow (cctk_complex a, cctk_real b); \
+PREFIX cctk_complex CCTK_Cmplx##IPow (cctk_complex a, int b); \
+EXTERN_C_END \
+DECLARE_CMPLX_CXX_OPERATORS(CCTK_Cmplx, cctk_real, cctk_complex)
/* declare complex functions for all available precisions */
@@ -48,57 +113,79 @@ DECLARE_CMPLX_FUNCTIONS (CCTK_Cmplx16, CCTK_REAL8, CCTK_COMPLEX16)
DECLARE_CMPLX_FUNCTIONS (CCTK_Cmplx32, CCTK_REAL16, CCTK_COMPLEX32)
#endif
+#undef PREFIX
+
/* declare the default precision complex functions as #define'd macros */
#ifdef CCTK_REAL_PRECISION_4
#define CCTK_Cmplx CCTK_Cmplx8
#define CCTK_CmplxReal CCTK_Cmplx8Real
#define CCTK_CmplxImag CCTK_Cmplx8Imag
+#define CCTK_CmplxNeg CCTK_Cmplx8Neg
#define CCTK_CmplxConjg CCTK_Cmplx8Conjg
#define CCTK_CmplxAbs CCTK_Cmplx8Abs
+#define CCTK_CmplxArg CCTK_Cmplx8Arg
+#define CCTK_CmplxNorm CCTK_Cmplx8Norm
#define CCTK_CmplxAdd CCTK_Cmplx8Add
#define CCTK_CmplxSub CCTK_Cmplx8Sub
#define CCTK_CmplxMul CCTK_Cmplx8Mul
#define CCTK_CmplxDiv CCTK_Cmplx8Div
+#define CCTK_CmplxCPow CCTK_Cmplx8CPow
#define CCTK_CmplxSin CCTK_Cmplx8Sin
#define CCTK_CmplxCos CCTK_Cmplx8Cos
#define CCTK_CmplxExp CCTK_Cmplx8Exp
+#define CCTK_CmplxLog CCTK_Cmplx8Log
#define CCTK_CmplxSqrt CCTK_Cmplx8Sqrt
#define CCTK_CmplxPow CCTK_Cmplx8Pow
+#define CCTK_CmplxIPow CCTK_Cmplx8IPow
#elif CCTK_REAL_PRECISION_8
#define CCTK_Cmplx CCTK_Cmplx16
#define CCTK_CmplxReal CCTK_Cmplx16Real
#define CCTK_CmplxImag CCTK_Cmplx16Imag
+#define CCTK_CmplxNeg CCTK_Cmplx16Neg
#define CCTK_CmplxConjg CCTK_Cmplx16Conjg
#define CCTK_CmplxAbs CCTK_Cmplx16Abs
+#define CCTK_CmplxArg CCTK_Cmplx16Arg
+#define CCTK_CmplxNorm CCTK_Cmplx16Norm
#define CCTK_CmplxAdd CCTK_Cmplx16Add
#define CCTK_CmplxSub CCTK_Cmplx16Sub
#define CCTK_CmplxMul CCTK_Cmplx16Mul
#define CCTK_CmplxDiv CCTK_Cmplx16Div
+#define CCTK_CmplxCPow CCTK_Cmplx16CPow
#define CCTK_CmplxSin CCTK_Cmplx16Sin
#define CCTK_CmplxCos CCTK_Cmplx16Cos
#define CCTK_CmplxExp CCTK_Cmplx16Exp
+#define CCTK_CmplxLog CCTK_Cmplx16Log
#define CCTK_CmplxSqrt CCTK_Cmplx16Sqrt
#define CCTK_CmplxPow CCTK_Cmplx16Pow
+#define CCTK_CmplxIPow CCTK_Cmplx16IPow
#elif CCTK_REAL_PRECISION_16
#define CCTK_Cmplx CCTK_Cmplx32
#define CCTK_CmplxReal CCTK_Cmplx32Real
#define CCTK_CmplxImag CCTK_Cmplx32Imag
+#define CCTK_CmplxNeg CCTK_Cmplx32Neg
#define CCTK_CmplxConjg CCTK_Cmplx32Conjg
#define CCTK_CmplxAbs CCTK_Cmplx32Abs
+#define CCTK_CmplxArg CCTK_Cmplx32Arg
+#define CCTK_CmplxNorm CCTK_Cmplx32Norm
#define CCTK_CmplxAdd CCTK_Cmplx32Add
#define CCTK_CmplxSub CCTK_Cmplx32Sub
#define CCTK_CmplxMul CCTK_Cmplx32Mul
#define CCTK_CmplxDiv CCTK_Cmplx32Div
+#define CCTK_CmplxCPow CCTK_Cmplx32CPow
#define CCTK_CmplxSin CCTK_Cmplx32Sin
#define CCTK_CmplxCos CCTK_Cmplx32Cos
#define CCTK_CmplxExp CCTK_Cmplx32Exp
+#define CCTK_CmplxLog CCTK_Cmplx32Log
#define CCTK_CmplxSqrt CCTK_Cmplx32Sqrt
#define CCTK_CmplxPow CCTK_Cmplx32Pow
+#define CCTK_CmplxIPow CCTK_Cmplx32IPow
#endif
-#ifdef __cplusplus
-}
+#ifndef DEFINE_CCTK_COMPLEX_EXTERN_FUNCTIONS
+# define DEFINE_CCTK_COMPLEX_INLINE_FUNCTIONS
+# include "../main/Complex.c"
+# undef DEFINE_CCTK_COMPLEX_INLINE_FUNCTIONS
#endif
#endif /* _CCTK_COMPLEX_H_ */
diff --git a/src/include/cctk_Types.h b/src/include/cctk_Types.h
index b6d20771..734a0c17 100644
--- a/src/include/cctk_Types.h
+++ b/src/include/cctk_Types.h
@@ -38,28 +38,40 @@ typedef const char * CCTK_STRING;
#ifdef HAVE_CCTK_REAL16
#define HAVE_CCTK_COMPLEX32 1
-typedef struct
+typedef struct CCTK_COMPLEX32
{
CCTK_REAL16 Re;
CCTK_REAL16 Im;
+#ifdef __cplusplus
+ CCTK_REAL16 real() const { return Re; }
+ CCTK_REAL16 imag() const { return Im; }
+#endif
} CCTK_COMPLEX32;
#endif
#ifdef HAVE_CCTK_REAL8
#define HAVE_CCTK_COMPLEX16 1
-typedef struct
+typedef struct CCTK_COMPLEX16
{
CCTK_REAL8 Re;
CCTK_REAL8 Im;
+#ifdef __cplusplus
+ CCTK_REAL8 real() const { return Re; }
+ CCTK_REAL8 imag() const { return Im; }
+#endif
} CCTK_COMPLEX16;
#endif
#ifdef HAVE_CCTK_REAL4
#define HAVE_CCTK_COMPLEX8 1
-typedef struct
+typedef struct CCTK_COMPLEX8
{
CCTK_REAL4 Re;
CCTK_REAL4 Im;
+#ifdef __cplusplus
+ CCTK_REAL4 real() const { return Re; }
+ CCTK_REAL4 imag() const { return Im; }
+#endif
} CCTK_COMPLEX8;
#endif
diff --git a/src/main/Complex.c b/src/main/Complex.c
index 148ce60b..47e3a888 100644
--- a/src/main/Complex.c
+++ b/src/main/Complex.c
@@ -11,12 +11,22 @@
#include <math.h>
#include "cctk_Flesh.h"
-#include "cctk_Complex.h"
+
+#ifndef DEFINE_CCTK_COMPLEX_INLINE_FUNCTIONS
+# define DEFINE_CCTK_COMPLEX_EXTERN_FUNCTIONS
+# include "cctk_Complex.h"
+# undef DEFINE_CCTK_COMPLEX_EXTERN_FUNCTIONS
+#endif
+
+
+#ifndef DEFINE_CCTK_COMPLEX_INLINE_FUNCTIONS
static const char *rcsid = "$Header$";
CCTK_FILEVERSION(main_Complex_c);
+#endif
+
/********************************************************************
********************* Local Data Types ***********************
@@ -26,7 +36,6 @@ CCTK_FILEVERSION(main_Complex_c);
********************* Local Routine Prototypes *********************
********************************************************************/
-#define NORM(cx) sqrt((cx.Re) * (cx.Re) + (cx.Im) * (cx.Im))
/********************************************************************
********************* Other Routine Prototypes *********************
********************************************************************/
@@ -129,6 +138,37 @@ cctk_real CCTK_Cmplx##Imag (cctk_complex complex_number) \
/*@@
+ @routine CCTK_CmplxNeg
+ @date 2006-06-07
+ @author Erik Schnetter
+ @desc
+ Returns the negative of a complex number.
+ @enddesc
+
+ @var in
+ @vdesc The complex number
+ @vtype CCTK_COMPLEX
+ @vio in
+ @endvar
+
+ @returntype CCTK_COMPLEX
+ @returndesc
+ The negative
+ @endreturndesc
+@@*/
+#define DEFINE_CCTK_CMPLX_NEG(CCTK_Cmplx, cctk_real, cctk_complex) \
+cctk_complex CCTK_Cmplx##Neg (cctk_complex complex_number) \
+{ \
+ cctk_complex result; \
+ \
+ \
+ result.Re = -complex_number.Re; \
+ result.Im = -complex_number.Im; \
+ return (result); \
+}
+
+
+ /*@@
@routine CCTK_CmplxConjg
@date Tue Dec 14 12:22:41 1999
@author Tom Goodale
@@ -178,10 +218,62 @@ cctk_complex CCTK_Cmplx##Conjg (cctk_complex complex_number) \
The absolute value of the complex number
@endreturndesc
@@*/
-#define DEFINE_CCTK_CMPLX_ABS(CCTK_Cmplx, cctk_real, cctk_complex) \
+#define DEFINE_CCTK_CMPLX_ABS(CCTK_Cmplx, cctk_real, cctk_complex, type) \
cctk_real CCTK_Cmplx##Abs (cctk_complex complex_number) \
{ \
- return ((cctk_real)hypot (complex_number.Re, complex_number.Im)); \
+ return ((cctk_real)hypot##type (complex_number.Re, complex_number.Im)); \
+}
+
+
+ /*@@
+ @routine CCTK_CmplxNorm
+ @date 2006-07-06
+ @author Erik Schnetter
+ @desc
+ Return the absolute value squared of a complex number.
+ @enddesc
+
+ @var in
+ @vdesc The complex number
+ @vtype CCTK_COMPLEX
+ @vio in
+ @endvar
+
+ @returntype CCTK_REAL
+ @returndesc
+ The absolute value squared of the complex number
+ @endreturndesc
+@@*/
+#define DEFINE_CCTK_CMPLX_NORM(CCTK_Cmplx, cctk_real, cctk_complex, type) \
+cctk_real CCTK_Cmplx##Norm (cctk_complex complex_number) \
+{ \
+ return (pow##type (complex_number.Re, 2) + pow##type (complex_number.Im, 2));\
+}
+
+
+ /*@@
+ @routine CCTK_CmplxArg
+ @date 2005-11-17
+ @author Erik Schnetter
+ @desc
+ Return the argument of a complex number.
+ @enddesc
+
+ @var in
+ @vdesc The complex number
+ @vtype CCTK_COMPLEX
+ @vio in
+ @endvar
+
+ @returntype CCTK_REAL
+ @returndesc
+ The argument of the complex number
+ @endreturndesc
+@@*/
+#define DEFINE_CCTK_CMPLX_ARG(CCTK_Cmplx, cctk_real, cctk_complex, type) \
+cctk_real CCTK_Cmplx##Arg (cctk_complex complex_number) \
+{ \
+ return (atan2##type (complex_number.Im, complex_number.Re)); \
}
@@ -321,17 +413,17 @@ cctk_complex CCTK_Cmplx##Mul (cctk_complex a, cctk_complex b) \
The quotient
@endreturndesc
@@*/
-#define DEFINE_CCTK_CMPLX_DIV(CCTK_Cmplx, cctk_real, cctk_complex) \
+#define DEFINE_CCTK_CMPLX_DIV(CCTK_Cmplx, cctk_real, cctk_complex, type) \
cctk_complex CCTK_Cmplx##Div (cctk_complex a, cctk_complex b) \
{ \
cctk_real afact, bfact, factor; \
cctk_complex result; \
\
\
- afact = (cctk_real)(fabs(a.Re) + fabs(a.Im)); \
+ afact = (cctk_real)(fabs##type(a.Re) + fabs##type(a.Im)); \
a.Re /= afact; \
a.Im /= afact; \
- bfact = (cctk_real)(fabs(b.Re) + fabs(b.Im)); \
+ bfact = (cctk_real)(fabs##type(b.Re) + fabs##type(b.Im)); \
b.Re /= bfact; \
b.Im /= bfact; \
factor = b.Re*b.Re + b.Im*b.Im; \
@@ -343,6 +435,39 @@ cctk_complex CCTK_Cmplx##Div (cctk_complex a, cctk_complex b) \
/*@@
+ @routine CCTK_CmplxCPow
+ @date 2005-11-17
+ @author Erik Schnetter
+ @desc
+ Raises a complex number to a given power
+ This algorithm was taken from glibc 2.3.5,
+ file sysdeps/generic/s_cpow.c.
+ @enddesc
+
+ @var a
+ @vdesc The base
+ @vtype CCTK_COMPLEX
+ @vio in
+ @endvar
+ @var b
+ @vdesc The exponent
+ @vtype CCTK_COMPLEX
+ @vio in
+ @endvar
+
+ @returntype CCTK_COMPLEX
+ @returndesc
+ a to the power of b
+ @endreturndesc
+@@*/
+#define DEFINE_CCTK_CMPLX_CPOW(CCTK_Cmplx, cctk_real, cctk_complex) \
+cctk_complex CCTK_Cmplx##CPow (cctk_complex a, cctk_complex b) \
+{ \
+ return CCTK_Cmplx##Exp (CCTK_Cmplx##Mul (b, CCTK_Cmplx##Log (a))); \
+}
+
+
+ /*@@
@routine CCTK_CmplxSin
@date Wed 12 Dec 2001
@author Thomas Radke
@@ -361,24 +486,14 @@ cctk_complex CCTK_Cmplx##Div (cctk_complex a, cctk_complex b) \
The sine
@endreturndesc
@@*/
-#define DEFINE_CCTK_CMPLX_SIN(CCTK_Cmplx, cctk_real, cctk_complex) \
+#define DEFINE_CCTK_CMPLX_SIN(CCTK_Cmplx, cctk_real, cctk_complex, type) \
cctk_complex CCTK_Cmplx##Sin (cctk_complex complex_number) \
{ \
cctk_complex result; \
\
\
- if (complex_number.Im == (cctk_real)0.0) \
- { \
- result.Re = (cctk_real)sin (complex_number.Re); \
- result.Im = (cctk_real)0.0; \
- } \
- else \
- { \
- result.Re = (cctk_real)(sin (complex_number.Re) \
- * cosh (complex_number.Im)); \
- result.Im = (cctk_real)(cos (complex_number.Re) \
- * sinh (complex_number.Im)); \
- } \
+ result.Re = sin##type (complex_number.Re) * cosh##type (complex_number.Im); \
+ result.Im = cos##type (complex_number.Re) * sinh##type (complex_number.Im); \
\
return (result); \
}
@@ -403,24 +518,14 @@ cctk_complex CCTK_Cmplx##Sin (cctk_complex complex_number) \
The cosine
@endreturndesc
@@*/
-#define DEFINE_CCTK_CMPLX_COS(CCTK_Cmplx, cctk_real, cctk_complex) \
+#define DEFINE_CCTK_CMPLX_COS(CCTK_Cmplx, cctk_real, cctk_complex, type) \
cctk_complex CCTK_Cmplx##Cos (cctk_complex complex_number) \
{ \
cctk_complex result; \
\
\
- if (complex_number.Im == (cctk_real)0.0) \
- { \
- result.Re = (cctk_real)cos (complex_number.Re); \
- result.Im = (cctk_real)0.0; \
- } \
- else \
- { \
- result.Re = (cctk_real)(cos (complex_number.Re) \
- * cosh (complex_number.Im)); \
- result.Im = (cctk_real)(sin (complex_number.Re) \
- * sinh (complex_number.Im)); \
- } \
+ result.Re = cos##type (complex_number.Re) * cosh##type (complex_number.Im); \
+ result.Im = sin##type (complex_number.Re) * sinh##type (complex_number.Im); \
\
return (result); \
}
@@ -445,17 +550,51 @@ cctk_complex CCTK_Cmplx##Cos (cctk_complex complex_number) \
The exponential
@endreturndesc
@@*/
-#define DEFINE_CCTK_CMPLX_EXP(CCTK_Cmplx, cctk_real, cctk_complex) \
+#define DEFINE_CCTK_CMPLX_EXP(CCTK_Cmplx, cctk_real, cctk_complex, type) \
cctk_complex CCTK_Cmplx##Exp (cctk_complex complex_number) \
{ \
cctk_real rho, theta; \
cctk_complex result; \
\
\
- rho = (cctk_real)exp (complex_number.Re); \
+ rho = (cctk_real)exp##type (complex_number.Re); \
theta = complex_number.Im; \
- result.Re = rho * (cctk_real)cos (theta); \
- result.Im = rho * (cctk_real)sin (theta); \
+ result.Re = rho * (cctk_real)cos##type (theta); \
+ result.Im = rho * (cctk_real)sin##type (theta); \
+ \
+ return (result); \
+}
+
+
+ /*@@
+ @routine CCTK_CmplxLog
+ @date 2005-11-17
+ @author Erik Schnetter
+ @desc
+ Returns the logarithm of a complex number.
+ This algorithm was taken from glibc 2.3.5,
+ file sysdeps/generic/s_clog.c.
+ @enddesc
+
+ @var complex_number
+ @vdesc The complex number
+ @vtype CCTK_COMPLEX
+ @vio in
+ @endvar
+
+ @returntype CCTK_COMPLEX
+ @returndesc
+ The logarithm
+ @endreturndesc
+@@*/
+#define DEFINE_CCTK_CMPLX_LOG(CCTK_Cmplx, cctk_real, cctk_complex, type) \
+cctk_complex CCTK_Cmplx##Log (cctk_complex complex_number) \
+{ \
+ cctk_complex result; \
+ \
+ \
+ result.Re = log##type (hypot##type (complex_number.Re, complex_number.Im)); \
+ result.Im = atan2##type (complex_number.Im, complex_number.Re); \
\
return (result); \
}
@@ -467,6 +606,8 @@ cctk_complex CCTK_Cmplx##Exp (cctk_complex complex_number) \
@author Thomas Radke
@desc
Returns the square root of a complex number.
+ This algorithm was taken from glibc 2.3.5,
+ file sysdeps/generic/s_csqrt.c.
@enddesc
@history
This code bears a resemblance to that in the GSL. That
@@ -484,55 +625,42 @@ cctk_complex CCTK_Cmplx##Exp (cctk_complex complex_number) \
The square root
@endreturndesc
@@*/
-#define DEFINE_CCTK_CMPLX_SQRT(CCTK_Cmplx, cctk_real, cctk_complex) \
+#define DEFINE_CCTK_CMPLX_SQRT(CCTK_Cmplx, cctk_real, cctk_complex, type) \
cctk_complex CCTK_Cmplx##Sqrt (cctk_complex complex_number) \
{ \
- cctk_real x, y, w, t; \
+ cctk_real d, r, s; \
cctk_complex result; \
\
\
- if (complex_number.Re == (cctk_real)0.0 \
- && complex_number.Im == (cctk_real)0.0) \
+ d = hypot##type (complex_number.Re, complex_number.Im); \
+ /* Use the identity 2 Re res Im res = Im x \
+ to avoid cancellation error in d +/- Re x. */ \
+ if (complex_number.Re > 0) \
{ \
- result.Re = result.Im = (cctk_real)0.0; \
+ r = sqrt##type (0.5##type * d + 0.5##type * complex_number.Re); \
+ s = (0.5##type * complex_number.Im) / r; \
} \
else \
{ \
- x = (cctk_real)fabs (complex_number.Re); \
- y = (cctk_real)fabs (complex_number.Im); \
- if (x >= y) \
- { \
- t = y / x; \
- w = (cctk_real)(sqrt (x) * sqrt (0.5 * (1.0 + sqrt (1.0 * t * t)))); \
- } \
- else \
- { \
- t = x / y; \
- w = (cctk_real)(sqrt (y) * sqrt (0.5 * (t + sqrt (1.0 * t * t)))); \
- } \
- \
- if (complex_number.Re >= (cctk_real)0.0) \
- { \
- result.Re = w; \
- result.Im = complex_number.Im / ((cctk_real)2.0 * w); \
- } \
- else \
- { \
- x = complex_number.Im >= (cctk_real)0.0 ? w : -w; \
- result.Re = complex_number.Im / ((cctk_real)2.0 * x); \
- result.Im = x; \
- } \
+ s = sqrt##type (0.5##type * d - 0.5##type * complex_number.Re); \
+ r = fabs##type ((0.5##type * complex_number.Im) / s); \
} \
\
+ result.Re = r; \
+ result.Im = copysign##type (s, complex_number.Im); \
+ \
return (result); \
}
+
/*@@
@routine CCTK_CmplxPow
@date
@author Yaakoub Y El Khamra
@desc
Raises a complex number to a given power
+ This algorithm was taken from glibc 2.3.5,
+ file sysdeps/generic/s_cpow.c.
@enddesc
@var cx_number
@@ -540,61 +668,184 @@ cctk_complex CCTK_Cmplx##Sqrt (cctk_complex complex_number) \
@vtype CCTK_COMPLEX
@vio in
@endvar
-
@var w
- @vdesc The power
+ @vdesc The exponent
@vtype CCTK_REAL
@vio in
@endvar
@returntype CCTK_COMPLEX
@returndesc
- The given power of the number
+ The power of the complex number
@endreturndesc
@@*/
#define DEFINE_CCTK_CMPLX_POW(CCTK_Cmplx, cctk_real, cctk_complex) \
-cctk_complex CCTK_Cmplx##Pow (cctk_complex cx_number, cctk_real w) \
+cctk_complex CCTK_Cmplx##Pow (cctk_complex complex_number, cctk_real w) \
{ \
cctk_complex result; \
\
- cctk_real phi = (cctk_real)atan2 (cx_number.Im, cx_number.Re) * w; \
- cctk_real r = (cctk_real)pow (NORM (cx_number), w); \
\
- result.Re = r * (cctk_real)cos (phi); \
- result.Im = r * (cctk_real)sin (phi); \
+ result = CCTK_Cmplx##Log (complex_number); \
+ result.Re *= w; \
+ result.Im *= w; \
+ result = CCTK_Cmplx##Exp (result); \
\
return (result); \
}
+
+ /*@@
+ @routine CCTK_CmplxIPow
+ @date
+ @author Erik Schnetter
+ @desc
+ Raises a complex number to a given power
+ This algorithm was taken from glibc 2.3.5,
+ file sysdeps/generic/s_cpow.c.
+ @enddesc
+
+ @var cx_number
+ @vdesc The complex number
+ @vtype CCTK_COMPLEX
+ @vio in
+ @endvar
+ @var w
+ @vdesc The exponent
+ @vtype int
+ @vio in
+ @endvar
+
+ @returntype CCTK_COMPLEX
+ @returndesc
+ The power of the complex number
+ @endreturndesc
+@@*/
+#define DEFINE_CCTK_CMPLX_IPOW(CCTK_Cmplx, cctk_real, cctk_complex) \
+cctk_complex CCTK_Cmplx##IPow (cctk_complex complex_number, int w) \
+{ \
+ cctk_complex result; \
+ \
+ \
+ result = CCTK_Cmplx (1, 0); \
+ if (w < 0) \
+ { \
+ result = CCTK_Cmplx##Div (result, \
+ CCTK_Cmplx##IPow (complex_number, - w)); \
+ } \
+ else \
+ { \
+ while (w) \
+ { \
+ if (w % 2) \
+ { \
+ result = CCTK_Cmplx##Mul (result, complex_number); \
+ } \
+ w /= 2; \
+ complex_number = CCTK_Cmplx##Mul (complex_number, complex_number); \
+ } \
+ } \
+ \
+ return (result); \
+}
+
+
+
+#ifdef DEFINE_CCTK_COMPLEX_INLINE_FUNCTIONS
+# define PREFIX static inline
+#else
+# define PREFIX
+#endif
+
+
+#ifdef __cplusplus
+ /* define C++ overloaded operators */
+# define DEFINE_CMPLX_CXX_OPERATORS(CCTK_Cmplx, cctk_real, cctk_complex) \
+PREFIX cctk_complex operator+ (cctk_complex const & a) { return a; } \
+PREFIX cctk_complex operator- (cctk_complex const & a) { return CCTK_Cmplx##Neg (a); } \
+PREFIX cctk_complex conj (cctk_complex const & a) { return CCTK_Cmplx##Conjg (a); } \
+PREFIX cctk_real abs (cctk_complex const & a) { return CCTK_Cmplx##Abs (a); } \
+PREFIX cctk_real arg (cctk_complex const & a) { return CCTK_Cmplx##Arg (a); } \
+PREFIX cctk_real norm (cctk_complex const & a) { return CCTK_Cmplx##Norm (a); } \
+PREFIX cctk_complex operator+ (cctk_complex const & a, cctk_complex const & b) { return CCTK_Cmplx##Add (a, b); } \
+PREFIX cctk_complex operator+ (cctk_real const & a, cctk_complex const & b) { return CCTK_Cmplx##Add (CCTK_Cmplx (a, 0), b); } \
+PREFIX cctk_complex operator+ (cctk_complex const & a, cctk_real const & b) { return CCTK_Cmplx##Add (a, CCTK_Cmplx (b, 0)); } \
+PREFIX cctk_complex operator- (cctk_complex const & a, cctk_complex const & b) { return CCTK_Cmplx##Sub (a, b); } \
+PREFIX cctk_complex operator- (cctk_real const & a, cctk_complex const & b) { return CCTK_Cmplx##Sub (CCTK_Cmplx (a, 0), b); } \
+PREFIX cctk_complex operator- (cctk_complex const & a, cctk_real const & b) { return CCTK_Cmplx##Sub (a, CCTK_Cmplx (b, 0)); } \
+PREFIX cctk_complex operator* (cctk_complex const & a, cctk_complex const & b) { return CCTK_Cmplx##Mul (a, b); } \
+PREFIX cctk_complex operator* (cctk_real const & a, cctk_complex const & b) { return CCTK_Cmplx##Mul (CCTK_Cmplx (a, 0), b); } \
+PREFIX cctk_complex operator* (cctk_complex const & a, cctk_real const & b) { return CCTK_Cmplx##Mul (a, CCTK_Cmplx (b, 0)); } \
+PREFIX cctk_complex operator/ (cctk_complex const & a, cctk_complex const & b) { return CCTK_Cmplx##Div (a, b); } \
+PREFIX cctk_complex operator/ (cctk_real const & a, cctk_complex const & b) { return CCTK_Cmplx##Div (CCTK_Cmplx (a, 0), b); } \
+PREFIX cctk_complex operator/ (cctk_complex const & a, cctk_real const & b) { return CCTK_Cmplx##Div (a, CCTK_Cmplx (b, 0)); } \
+PREFIX cctk_complex operator+= (cctk_complex & a, cctk_complex const & b) { return a = CCTK_Cmplx##Add (a, b); } \
+PREFIX cctk_complex operator+= (cctk_complex & a, cctk_real const & b) { return a = CCTK_Cmplx##Add (a, CCTK_Cmplx (b, 0)); } \
+PREFIX cctk_complex operator-= (cctk_complex & a, cctk_complex const & b) { return a = CCTK_Cmplx##Sub (a, b); } \
+PREFIX cctk_complex operator-= (cctk_complex & a, cctk_real const & b) { return a = CCTK_Cmplx##Sub (a, CCTK_Cmplx (b, 0)); } \
+PREFIX cctk_complex operator*= (cctk_complex & a, cctk_complex const & b) { return a = CCTK_Cmplx##Mul (a, b); } \
+PREFIX cctk_complex operator*= (cctk_complex & a, cctk_real const & b) { return a = CCTK_Cmplx##Mul (a, CCTK_Cmplx (b, 0)); } \
+PREFIX cctk_complex operator/= (cctk_complex & a, cctk_complex const & b) { return a = CCTK_Cmplx##Div (a, b); } \
+PREFIX cctk_complex operator/= (cctk_complex & a, cctk_real const & b) { return a = CCTK_Cmplx##Div (a, CCTK_Cmplx (b, 0)); } \
+PREFIX cctk_complex pow (cctk_complex const & a, int i) { return CCTK_Cmplx##IPow (a, i); } \
+PREFIX cctk_complex pow (cctk_complex const & a, cctk_real const & r) { return CCTK_Cmplx##Pow (a, r); } \
+PREFIX cctk_complex pow (cctk_complex const & a, cctk_complex const & b) { return CCTK_Cmplx##CPow (a, b); } \
+PREFIX cctk_complex sin (cctk_complex const & a) { return CCTK_Cmplx##Sin (a); } \
+PREFIX cctk_complex cos (cctk_complex const & a) { return CCTK_Cmplx##Cos (a); } \
+PREFIX cctk_complex exp (cctk_complex const & a) { return CCTK_Cmplx##Exp (a); } \
+PREFIX cctk_complex log (cctk_complex const & a) { return CCTK_Cmplx##Log (a); } \
+PREFIX cctk_complex sqrt (cctk_complex const & a) { return CCTK_Cmplx##Sqrt (a); } \
+PREFIX std::ostream & operator << (std::ostream & os, cctk_complex const & a) { return os << CCTK_Cmplx##Real (a) << " " << CCTK_Cmplx##Imag (a); }
+#else
+ /* define no C++ overloaded operators */
+# define DEFINE_CMPLX_CXX_OPERATORS(CCTK_Cmplx, cctk_real, cctk_complex)
+#endif
+
+
/* macro to define a set of complex functions for a given precision */
-#define DEFINE_CMPLX_FUNCTIONS(CCTK_Cmplx, cctk_real, cctk_complex) \
- DEFINE_CCTK_CMPLX (CCTK_Cmplx, cctk_real, cctk_complex) \
- DEFINE_CCTK_CMPLX_REAL (CCTK_Cmplx, cctk_real, cctk_complex) \
- DEFINE_CCTK_CMPLX_IMAG (CCTK_Cmplx, cctk_real, cctk_complex) \
- DEFINE_CCTK_CMPLX_CONJG (CCTK_Cmplx, cctk_real, cctk_complex) \
- DEFINE_CCTK_CMPLX_ABS (CCTK_Cmplx, cctk_real, cctk_complex) \
- DEFINE_CCTK_CMPLX_ADD (CCTK_Cmplx, cctk_real, cctk_complex) \
- DEFINE_CCTK_CMPLX_SUB (CCTK_Cmplx, cctk_real, cctk_complex) \
- DEFINE_CCTK_CMPLX_MUL (CCTK_Cmplx, cctk_real, cctk_complex) \
- DEFINE_CCTK_CMPLX_DIV (CCTK_Cmplx, cctk_real, cctk_complex) \
- DEFINE_CCTK_CMPLX_SIN (CCTK_Cmplx, cctk_real, cctk_complex) \
- DEFINE_CCTK_CMPLX_COS (CCTK_Cmplx, cctk_real, cctk_complex) \
- DEFINE_CCTK_CMPLX_EXP (CCTK_Cmplx, cctk_real, cctk_complex) \
- DEFINE_CCTK_CMPLX_SQRT (CCTK_Cmplx, cctk_real, cctk_complex) \
- DEFINE_CCTK_CMPLX_POW (CCTK_Cmplx, cctk_real, cctk_complex) \
+#define DEFINE_CMPLX_FUNCTIONS(CCTK_Cmplx, cctk_real, cctk_complex, type) \
+PREFIX DEFINE_CCTK_CMPLX (CCTK_Cmplx, cctk_real, cctk_complex) \
+PREFIX DEFINE_CCTK_CMPLX_REAL (CCTK_Cmplx, cctk_real, cctk_complex) \
+PREFIX DEFINE_CCTK_CMPLX_IMAG (CCTK_Cmplx, cctk_real, cctk_complex) \
+PREFIX DEFINE_CCTK_CMPLX_NEG (CCTK_Cmplx, cctk_real, cctk_complex) \
+PREFIX DEFINE_CCTK_CMPLX_CONJG (CCTK_Cmplx, cctk_real, cctk_complex) \
+PREFIX DEFINE_CCTK_CMPLX_ABS (CCTK_Cmplx, cctk_real, cctk_complex, type) \
+PREFIX DEFINE_CCTK_CMPLX_NORM (CCTK_Cmplx, cctk_real, cctk_complex, type) \
+PREFIX DEFINE_CCTK_CMPLX_ARG (CCTK_Cmplx, cctk_real, cctk_complex, type) \
+PREFIX DEFINE_CCTK_CMPLX_ADD (CCTK_Cmplx, cctk_real, cctk_complex) \
+PREFIX DEFINE_CCTK_CMPLX_SUB (CCTK_Cmplx, cctk_real, cctk_complex) \
+PREFIX DEFINE_CCTK_CMPLX_MUL (CCTK_Cmplx, cctk_real, cctk_complex) \
+PREFIX DEFINE_CCTK_CMPLX_DIV (CCTK_Cmplx, cctk_real, cctk_complex, type) \
+PREFIX DEFINE_CCTK_CMPLX_CPOW (CCTK_Cmplx, cctk_real, cctk_complex) \
+PREFIX DEFINE_CCTK_CMPLX_SIN (CCTK_Cmplx, cctk_real, cctk_complex, type) \
+PREFIX DEFINE_CCTK_CMPLX_COS (CCTK_Cmplx, cctk_real, cctk_complex, type) \
+PREFIX DEFINE_CCTK_CMPLX_EXP (CCTK_Cmplx, cctk_real, cctk_complex, type) \
+PREFIX DEFINE_CCTK_CMPLX_LOG (CCTK_Cmplx, cctk_real, cctk_complex, type) \
+PREFIX DEFINE_CCTK_CMPLX_SQRT (CCTK_Cmplx, cctk_real, cctk_complex, type) \
+PREFIX DEFINE_CCTK_CMPLX_POW (CCTK_Cmplx, cctk_real, cctk_complex) \
+PREFIX DEFINE_CCTK_CMPLX_IPOW (CCTK_Cmplx, cctk_real, cctk_complex) \
+DEFINE_CMPLX_CXX_OPERATORS(CCTK_Cmplx, cctk_real, cctk_complex)
+
/* define complex functions for all available precisions */
/* NOTE: The code in this file will function correctly with C float
* and double types, but not with long double.
*/
-#ifdef CCTK_REAL4
- DEFINE_CMPLX_FUNCTIONS (CCTK_Cmplx8, CCTK_REAL4, CCTK_COMPLEX8)
+#ifdef HAVE_CCTK_REAL4
+ #define _REAL_TYPE f
+ DEFINE_CMPLX_FUNCTIONS (CCTK_Cmplx8, CCTK_REAL4, CCTK_COMPLEX8, _REAL_TYPE)
+ #undef _REAL_TYPE
#endif
-#ifdef CCTK_REAL8
- DEFINE_CMPLX_FUNCTIONS (CCTK_Cmplx16, CCTK_REAL8, CCTK_COMPLEX16)
+#ifdef HAVE_CCTK_REAL8
+ #define _REAL_TYPE
+ DEFINE_CMPLX_FUNCTIONS (CCTK_Cmplx16, CCTK_REAL8, CCTK_COMPLEX16, _REAL_TYPE)
+ #undef _REAL_TYPE
#endif
-#ifdef CCTK_REAL16
- DEFINE_CMPLX_FUNCTIONS (CCTK_Cmplx32, CCTK_REAL16, CCTK_COMPLEX32)
+#ifdef HAVE_CCTK_REAL16
+ #define _REAL_TYPE l
+ DEFINE_CMPLX_FUNCTIONS (CCTK_Cmplx32, CCTK_REAL16, CCTK_COMPLEX32, _REAL_TYPE)
+ #undef _REAL_TYPE
#endif
+
+#undef PREFIX
generated by cgit v1.2.3 (git 2.39.1) at 2024-07-25 23:13:17 +0000