正整数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見ても同じようなやり方がほとんどだった。