ダメ出し：アトキンの篩

2011-05-11　R でエラトステネスの篩　の後半に，アトキンの篩のプログラムが掲載されている。

ご本人も「多分条件判定減らすとかしてRに合った実装しないとダメなんだと思う」というとおり，ベクトル化をはかる。ちなみに，元のプログラムでは for を避けるためか while で書いているがこれはほとんど意味がない。

limit = 10000000 で 23 秒かかっていたものが，書き直すと 0.4 秒ほどで計算終了ということになった。

プログラムは以下の通り。

atkin <- function(limit = 1e+06, return = FALSE) {
    sqrt.limit <- sqrt(limit)
    isprime <- logical(limit)
    for (z in c(1, 5)) {
        for (y in seq.int(z, sqrt.limit, by = 6)) {
            y2 <- y * y
            max.x <- floor(sqrt((limit - y2) / 4))
            if (1 <= max.x) {
                n <- 4 * (1:max.x) ^ 2 + y2
                isprime[n] <- !isprime[n]
            }
            x <- y + 1
            max.x <- floor(sqrt((limit + y2) / 3))
            if (x <= max.x) {
                n <- 3 * seq.int(x, max.x, by = 2) ^ 2 - y2
                isprime[n] <- !isprime[n]
            }
        }
    }
    for (z in c(2, 4)) {
        for (y in seq.int(z, sqrt.limit, by = 6)) {
            y2 <- y * y
            max.x <- floor(sqrt((limit - y2) / 3))
            if (1 <= max.x) {
                n <- 3 * seq.int(1, max.x, by = 2) ^ 2 + y2
                isprime[n] <- !isprime[n]
            }
            x <- y + 1
            max.x <- floor(sqrt((limit + y2) / 3))
            if (x <= max.x) {
                n <- 3 * seq.int(x, max.x, by = 2) ^ 2 - y2
                isprime[n] <- !isprime[n]
            }
        }
    }
    for (y in seq.int(3, sqrt.limit, by = 6)) {
        y2 <- y * y
        for (x in 1:2) {
            max.x <- floor(sqrt((limit - y2) / 4))
            if (x <= max.x) {
                n <- 4 * seq.int(x, max.x, by = 3) ^ 2 + y2
                isprime[n] <- !isprime[n]
            }
        }
    }
    for (n in seq.int(5, sqrt.limit, by = 2)) {
        if (isprime[n]) {
            isprime[1:(limit %/% (n * n)) * n * n] <- FALSE
        }
    }
    isprime[2] <- isprime[3] <- TRUE
    if (return) {
        which(isprime)
    }
}