[Rcpp-commits] r1999 - in pkg/Rcpp/inst: include/Rcpp/sugar/functions unitTests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Aug 13 13:25:06 CEST 2010
Author: romain
Date: 2010-08-13 13:25:05 +0200 (Fri, 13 Aug 2010)
New Revision: 1999
Modified:
pkg/Rcpp/inst/include/Rcpp/sugar/functions/complex.h
pkg/Rcpp/inst/unitTests/runit.sugar.R
Log:
some more complex functions
Modified: pkg/Rcpp/inst/include/Rcpp/sugar/functions/complex.h
===================================================================
--- pkg/Rcpp/inst/include/Rcpp/sugar/functions/complex.h 2010-08-13 09:48:25 UTC (rev 1998)
+++ pkg/Rcpp/inst/include/Rcpp/sugar/functions/complex.h 2010-08-13 11:25:05 UTC (rev 1999)
@@ -112,8 +112,106 @@
r.i = - ::sin(-z.i) * ::sinh(z.r);
return r ;
}
+ inline Rcomplex complex__sin(Rcomplex z){
+ Rcomplex r ;
+ r.r = ::sin(z.r) * ::cosh(z.i);
+ r.i = ::cos(z.r) * ::sinh(z.i);
+ return r;
+ }
+ inline Rcomplex complex__tan(Rcomplex z){
+ Rcomplex r ;
+ double x2, y2, den;
+ x2 = 2.0 * z.r;
+ y2 = 2.0 * z.i;
+ den = ::cos(x2) + ::cosh(y2);
+ r.r = ::sin(x2)/den;
+ /* any threshold between -log(DBL_EPSILON)
+ and log(DBL_XMAX) will do*/
+ if (ISNAN(y2) || ::fabs(y2) < 50.0)
+ r.i = ::sinh(y2)/den;
+ else
+ r.i = (y2 <0 ? -1.0 : 1.0);
+ return r ;
+ }
+inline Rcomplex complex__asin(Rcomplex z)
+{
+ Rcomplex r ;
+ double alpha, bet, t1, t2, x, y;
+ x = z.r;
+ y = z.i;
+ t1 = 0.5 * RCPP_HYPOT(x + 1, y);
+ t2 = 0.5 * RCPP_HYPOT(x - 1, y);
+ alpha = t1 + t2;
+ bet = t1 - t2;
+ r.r = ::asin(bet);
+ r.i = ::log(alpha + ::sqrt(alpha*alpha - 1));
+ if(y < 0 || (y == 0 && x > 1)) r.i *= -1;
+ return r ;
+}
+inline Rcomplex complex__acos(Rcomplex z)
+{
+ Rcomplex r, Asin = complex__asin(z);
+ r.r = M_PI_2 - Asin.r;
+ r.i = - Asin.i;
+ return r ;
+}
+
+ /* Complex Arctangent Function */
+ /* Equation (4.4.39) Abramowitz and Stegun */
+ /* with additional terms to force the branch cuts */
+ /* to agree with figure 4.4, p79. Continuity */
+ /* on the branch cuts (pure imaginary axis; x==0, |y|>1) */
+ /* is standard: z_asin() is continuous from the right */
+ /* if y >= 1, and continuous from the left if y <= -1. */
+
+inline Rcomplex complex__atan(Rcomplex z)
+{
+ Rcomplex r;
+ double x, y;
+ x = z.r;
+ y = z.i;
+ r.r = 0.5 * ::atan(2 * x / ( 1 - x * x - y * y));
+ r.i = 0.25 * ::log((x * x + (y + 1) * (y + 1)) /
+ (x * x + (y - 1) * (y - 1)));
+ if(x*x + y*y > 1) {
+ r.r += M_PI_2;
+ if(x < 0 || (x == 0 && y < 0)) r.r -= M_PI;
+ }
+ return r ;
+}
+
+
+ inline Rcomplex complex__acosh(Rcomplex z){
+ Rcomplex r, a = complex__acos(z);
+ r.r = -a.i;
+ r.i = a.r;
+ return r ;
+ }
+
+ inline Rcomplex complex__asinh(Rcomplex z){
+ Rcomplex r, b;
+ b.r = -z.i;
+ b.i = z.r;
+ Rcomplex a = complex__asin(b);
+ r.r = a.i;
+ r.i = -a.r;
+ return r ;
+ }
+
+ inline Rcomplex complex__atanh(Rcomplex z){
+ Rcomplex r, b;
+ b.r = -z.i;
+ b.i = z.r;
+ Rcomplex a = complex__atan(b);
+ r.r = a.i;
+ r.i = -a.r;
+ return r ;
+ }
+
+
+
} // internal
#define RCPP_SUGAR_COMPLEX(__NAME__,__OUT__) \
@@ -135,8 +233,18 @@
RCPP_SUGAR_COMPLEX( log, Rcomplex )
RCPP_SUGAR_COMPLEX( sqrt, Rcomplex )
RCPP_SUGAR_COMPLEX( cos, Rcomplex )
+RCPP_SUGAR_COMPLEX( sin, Rcomplex )
+RCPP_SUGAR_COMPLEX( tan, Rcomplex )
+RCPP_SUGAR_COMPLEX( acos, Rcomplex )
+RCPP_SUGAR_COMPLEX( asin, Rcomplex )
+RCPP_SUGAR_COMPLEX( atan, Rcomplex )
+RCPP_SUGAR_COMPLEX( acosh, Rcomplex )
+RCPP_SUGAR_COMPLEX( asinh, Rcomplex )
+RCPP_SUGAR_COMPLEX( atanh, Rcomplex )
+
RCPP_SUGAR_COMPLEX( cosh, Rcomplex )
+
#undef RCPP_SUGAR_COMPLEX
} // Rcpp
Modified: pkg/Rcpp/inst/unitTests/runit.sugar.R
===================================================================
--- pkg/Rcpp/inst/unitTests/runit.sugar.R 2010-08-13 09:48:25 UTC (rev 1998)
+++ pkg/Rcpp/inst/unitTests/runit.sugar.R 2010-08-13 11:25:05 UTC (rev 1999)
@@ -473,15 +473,23 @@
'
ComplexVector cx( x );
return List::create(
- _["Re"] = Re( cx ),
- _["Im"] = Im( cx ),
- _["Conj"] = Conj( cx ),
- _["Mod"] = Mod( cx ),
- _["exp"] = exp( cx ),
- _["log"] = log( cx ),
- _["sqrt"] = sqrt( cx ),
- _["cos"] = cos( cx ),
- _["cosh"] = cosh( cx )
+ _["Re"] = Re( cx ),
+ _["Im"] = Im( cx ),
+ _["Conj"] = Conj( cx ),
+ _["Mod"] = Mod( cx ),
+ _["exp"] = exp( cx ),
+ _["log"] = log( cx ),
+ _["sqrt"] = sqrt( cx ),
+ _["cos"] = cos( cx ),
+ _["sin"] = sin( cx ),
+ _["tan"] = tan( cx ),
+ _["acos"] = acos( cx ),
+ _["asin"] = asin( cx ),
+ _["atan"] = atan( cx ),
+ // _["acosh"] = acosh( cx ),
+ _["asinh"] = asinh( cx ),
+ _["atanh"] = atanh( cx ),
+ _["cosh"] = cosh( cx )
) ;
'
),
@@ -1073,16 +1081,25 @@
x <- c( rnorm(10), NA ) + 1i*c( rnorm(10), NA )
fx <- .rcpp.sugar$runit_complex
checkEquals( fx(x), list(
- Re = Re(x),
- Im = Im(x),
- Conj = Conj(x),
- Mod = Mod(x),
- exp = exp(x),
- log = log(x),
- sqrt = sqrt(x),
- cos = cos(x),
- cosh = cosh(x)
- )
+ Re = Re(x),
+ Im = Im(x),
+ Conj = Conj(x),
+ Mod = Mod(x),
+ exp = exp(x),
+ log = log(x),
+ sqrt = sqrt(x),
+ cos = cos(x),
+ sin = sin(x),
+ tan = tan(x),
+ acos = acos(x),
+ asin = asin(x),
+ atan = atan(x),
+ # acosh = acosh(x),
+ asinh = asinh(x),
+ atanh = atanh(x),
+ cosh = cosh(x)
+
+ )
)
}
More information about the Rcpp-commits
mailing list