R: Rcppでマンデルブロー集合計算を高速化 [統計]
承前。やはりインタプリタではおそいので、Rcppをつかって高速化してみた。劇的に速くなった。
[追記] 画像とコードを更新しました(10/19 23:20)。
コード
[追記] 画像とコードを更新しました(10/19 23:20)。
コード
library(Rcpp) library(inline) ## C++ source mandelbrot.src <- " Rcpp::NumericVector x(xx); Rcpp::NumericVector y(yy); unsigned short t = Rcpp::as<unsigned short>(threshold); int nx = x.size(), ny = y.size(); Rcpp::NumericVector u(nx); if (nx != ny) { return Rcpp::wrap(-1); } else { for (int i = 0; i < nx; i++) { std::complex<double> z(0.0, 0.0); std::complex<double> c(x[i], y[i]); unsigned short k = 0; while (k < t) { z = z * z + c; if (abs(z) > 2.0) break; k++; } u[i] = k; } return Rcpp::wrap(u); } " ## Call C++ mandelbrot.c <- cfunction(signature(xx = "numeric", yy = "numeric", threshold = "numeric"), mandelbrot.src, includes = "#include <complex>", Rcpp = TRUE) ## R function mandelbrot <- function(min.x, max.x, min.y, max.y, px = 256, py = round((max.y - min.y) / (max.x - min.x) * px), threshold = 255) { if (min.x > max.x) { a <- min.x min.x <- max.x max.x <- a } if (min.y > max.y) { a <- min.y min.y <- max.y max.y <- a } if (px < 0) px <- -px if (py < 0) py <- -py rx <- seq(min.x, max.x, length = px) ry <- seq(min.y, max.y, length = py) xx <- rep(rx, py) yy <- rep(ry, each = px) z <- mandelbrot.c(xx, yy, threshold) matz <- matrix(z, nrow = px) } ## sample parameters px <- 768 py <- 768 threshold <- 1023 min.x <- -0.7625 max.x <- -0.7605 min.y <- -0.0860 max.y <- -0.0840 matz <- mandelbrot(min.x, max.x, min.y, max.y, px, py, threshold = threshold) png(file = "Mandelbrot%03d.png", width = 900, height = 900) image(seq(min.x, max.x, length = nrow(matz)), seq(min.y, max.y, length = ncol(matz)), matz, col = rep(rainbow(256), (threshold + 1) / 256), xlab = "x", ylab = "y") dev.off()
コメント 0