Problem 88

2つ以上の自然数の集合の積かつ和として表すことのできる数を積和数と呼ぶ。
集合の大きさkに対し、最小の積和数Nを最小積和数と呼ぶ。
2≦k≦12000について最小積和数の集合の和(集合なので重複は除く)を求めよ、という問題。
とりあえず答えは出たもののたっぷり30分はかかった。やはりRは再帰を使うと時間がかかる。かといってほかの方法も思いつかない。フォーラム覗けるようになったのでまた考えようと思う。

## create prime number list
limit <- 1e5
sieve <- logical(limit)                 # prime number = FALSE
sieve[1] <- TRUE
for(i in 2:sqrt(limit)){
  if(!sieve[i]){
    for(j in seq(i * 2, limit, by = i)){
      sieve[j] <- TRUE
    }
  }
}
PNlist <- (1:limit)[!sieve]
head(PNlist) 

# prime factorization
Fact <- function(n){
  if(!identical(intersect(PNlist, n), numeric(0))) return(n)
  ans <- numeric(0)
  i <- 2
  while(n > 1){
    while(n %% i == 0){
      n <- n / i
      ans <- c(ans, i)
    }
    i <- i + 1
  }
  ans
}

# all factor
fact <- function(n){
  if(!identical(intersect(PNlist, n), numeric(0))) return(n)
  (2:(n/2))[n %% (2:(n/2)) == 0]
}

mydiff   <- function(x, y){
  tmp    <- intersect(names(table(x)), names(table(y)))
  result <- c(table(x)[tmp] - table(y)[tmp],
              table(x)[setdiff(names(table(x)), tmp)])
  result <- rep(names(result), result)
  ifelse(is.numeric(c(x,y)),
         return(as.numeric(result)),
         return(result))
}

# generate factor list
factlist <- function(f, fl=paste(f,collapse=" ")){
  if(length(f) <= 2){
    return(fl)
  }
  comb <- t(unique(t(combn(sort(f), 2))))
  ansl <- numeric(0)
  for(i in 1:(length(comb[1,]))){
    newf <- sort(c(mydiff(f, comb[,i]),prod(comb[,i])))
    fc <- paste(newf, collapse=" ")
    if(is.element(fc,fl)) next
    fl <- factlist(newf, c(fl, fc))
  }
  return(fl)
}
n2flist <- function(n) factlist(Fact(n))
# compute k from factor list
factlist2k <- function(f){
  f <- as.numeric(unlist(strsplit(f, " ")))
  prod(f) - sum(f) + length(f)
}
# n2k
n2k <- function(n) mapply(factlist2k, n2flist(n))

library(gmp)
limit <- 12000
candidates <- numeric(limit)
for(i in 4:(limit*2)){
  if(isprime(i)) next
  tmp <- n2k(as.numeric(i))
  tmp <- tmp[tmp <= limit]
  candidates[tmp][candidates[tmp]==0] <- i
  if(length(candidates[candidates==0])==1) break
}
sum(unique(candidates[2:limit]))

考え方

  1. nの約数の組み合わせf = (a1, a2, a3,...)に対し、kはn - sum(f) + length(f)で計算できる。
  2. よって、nを約数に分解できれば簡単にkを計算できる。
  3. 小さいnから対応するkを計算していって、kとnの対応リストに早いもの順でnを埋めていく。
  4. リストが全て埋まったらストップ。重複を除いて和を計算。

nに対する約数への分解パターンの列挙に最大で2秒近くかかっているので、ここをどうにかできればもう少し短縮できるかもしれない。それでも再帰を使う以上多少時間かかるのは避けられないし、それを12000も繰り返したら数分はどうしてもかかるだろう。
「12000程度の数までなら約数の数はそれほど多くならない」という点(実際、2^14で約数は14個しかないが、そこから計算されるkは16370になる)に注目し、うまいことアレした例がフォーラムの一つ目に投稿されてるけど、再帰や繰り返しが苦手な言語だとああいう方法が正解なのかもしれない。