2013-07-06

math.h関連

cpp Rcpp

math.h

#include <RcppArmadilloExtensions/sample.h>

// [[Rcpp::depends(RcppArmadillo)]]

using namespace Rcpp ;

// [[Rcpp::export]]

NumericVector arith(double a, double b) {
	NumericVector x(28);
	
	x[0] = a+b;
	x[1] = a-b;
	x[2] = a*b;
	x[3] = a/b;
	x[4] = fmax(a,b);
	x[5] = round(a);
	x[6] = lgamma(a);
	x[7] = sin(a);
	x[8] = acos(a);
	x[9] = asin(a);
	x[10] = atan(a);
	x[11] = atan2(a,b);
	x[12] = ceil(a);
	x[13] = cos(a);
	x[14] = cosh(a);
	x[15] = exp(a);
	x[16] = fabs(a);
	x[17] = floor(a);
	x[18] = fmod(a,b);
	x[19] = log(a);
	x[20] = log10(a);
	x[21] = M_PI;
	x[22] = pow(a,b);
	x[23] = sin(a);
	x[24] = sinh(a);
	x[25] = sqrt(a);
	x[26] = tan(a);
	x[27] = tanh(a);
	return x;
}

arith.R <- function(a, b) {
	x = rep(0,28);
	
	x[1] = a+b;
	x[2] = a-b;
	x[3] = a*b;
	x[4] = a/b;
	x[5] = max(a,b);
	x[6] = round(a);
	x[7] = log(gamma(a));
	x[8] = sin(a);
	x[9] = acos(a);
	x[10] = asin(a);
	x[11] = atan(a);
	x[12] = atan2(a,b);
	x[13] = ceiling(a);
	x[14] = cos(a);
	x[15] = cosh(a);
	x[16] = exp(a);
	x[17] = abs(a);
	x[18] = floor(a);
	x[19] = a%%b;
	x[20] = log(a);
	x[21] = log10(a);
	x[22] = pi;
	x[23] = a^b;
	x[24] = sin(a);
	x[25] = sinh(a);
	x[26] = sqrt(a);
	x[27] = tan(a);
	x[28] = tanh(a);
	x
}

a <- runif(1)
b <- runif(1)
a.R <- arith.R(a,b)
a.cpp <- arith(a,b)
a.R-a.cpp

> sourceCpp("arith.cpp")
> a <- runif(1)
> b <- runif(1)
> a.R <- arith.R(a,b)
> a.cpp <- arith(a,b)
> a.R-a.cpp
 [1]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
 [7] -2.220446e-16  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
[13]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
[19]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
[25]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00

真偽判断

#include <RcppArmadilloExtensions/sample.h>

// [[Rcpp::depends(RcppArmadillo)]]

using namespace Rcpp ;

// [[Rcpp::export]]

LogicalVector logic() {
	LogicalVector x(4);
	if(M_PI == M_PI){
		x[0] = TRUE;
	}else{
		x[0] = FALSE;
	}
	if(M_PI != log(1)){
		x[1] = FALSE;
	}else{
		x[1] = TRUE;
	}
	if(abs(sin(M_PI)) < 0.000001){
		x[2] = TRUE;
	}else{
		x[2] = FALSE;
	}
	if(abs(sin(1/2*M_PI))-1 < 0.000001){
		x[3] = FALSE;
	}else{
		x[3] = TRUE;
	}
	return x;
}

> sourceCpp("logic.cpp")
> logic()
[1]  TRUE FALSE  TRUE FALSE

台形の面積計算関数

#include <RcppArmadilloExtensions/sample.h>

// [[Rcpp::depends(RcppArmadillo)]]

using namespace Rcpp ;

// [[Rcpp::export]]

NumericVector daikei(NumericVector a,NumericVector b,NumericVector h) {
	int n = a.size();
	NumericVector x(n);
	for(int i=0;i<n;++i){
		x[i] = (a[i]+b[i])*h[i]/2;
	}
	return x;
}

> sourceCpp("daikei.cpp")
> a <- runif(10)
> b <- runif(10)
> h <- runif(10)
> daikei(a,b,h)
 [1] 0.001660718 0.286582070 0.346243689 0.104847673 0.326411176 0.244554508 0.438464121
 [8] 0.581787228 0.020789185 0.064599295
> (a+b)*h/2
 [1] 0.001660718 0.286582070 0.346243689 0.104847673 0.326411176 0.244554508 0.438464121
 [8] 0.581787228 0.020789185 0.064599295

多次元球の表面積と体積

#include <RcppArmadilloExtensions/sample.h>

// [[Rcpp::depends(RcppArmadillo)]]

using namespace Rcpp ;

// [[Rcpp::export]]

NumericMatrix sphere_n(int n,double r) {
	NumericMatrix x(n,2);
	for(int i=0;i<n;++i){
		x(i,0) = 2*pow(M_PI,(i+1.0)/2)/exp(lgamma((i+1.0)/2))*pow(r,i);
		x(i,1) = pow(M_PI,(i+1.0)/2)/exp(lgamma((i+1.0)/2+1))*pow(r,i+1);
	}
	return x;
}

sphere_n.R <- function(n,r){
	ret <- matrix(0,n,2)
	for(i in 1:n){
		ret[i,1] <- 2*pi^(i/2)/gamma(i/2)*r^(i-1)
		ret[i,2] <- pi^(i/2)/gamma(i/2+1)*r^i
	}
	ret
}

sphere_n.R(5,2)
sphere_n(5,2)

> sourceCpp("sphere_n.cpp")
> sphere_n.R(5,2)
          [,1]      [,2]
[1,]   2.00000   4.00000
[2,]  12.56637  12.56637
[3,]  50.26548  33.51032
[4,] 157.91367  78.95684
[5,] 421.10312 168.44125
> sphere_n(5,2)
          [,1]      [,2]
[1,]   2.00000   4.00000
[2,]  12.56637  12.56637
[3,]  50.26548  33.51032
[4,] 157.91367  78.95684
[5,] 421.10312 168.44125