Rcpp Version 0.12.12
complex.h
Go to the documentation of this file.
1 // -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; tab-width: 8 -*-
2 //
3 // complex.h: Rcpp R/C++ interface class library -- complex
4 //
5 // Copyright (C) 2010 - 2011 Dirk Eddelbuettel and Romain Francois
6 //
7 // This file is part of Rcpp.
8 //
9 // Rcpp is free software: you can redistribute it and/or modify it
10 // under the terms of the GNU General Public License as published by
11 // the Free Software Foundation, either version 2 of the License, or
12 // (at your option) any later version.
13 //
14 // Rcpp is distributed in the hope that it will be useful, but
15 // WITHOUT ANY WARRANTY; without even the implied warranty of
16 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 // GNU General Public License for more details.
18 //
19 // You should have received a copy of the GNU General Public License
20 // along with Rcpp. If not, see <http://www.gnu.org/licenses/>.
21 
22 #ifndef Rcpp__sugar__complex_h
23 #define Rcpp__sugar__complex_h
24 
25 #ifdef HAVE_HYPOT
26 # define RCPP_HYPOT ::hypot
27 #else
28 # define RCPP_HYPOT ::Rf_pythag
29 #endif
30 
31 namespace Rcpp{
32 namespace sugar{
33 
34 
35 template <bool NA, typename RESULT_TYPE, typename T, typename FunPtr>
37  Rcpp::traits::r_sexptype_traits<RESULT_TYPE>::rtype ,
38  NA,
39  SugarComplex<NA,RESULT_TYPE,T,FunPtr>
40  > {
41 public:
42 
44 
45  SugarComplex( FunPtr ptr_, const VEC_TYPE & vec_) : ptr(ptr_), vec(vec_){}
46 
47  inline RESULT_TYPE operator[]( R_xlen_t i) const {
48  Rcomplex x = vec[i] ;
50  return Rcpp::traits::get_na< Rcpp::traits::r_sexptype_traits<RESULT_TYPE>::rtype >() ;
51  return ptr( x );
52  }
53  inline R_xlen_t size() const { return vec.size() ; }
54 
55 private:
56  FunPtr ptr ;
57  const VEC_TYPE& vec ;
58 };
59 } // sugar
60 
61 namespace internal{
62 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;
89 
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  }
137 
138 inline Rcomplex complex__asin(Rcomplex z)
139 {
140  Rcomplex r ;
141  double alpha, bet, t1, t2, x, y;
142  x = z.r;
143  y = z.i;
144  t1 = 0.5 * RCPP_HYPOT(x + 1, y);
145  t2 = 0.5 * RCPP_HYPOT(x - 1, y);
146  alpha = t1 + t2;
147  bet = t1 - t2;
148  r.r = ::asin(bet);
149  r.i = ::log(alpha + ::sqrt(alpha*alpha - 1));
150  if(y < 0 || (y == 0 && x > 1)) r.i *= -1;
151  return r ;
152 }
153 
154 inline Rcomplex complex__acos(Rcomplex z)
155 {
156  Rcomplex r, Asin = complex__asin(z);
157  r.r = M_PI_2 - Asin.r;
158  r.i = - Asin.i;
159  return r ;
160 }
161 
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. */
169 
170 inline Rcomplex complex__atan(Rcomplex z)
171 {
172  Rcomplex r;
173  double x, y;
174  x = z.r;
175  y = z.i;
176  r.r = 0.5 * ::atan(2 * x / ( 1 - x * x - y * y));
177  r.i = 0.25 * ::log((x * x + (y + 1) * (y + 1)) /
178  (x * x + (y - 1) * (y - 1)));
179  if(x*x + y*y > 1) {
180  r.r += M_PI_2;
181  if(x < 0 || (x == 0 && y < 0)) r.r -= M_PI;
182  }
183  return r ;
184 }
185 
186 
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  }
193 
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  }
203 
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  }
213 inline Rcomplex complex__sinh(Rcomplex z)
214 {
215  Rcomplex r, b;
216  b.r = -z.i;
217  b.i = z.r;
218  Rcomplex a = complex__sin(b);
219  r.r = a.i;
220  r.i = -a.r;
221  return r ;
222 }
223 
224 inline Rcomplex complex__tanh(Rcomplex z)
225 {
226  Rcomplex r, b;
227  b.r = -z.i;
228  b.i = z.r;
229  Rcomplex a = complex__tan(b);
230  r.r = a.i;
231  r.i = -a.r;
232  return r ;
233 }
234 
235 
236 
237 } // internal
238 
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  }
249 
250 RCPP_SUGAR_COMPLEX( Re, double )
251 RCPP_SUGAR_COMPLEX( Im, double )
252 RCPP_SUGAR_COMPLEX( Mod, double )
253 RCPP_SUGAR_COMPLEX( Arg, double )
254 RCPP_SUGAR_COMPLEX( Conj, Rcomplex )
255 RCPP_SUGAR_COMPLEX( exp, Rcomplex )
256 RCPP_SUGAR_COMPLEX( log, Rcomplex )
257 RCPP_SUGAR_COMPLEX( sqrt, Rcomplex )
258 RCPP_SUGAR_COMPLEX( cos, Rcomplex )
259 RCPP_SUGAR_COMPLEX( sin, Rcomplex )
260 RCPP_SUGAR_COMPLEX( tan, Rcomplex )
261 RCPP_SUGAR_COMPLEX( acos, Rcomplex )
262 RCPP_SUGAR_COMPLEX( asin, Rcomplex )
263 RCPP_SUGAR_COMPLEX( atan, Rcomplex )
264 RCPP_SUGAR_COMPLEX( acosh, Rcomplex )
265 RCPP_SUGAR_COMPLEX( asinh, Rcomplex )
266 RCPP_SUGAR_COMPLEX( atanh, Rcomplex )
267 RCPP_SUGAR_COMPLEX( cosh, Rcomplex )
268 RCPP_SUGAR_COMPLEX( sinh, Rcomplex )
269 RCPP_SUGAR_COMPLEX( tanh, Rcomplex )
270 
271 #undef RCPP_SUGAR_COMPLEX
272 
273 } // Rcpp
274 #endif
275 
void sqrt(InputIterator begin, InputIterator end, OutputIterator out)
Definition: algorithm.h:479
SugarComplex(FunPtr ptr_, const VEC_TYPE &vec_)
Definition: complex.h:45
Rcomplex complex__cosh(Rcomplex z)
Definition: complex.h:110
R_xlen_t size() const
Definition: VectorBase.h:49
Rcomplex complex__log(Rcomplex x)
Definition: complex.h:80
bool is_na< CPLXSXP >(Rcomplex x)
Definition: is_na.h:46
void exp(InputIterator begin, InputIterator end, OutputIterator out)
Definition: algorithm.h:474
Rcomplex complex__cos(Rcomplex z)
Definition: complex.h:104
Rcomplex complex__atan(Rcomplex z)
Definition: complex.h:170
double complex__Im(Rcomplex x)
Definition: complex.h:63
double complex__Arg(Rcomplex x)
Definition: complex.h:71
Rcomplex complex__sin(Rcomplex z)
Definition: complex.h:116
Rcomplex complex__acosh(Rcomplex z)
Definition: complex.h:187
double complex__Re(Rcomplex x)
Definition: complex.h:62
Rcomplex complex__Conj(Rcomplex x)
Definition: complex.h:65
RESULT_TYPE operator[](R_xlen_t i) const
Definition: complex.h:47
Rcpp::VectorBase< CPLXSXP, NA, T > VEC_TYPE
Definition: complex.h:43
Rcomplex complex__sinh(Rcomplex z)
Definition: complex.h:213
Rcomplex complex__asinh(Rcomplex z)
Definition: complex.h:194
#define RCPP_SUGAR_COMPLEX(__NAME__, __OUT__)
Definition: complex.h:239
double complex__Mod(Rcomplex x)
Definition: complex.h:64
R_xlen_t size() const
Definition: complex.h:53
Rcomplex complex__exp(Rcomplex x)
Definition: complex.h:73
const VEC_TYPE & vec
Definition: complex.h:57
void log(InputIterator begin, InputIterator end, OutputIterator out)
Definition: algorithm.h:469
Rcomplex complex__acos(Rcomplex z)
Definition: complex.h:154
#define RCPP_HYPOT
Definition: complex.h:28
Rcpp API.
Definition: algo.h:28
Rcomplex complex__tan(Rcomplex z)
Definition: complex.h:122
Rcomplex complex__sqrt(Rcomplex z)
Definition: complex.h:86
Rcomplex complex__atanh(Rcomplex z)
Definition: complex.h:204
Rcomplex complex__asin(Rcomplex z)
Definition: complex.h:138
Rcomplex complex__tanh(Rcomplex z)
Definition: complex.h:224