R 中的 Collatz 猜想
Collatz conjecture in R
我仍在教一些 R,主要是给自己(和我的学生)。
这是 Collatz 序列在 R 中的实现:
f <- function(n)
{
# construct the entire Collatz path starting from n
if (n==1) return(1)
if (n %% 2 == 0) return(c(n, f(n/2)))
return(c(n, f(3*n + 1)))
}
调用 f(13) 我得到
13, 40, 20, 10, 5, 16, 8, 4, 2, 1
但是请注意,这里矢量的大小是动态增长的。此类移动往往会导致代码效率低下。有没有更高效的版本?
在Python我会用
def collatz(n):
assert isinstance(n, int)
assert n >= 1
def __colla(n):
while n > 1:
yield n
if n % 2 == 0:
n = int(n / 2)
else:
n = int(3 * n + 1)
yield 1
return list([x for x in __colla(n)])
我找到了一种无需先验指定维度即可写入向量的方法。因此,解决方案可能是
collatz <-function(n)
{
stopifnot(n >= 1)
# define a vector without specifying the length
x = c()
i = 1
while (n > 1)
{
x[i] = n
i = i + 1
n = ifelse(n %% 2, 3*n + 1, n/2)
}
x[i] = 1
# now "cut" the vector
dim(x) = c(i)
return(x)
}
我很想知道通过 Rcpp
实现的 C++ 与您的两种基本 R 方法相比如何。这是我的结果。
首先让我们定义一个函数collatz_Rcpp
,它returns 给定整数n
的冰雹序列。 (非递归)实现改编自 Rosetta Code.
library(Rcpp)
cppFunction("
std::vector<int> collatz_Rcpp(int i) {
std::vector<int> v;
while(true) {
v.push_back(i);
if (i == 1) break;
i = (i % 2) ? (3 * i + 1) : (i / 2);
}
return v;
}
")
我们现在 运行 使用您的基础 R 和 Rcpp
实现进行 microbenchmark
分析。我们计算前 10000 个整数的 Hailstone 序列
# base R implementation
collatz_R <- function(n) {
# construct the entire Collatz path starting from n
if (n==1) return(1)
if (n %% 2 == 0) return(c(n, collatz(n/2)))
return(c(n, collatz(3*n + 1)))
}
# "updated" base R implementation
collatz_R_updated <-function(n) {
stopifnot(n >= 1)
# define a vector without specifying the length
x = c()
i = 1
while (n > 1) {
x[i] = n
i = i + 1
n = ifelse(n %% 2, 3*n + 1, n/2)
}
x[i] = 1
# now "cut" the vector
dim(x) = c(i)
return(x)
}
library(microbenchmark)
n <- 10000
res <- microbenchmark(
baseR = sapply(1:n, collatz_R),
baseR_updated = sapply(1:n, collatz_R_updated),
Rcpp = sapply(1:n, collatz_Rcpp))
res
# expr min lq mean median uq max
# baseR 65.68623 73.56471 81.42989 77.46592 83.87024 193.2609
#baseR_updated 3861.99336 3997.45091 4240.30315 4122.88577 4348.97153 5463.7787
# Rcpp 36.52132 46.06178 51.61129 49.27667 53.10080 168.9824
library(ggplot2)
autoplot(res)
(非递归)Rcpp
实现似乎比原始(递归)基础 R 实现快 30% 左右。 "updated"(非递归)base R 实现比原始(递归)base R 方法慢得多(由于 baseR_updated
,microbenchmark
在我的 MacBook Air 上需要大约 10 分钟才能完成).
我仍在教一些 R,主要是给自己(和我的学生)。
这是 Collatz 序列在 R 中的实现:
f <- function(n)
{
# construct the entire Collatz path starting from n
if (n==1) return(1)
if (n %% 2 == 0) return(c(n, f(n/2)))
return(c(n, f(3*n + 1)))
}
调用 f(13) 我得到 13, 40, 20, 10, 5, 16, 8, 4, 2, 1
但是请注意,这里矢量的大小是动态增长的。此类移动往往会导致代码效率低下。有没有更高效的版本?
在Python我会用
def collatz(n):
assert isinstance(n, int)
assert n >= 1
def __colla(n):
while n > 1:
yield n
if n % 2 == 0:
n = int(n / 2)
else:
n = int(3 * n + 1)
yield 1
return list([x for x in __colla(n)])
我找到了一种无需先验指定维度即可写入向量的方法。因此,解决方案可能是
collatz <-function(n)
{
stopifnot(n >= 1)
# define a vector without specifying the length
x = c()
i = 1
while (n > 1)
{
x[i] = n
i = i + 1
n = ifelse(n %% 2, 3*n + 1, n/2)
}
x[i] = 1
# now "cut" the vector
dim(x) = c(i)
return(x)
}
我很想知道通过 Rcpp
实现的 C++ 与您的两种基本 R 方法相比如何。这是我的结果。
首先让我们定义一个函数collatz_Rcpp
,它returns 给定整数n
的冰雹序列。 (非递归)实现改编自 Rosetta Code.
library(Rcpp)
cppFunction("
std::vector<int> collatz_Rcpp(int i) {
std::vector<int> v;
while(true) {
v.push_back(i);
if (i == 1) break;
i = (i % 2) ? (3 * i + 1) : (i / 2);
}
return v;
}
")
我们现在 运行 使用您的基础 R 和 Rcpp
实现进行 microbenchmark
分析。我们计算前 10000 个整数的 Hailstone 序列
# base R implementation
collatz_R <- function(n) {
# construct the entire Collatz path starting from n
if (n==1) return(1)
if (n %% 2 == 0) return(c(n, collatz(n/2)))
return(c(n, collatz(3*n + 1)))
}
# "updated" base R implementation
collatz_R_updated <-function(n) {
stopifnot(n >= 1)
# define a vector without specifying the length
x = c()
i = 1
while (n > 1) {
x[i] = n
i = i + 1
n = ifelse(n %% 2, 3*n + 1, n/2)
}
x[i] = 1
# now "cut" the vector
dim(x) = c(i)
return(x)
}
library(microbenchmark)
n <- 10000
res <- microbenchmark(
baseR = sapply(1:n, collatz_R),
baseR_updated = sapply(1:n, collatz_R_updated),
Rcpp = sapply(1:n, collatz_Rcpp))
res
# expr min lq mean median uq max
# baseR 65.68623 73.56471 81.42989 77.46592 83.87024 193.2609
#baseR_updated 3861.99336 3997.45091 4240.30315 4122.88577 4348.97153 5463.7787
# Rcpp 36.52132 46.06178 51.61129 49.27667 53.10080 168.9824
library(ggplot2)
autoplot(res)
(非递归)Rcpp
实现似乎比原始(递归)基础 R 实现快 30% 左右。 "updated"(非递归)base R 实现比原始(递归)base R 方法慢得多(由于 baseR_updated
,microbenchmark
在我的 MacBook Air 上需要大约 10 分钟才能完成).