下面是一个比你的版本快20倍的版本，并附有对更改的评论：

# Set a to the final length from the start.
a <- c(1, 2, 3, rep(NA, 997))
p <- 3

# Define a vectorized gcd() function.  We'll be testing
# lots of gcds at once.  This uses the Euclidean algorithm.
gcd <- function(x, y) { # vectorized gcd

  while (any(y != 0)) {
    x1 <- ifelse(y == 0, x, y)
    y <- ifelse(y == 0, 0, x %% y)
    x <- x1
  }
  x
}

# Guess at a reasonably large vector to work from,
# but we'll grow it later if not big enough.
allnum <- 1:1000

# Keep a logical record of what has been used
used <- c(rep(TRUE, 3), rep(FALSE, length(allnum) - 3))

for (n in p:1000) {
  # rule 1 - remove all number that are in sequence already
  # nothing to do -- used already records that.
  
  repeat {
    # rule 3 - search the remaining set for numbers that have gcd == 1

    keep <- !used & gcd(a[n], allnum) == 1

    # rule 4 - search the remaining number for gcd > 1
    keep <- keep & gcd(a[n-1], allnum) > 1
  
    # If we found anything, break out of this loop
    if (any(keep))
      break
    
    # Otherwise, make the set of possible values twice as big,
    # and try again
    allnum <- seq_len(2*length(allnum))
    used <- c(used, rep(FALSE, length(used)))
  }
  # select the lowest
  newval <- which.max(keep)
  # Assign into the appropriate place
  a[n+1] <- newval
  # Record that it has been used
  used[newval] <- TRUE  
}

如果你分析它，你会发现它大部分时间都花在gcd()函数上，你可以用C或C++重写它，这样做会更快。

这里最大的变化是预分配和限制搜索尚未使用的号码。

library(numbers)
N <- 5e3
a <- integer(N)
a[1:3] <- 1:3
b <- logical(N) # which numbers have been used already?
b[1:3] <- TRUE
NN <- 1:N

system.time({
  for (n in 4:N) {
    a1 <- a[n - 1L]
    a2 <- a[n - 2L]
    for (k in NN[!b]) {
      if (GCD(k, a1) == 1L & GCD(k, a2) > 1L) {
        a[n] <- k
        b[k] <- TRUE
        break
      }
    }
    if (!a[n]) {
      a <- a[1:(n - 1L)]
      break
    }
  }
})
#>    user  system elapsed 
#>    1.28    0.00    1.28
length(a)
#> [1] 1137

有关快速C++算法，请参见here。