考拉兹猜想在R中的实现

Question

考拉兹猜想在R中的实现

3

我仍在自学R语言并向我的学生教授。

以下是R语言中Collatz数列的实现：

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)
}

- tschm

2

a-look-a-r-vectorization-through-the-collatz-conjecture - Patrick Artner

2

r-and-the-collatz-conjecture-part-2.html - Patrick Artner

2

谷歌 R Collatz 高效 - 前两个结果 - 应该会有帮助 :) 对 R 没有任何概念 - Patrick Artner

很遗憾，不是真的。我之前看过那些链接。它们甚至没有试图明确构造猜想中数字的序列。它们只是更新计数器以获取序列长度。它们运行于参数向量上。尽管如此，非常有趣。谢谢... - tschm

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Maurits Evers · Accepted Answer

我很好奇通过Rcpp实现的C++方法与您的两种基本R方法相比如何。这是我的结果。

首先，让我们定义一个函数collatz_Rcpp，它返回给定整数n的Hailstone序列。这个（非递归）实现是从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实现运行微基准分析，计算前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％。 "更新"的（非递归）基础R实现比原始（递归）基础R方法慢得多（由于baseR_updated，microbenchmark在我的MacBook Air上需要大约10分钟才能完成)。