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: [email protected] 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://webptt.com/m.aspx?n=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