Problem 127

R Project Euler

正整数nの素因数の積をrad(n)とする．また，二つの正整数x, yの最大公約数(greatest common divisor)をGCD(x, y)とする．
ここで，次の条件を満たす3つの正整数(a, b, c)をabc-hitであると呼ぶ．

  1. GCD(a, b) = GCD(b, c) = GCD(a, c) = 1
  2. a < b
  3. a + b = c
  4. rad(abc) < c

c < 120,000について，abc-hitであるときのcの和Σcを求めよ，という問題．
15分くらい．話にならないですね．本当に色々と話にならない．というかもう何週間か考えてて完全に諦めてたので，アルゴリズムはwonderfl build flash online | 面白法人カヤックを参考にしました．それで15分ってことはもうRが遅い(or僕の実装がクソ)ということなのでもう何だかアレがアレでアレじゃないですか．
ただ副作用としてProblem 124は2秒程度に短縮できるようになった．エラトステネスの篩に少し手を加えるだけでradのリストが求められるので．
一応解いた問題の数だけならRユーザの中でトップタイになった(Project Euler Language: R)．ただProject Eulerのスコアには解いた問題数の他にCurrent Performanceというものがあって，これは最新の問題から25問のうちどれだけを解いたかに基づいてスコアが付く．その関係でまだ1位ではない．というか新しい方の問題は僕の力では歯が立たないものばかり．というかすでにかなり歯が立たなくなってきている．ボチボチ番号順に解くのは止めて正解者の多い順にやってこうか．

library(gmp)

limit <- 120000
## gen. prime list
sieve <- logical(limit)
sievelim <- sqrt(limit)
sieve[1] <- TRUE
sieve[seq(4, limit, by=2)] <- TRUE
for(i in seq(3, sievelim, by=2)) if(!sieve[i]) sieve[seq(i*2, limit, by=i)] <- TRUE
prime <- (1:limit)[!sieve]
## gen. rad list
rad <- rep(1, limit)
for(p in prime){
  rad[seq(p, limit, by=p)] <- rad[seq(p, limit, by=p)] * p
}
count <- 0
ans <- 0
for(c in 116874:limit){
  cr <- c/rad[c]
  if(cr < 3) next
  fl <- unique(as.numeric(factorize(c)))
  ch <- logical(c)
  ch[(1:c) > (ceiling(c/2)-1)] <- TRUE
  for(f in fl) ch[seq(f, c, by=f)] <- TRUE
  a <- (1:c)[!ch]
  a <- a[-1]
  ch <- logical(length(a))
  ch[cr <= rad[a]*rad[c-a]] <- TRUE
  a <- a[!ch]
  if(cr > rad[c-1]){
    ans <- ans + c
    count <- count + 1
  }
  for(a in a){
    if(sum(unique(gcd(a, c-a))) == 1 &&
       sum(unique(gcd(c, c-a))) == 1){
      ans <- ans + c
      count <- count + 1
    }
  }
}
ans