r - 关于摆脱循环的建议

Question

我编写了一个程序来解决 3n + 1 问题(又名“奇妙的数字”和其他各种问题)。但它有一个双循环。我如何对其进行矢量化？

代码是

count <- vector("numeric", 100000)
L <- length(count)

for (i in 1:L)
{
x <- i
   while (x > 1)
   {
   if (round(x/2) == x/2) 
     {
     x <- x/2
     count[i] <- count[i] + 1 
     } else
     {
     x <- 3*x + 1
     count[i] <- count[i] + 1
     }
   }
}

谢谢!

Answer 1

我通过创建一个向量 x 来彻底改变这个过程，其中第 i 个元素是算法每次迭代后的值。结果比较好理解为

f1 <- function(L) {
    x <- seq_len(L)
    count <- integer(L)
    while (any(i <- x > 1)) {
        count[i] <- count[i] + 1L
        x <- ifelse(round(x/2) == x/2, x / 2, 3 * x + 1) * i
    }
    count
}

这可以优化为 (a) 仅跟踪那些仍在发挥作用的值(通过 idx)和 (b) 避免不必要的操作，例如，ifelse 计算所有 x 值的两个参数，x/2 计算两次。

f2 <- function(L) {
    idx <- x <- seq_len(L)
    count <- integer(L)
    while (length(x)) {
        ix <- x > 1
        x <- x[ix]
        idx <- idx[ix]
        count[idx] <- count[idx] + 1L
        i <- as.logical(x %% 2)
        x[i] <- 3 * x[i] + 1
        i <- !i
        x[i] <- x[i] / 2
    }
    count
}

有了f0这个原函数，我有

> L <- 10000
> system.time(ans0 <- f0(L))
   user  system elapsed 
  7.785   0.000   7.812 
> system.time(ans1 <- f1(L))
   user  system elapsed 
  1.738   0.000   1.741 
> identical(ans0, ans1)
[1] TRUE
> system.time(ans2 <- f2(L))
   user  system elapsed 
  0.301   0.000   0.301 
> identical(ans1, ans2)
[1] TRUE

一个调整是将奇数更新为 3 * x[i] + 1 然后无条件地除以二

x[i] <- 3 * x[i] + 1
count[idx[i]] <- count[idx[i]] + 1L
x <- x / 2
count[idx] <- count[idx] + 1

用这个作为 f3(不知道为什么今天早上 f2 变慢了!)我明白了

> system.time(ans2 <- f2(L))
   user  system elapsed 
   0.36    0.00    0.36 
> system.time(ans3 <- f3(L))
   user  system elapsed 
  0.201   0.003   0.206 
> identical(ans2, ans3)
[1] TRUE

似乎在除以二阶段可以采取更大的步骤，例如，8 是 2^3，所以我们可以采取 3 步(加 3 计数)并完成，20 是 2^2 * 5所以我们可以采取两个步骤并在 5 处进入下一个迭代。实现？