Re: [問題] 質因數分解
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
--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 140.109.73.190
※ 文章網址: https://www.ptt.cc/bbs/R_Language/M.1446631193.A.0DC.html
※ 編輯: celestialgod (140.109.73.190), 11/04/2015 18:43:36
※ 編輯: celestialgod (180.218.154.163), 11/04/2015 19:51:09
推
11/05 09:43, , 1F
11/05 09:43, 1F
討論串 (同標題文章)