作者showfeb ()
看板R_Language
标题Re: [问题] 质因数分解
时间Thu Nov 5 09:17:35 2015
利用recursive function求解
fn=function(k){
if(k%%1!=0|k<=1) return(k)
i=2
repeat{
if(k%%i==0) {
return(c(i,fn(k%/%i)))
}
if(i==k) return(k)
i=i+1
}
}
fn(k=12345678)
system.time(fn(k=2^11*3^6*47))
※ 引述《taurusrob (taurusrob)》之铭言：
: 供你参考：
: 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秒
--
※ 发信站: 批踢踢实业坊(ptt.cc), 来自: 140.128.153.200
※ 文章网址: https://webptt.com/cn.aspx?n=bbs/R_Language/M.1446686257.A.B72.html
1F：→ celestialgod: 从来没写过递回(蒙脸 11/05 09:47
2F：→ celestialgod: 漂亮的解法!! 11/05 09:47
3F：推 haolihy: 谢谢! 11/05 13:32
4F：推 Edster: 算阶乘, 碎形的时候就会用到递回了. 没跟到这波, 高手真多 11/05 14:10