RでProject Euler 11〜15

R

Problem 11

文字列のベクトルを作っておいて、各文字列をsplitしたかったんですが、うまくいかず。しかたなくファイルにしてそれを読み込みました。それだと簡単に行列になるんですよね。

そうすれば、あとは始点と方向を列挙するだけ。これはデータフレームにして、あとはapplyで1レコードずつ処理します。データフレームの使い方もやっと少しわかりました。
あと、applyすると

pos_vec <- function(M) {
    xs <- c()
    ys <- c()
    vxs <- c()
    vys <- c()
    H <- dim(M)[1]
    W <- dim(M)[2]
    vs <- matrix(c(c(1, 0), c(1, 1), c(0, 1), c(-1, 1)), nrow=4, ncol=2)
    print(vs)
    for(x0 in 1:H) {
        for(y0 in 1:W) {
            for(k in 1:4) {
                v <- vs[k,]
                x1 <- x0 + v[1] * 3
                y1 <- y0 + v[2] * 3
                if(1 <= x1 && x1 <= W && y1 <= H) {
                    xs <- c(xs, x0)
                    ys <- c(ys, y0)
                    vxs <- c(vxs, v[1])
                    vys <- c(vys, v[2])
                }
            }
        }
    }
    
    return(data.frame(x0 = xs, y0 = ys, vx = vxs, vy = vys))
}

mul <- function(pt) {
    x0 <- pt[1]
    y0 <- pt[2]
    vx <- pt[3]
    vy <- pt[4]
    return(M[x0+vx*0,y0+vy*0] * M[x0+vx*1,y0+vy*1] *
           M[x0+vx*2,y0+vy*2] * M[x0+vx*3,y0+vy*3])
}

M <- as.matrix(read.table("e011_table.txt", sep=" "))
starts <- pos_vec(M)
v <- apply(starts, 1, mul)
print(max(unlist(v)))

あとvはリストなんですが、リストのままだとmaxとかが適用できないので、unlistでリストでなくする必要があります（よくわかっていない）。

Problem 12

普通に試し割りで素因数分解しています。素因数分解した結果はデータフレームにしています。
pとeの列があるので、そこからeの列だけ選んで、約数の個数を数えています。

factorize <- function(n) {
    ps <- c()
    es <- c()
    p = 2
    while(n >= p * p) {
        if(n %% p == 0) {
            e <- 1
            n <- n / p
            while(n %% p == 0) {
                e <- e + 1
                n <- n / p
            }
            ps <- c(ps, p)
            es <- c(es, e)
        }
        p <- p + 1
    }
    if(n > 1) {
        ps <- c(ps, n)
        es <- c(es, 1)
    }
    return(data.frame(p = ps, e = es))
}

num_divisors <- function(n) {
    f <- factorize(n)
    m <- 1
    for(e in f[,colnames(f) == "e"]) {
        m <- m * (e + 1)
    }
    return(m)
}

N <- 500
n <- 1
while(TRUE) {
    m <- n * (n + 1) / 2
    nd <- num_divisors(m)
    if(nd > N)
        break
    n <- n + 1
}
print(m)

Problem 13

ふつうに1桁をベクトルの要素として和を関数にします。

s <- c("37107287533902102798797998220837590246510135740250",
       "46376937677490009712648124896970078050417018260538",
        ...
)

long <- function(s) {
    return(sapply(nchar(s):1, function(k) as.integer(substring(s, k, k))))
}

add <- function(x, y) {
    if(length(x) < length(y)) {
        return(add(y, x))
    }
    
    z <- c()
    c <- 0
    for(k in 1:length(y)) {
        a <- x[k] + y[k] + c
        z <- c(z, a %% 10)
        c <- a %/% 10
    }
    k <- k + 1
    while(k <= length(x)) {
        a <- x[k] + c
        z <- c(z, a %% 10)
        c <- a %/% 10
        k <- k + 1
    }
    if(c > 0) {
        z <- c(z, c)
    }
    return(z)
}

z <- long(s[1])
for(k in 2:length(s)) {
    z <- add(z, long(s[k]))
}
print(as.numeric(paste(z[length(z):(length(z)-9)], collapse="")))

Problem 14

ただのコラッツなので、100万までのベクトルを用意して、メモ化します。
なのですが、ここで注意点が。グローバルな変数に代入するには<-ではなく<<-を使わなければならないそうです。

step <- function(n) {
    if(n %% 2 == 0) {
        return(n %/% 2)
    }
    else {
        return(n * 3 + 1)
    }
}

num_steps <- function(n) {
    if(!is.na(a[n])) {
        return(a[n])
    }
    else {
        m <- 1 + num_steps(step(n))
        if(n <= N) {
            a[n] <<- m
        }
        return(m)
    }
}

N = 1000000
a <- rep(NA, N)
a[1] <- 1
for(n in 2:N) {
    num_steps(n)
}
max_n <- 0
max_len <- 0
for(n in 1:N) {
    if(a[n] > max_len) {
        max_len <- a[n]
        max_n <- n
    }
}
print(max_n)

Problem 15

Cの計算は再帰で計算できるから簡単ですよね。本当はreduceが使えるといいんですが。

C <- function(n, m) {
    if(m == 0) {
        return(1)
    }
    else {
        return(C(n, m - 1) * (n - m + 1) / m)
    }
}

N = 20
print(C(N * 2, N))