如何在R中高效生成黄石序列（A098550）？

我通过Numberphile的YouTube视频了解到黄石序列（OEIS编号A098550），该序列以1和2起始，后续项需遵循以下规则生成：

• 无重复项
• 始终选取最小整数
• gcd(aₙ, aₙ₋₁) = 1
• gcd(aₙ, aₙ₋₂) > 1

其前15项为：1 2 3 4 9 8 15 14 5 6 25 12 35 16 7

我用R写了一个快速粗糙（Q&D）的实现，但生成更长序列时速度极慢，而且代码还假设了序列的最大数值（注：10000项的序列数值不超过5000）。原始代码如下：

library(DescTools)
a <- c(1, 2, 3)
p <- length(a)

# all natural numbers
all_ints <- 1:5000

for (n in p:1000) {
    # rule 1 - remove all number that are in sequence already
    next_a_set <- all_ints[which(!all_ints %in% a)]

    # rule 3 - search the remaining set for numbers that have gcd == 1
    next_a_option <- next_a_set[which(
        sapply(
            next_a_set,
            function(x) GCD(a[n], x)
        ) == 1
    )]

    # rule 4 - search the remaining number for gcd > 1
    next_a <- next_a_option[which(
        sapply(
            next_a_option,
            function(x) GCD(a[n - 1], x)
        ) > 1
    )]

    # select the lowest
    a <- c(a, min(next_a))
    n <- n + 1
}

核心优化思路 #

原始代码效率低下的核心问题：频繁复制向量、低效的已使用数查询、逐元素sapply计算。针对这些问题，给出以下优化方案：

1. 预分配向量，避免频繁内存复制 #

每次用c(a, ...)扩展向量会触发内存复制，预分配足够长度的向量直接赋值能大幅减少开销：

target_length <- 10000
a <- integer(target_length)
a[1:3] <- c(1, 2, 3)

2. 用逻辑向量跟踪已使用数，实现O(1)查询 #

用logical向量标记已使用的整数，替代!all_ints %in% a这种O(n)级别的低效查询：

max_val <- 5000
used <- logical(max_val)
used[a[1:3]] <- TRUE

3. 内置gcd替代第三方包，提前终止判断 #

R 4.0+版本内置了base::gcd函数，无需加载DescTools；同时从最小未使用数开始检查，找到符合条件的就停止遍历，不用计算所有候选。

优化后的完整代码 #

target_length <- 10000
max_val <- 5000

# 预分配序列存储向量
a <- integer(target_length)
a[1:3] <- c(1, 2, 3)

# 标记已使用的数值
used <- logical(max_val)
used[a[1:3]] <- TRUE

for (n in 3:(target_length - 1)) {
    prev1 <- a[n]
    prev2 <- a[n - 1]
    
    # 从最小整数开始遍历，找到第一个符合所有规则的数
    for (x in seq_len(max_val)) {
        if (used[x]) next
        if (gcd(x, prev1) == 1 && gcd(x, prev2) > 1) {
            a[n + 1] <- x
            used[x] <- TRUE
            break
        }
    }
}

进阶优化建议 #

• 动态扩展数值范围：如果不确定max_val，可以在循环中检查是否所有数都已使用，动态扩展used向量；
• 预计算质因数：提前为每个数计算质因数集合，判断gcd条件时只需检查质因数交集，进一步提升效率；
• Rcpp实现核心循环：将内层判断逻辑用C++实现，能把速度提升10-100倍，适合生成十万级以上的序列。

内容的提问来源于stack exchange，提问作者Rob Hanssen