From 65b32c63b24b094d6a0457e801bd871958dd9ab6 Mon Sep 17 00:00:00 2001 From: eschnett Date: Fri, 8 Mar 2013 20:29:17 +0000 Subject: Replace Cactus complex number type with C/C++ complex numbers Map CCTK_COMPLEX to "double complex" in C, and "complex" in C++. (It is already mapped to "double complex" in Fortran.) Update type definitions. Re-implement Cactus complex number math functions by calling the respective C functions. Update thorn that access real and imaginary parts of complex numbers to use standard-conforming methods instead. git-svn-id: http://svn.cactuscode.org/flesh/trunk@4979 17b73243-c579-4c4c-a9d2-2d5706c11dac --- src/include/cctk_Complex.h | 110 +++-------- src/include/cctk_Types.h | 52 +++--- src/main/Complex.c | 452 ++++++++++++++------------------------------- src/util/Table.c | 9 +- 4 files changed, 185 insertions(+), 438 deletions(-) diff --git a/src/include/cctk_Complex.h b/src/include/cctk_Complex.h index 977896d3..96eaf17f 100644 --- a/src/include/cctk_Complex.h +++ b/src/include/cctk_Complex.h @@ -12,92 +12,32 @@ #define _CCTK_COMPLEX_H_ #ifdef __cplusplus -# 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 -# include -# 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) +extern "C" { #endif /* macro to declare a set of complex functions for a given precision */ #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) +cctk_complex CCTK_Cmplx (cctk_real Re, cctk_real Im); \ +cctk_real CCTK_Cmplx##Real (cctk_complex a); \ +cctk_real CCTK_Cmplx##Imag (cctk_complex a); \ +cctk_complex CCTK_Cmplx##Neg (cctk_complex a); \ +cctk_complex CCTK_Cmplx##Conjg (cctk_complex a); \ +cctk_real CCTK_Cmplx##Abs (cctk_complex a); \ +cctk_real CCTK_Cmplx##Arg (cctk_complex a); \ +cctk_real CCTK_Cmplx##Norm (cctk_complex a); \ +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##CPow (cctk_complex a, cctk_complex b); \ +cctk_complex CCTK_Cmplx##Sin (cctk_complex a); \ +cctk_complex CCTK_Cmplx##Cos (cctk_complex a); \ +cctk_complex CCTK_Cmplx##Exp (cctk_complex a); \ +cctk_complex CCTK_Cmplx##Log (cctk_complex a); \ +cctk_complex CCTK_Cmplx##Sqrt (cctk_complex a); \ +cctk_complex CCTK_Cmplx##Pow (cctk_complex a, cctk_real b); \ +cctk_complex CCTK_Cmplx##IPow (cctk_complex a, int b); /* declare complex functions for all available precisions */ @@ -113,8 +53,6 @@ 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 @@ -182,10 +120,8 @@ DECLARE_CMPLX_FUNCTIONS (CCTK_Cmplx32, CCTK_REAL16, CCTK_COMPLEX32) #define CCTK_CmplxIPow CCTK_Cmplx32IPow #endif -#ifndef DEFINE_CCTK_COMPLEX_EXTERN_FUNCTIONS -# define DEFINE_CCTK_COMPLEX_INLINE_FUNCTIONS -# include "../main/Complex.c" -# undef DEFINE_CCTK_COMPLEX_INLINE_FUNCTIONS +#ifdef __cplusplus +} #endif #endif /* _CCTK_COMPLEX_H_ */ diff --git a/src/include/cctk_Types.h b/src/include/cctk_Types.h index e3fb55a6..eeda7d77 100644 --- a/src/include/cctk_Types.h +++ b/src/include/cctk_Types.h @@ -34,45 +34,37 @@ typedef const char * CCTK_STRING; #define HAVE_CCTK_CHAR 1 #define HAVE_CCTK_STRING 1 -/* Structures for complex types */ +/* Declarations for complex types */ -#ifdef HAVE_CCTK_REAL16 -#define HAVE_CCTK_COMPLEX32 1 -typedef struct CCTK_COMPLEX32 -{ - CCTK_REAL16 Re; - CCTK_REAL16 Im; #ifdef __cplusplus - CCTK_REAL16 real() const { return Re; } - CCTK_REAL16 imag() const { return Im; } +# include #endif -} CCTK_COMPLEX32; + +#ifdef HAVE_CCTK_REAL16 +# define HAVE_CCTK_COMPLEX32 1 +# ifdef __cplusplus +typedef std::complex CCTK_COMPLEX32; +# else +typedef long double _Complex CCTK_COMPLEX32; +# endif #endif #ifdef HAVE_CCTK_REAL8 #define HAVE_CCTK_COMPLEX16 1 -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; +# ifdef __cplusplus +typedef std::complex CCTK_COMPLEX16; +# else +typedef double _Complex CCTK_COMPLEX16; +# endif #endif #ifdef HAVE_CCTK_REAL4 -#define HAVE_CCTK_COMPLEX8 1 -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; +# define HAVE_CCTK_COMPLEX8 1 +# ifdef __cplusplus +typedef std::complex CCTK_COMPLEX8; +# else +typedef float _Complex CCTK_COMPLEX8; +# endif #endif /* Small positive integer type */ @@ -262,7 +254,7 @@ typedef unsigned char CCTK_BYTE; /* We have __attribute__((unused)), so use it */ #define CCTK_DECLARE_INIT(typ,nam,val) \ - typ nam __attribute__((unused)) = (val); + typ nam CCTK_ATTRIBUTE_UNUSED = (val); #else diff --git a/src/main/Complex.c b/src/main/Complex.c index 47e3a888..84bcf7c6 100644 --- a/src/main/Complex.c +++ b/src/main/Complex.c @@ -8,6 +8,7 @@ @version $Id$ @@*/ +#include #include #include "cctk_Flesh.h" @@ -72,16 +73,10 @@ CCTK_FILEVERSION(main_Complex_c); The complex number @endreturndesc @@*/ -#define DEFINE_CCTK_CMPLX(CCTK_Cmplx, cctk_real, cctk_complex) \ -cctk_complex CCTK_Cmplx (cctk_real Re, cctk_real Im) \ -{ \ - cctk_complex result; \ - \ - \ - result.Re = Re; \ - result.Im = Im; \ - \ - return (result); \ +#define DEFINE_CCTK_CMPLX(CCTK_Cmplx, cctk_real, cctk_complex, type) \ +cctk_complex CCTK_Cmplx (cctk_real Re, cctk_real Im) \ +{ \ + return Re + _Complex_I * Im; \ } @@ -93,7 +88,7 @@ cctk_complex CCTK_Cmplx (cctk_real Re, cctk_real Im) \ Returns the real part of a complex number. @enddesc - @var complex_number + @var x @vdesc The complex number @vtype CCTK_COMPLEX @vio in @@ -104,10 +99,10 @@ cctk_complex CCTK_Cmplx (cctk_real Re, cctk_real Im) \ The real part @endreturndesc @@*/ -#define DEFINE_CCTK_CMPLX_REAL(CCTK_Cmplx, cctk_real, cctk_complex) \ -cctk_real CCTK_Cmplx##Real (cctk_complex complex_number) \ -{ \ - return (complex_number.Re); \ +#define DEFINE_CCTK_CMPLX_REAL(CCTK_Cmplx, cctk_real, cctk_complex, type) \ +cctk_real CCTK_Cmplx##Real (cctk_complex x) \ +{ \ + return creal##type(x); \ } @@ -119,7 +114,7 @@ cctk_real CCTK_Cmplx##Real (cctk_complex complex_number) \ Returns the imaginary part of a complex number. @enddesc - @var complex_number + @var x @vdesc The complex number @vtype CCTK_COMPLEX @vio in @@ -130,10 +125,10 @@ cctk_real CCTK_Cmplx##Real (cctk_complex complex_number) \ The imaginary part @endreturndesc @@*/ -#define DEFINE_CCTK_CMPLX_IMAG(CCTK_Cmplx, cctk_real, cctk_complex) \ -cctk_real CCTK_Cmplx##Imag (cctk_complex complex_number) \ -{ \ - return (complex_number.Im); \ +#define DEFINE_CCTK_CMPLX_IMAG(CCTK_Cmplx, cctk_real, cctk_complex, type) \ +cctk_real CCTK_Cmplx##Imag (cctk_complex x) \ +{ \ + return cimag##type(x); \ } @@ -145,7 +140,7 @@ cctk_real CCTK_Cmplx##Imag (cctk_complex complex_number) \ Returns the negative of a complex number. @enddesc - @var in + @var x @vdesc The complex number @vtype CCTK_COMPLEX @vio in @@ -156,15 +151,10 @@ cctk_real CCTK_Cmplx##Imag (cctk_complex complex_number) \ 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); \ +#define DEFINE_CCTK_CMPLX_NEG(CCTK_Cmplx, cctk_real, cctk_complex, type) \ +cctk_complex CCTK_Cmplx##Neg (cctk_complex x) \ +{ \ + return -x; \ } @@ -176,7 +166,7 @@ cctk_complex CCTK_Cmplx##Neg (cctk_complex complex_number) \ Returns the complex conjugate of a complex number. @enddesc - @var in + @var x @vdesc The complex number @vtype CCTK_COMPLEX @vio in @@ -187,15 +177,10 @@ cctk_complex CCTK_Cmplx##Neg (cctk_complex complex_number) \ The complex conjugate @endreturndesc @@*/ -#define DEFINE_CCTK_CMPLX_CONJG(CCTK_Cmplx, cctk_real, cctk_complex) \ -cctk_complex CCTK_Cmplx##Conjg (cctk_complex complex_number) \ -{ \ - cctk_complex result; \ - \ - \ - result.Re = complex_number.Re; \ - result.Im = -complex_number.Im; \ - return (result); \ +#define DEFINE_CCTK_CMPLX_CONJG(CCTK_Cmplx, cctk_real, cctk_complex, type) \ +cctk_complex CCTK_Cmplx##Conjg (cctk_complex x) \ +{ \ + return conj##type(x); \ } @@ -207,7 +192,7 @@ cctk_complex CCTK_Cmplx##Conjg (cctk_complex complex_number) \ Return the absolute value of a complex number. @enddesc - @var in + @var x @vdesc The complex number @vtype CCTK_COMPLEX @vio in @@ -218,10 +203,10 @@ 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, type) \ -cctk_real CCTK_Cmplx##Abs (cctk_complex complex_number) \ -{ \ - return ((cctk_real)hypot##type (complex_number.Re, complex_number.Im)); \ +#define DEFINE_CCTK_CMPLX_ABS(CCTK_Cmplx, cctk_real, cctk_complex, type) \ +cctk_real CCTK_Cmplx##Abs (cctk_complex x) \ +{ \ + return cabs##type(x); \ } @@ -233,7 +218,7 @@ cctk_real CCTK_Cmplx##Abs (cctk_complex complex_number) \ Return the absolute value squared of a complex number. @enddesc - @var in + @var x @vdesc The complex number @vtype CCTK_COMPLEX @vio in @@ -244,10 +229,10 @@ cctk_real CCTK_Cmplx##Abs (cctk_complex complex_number) \ 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));\ +#define DEFINE_CCTK_CMPLX_NORM(CCTK_Cmplx, cctk_real, cctk_complex, type) \ +cctk_real CCTK_Cmplx##Norm (cctk_complex x) \ +{ \ + return creal##type(x)*creal##type(x) + cimag##type(x)*cimag##type(x); \ } @@ -259,7 +244,7 @@ cctk_real CCTK_Cmplx##Norm (cctk_complex complex_number) \ Return the argument of a complex number. @enddesc - @var in + @var x @vdesc The complex number @vtype CCTK_COMPLEX @vio in @@ -270,10 +255,10 @@ cctk_real CCTK_Cmplx##Norm (cctk_complex complex_number) \ 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)); \ +#define DEFINE_CCTK_CMPLX_ARG(CCTK_Cmplx, cctk_real, cctk_complex, type) \ +cctk_real CCTK_Cmplx##Arg (cctk_complex x) \ +{ \ + return carg##type(x); \ } @@ -301,15 +286,10 @@ cctk_real CCTK_Cmplx##Arg (cctk_complex complex_number) \ The sum of a and b @endreturndesc @@*/ -#define DEFINE_CCTK_CMPLX_ADD(CCTK_Cmplx, cctk_real, cctk_complex) \ -cctk_complex CCTK_Cmplx##Add (cctk_complex a, cctk_complex b) \ -{ \ - cctk_complex result; \ - \ - \ - result.Re = a.Re + b.Re; \ - result.Im = a.Im + b.Im; \ - return (result); \ +#define DEFINE_CCTK_CMPLX_ADD(CCTK_Cmplx, cctk_real, cctk_complex, type) \ +cctk_complex CCTK_Cmplx##Add (cctk_complex a, cctk_complex b) \ +{ \ + return a+b; \ } @@ -337,15 +317,10 @@ cctk_complex CCTK_Cmplx##Add (cctk_complex a, cctk_complex b) \ The difference @endreturndesc @@*/ -#define DEFINE_CCTK_CMPLX_SUB(CCTK_Cmplx, cctk_real, cctk_complex) \ -cctk_complex CCTK_Cmplx##Sub (cctk_complex a, cctk_complex b) \ -{ \ - cctk_complex result; \ - \ - \ - result.Re = a.Re - b.Re; \ - result.Im = a.Im - b.Im; \ - return (result); \ +#define DEFINE_CCTK_CMPLX_SUB(CCTK_Cmplx, cctk_real, cctk_complex, type) \ +cctk_complex CCTK_Cmplx##Sub (cctk_complex a, cctk_complex b) \ +{ \ + return a-b; \ } @@ -373,15 +348,10 @@ cctk_complex CCTK_Cmplx##Sub (cctk_complex a, cctk_complex b) \ The product @endreturndesc @@*/ -#define DEFINE_CCTK_CMPLX_MUL(CCTK_Cmplx, cctk_real, cctk_complex) \ -cctk_complex CCTK_Cmplx##Mul (cctk_complex a, cctk_complex b) \ -{ \ - cctk_complex result; \ - \ - \ - result.Re = a.Re*b.Re - a.Im*b.Im; \ - result.Im = a.Im*b.Re + a.Re*b.Im; \ - return (result); \ +#define DEFINE_CCTK_CMPLX_MUL(CCTK_Cmplx, cctk_real, cctk_complex, type) \ +cctk_complex CCTK_Cmplx##Mul (cctk_complex a, cctk_complex b) \ +{ \ + return a*b; \ } @@ -391,10 +361,6 @@ cctk_complex CCTK_Cmplx##Mul (cctk_complex a, cctk_complex b) \ @author Gabrielle Allen @desc Divide two complex numbers. - Make sure that the intermediate values do not overflow. - See e.g. the C99 compliant implementation in libgcc - (which ships with gcc). - This implementation here does not handle nan and inf. @enddesc @var a @@ -413,24 +379,10 @@ 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, type) \ -cctk_complex CCTK_Cmplx##Div (cctk_complex a, cctk_complex b) \ -{ \ - cctk_real afact, bfact, factor; \ - cctk_complex result; \ - \ - \ - afact = (cctk_real)(fabs##type(a.Re) + fabs##type(a.Im)); \ - a.Re /= afact; \ - a.Im /= afact; \ - 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; \ - result.Re = afact / bfact * (a.Re*b.Re + a.Im*b.Im) / factor; \ - result.Im = afact / bfact * (a.Im*b.Re - a.Re*b.Im) / factor; \ - \ - return (result); \ +#define DEFINE_CCTK_CMPLX_DIV(CCTK_Cmplx, cctk_real, cctk_complex, type) \ +cctk_complex CCTK_Cmplx##Div (cctk_complex a, cctk_complex b) \ +{ \ + return a/b; \ } @@ -440,8 +392,6 @@ cctk_complex CCTK_Cmplx##Div (cctk_complex a, cctk_complex b) \ @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 @@ -460,10 +410,10 @@ cctk_complex CCTK_Cmplx##Div (cctk_complex a, cctk_complex b) \ 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))); \ +#define DEFINE_CCTK_CMPLX_CPOW(CCTK_Cmplx, cctk_real, cctk_complex, type) \ +cctk_complex CCTK_Cmplx##CPow (cctk_complex a, cctk_complex b) \ +{ \ + return cpow##type(a, b); \ } @@ -475,7 +425,7 @@ cctk_complex CCTK_Cmplx##CPow (cctk_complex a, cctk_complex b) \ Returns the sine of a complex number. @enddesc - @var complex_number + @var x @vdesc The complex number @vtype CCTK_COMPLEX @vio in @@ -486,16 +436,10 @@ cctk_complex CCTK_Cmplx##CPow (cctk_complex a, cctk_complex b) \ The sine @endreturndesc @@*/ -#define DEFINE_CCTK_CMPLX_SIN(CCTK_Cmplx, cctk_real, cctk_complex, type) \ -cctk_complex CCTK_Cmplx##Sin (cctk_complex complex_number) \ -{ \ - cctk_complex result; \ - \ - \ - 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); \ +#define DEFINE_CCTK_CMPLX_SIN(CCTK_Cmplx, cctk_real, cctk_complex, type) \ +cctk_complex CCTK_Cmplx##Sin (cctk_complex x) \ +{ \ + return csin##type(x); \ } @@ -507,7 +451,7 @@ cctk_complex CCTK_Cmplx##Sin (cctk_complex complex_number) \ Returns the cosine of a complex number. @enddesc - @var complex_number + @var x @vdesc The complex number @vtype CCTK_COMPLEX @vio in @@ -518,16 +462,10 @@ cctk_complex CCTK_Cmplx##Sin (cctk_complex complex_number) \ The cosine @endreturndesc @@*/ -#define DEFINE_CCTK_CMPLX_COS(CCTK_Cmplx, cctk_real, cctk_complex, type) \ -cctk_complex CCTK_Cmplx##Cos (cctk_complex complex_number) \ -{ \ - cctk_complex result; \ - \ - \ - 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); \ +#define DEFINE_CCTK_CMPLX_COS(CCTK_Cmplx, cctk_real, cctk_complex, type) \ +cctk_complex CCTK_Cmplx##Cos (cctk_complex x) \ +{ \ + return ccos##type(x); \ } @@ -539,7 +477,7 @@ cctk_complex CCTK_Cmplx##Cos (cctk_complex complex_number) \ Returns the exponential of a complex number. @enddesc - @var complex_number + @var x @vdesc The complex number @vtype CCTK_COMPLEX @vio in @@ -550,19 +488,10 @@ cctk_complex CCTK_Cmplx##Cos (cctk_complex complex_number) \ The exponential @endreturndesc @@*/ -#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##type (complex_number.Re); \ - theta = complex_number.Im; \ - result.Re = rho * (cctk_real)cos##type (theta); \ - result.Im = rho * (cctk_real)sin##type (theta); \ - \ - return (result); \ +#define DEFINE_CCTK_CMPLX_EXP(CCTK_Cmplx, cctk_real, cctk_complex, type) \ +cctk_complex CCTK_Cmplx##Exp (cctk_complex x) \ +{ \ + return cexp##type(x); \ } @@ -572,11 +501,9 @@ cctk_complex CCTK_Cmplx##Exp (cctk_complex complex_number) \ @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 + @var x @vdesc The complex number @vtype CCTK_COMPLEX @vio in @@ -587,16 +514,10 @@ cctk_complex CCTK_Cmplx##Exp (cctk_complex complex_number) \ 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); \ +#define DEFINE_CCTK_CMPLX_LOG(CCTK_Cmplx, cctk_real, cctk_complex, type) \ +cctk_complex CCTK_Cmplx##Log (cctk_complex x) \ +{ \ + return clog##type(x); \ } @@ -606,15 +527,9 @@ cctk_complex CCTK_Cmplx##Log (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 - code contains references. - @endhistory - @var complex_number + @var x @vdesc The complex number @vtype CCTK_COMPLEX @vio in @@ -625,31 +540,10 @@ cctk_complex CCTK_Cmplx##Log (cctk_complex complex_number) \ The square root @endreturndesc @@*/ -#define DEFINE_CCTK_CMPLX_SQRT(CCTK_Cmplx, cctk_real, cctk_complex, type) \ -cctk_complex CCTK_Cmplx##Sqrt (cctk_complex complex_number) \ -{ \ - cctk_real d, r, s; \ - cctk_complex result; \ - \ - \ - 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) \ - { \ - r = sqrt##type (0.5##type * d + 0.5##type * complex_number.Re); \ - s = (0.5##type * complex_number.Im) / r; \ - } \ - else \ - { \ - 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); \ +#define DEFINE_CCTK_CMPLX_SQRT(CCTK_Cmplx, cctk_real, cctk_complex, type) \ +cctk_complex CCTK_Cmplx##Sqrt (cctk_complex x) \ +{ \ + return csqrt##type(x); \ } @@ -659,11 +553,9 @@ cctk_complex CCTK_Cmplx##Sqrt (cctk_complex complex_number) \ @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 + @var x @vdesc The complex number @vtype CCTK_COMPLEX @vio in @@ -679,18 +571,10 @@ cctk_complex CCTK_Cmplx##Sqrt (cctk_complex complex_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 complex_number, cctk_real w) \ -{ \ - cctk_complex result; \ - \ - \ - result = CCTK_Cmplx##Log (complex_number); \ - result.Re *= w; \ - result.Im *= w; \ - result = CCTK_Cmplx##Exp (result); \ - \ - return (result); \ +#define DEFINE_CCTK_CMPLX_POW(CCTK_Cmplx, cctk_real, cctk_complex, type) \ +cctk_complex CCTK_Cmplx##Pow (cctk_complex x, cctk_real w) \ +{ \ + return cpow##type(x, w); \ } @@ -700,11 +584,9 @@ cctk_complex CCTK_Cmplx##Pow (cctk_complex complex_number, cctk_real w) \ @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 + @var x @vdesc The complex number @vtype CCTK_COMPLEX @vio in @@ -720,132 +602,68 @@ cctk_complex CCTK_Cmplx##Pow (cctk_complex complex_number, cctk_real w) \ 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) \ +#define DEFINE_CCTK_CMPLX_IPOW(CCTK_Cmplx, cctk_real, cctk_complex, type) \ +cctk_complex CCTK_Cmplx##IPow (cctk_complex x, int w) \ { \ - cctk_complex result; \ - \ - \ - result = CCTK_Cmplx (1, 0); \ - if (w < 0) \ + int w0 = w; \ + cctk_complex result = 1; \ + while (w) \ { \ - result = CCTK_Cmplx##Div (result, \ - CCTK_Cmplx##IPow (complex_number, - w)); \ - } \ - else \ - { \ - while (w) \ + if (w % 2) \ { \ - if (w % 2) \ - { \ - result = CCTK_Cmplx##Mul (result, complex_number); \ - } \ - w /= 2; \ - complex_number = CCTK_Cmplx##Mul (complex_number, complex_number); \ + result *= x; \ } \ + w /= 2; \ + x *= x; \ + } \ + if (w0 < 0) \ + { \ + result = 1 / result; \ } \ - \ - return (result); \ + 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, 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 DEFINE_CMPLX_FUNCTIONS(CCTK_Cmplx, cctk_real, cctk_complex, type) \ +DEFINE_CCTK_CMPLX (CCTK_Cmplx, cctk_real, cctk_complex, type) \ +DEFINE_CCTK_CMPLX_REAL (CCTK_Cmplx, cctk_real, cctk_complex, type) \ +DEFINE_CCTK_CMPLX_IMAG (CCTK_Cmplx, cctk_real, cctk_complex, type) \ +DEFINE_CCTK_CMPLX_NEG (CCTK_Cmplx, cctk_real, cctk_complex, type) \ +DEFINE_CCTK_CMPLX_CONJG (CCTK_Cmplx, cctk_real, cctk_complex, type) \ +DEFINE_CCTK_CMPLX_ABS (CCTK_Cmplx, cctk_real, cctk_complex, type) \ +DEFINE_CCTK_CMPLX_NORM (CCTK_Cmplx, cctk_real, cctk_complex, type) \ +DEFINE_CCTK_CMPLX_ARG (CCTK_Cmplx, cctk_real, cctk_complex, type) \ +DEFINE_CCTK_CMPLX_ADD (CCTK_Cmplx, cctk_real, cctk_complex, type) \ +DEFINE_CCTK_CMPLX_SUB (CCTK_Cmplx, cctk_real, cctk_complex, type) \ +DEFINE_CCTK_CMPLX_MUL (CCTK_Cmplx, cctk_real, cctk_complex, type) \ +DEFINE_CCTK_CMPLX_DIV (CCTK_Cmplx, cctk_real, cctk_complex, type) \ +DEFINE_CCTK_CMPLX_CPOW (CCTK_Cmplx, cctk_real, cctk_complex, type) \ +DEFINE_CCTK_CMPLX_SIN (CCTK_Cmplx, cctk_real, cctk_complex, type) \ +DEFINE_CCTK_CMPLX_COS (CCTK_Cmplx, cctk_real, cctk_complex, type) \ +DEFINE_CCTK_CMPLX_EXP (CCTK_Cmplx, cctk_real, cctk_complex, type) \ +DEFINE_CCTK_CMPLX_LOG (CCTK_Cmplx, cctk_real, cctk_complex, type) \ +DEFINE_CCTK_CMPLX_SQRT (CCTK_Cmplx, cctk_real, cctk_complex, type) \ +DEFINE_CCTK_CMPLX_POW (CCTK_Cmplx, cctk_real, cctk_complex, type) \ +DEFINE_CCTK_CMPLX_IPOW (CCTK_Cmplx, cctk_real, cctk_complex, type) /* 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 HAVE_CCTK_REAL4 - #define _REAL_TYPE f - DEFINE_CMPLX_FUNCTIONS (CCTK_Cmplx8, CCTK_REAL4, CCTK_COMPLEX8, _REAL_TYPE) - #undef _REAL_TYPE + #define KIND_SUFFIX f + DEFINE_CMPLX_FUNCTIONS (CCTK_Cmplx8, CCTK_REAL4, CCTK_COMPLEX8, KIND_SUFFIX) + #undef KIND_SUFFIX #endif #ifdef HAVE_CCTK_REAL8 - #define _REAL_TYPE - DEFINE_CMPLX_FUNCTIONS (CCTK_Cmplx16, CCTK_REAL8, CCTK_COMPLEX16, _REAL_TYPE) - #undef _REAL_TYPE + #define KIND_SUFFIX + DEFINE_CMPLX_FUNCTIONS (CCTK_Cmplx16, CCTK_REAL8, CCTK_COMPLEX16, KIND_SUFFIX) + #undef KIND_SUFFIX #endif #ifdef HAVE_CCTK_REAL16 - #define _REAL_TYPE l - DEFINE_CMPLX_FUNCTIONS (CCTK_Cmplx32, CCTK_REAL16, CCTK_COMPLEX32, _REAL_TYPE) - #undef _REAL_TYPE + #define KIND_SUFFIX l + DEFINE_CMPLX_FUNCTIONS (CCTK_Cmplx32, CCTK_REAL16, CCTK_COMPLEX32, KIND_SUFFIX) + #undef KIND_SUFFIX #endif - -#undef PREFIX diff --git a/src/util/Table.c b/src/util/Table.c index 1c9a3db8..9a613a41 100644 --- a/src/util/Table.c +++ b/src/util/Table.c @@ -94,6 +94,7 @@ typedef int bool; #endif #include "cctk_Types.h" +#include "cctk_Complex.h" #include "cctk_Constants.h" #include "cctk_Groups.h" #include "cctk_Flesh.h" @@ -5591,8 +5592,8 @@ int Util_TablePrint(FILE *stream, int handle) for (i = 0 ; i < tep->N_elements ; ++i) { fprintf(stream, "\t(%g,%g)", - (double) value_ptr_complex[i].Re, - (double) value_ptr_complex[i].Im); + (double) CCTK_CmplxReal(value_ptr_complex[i]), + (double) CCTK_CmplxImag(value_ptr_complex[i])); } } break; @@ -5729,8 +5730,8 @@ int Util_TablePrintPretty(FILE *stream, int handle) fprintf(stream, " "); } fprintf(stream, "(%#.17g,%#.17g)", - (double)value_ptr_complex[i].Re, - (double)value_ptr_complex[i].Im); + (double)CCTK_CmplxReal(value_ptr_complex[i]), + (double)CCTK_CmplxImag(value_ptr_complex[i])); } if (tep->N_elements != 1) { -- cgit v1.2.3