PS: 希望你这个不是作业，作业要自己做XD
primefactor <- function(num){
p <- which(num %% 1:num == 0)
isPrime <- function(num){
return(sum(num %% 1:num == 0) == 2)
}
p <- p[sapply(p, isPrime)]
e <- vector('numeric', length(p))
for (i in 1:length(p))
{
repeat{
if (num %% p[i]^e[i] != 0)
break
e[i] = e[i] + 1
}
}
return(list('pfactor'=p, 'multi'=e - 1))
}
st = proc.time()
primefactor(12345678)
# $pfactor
# [1] 2 3 47 14593
#
# $multi
# [1] 1 2 1 1
proc.time() - st
# user system elapsed
# 20.53 5.54 26.29
env: i5-750@2.67GHz machine with RRO-3.2.2 in windows 7 SP1
在p <- which(num %% 1:num == 0)跟isPrime这两步时要小心内存会被用完
如果担心的话可以拆分成repeat去做
例如：
num = 12345678
p = vector('numeric', 1000)
k = 1
step_p = 100000
repeat{
tmp = which(num %% k:(k+step_p) == 0)
if (length(tmp) > 0)
p[(sum(p!=0)+1):(sum(p!=0)+length(tmp))] = tmp + k - 1
if (k > num) break
k = k + step_p
}
p = p[p!=0]
# 1 2 3 6 9 18 47
# 94 141 282 423 846 14593 29186
# 43779 87558 131337 262674 685871 1371742 2057613
# 4115226 6172839 12345678
给一个更快的方法(直接用向量化去做)：
primefactorFast <- function(num){
p <- which(num %% 1:num == 0)
isPrime <- function(num){
return(sum(num %% 1:num == 0) == 2)
}
p <- p[sapply(p, isPrime)]
e <- vector('numeric', length(p))
e_basic <- 1
step_e <- 3
subset_e <- 1:length(p)
repeat{
m = num %% sweep(matrix(replicate(step_e, p[subset_e]),
length(p[subset_e])), 2, e_basic:(e_basic + step_e - 1), '^') == 0
e[subset_e] = e[subset_e] + rowSums(m)
e_basic = e_basic + step_e
subset_e = which(apply(m,1,all))
if (is.null(subset_e) || length(subset_e) == 0)
break
}
return(list('pfactor'=p, 'multi'=e))
}
st = proc.time()
primefactorFast(12345678)
# $pfactor
# [1] 2 3 47 14593
#
# $multi
# [1] 1 2 1 1
proc.time() - st
# user system elapsed
# 0.06 0.00 0.06
PS: 更动step_e可以得到更好的效能
final version: http://pastebin.com/zkD3N2kp
test 2^11*3^6*47 (10^7.8) 只需要6.5秒
※ 引述《haolihy (好厉害)》之铭言：
: [问题类型]:
: 效能咨询(我想让R 跑更快)
: [软件熟悉度]:
: 入门(写过其他程式，只是对语法不熟悉)
: [问题叙述]:
: 质因子分解
: [程式范例]:
: 想写一个程式对正整数n作质因子分解
: 目标是n<10^8 可以1min跑完
: 我的程式大概要6~7min 改进空间应该很大
: 而太慢的主因应该是列出<n质数表的程式太慢
: 请大家赐教
: 质数表的程式:
: primeY <- function(num){
: sqnum <- floor(sqrt(num))
: x <- c(2:num)
: y <- c()
: repeat{
: if(x[1] > sqnum) break
: y <- c(y, x[1])
: x <- x[x %% x[1] != 0]
: }
: return(c(y,x))
: }
: 质因子分解的程式:
: primefactor <- function(num){
: p <- primeY(num)[which(num %% primeY(num) == 0)]
: q <- c()
: e <- c()
: for(i in 1:length(p)){
: e[i] <- 1
: q[i] <- num / p[i]
: repeat{
: if(q[i] %% p[i] != 0) break
: e[i] <- e[i] + 1
: q[i] <- q[i] / p[i]
: }
: }
: return(list('pfactor'=p,
: 'multi'=e))
: }
: [环境叙述]:
: R 3.2.2