2
2
//
3
3
// complex.h: Rcpp R/C++ interface class library -- complex
4
4
//
5
- // Copyright (C) 2010 - 2011 Dirk Eddelbuettel and Romain Francois
5
+ // Copyright (C) 2010 - 2018 Dirk Eddelbuettel and Romain Francois
6
6
//
7
7
// This file is part of Rcpp.
8
8
//
22
22
#ifndef Rcpp__sugar__complex_h
23
23
#define Rcpp__sugar__complex_h
24
24
25
- #ifdef HAVE_HYPOT
26
- # define RCPP_HYPOT ::hypot
27
- #else
28
- # define RCPP_HYPOT ::Rf_pythag
29
- #endif
30
-
31
25
namespace Rcpp {
32
26
namespace sugar {
33
27
@@ -60,89 +54,89 @@ class SugarComplex : public Rcpp::VectorBase<
60
54
61
55
namespace internal {
62
56
inline 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;
89
83
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
+ }
137
131
138
132
inline Rcomplex complex__asin (Rcomplex z)
139
133
{
140
134
Rcomplex r ;
141
135
double alpha, bet, t1, t2, x, y;
142
136
x = z.r ;
143
137
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);
146
140
alpha = t1 + t2;
147
141
bet = t1 - t2;
148
142
r.r = ::asin (bet);
@@ -159,13 +153,13 @@ inline Rcomplex complex__acos(Rcomplex z)
159
153
return r ;
160
154
}
161
155
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. */
169
163
170
164
inline Rcomplex complex__atan (Rcomplex z)
171
165
{
@@ -175,7 +169,7 @@ inline Rcomplex complex__atan(Rcomplex z)
175
169
y = z.i ;
176
170
r.r = 0.5 * ::atan (2 * x / ( 1 - x * x - y * y));
177
171
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 )));
179
173
if (x*x + y*y > 1 ) {
180
174
r.r += M_PI_2;
181
175
if (x < 0 || (x == 0 && y < 0 )) r.r -= M_PI;
@@ -184,32 +178,32 @@ inline Rcomplex complex__atan(Rcomplex z)
184
178
}
185
179
186
180
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
+ }
193
187
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
+ }
203
197
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
+ }
213
207
inline Rcomplex complex__sinh (Rcomplex z)
214
208
{
215
209
Rcomplex r, b;
@@ -232,20 +226,15 @@ inline Rcomplex complex__tanh(Rcomplex z)
232
226
return r ;
233
227
}
234
228
235
-
236
-
237
229
} // internal
238
230
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
+ }
249
238
250
239
RCPP_SUGAR_COMPLEX ( Re, double )
251
240
RCPP_SUGAR_COMPLEX ( Im, double )
0 commit comments