星期一, 八月 20, 2012

笨办法学R编程(3)

经历了前面两个小挑战,你应该对R有点理解了。我们继续推进,今天的问题有点点复杂,复杂的不是R,而是一个数学概念:质数和质因子。任何一个合数都可以被几个质数所分解,这个性质很重要,我们将用它来解决Project Euler的第三个问题。还是和之前一样的,你需要自己在R控制台中敲打下面这些命令,根据结果自行揣摩其用处。
# 预备练习,学习for循环、建立自定义函数和其它一些函数
 
for (n in 1:10) {
    print(sqrt(n))
}
 
x <- c('hello','world','I','love','R')
for (n in x) {
    print(n)
}

x <- seq(from=1,to=10,by=1)
print(x)
x <- seq(from=1,to=10,by=2)
print(x)
x <- seq(from=1,to=2,length.out=10)
print(x)
round(x) 
x > 1.5
all(x>1.5)
any(x>1.5)
 
# 如何自定义一个求圆面积的函数
myfunc <- function(r) {
    area <- pi*r^2
    return(area)
}
print(myfunc(4))
 
# 同时求四个不同半径圆的面积
r <- c(2,2,4,3)
sapply(X=r,FUN=myfunc)
 
# Project Euler 3
# 找到600851475143这个数的最大质因子
# 先建立一个函数以判断某个数是否为质数
 
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到100的质数,看函数对不对
x = 1:100
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。本例中除了使用for循环外,还见到了sapply函数,这是R语言中非常重要的一类向量化计算函数。求质数的方法可以参考这个文章,本例使用的是其中的境界4。实际上根据质因子的性质,本例不一定非要建立判断质数的函数,不过这个函数我们在后面会用到的。另外如果你想用其它软件找这个数字的质因子,也可以看看这里

10 条评论:

  1. 贴一个我的代码(http://ruready.sinaapp.com)

    # 构建一个函数,判断任何整型数字是否质数
    if.prime = function(n){
    if (n <= 1) return (FALSE)
    if (n == 2) return (TRUE)
    for (i in 2:sqrt(n))
    if (n %% i == 0) return (FALSE)
    else return (TRUE)
    }

    # 用for循环对x进行试除,输出最大的质因数
    x = 600851475143
    for (i in 2:sqrt(x)){
    if (if.prime(i) & x %% i == 0){
    x = x/i
    max.prime = i
    if (i >= x) break
    }
    }
    print(max.prime)

    回复删除
  2. 多谢,收益良多!!

    回复删除
  3. 博主,x %in% c(2,3,5,7)那个符号%in%我怎么查不到什么意思。

    回复删除
    回复
    1. 那表示x 是否为2,3,5,7中的任意一个。

      删除
  4. 和除法没有任何关系吧?其余的带%%的都和整除有点关系的,这个属于例外,能这么理解吗?

    回复删除
  5. 此评论已被作者删除。

    回复删除
  6. 不用for循环就可以的
    nseq<-seq(1,sqrt(600851475143))
    x<-nseq[600851475143%%xseq==0]
    y<-600851475143/x
    z<-c(x,y)
    max(z[sapply(z,findprime)])

    回复删除