和果子一起来做题-Project Euler-03-R语言版本

这是第3题

The prime factors of 13195 are 5, 7, 13 and 29.
What is the largest prime factor of the number 600851475143 ?

就是找出600851475143中最大的质数。
我的思路是，找出小于600851475143的质数，得到质数表，然后用600851475143除以这些质数再判断
首先找出小于600851475143的质数
思路是，如果一个数不能被其他除了自己和1的数整除，那么他就是质数，很明显我是按照定义来的

k <- 2
x <-c(2,3)
for (j in 5:600851475143){
  a <-c()
  for (i in 2:(j-1)){
    a[i-1] = sum( j%%i == 0)
  }
  if (sum(a)==0){
    k = k + 1
    x[k] = j
  }
}

运行后提示：

Error: cannot allocate vector of size 4476.7 Gb

我想了一下，不需要数到600851475143，到他的平方根即可

k <- 2
x <-c(2,3)
for (j in 5:sqrt(600851475143)){
  a <-c()
  for (i in 2:(j-1)){
    a[i-1] = sum( j%%i == 0)
  }
  if (sum(a)==0){
    k = k + 1
    x[k] = j
  }
}

这一次应该可以，但是运行时间太久，我等了半天都不行，我先测试一下代码是否正确

k <- 2
x <-c(2,3)
for (j in 5:10000){
  a <-c()
  for (i in 2:(j-1)){
    a[i-1] = sum( j%%i == 0)
  }
  if (sum(a)==0){
    k = k + 1
    x[k] = j
  }
}

运行大概需要20s，k是1229,说明找到了1229个质数，最大的是max(x),9973
那么用600851475143除以x[k],先看看这一群中能除尽的最大数

y <- 600851475143
x[which(y %% x==0)]

输出结果是

x[which(y %% x==0)]
[1] 71 839 1471 6857

最大是6857，恰巧这四个数的乘积就是600851475143，说明蒙对了。
我在想有没有什么通用的方法呢？在网上检索了一下质数的查找方法
找质数算法之埃拉托色尼筛选法-Sieve of Eratosthenes算法
在R语言中使用这个方法获取质数

x <- 600851475143
num <- 2:sqrt(x)
i <- 0
wt <- c(2) #用于储存质数
repeat{
  i <- i +1
  num <- num[which(num %% num[1] !=0)] #每次把第一个数以及他的倍数去除
  wt[i+1] <- num[1] #把第一个数传入向量中
  if (length(num) == 0) break #repeat需要设定终止条件
}
wt <- wt[-length(wt)] #去掉最后一个NA

最终发现这个方法是可行的

y <- 600851475143
wt[which(y %% wt==0)]

wt[which(y %% wt==0)]
[1] 71 839 1471 6857

这时候确定无疑就是6857，并且很高兴解决了这么多问题，然后就去参考一下别人博客的答案

这个先写了一个能判断质数的函数，其中all用的很传神，我在之前不知道有这个函数，就用了迂回的方法代替

findprime <- function(x) {
  if (x %in% c(2,3,5,7)) return(TRUE)
  if (x%%2 == 0 | x==1) return(FALSE)
  xsqrt <- round(sqrt(x))
  xseq <- seq(from=3,to=xsqrt,by=2)
  if (all(x %% xseq !=0)) return(TRUE)
  else return(FALSE)
}

测试一下效果

x = 1:700000
x[sapply(x,findprime)]

发现速度特别快,下面找出最大的质数

n <- 600851475143
for (i in seq(from=3, to=round(sqrt(n)), by=2)) {
  if (findprime(i) & n %% i == 0) {
    n <- n / i
    prime.factor <- i       
    if (i >= n) break
  }
}
print(prime.factor)

得到的答案是6857，而且速度很快，几乎是秒出，太惊讶了！！！
修改一下代码，得到所有的质数

n <- 600851475143
j <- 0
factor <- c()
for (i in seq(from=3, to=round(sqrt(n)), by=2)) {
  if (findprime(i) & n %% i == 0) {
    n <- n / i
    j <- j + 1
    factor[j] <- i       
    if (i >= n) break
  }
}
print(factor)

print(factor)
[1] 71 839 1471 6857

总结：
1.使用字母时为什么常用i和j呢？因为他们大小写不容易混淆，我今天使用的是k，大小写分不清，会出现object "k" not found 的提示
2.要学会写函数
3.尽量不要循环里面套循环，速度会很慢
4.all()和any()的用法
5.学习Adcanced R里面的函数章节
函数Functions
泛函数Functionals