作者taurusrob (taurusrob)
看板R_Language
標題Re: [問題] 質因數分解
時間Thu Nov 5 02:11:37 2015
供你參考：
factorizePrime <- function(num){
isPrime <-function(num){
return(sum(num %% 1:floor(sqrt(num)) == 0) <= 1)
}
FindMidTwo <- function(num){
left <- max(which(num %% 1:floor(sqrt(num)) == 0))
return(c(left, num/left))
}
primeGroup <- FindMidTwo(num)
repeat{
if(all(sapply(primeGroup, isPrime))) break
keep <- sapply(primeGroup, isPrime)
primeGroup <- c(primeGroup[keep], as.vector(sapply(primeGroup[!keep],
FindMidTwo)))
}
return(table(primeGroup))
}
※ 引述《celestialgod (天)》之銘言：
: 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), 來自: 1.163.29.247
※ 文章網址: https://webptt.com/m.aspx?n=bbs/R_Language/M.1446660699.A.794.html
1F：→ haolihy: 謝謝! 11/05 13:32