Problem 132 - もうカツ丼はいいよな

Problem 132 - Project Euler

1のみからなる数(1, 11, 111, ...)をrepunitと呼ぶ。k桁のrepunitをR(k)で表すことにしよう。
さて、R(10) = 1111111111を素因数分解すると素因数は11, 41, 271, 9091であるから、素因数の総和は9414である。
R( $10^9$ )の素因数のうち小さいものから40個の総和を求めよ。

R(A(n))はnで割り切れるが、剰余が循環するという性質から、mを1以上の整数としてR(m * A(n))もnで割り切れる(c.f. Problem 129 - もうカツ丼でいいよな)。
つまり、10^9がA(n)で割り切れるかどうかを調べればR(10^9)がnで割り切れるかどうかを調べることができる。ただ、このやりかたで素因数を網羅的に調べようとすると時間がかかりすぎるので別のアプローチをする。
R(k)を9倍して1を足すと10^kになる。よって
$10^k \bmod n = 1$
であるならば、nはR(k)の9倍を割り切れる。
冪剰余は高速に求められる(c.f.冪剰余 - Wikipedia)ので、
$10^9 \bmod n = 1$
を満たしかつ素数であるようなnを小さいものから40個探せばよい。
ただ、R(k)の9倍はkの値に関わらず3で割りきれてしまうので、3は最初のやり方で別にチェックする。9も割り切れるが、9は素数でないのでそもそも対象外。

library(gmp)

lim <- 1e6
sieve <- logical(lim)
sieve[c(1, seq(4, lim, by=2))] <- TRUE
for(i in 3:(sqrt(lim))){
  if(!sieve[i]) sieve[seq(i*2, lim, by=i)] <- TRUE
}

prm.list <- (1:lim)[!sieve]

A <- function(n){
  if(n %% 2 == 0 || n %% 5 == 0) return(0)
  k <- 1
  m <- 1
  while(m != 0){
    k <- k + 1
    m <- (m * 10 + 1) %% n
  }
  k
}
lim <- 1e9
f <- lim %% A(3)

cnt <- 0
ans <- 0
if(f){
  i <- 3
}else{
  i <- 2
}
while(cnt < 40){
  if(powm(10, lim, prm.list[i]) == 1){
    ans <- ans + prm.list[i]
    cnt <- cnt + 1
  }
  i <- i + 1
}
ans

610ms。