想要摆脱循环的建议

Pet*_*ica 2 loops r vectorization

我编写了一个程序,可以解决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
     }
   }
}
Run Code Online (Sandbox Code Playgroud)

谢谢!

Mar*_*gan 9

我通过创建一个向量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
}
Run Code Online (Sandbox Code Playgroud)

这可以被优化为(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
}
Run Code Online (Sandbox Code Playgroud)

用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
Run Code Online (Sandbox Code Playgroud)

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

x[i] <- 3 * x[i] + 1
count[idx[i]] <- count[idx[i]] + 1L
x <- x / 2
count[idx] <- count[idx] + 1
Run Code Online (Sandbox Code Playgroud)

这就是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
Run Code Online (Sandbox Code Playgroud)

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