22//
33// complex.h: Rcpp R/C++ interface class library -- complex
44//
5- // Copyright (C) 2010 - 2011 Dirk Eddelbuettel and Romain Francois
5+ // Copyright (C) 2010 - 2018 Dirk Eddelbuettel and Romain Francois
66//
77// This file is part of Rcpp.
88//
2222#ifndef Rcpp__sugar__complex_h
2323#define Rcpp__sugar__complex_h
2424
25- #ifdef HAVE_HYPOT
26- # define RCPP_HYPOT ::hypot
27- #else
28- # define RCPP_HYPOT ::Rf_pythag
29- #endif
30-
3125namespace Rcpp {
3226namespace sugar {
3327
@@ -60,89 +54,89 @@ class SugarComplex : public Rcpp::VectorBase<
6054
6155namespace internal {
6256inline double complex__Re ( Rcomplex x){ return x.r ; }
63- inline double complex__Im ( Rcomplex x){ return x.i ; }
64- inline double complex__Mod ( Rcomplex x){ return ::sqrt ( x.i * x.i + x.r * x.r ) ; }
65- inline Rcomplex complex__Conj ( Rcomplex x){
66- Rcomplex y ;
67- y.r = x.r ;
68- y.i = -x.i ;
69- return y ;
70- }
71- inline double complex__Arg ( Rcomplex x ){ return ::atan2 (x.i , x.r ); }
72- // TODO: this does not use HAVE_C99_COMPLEX as in R, perhaps it should
73- inline Rcomplex complex__exp ( Rcomplex x){
74- Rcomplex y ;
75- double expx = ::exp (x.r );
76- y.r = expx * ::cos (x.i );
77- y.i = expx * ::sin (x.i );
78- return y ;
79- }
80- inline Rcomplex complex__log ( Rcomplex x){
81- Rcomplex y ;
82- y.i = ::atan2 (x.i , x.r );
83- y.r = ::log ( RCPP_HYPOT ( x.r , x.i ) );
84- return y ;
85- }
86- inline Rcomplex complex__sqrt (Rcomplex z){
87- Rcomplex r ;
88- double mag;
57+ inline double complex__Im ( Rcomplex x){ return x.i ; }
58+ inline double complex__Mod ( Rcomplex x){ return ::sqrt ( x.i * x.i + x.r * x.r ) ; }
59+ inline Rcomplex complex__Conj ( Rcomplex x){
60+ Rcomplex y ;
61+ y.r = x.r ;
62+ y.i = -x.i ;
63+ return y ;
64+ }
65+ inline double complex__Arg ( Rcomplex x ){ return ::atan2 (x.i , x.r ); }
66+ // TODO: this does not use HAVE_C99_COMPLEX as in R, perhaps it should
67+ inline Rcomplex complex__exp ( Rcomplex x){
68+ Rcomplex y ;
69+ double expx = ::exp (x.r );
70+ y.r = expx * ::cos (x.i );
71+ y.i = expx * ::sin (x.i );
72+ return y ;
73+ }
74+ inline Rcomplex complex__log ( Rcomplex x){
75+ Rcomplex y ;
76+ y.i = ::atan2 (x.i , x.r );
77+ y.r = ::log (:: hypot ( x.r , x.i ) );
78+ return y ;
79+ }
80+ inline Rcomplex complex__sqrt (Rcomplex z){
81+ Rcomplex r ;
82+ double mag;
8983
90- if ( (mag = RCPP_HYPOT (z.r , z.i )) == 0.0 )
91- r.r = r.i = 0.0 ;
92- else if (z.r > 0 ) {
93- r.r = ::sqrt (0.5 * (mag + z.r ) );
94- r.i = z.i / r.r / 2 ;
95- }
96- else {
97- r.i = ::sqrt (0.5 * (mag - z.r ) );
98- if (z.i < 0 )
99- r.i = - r.i ;
100- r.r = z.i / r.i / 2 ;
101- }
102- return r ;
103- }
104- inline Rcomplex complex__cos (Rcomplex z){
105- Rcomplex r ;
106- r.r = ::cos (z.r ) * ::cosh (z.i );
107- r.i = - ::sin (z.r ) * ::sinh (z.i );
108- return r ;
109- }
110- inline Rcomplex complex__cosh (Rcomplex z){
111- Rcomplex r;
112- r.r = ::cos (-z.i ) * ::cosh ( z.r );
113- r.i = - ::sin (-z.i ) * ::sinh (z.r );
114- return r ;
115- }
116- inline Rcomplex complex__sin (Rcomplex z){
117- Rcomplex r ;
118- r.r = ::sin (z.r ) * ::cosh (z.i );
119- r.i = ::cos (z.r ) * ::sinh (z.i );
120- return r;
121- }
122- inline Rcomplex complex__tan (Rcomplex z){
123- Rcomplex r ;
124- double x2, y2, den;
125- x2 = 2.0 * z.r ;
126- y2 = 2.0 * z.i ;
127- den = ::cos (x2) + ::cosh (y2);
128- r.r = ::sin (x2)/den;
129- /* any threshold between -log(DBL_EPSILON)
130- and log(DBL_XMAX) will do*/
131- if (ISNAN (y2) || ::fabs (y2) < 50.0 )
132- r.i = ::sinh (y2)/den;
133- else
134- r.i = (y2 <0 ? -1.0 : 1.0 );
135- return r ;
136- }
84+ if ( (mag = :: hypot (z.r , z.i )) == 0.0 )
85+ r.r = r.i = 0.0 ;
86+ else if (z.r > 0 ) {
87+ r.r = ::sqrt (0.5 * (mag + z.r ) );
88+ r.i = z.i / r.r / 2 ;
89+ }
90+ else {
91+ r.i = ::sqrt (0.5 * (mag - z.r ) );
92+ if (z.i < 0 )
93+ r.i = - r.i ;
94+ r.r = z.i / r.i / 2 ;
95+ }
96+ return r ;
97+ }
98+ inline Rcomplex complex__cos (Rcomplex z){
99+ Rcomplex r ;
100+ r.r = ::cos (z.r ) * ::cosh (z.i );
101+ r.i = - ::sin (z.r ) * ::sinh (z.i );
102+ return r ;
103+ }
104+ inline Rcomplex complex__cosh (Rcomplex z){
105+ Rcomplex r;
106+ r.r = ::cos (-z.i ) * ::cosh ( z.r );
107+ r.i = - ::sin (-z.i ) * ::sinh (z.r );
108+ return r ;
109+ }
110+ inline Rcomplex complex__sin (Rcomplex z){
111+ Rcomplex r ;
112+ r.r = ::sin (z.r ) * ::cosh (z.i );
113+ r.i = ::cos (z.r ) * ::sinh (z.i );
114+ return r;
115+ }
116+ inline Rcomplex complex__tan (Rcomplex z){
117+ Rcomplex r ;
118+ double x2, y2, den;
119+ x2 = 2.0 * z.r ;
120+ y2 = 2.0 * z.i ;
121+ den = ::cos (x2) + ::cosh (y2);
122+ r.r = ::sin (x2)/den;
123+ /* any threshold between -log(DBL_EPSILON)
124+ and log(DBL_XMAX) will do*/
125+ if (ISNAN (y2) || ::fabs (y2) < 50.0 )
126+ r.i = ::sinh (y2)/den;
127+ else
128+ r.i = (y2 <0 ? -1.0 : 1.0 );
129+ return r ;
130+ }
137131
138132inline Rcomplex complex__asin (Rcomplex z)
139133{
140134 Rcomplex r ;
141135 double alpha, bet, t1, t2, x, y;
142136 x = z.r ;
143137 y = z.i ;
144- t1 = 0.5 * RCPP_HYPOT (x + 1 , y);
145- t2 = 0.5 * RCPP_HYPOT (x - 1 , y);
138+ t1 = 0.5 * :: hypot (x + 1 , y);
139+ t2 = 0.5 * :: hypot (x - 1 , y);
146140 alpha = t1 + t2;
147141 bet = t1 - t2;
148142 r.r = ::asin (bet);
@@ -159,13 +153,13 @@ inline Rcomplex complex__acos(Rcomplex z)
159153 return r ;
160154}
161155
162- /* Complex Arctangent Function */
163- /* Equation (4.4.39) Abramowitz and Stegun */
164- /* with additional terms to force the branch cuts */
165- /* to agree with figure 4.4, p79. Continuity */
166- /* on the branch cuts (pure imaginary axis; x==0, |y|>1) */
167- /* is standard: z_asin() is continuous from the right */
168- /* if y >= 1, and continuous from the left if y <= -1. */
156+ /* Complex Arctangent Function */
157+ /* Equation (4.4.39) Abramowitz and Stegun */
158+ /* with additional terms to force the branch cuts */
159+ /* to agree with figure 4.4, p79. Continuity */
160+ /* on the branch cuts (pure imaginary axis; x==0, |y|>1) */
161+ /* is standard: z_asin() is continuous from the right */
162+ /* if y >= 1, and continuous from the left if y <= -1. */
169163
170164inline Rcomplex complex__atan (Rcomplex z)
171165{
@@ -175,7 +169,7 @@ inline Rcomplex complex__atan(Rcomplex z)
175169 y = z.i ;
176170 r.r = 0.5 * ::atan (2 * x / ( 1 - x * x - y * y));
177171 r.i = 0.25 * ::log ((x * x + (y + 1 ) * (y + 1 )) /
178- (x * x + (y - 1 ) * (y - 1 )));
172+ (x * x + (y - 1 ) * (y - 1 )));
179173 if (x*x + y*y > 1 ) {
180174 r.r += M_PI_2;
181175 if (x < 0 || (x == 0 && y < 0 )) r.r -= M_PI;
@@ -184,32 +178,32 @@ inline Rcomplex complex__atan(Rcomplex z)
184178}
185179
186180
187- inline Rcomplex complex__acosh (Rcomplex z){
188- Rcomplex r, a = complex__acos (z);
189- r.r = -a.i ;
190- r.i = a.r ;
191- return r ;
192- }
181+ inline Rcomplex complex__acosh (Rcomplex z){
182+ Rcomplex r, a = complex__acos (z);
183+ r.r = -a.i ;
184+ r.i = a.r ;
185+ return r ;
186+ }
193187
194- inline Rcomplex complex__asinh (Rcomplex z){
195- Rcomplex r, b;
196- b.r = -z.i ;
197- b.i = z.r ;
198- Rcomplex a = complex__asin (b);
199- r.r = a.i ;
200- r.i = -a.r ;
201- return r ;
202- }
188+ inline Rcomplex complex__asinh (Rcomplex z){
189+ Rcomplex r, b;
190+ b.r = -z.i ;
191+ b.i = z.r ;
192+ Rcomplex a = complex__asin (b);
193+ r.r = a.i ;
194+ r.i = -a.r ;
195+ return r ;
196+ }
203197
204- inline Rcomplex complex__atanh (Rcomplex z){
205- Rcomplex r, b;
206- b.r = -z.i ;
207- b.i = z.r ;
208- Rcomplex a = complex__atan (b);
209- r.r = a.i ;
210- r.i = -a.r ;
211- return r ;
212- }
198+ inline Rcomplex complex__atanh (Rcomplex z){
199+ Rcomplex r, b;
200+ b.r = -z.i ;
201+ b.i = z.r ;
202+ Rcomplex a = complex__atan (b);
203+ r.r = a.i ;
204+ r.i = -a.r ;
205+ return r ;
206+ }
213207inline Rcomplex complex__sinh (Rcomplex z)
214208{
215209 Rcomplex r, b;
@@ -232,20 +226,15 @@ inline Rcomplex complex__tanh(Rcomplex z)
232226 return r ;
233227}
234228
235-
236-
237229} // internal
238230
239- #define RCPP_SUGAR_COMPLEX (__NAME__,__OUT__ ) \
240- template <bool NA, typename T> \
241- inline sugar::SugarComplex<NA,__OUT__,T, __OUT__ (*)(Rcomplex) > \
242- __NAME__ ( \
243- const VectorBase<CPLXSXP,NA,T>& t \
244- ){ \
245- return sugar::SugarComplex<NA,__OUT__,T, __OUT__ (*)(Rcomplex) >( \
246- internal::complex__##__NAME__, t \
247- ) ; \
248- }
231+ #define RCPP_SUGAR_COMPLEX (__NAME__,__OUT__ ) \
232+ template <bool NA, typename T> \
233+ inline sugar::SugarComplex<NA,__OUT__,T, __OUT__ (*)(Rcomplex) > \
234+ __NAME__ (const VectorBase<CPLXSXP,NA,T>& t) { \
235+ return sugar::SugarComplex<NA,__OUT__,T, __OUT__ (*)(Rcomplex) >( \
236+ internal::complex__##__NAME__, t); \
237+ }
249238
250239RCPP_SUGAR_COMPLEX ( Re, double )
251240RCPP_SUGAR_COMPLEX ( Im, double )
0 commit comments