primes.R

# sieve
sieve <- function(N) {
  primes <- seq(1,N)
  sprimes <- numeric((N/2)-2)
  
  for(i in 2:(N/2)) {
    if(primes[i] != 0) { # prime
      primes[seq(2*i,N,i)] <- 0
      sprimes[i] <- sum(primes)
    }
  }
  mprimes <- tail(primes[which(primes!=0)],-2)
  dprimes <- diff(sprimes)
  sprimes <- sprimes[which(sprimes!=0)]
  nsieve <- length(unique(sprimes))
  nsieve
}


# par(mfrow=c(3,3))
# N <- seq(1000,18000,1000)
# for(n in N) {
#   D <- density(sieve(n))
#   plot(D, main=paste(D$n/n))
# }

N <- seq(1000, 50000, 1000) 
n <- numeric(50)

for(i in 1:50) {
  n[i] <- sieve(N[i])
}
par(mfrow=c(2,2))
plot(N, n, type='l')
# intercept 52, slope = 1.876e-04
# abline(model)