Problem 152

http://projecteuler.net/problem=152

平方数の逆数の和として1/2を表す方法はいくつかある。
例えば、35以下の数を使うとすると、{2, 3, 4, 5, 7, 12, 15, 20, 28, 35}を用いて、

2から45までの数を使うとすると、他に2つ方法がある。
2から80までの数を使うとして、何通りの方法があるか。

分母の因子から素数を消すには、同じ因子を分母に持つ分数を足さないといけない。例えば13^2なら13*1から13*6まで(80を超えない範囲)の数の平方を分母に持つ分数をいくつか選んで足す。足した結果、分母の因子に13^2が現れれば13^2を分母から消せる。
消したい因子の倍数が使えないと駄目なので、40以上の素数は消せない。よって候補から除外できる。また、少し計算すると17以上の素数は全て消す方法が無いということが分かる(ついでに11も消せない)。
また、5を消したい時に(5*7)^2を足してしまうと、因子に7が増えてしまう。大きい因子から見ていくと約束した上で、係数から今見ている素数より大きな素数とその倍数を除外することができる。
というところまではなんとかなったけどその先ゴニョゴニョするだけの能力が無いので諦めた。brute forceに持ち込むだけでいっぱいいっぱい。難しい。その上forumを見てもよく分からない。
とりあえず答えは出たけど30分くらいかかってしまう。

library(gmp)
nums <-
  list(
       n13 = 13 * rev(c(1, 2, 3, 4, 5, 6)),
       n11 = 11 * rev(c(1, 2, 3, 4, 5, 6, 7)),
       n7  = 7 * rev(c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)),
       n5  = 5 * rev(c(1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16)),
       n3  = 3 * rev(c(1, 2, 3, 4, 6, 8, 9, 12, 16, 18, 24)),
       n2  = 2 * rev(c(1, 2, 4, 8, 16, 32, 64))
       )
reduce <- function(now, num, fac, flag = FALSE){
  if(now == 1/2) return(now)
  if(flag && now != 0 &&
     max(as.numeric(factorize(denominator(now)))) < fac){
    return(c(now,
             Recall(now, num, fac, flag = FALSE)) )
  }
  if(length(num) == 0) return(NULL)
  if(now > 1/2) return(NULL)
  return(c(as.bigq(Recall(now + as.bigq(1, num[1]^2), num[-1], fac,
                          flag = TRUE)),
           as.bigq(Recall(now, num[-1], fac, flag = flag))))
}
fracs <- as.bigq(0)
ps <- c(13, 11, 7, 5, 3, 2)
for(i in 1:6){
  fracs.new <- as.bigq(NULL)
  for(j in 1:length(fracs)){
    fracs.new <- c(fracs.new, reduce(fracs[j], nums[[i]], ps[i]))
  }
  fracs <- c(fracs,fracs.new)
}