万万使不得!-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的数整除,那么他就是质数,很明显我是按照定义来的

1
2
3
4
5
6
7
8
9
10
11
12
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,到他的平方根即可

1
2
3
4
5
6
7
8
9
10
11
12
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
}
}

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

1
2
3
4
5
6
7
8
9
10
11
12
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],先看看这一群中能除尽的最大数

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

输出结果是

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

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

1
2
3
4
5
6
7
8
9
10
11
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

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

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

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

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

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

1
2
3
4
5
6
7
8
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)
}

测试一下效果

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

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

1
2
3
4
5
6
7
8
9
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,而且速度很快,几乎是秒出,太惊讶了!!!
修改一下代码,得到所有的质数

1
2
3
4
5
6
7
8
9
10
11
12
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

------ 本文结束------