Rcpp Version 1.0.9
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 - 2018 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 namespace Rcpp{
26 namespace sugar{
27 
28 
29 template <bool NA, typename RESULT_TYPE, typename T, typename FunPtr>
31  Rcpp::traits::r_sexptype_traits<RESULT_TYPE>::rtype ,
32  NA,
33  SugarComplex<NA,RESULT_TYPE,T,FunPtr>
34  > {
35 public:
36 
38 
39  SugarComplex( FunPtr ptr_, const VEC_TYPE & vec_) : ptr(ptr_), vec(vec_){}
40 
41  inline RESULT_TYPE operator[]( R_xlen_t i) const {
42  Rcomplex x = vec[i] ;
44  return Rcpp::traits::get_na< Rcpp::traits::r_sexptype_traits<RESULT_TYPE>::rtype >() ;
45  return ptr( x );
46  }
47  inline R_xlen_t size() const { return vec.size() ; }
48 
49 private:
50  FunPtr ptr ;
51  const VEC_TYPE& vec ;
52 };
53 } // sugar
54 
55 namespace internal{
56 inline double complex__Re( Rcomplex x){ return x.r ; }
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;
83 
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 }
131 
132 inline Rcomplex complex__asin(Rcomplex z)
133 {
134  Rcomplex r ;
135  double alpha, bet, t1, t2, x, y;
136  x = z.r;
137  y = z.i;
138  t1 = 0.5 * ::hypot(x + 1, y);
139  t2 = 0.5 * ::hypot(x - 1, y);
140  alpha = t1 + t2;
141  bet = t1 - t2;
142  r.r = ::asin(bet);
143  r.i = ::log(alpha + ::sqrt(alpha*alpha - 1));
144  if(y < 0 || (y == 0 && x > 1)) r.i *= -1;
145  return r ;
146 }
147 
148 inline Rcomplex complex__acos(Rcomplex z)
149 {
150  Rcomplex r, Asin = complex__asin(z);
151  r.r = M_PI_2 - Asin.r;
152  r.i = - Asin.i;
153  return r ;
154 }
155 
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. */
163 
164 inline Rcomplex complex__atan(Rcomplex z)
165 {
166  Rcomplex r;
167  double x, y;
168  x = z.r;
169  y = z.i;
170  r.r = 0.5 * ::atan(2 * x / ( 1 - x * x - y * y));
171  r.i = 0.25 * ::log((x * x + (y + 1) * (y + 1)) /
172  (x * x + (y - 1) * (y - 1)));
173  if(x*x + y*y > 1) {
174  r.r += M_PI_2;
175  if(x < 0 || (x == 0 && y < 0)) r.r -= M_PI;
176  }
177  return r ;
178 }
179 
180 
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 }
187 
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 }
197 
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 }
207 inline Rcomplex complex__sinh(Rcomplex z)
208 {
209  Rcomplex r, b;
210  b.r = -z.i;
211  b.i = z.r;
212  Rcomplex a = complex__sin(b);
213  r.r = a.i;
214  r.i = -a.r;
215  return r ;
216 }
217 
218 inline Rcomplex complex__tanh(Rcomplex z)
219 {
220  Rcomplex r, b;
221  b.r = -z.i;
222  b.i = z.r;
223  Rcomplex a = complex__tan(b);
224  r.r = a.i;
225  r.i = -a.r;
226  return r ;
227 }
228 
229 } // internal
230 
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  }
238 
239 RCPP_SUGAR_COMPLEX( Re, double )
240 RCPP_SUGAR_COMPLEX( Im, double )
241 RCPP_SUGAR_COMPLEX( Mod, double )
242 RCPP_SUGAR_COMPLEX( Arg, double )
243 RCPP_SUGAR_COMPLEX( Conj, Rcomplex )
244 RCPP_SUGAR_COMPLEX( exp, Rcomplex )
245 RCPP_SUGAR_COMPLEX( log, Rcomplex )
246 RCPP_SUGAR_COMPLEX( sqrt, Rcomplex )
247 RCPP_SUGAR_COMPLEX( cos, Rcomplex )
248 RCPP_SUGAR_COMPLEX( sin, Rcomplex )
249 RCPP_SUGAR_COMPLEX( tan, Rcomplex )
250 RCPP_SUGAR_COMPLEX( acos, Rcomplex )
251 RCPP_SUGAR_COMPLEX( asin, Rcomplex )
252 RCPP_SUGAR_COMPLEX( atan, Rcomplex )
253 RCPP_SUGAR_COMPLEX( acosh, Rcomplex )
254 RCPP_SUGAR_COMPLEX( asinh, Rcomplex )
255 RCPP_SUGAR_COMPLEX( atanh, Rcomplex )
256 RCPP_SUGAR_COMPLEX( cosh, Rcomplex )
257 RCPP_SUGAR_COMPLEX( sinh, Rcomplex )
258 RCPP_SUGAR_COMPLEX( tanh, Rcomplex )
259 
260 #undef RCPP_SUGAR_COMPLEX
261 
262 } // Rcpp
263 #endif
264 
R_xlen_t size() const
Definition: VectorBase.h:49
R_xlen_t size() const
Definition: complex.h:47
const VEC_TYPE & vec
Definition: complex.h:51
RESULT_TYPE operator[](R_xlen_t i) const
Definition: complex.h:41
SugarComplex(FunPtr ptr_, const VEC_TYPE &vec_)
Definition: complex.h:39
Rcpp::VectorBase< CPLXSXP, NA, T > VEC_TYPE
Definition: complex.h:37
double hypot(double a, double b)
Definition: Rmath.h:220
void exp(InputIterator begin, InputIterator end, OutputIterator out)
Definition: algorithm.h:474
void log(InputIterator begin, InputIterator end, OutputIterator out)
Definition: algorithm.h:469
void sqrt(InputIterator begin, InputIterator end, OutputIterator out)
Definition: algorithm.h:479
Rcomplex complex__Conj(Rcomplex x)
Definition: complex.h:59
Rcomplex complex__cosh(Rcomplex z)
Definition: complex.h:104
Rcomplex complex__exp(Rcomplex x)
Definition: complex.h:67
Rcomplex complex__atan(Rcomplex z)
Definition: complex.h:164
double complex__Im(Rcomplex x)
Definition: complex.h:57
double complex__Re(Rcomplex x)
Definition: complex.h:56
Rcomplex complex__tan(Rcomplex z)
Definition: complex.h:116
Rcomplex complex__atanh(Rcomplex z)
Definition: complex.h:198
Rcomplex complex__cos(Rcomplex z)
Definition: complex.h:98
Rcomplex complex__asin(Rcomplex z)
Definition: complex.h:132
Rcomplex complex__tanh(Rcomplex z)
Definition: complex.h:218
double complex__Mod(Rcomplex x)
Definition: complex.h:58
Rcomplex complex__acos(Rcomplex z)
Definition: complex.h:148
Rcomplex complex__sinh(Rcomplex z)
Definition: complex.h:207
Rcomplex complex__sqrt(Rcomplex z)
Definition: complex.h:80
Rcomplex complex__acosh(Rcomplex z)
Definition: complex.h:181
Rcomplex complex__log(Rcomplex x)
Definition: complex.h:74
Rcomplex complex__asinh(Rcomplex z)
Definition: complex.h:188
double complex__Arg(Rcomplex x)
Definition: complex.h:65
Rcomplex complex__sin(Rcomplex z)
Definition: complex.h:110
bool is_na< CPLXSXP >(Rcomplex x)
Definition: is_na.h:47
Rcpp API.
Definition: algo.h:28
#define RCPP_SUGAR_COMPLEX(__NAME__, __OUT__)
Definition: complex.h:231