裏 RjpWiki

Julia ときどき R, Python によるコンピュータプログラム,コンピュータ・サイエンス,統計学

そんなものじゃないだろう

2012年02月13日 | 雑感

>  記者会見した長沢寛道・同研究科長は「高い値の放射性セシウムが検出された水田で作付けし、データを正確に把握することは、福島県の農業復興に必要だ」と話した。(東京大大学院農学生命科学研究科 長沢寛道研究科長)

作付けして,データを取ってもらえれば科学的知見を得るのに役立つという考えなのだろう。

しかし,農業に携わる人は,金銭的保証をしてくれればよいなどとは思っていないはず。食べてもらえない米を作るなどと言うのは耐えられない苦痛であるはず。科学者は,その苦痛を思いやれないのだろうなあと,やりきれない思いがする。

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

ダメ出し:アトキンの篩

2012年02月13日 | ブログラミング

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)
    }
}

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

PVアクセスランキング にほんブログ村

PVアクセスランキング にほんブログ村