素数の日付を含む最長期間

素数の日付を含む最長期間

締め切りが 2017/05/30 10:00 AM なので，その 1 分後に投稿されるように予約

設問

日付をYYYYMMDD形式で表現し、8桁の数値としてみたとき、その値が素数かどうかを判定します。
1970年1月1日～2019年12月31日までの50年間のうち、素数がちょうど n 個含まれる期間で最長のものを考え、その日数を求めます。
なお、日数は両端の日付を含んで数えるものとします。
また、閏年は考慮するものとします。

例えば、n = 3 のとき、以下の青色で塗りつぶした範囲が最長になり、その日数は「202」です。

2015-04-11 ： 素数
2015-04-12 ： 開始日
2015-05-13 ： 素数
2015-08-21 ： 素数
2015-10-11 ： 素数
2015-10-30 ： 終了日
2015-10-31 ： 素数

標準入力から n が与えられたとき、その最長日数を求め、標準出力に出力してください。
なお、n は1000以下の自然数とします。

【入出力サンプル】
標準入力
3

標準出力
202

-------------------------------------

R で書くと実行時間が1秒以内に収まらない

is.prime = function(n) {
    if (n %% 2 == 0) return(FALSE)
    else if(n %% 3 == 0) return(FALSE)
    maxitr = as.integer(sqrt(n))
    i = 1
    repeat {
        i = i+4
        if (i > maxitr) return(TRUE)
        else if (n %% i == 0) return(n == i)
        i = i+2
        if (i > maxitr) return(TRUE)
        else if (n %% i == 0) return(n == i)
    }
}

J.day = function(iy, jm, kd) {
    tmp = -(jm < 3)
    kd - 32075 + (1461 * (iy + 4800 + tmp))%/%4 + (367 * (jm - 2 - tmp * 12))%/%12 - (3 * ((iy + 4900 + tmp)%/%100))%/%4
}

date2 = function(jul) {
    l = jul + 68569
    n = (4 * l)%/%146097
    l = l - (146097 * n + 3)%/%4
    iy = (4000 * (l + 1))%/%1461001
    l = l - (1461 * iy)%/%4 + 31
    jm = (80 * l)%/%2447
    kd = l - (2447 * jm)%/%80
    l = jm%/%11
    jm = jm + 2 - 12 * l
    iy = 100 * (n - 49) + iy + l
    iy*10000+jm*100+kd
}

f = function(n) {
    begin = J.day(1970, 1, 1)
    end = J.day(2019, 12, 31)
    count = 0
    p.date = integer(1100)
    for (i in begin:end) {
        date = date2(i)
        if (is.prime(date)) {
            count = count+1
            p.date[count] = i
        }
    }
    Max = 0
    for (i in 1:(count-n-1)) {
        Max = max(Max, p.date[i+n+1]-p.date[i]-1)
    }
    cat(Max)
}
# f(scan(file("stdin", "r")))

# f(3) # 202
# f(10) # 413
# f(333) # 5771
# f(876) # 14708
# f(999) # 16724

awk で書いて少し速くなって，パスした

function isPrime(n,     maxitr, i) {
    if (n % 2 == 0) return 0
    else if(n % 3 == 0) return 0
    maxitr = int(sqrt(n))
    i = 1
    for (;;) {
        i = i+4
        if (i > maxitr) return 1
        else if (n % i == 0) return n == i
        i = i+2
        if (i > maxitr) return 1
        else if (n % i == 0) return n == i
    }
}

function Jday(iy, jm, kd,     tmp) {
    tmp = -(3 > jm)
    return kd - 32075 + int((1461 * (iy + 4800 + tmp))/4) + int((367 * (jm - 2 - tmp * 12))/12) - int((3 * (int((iy + 4900 + tmp)/100)))/4)
}

function revJday(jul,     l, n, iy, jm, kd) {
    l = jul + 68569
    n = int((4 * l)/146097)
    l = l - int((146097 * n + 3)/4)
    iy = int((4000 * (l + 1))/1461001)
    l = l - int((1461 * iy)/4) + 31
    jm = int((80 * l)/2447)
    kd = l - int((2447 * jm)/80)
    l = int(jm/11)
    jm = jm + 2 - 12 * l
    iy = 100 * (n - 49) + iy + l
    return iy*10000+jm*100+kd
}

function max(x, y) {
    return x >= y ? x : y
}

function f(n,     i, begin, end, count, primeDate, date, Max) {
    begin = Jday(1970, 1, 1)
    end = Jday(2019, 12, 31)
    count = 0
    for (i = begin; end >= i; i++) {
        date = revJday(i)
        if (isPrime(date)) {
        primeDate[++count] = i
        }
    }
    Max = 0
    for (i  = 1; count-n-1 >= i; i++) {
        Max = max(Max, primeDate[i+n+1]-primeDate[i]-1)
    }
    print Max
}
BEGIN {
f(3)
f(10)
f(333)
f(876)
f(999)
}