Problem 141

Problem 141 - Project Euler

正整数nをdで割ったときの商をq、余りをrとおく。このとき、d, q, rが幾何級数の正整数の連続した項となるような場合がある(d, q, rの順序は入れ替えても良い)。
例えば、58を6で割ると9余り4となるが、4, 6, 9は公比3/2の幾何級数の連続した項である。このような数nをprogressive数と呼ぼう。
progressive数のうちいくつか、例えば9や10404(102^2)などは平方数でもある。10^5未満の平方数であるprogressive数の総和は124657である。
10^12未満の平方数かつprogressive数であるような数の総和を求めよ。

まず、n, d, q, rは次の関係にある。

ここで

だがqの大小関係は定まらない。qが(d, q, r)の中で1番小さいとすると、公比をaとして

となるため、nは平方数ではない。よってqは(d, q, r)のなかで1番大きいか2番目に大きいかだが、何れの場合も

という形になる。aをという形の既約分数に書き直すと、が整数にならなければならないことがわかる。このとき、rはcを因子として2つ以上持たなければならないので、

結局

であり、



という条件を満たす範囲でnが平方数となるようなものを探せばいいことになる。

library(gmp)
b <- 1
limit <- 10^12
ans <- 0
while(b^3 < limit){
  for(c in 1:b){
    e <- 1
    b3 <- b^3
    cb3 <- c*b3
    c2 <- c^2
    if(cb3 + c2 > limit) break
    if(max(gcd(c, b)) != 1) next
    ## while((tmp <- e^2*cb3 + c2*e) < limit){
    ##   if(round(sqrt(tmp))^2 == tmp) ans <- ans + tmp
    ##   e <- e + 1
    ## }
    e <- 1:(sqrt(limit/(cb3 + c2)))
    tmp <- e^2*cb3 + c2*e
    tmp <- tmp[tmp < limit]
    ans <- ans + sum(tmp[round(sqrt(tmp))^2 == tmp])
  }
  b <- b + 1
}
as.character(ans)

が、とても遅かった…
1分は切ったけど59.3秒なので切ってないようなものでなんともかんとも。Forum見ても同じようなやり方がほとんどだった。