factores - para que sirve factor en r

Función R para devolver TODOS los factores

Question

Mi búsqueda normal foo me está fallando. Estoy tratando de encontrar una función R que devuelva TODOS los factores de un entero. Hay al menos 2 paquetes con funciones factorize() : gmp y conf.design, sin embargo, estas funciones solo devuelven factores primos. Me gustaría una función que devuelva todos los factores.

Obviamente, la búsqueda de esto se hace difícil ya que R tiene un constructo llamado factores que hace mucho ruido en la búsqueda.

Answer 1

Actualización importante

A continuación se muestra mi último algoritmo de factorización R. Es mucho más rápido y rinde homenaje a la función rle .

Algoritmo 3 (Actualizado)

library(gmp)
MyFactors <- function(MyN) {
    myRle <- function (x1) {
        n1 <- length(x1)
        y1 <- x1[-1L] != x1[-n1]
        i <- c(which(y1), n1)
        list(lengths = diff(c(0L, i)), values = x1[i], uni = sum(y1)+1L)
    }

    if (MyN==1L) return(MyN)
    else {
        pfacs <- myRle(factorize(MyN))
        unip <- pfacs$values
        pv <- pfacs$lengths
        n <- pfacs$uni
        myf <- unip[1L]^(0L:pv[1L])
        if (n > 1L) {
            for (j in 2L:n) {
                myf <- c(myf, do.call(c,lapply(unip[j]^(1L:pv[j]), function(x) x*myf)))
            }
        }
    }
    myf[order(asNumeric(myf))]  ## ''order'' is faster than ''sort.list''
}

A continuación se muestran los nuevos puntos de referencia (como Dirk Eddelbuettel dice here , " No se puede discutir con empíricos "):

Caso 1 (factores primos grandes)

set.seed(100)
myList <- lapply(1:10^3, function(x) sample(10^6, 10^5))
benchmark(SortList=lapply(myList, function(x) sort.list(x)),
            OrderFun=lapply(myList, function(x) order(x)),
            replications=3,
            columns = c("test", "replications", "elapsed", "relative"))
      test replications elapsed relative
2 OrderFun            3   59.41    1.000
1 SortList            3   61.52    1.036

## The times are limited by "gmp::factorize" and since it relies on
## pseudo-random numbers, the times can vary (i.e. one pseudo random
## number may lead to a factorization faster than others). With this
## in mind, any differences less than a half of second
## (or so) should be viewed as the same. 
x <- pow.bigz(2,256)+1
system.time(z1 <- MyFactors(x))
user  system elapsed
14.94    0.00   14.94
system.time(z2 <- all_divisors(x))      ## system.time(factorize(x))
user  system elapsed                    ##  user  system elapsed
14.94    0.00   14.96                   ## 14.94    0.00   14.94 
all(z1==z2)
[1] TRUE

x <- as.bigz("12345678987654321321")
system.time(x1 <- MyFactors(x^2))
user  system elapsed 
20.66    0.02   20.71
system.time(x2 <- all_divisors(x^2))    ## system.time(factorize(x^2))
user  system elapsed                    ##  user  system elapsed
20.69    0.00   20.69                   ## 20.67    0.00   20.67
all(x1==x2)
[1] TRUE

Caso 2 (números más pequeños)

set.seed(199)
samp <- sample(10^9, 10^5)
benchmark(JosephDivs=sapply(samp, MyFactors),
            DontasDivs=sapply(samp, all_divisors),
            OldDontas=sapply(samp, Oldall_divisors),
            replications=10,
            columns = c("test", "replications", "elapsed", "relative"),
            order = "relative")
        test replications elapsed relative
1 JosephDivs           10  470.31    1.000
2 DontasDivs           10  567.10    1.206  ## with vapply(..., USE.NAMES = FALSE)
3  OldDontas           10  626.19    1.331  ## with sapply

Caso 3 (para una minuciosidad completa)

set.seed(97)
samp <- sample(10^6, 10^4)
benchmark(JosephDivs=sapply(samp, MyFactors),
            DontasDivs=sapply(samp, all_divisors),
            CottonDivs=sapply(samp, get_all_factors),
            ChaseDivs=sapply(samp, FUN),
            replications=5,
            columns = c("test", "replications", "elapsed", "relative"),
            order = "relative")
        test replications elapsed relative
1 JosephDivs            5   22.68    1.000
2 DontasDivs            5   27.66    1.220
3 CottonDivs            5  126.66    5.585
4  ChaseDivs            5  554.25   24.438

Mensaje original

El algoritmo del Sr. Cotton es una muy buena implementación de R El método de fuerza bruta solo te llevará hasta el momento y falla con grandes números (también es increíblemente lento). He proporcionado tres algoritmos que satisfarán diferentes necesidades. El primero (es el algoritmo original que publiqué el 15 de enero y se ha actualizado ligeramente), es un algoritmo de factorización independiente que ofrece un enfoque combinatorio eficiente, preciso y que puede traducirse fácilmente a otros idiomas. El segundo algoritmo es más bien un tamiz que es muy rápido y extremadamente útil cuando necesita la factorización de miles de números rápidamente. El tercero es un algoritmo autónomo corto (publicado arriba), pero potente, que es superior para cualquier número menor a 2 ^ 70 (eliminé casi todo de mi código original). Me plyr::count el uso de la función plyr::count de Richie Cotton (me inspiró a escribir mi propia función rle que tiene un retorno muy similar al de plyr::count ), la forma limpia de George Dontas de manejar el caso trivial (es decir, if (n==1) return(1) ), y la solución provista por Zelazny7 a una question que tuve con respecto a los vectores bigz.

Algoritmo 1 (original)

library(gmp)
factor2 <- function(MyN) {
    if (MyN == 1) return(1L)
    else {
        max_p_div <- factorize(MyN)
        prime_vec <- max_p_div <- max_p_div[sort.list(asNumeric(max_p_div))]
        my_factors <- powers <- as.bigz(vector())
        uni_p <- unique(prime_vec); maxp <- max(prime_vec)
        for (i in 1:length(uni_p)) {
            temp_size <- length(which(prime_vec == uni_p[i]))
            powers <- c(powers, pow.bigz(uni_p[i], 1:temp_size))
        }
        my_factors <- c(as.bigz(1L), my_factors, powers)
        temp_facs <- powers; r <- 2L
        temp_facs2 <- max_p_div2 <- as.bigz(vector())
        while (r <= length(uni_p)) {
            for (i in 1:length(temp_facs)) {
                a <- which(prime_vec >  max_p_div[i])
                temp <- mul.bigz(temp_facs[i], powers[a])
                temp_facs2 <- c(temp_facs2, temp)
                max_p_div2 <- c(max_p_div2, prime_vec[a])
            }
            my_sort <- sort.list(asNumeric(max_p_div2))
            temp_facs <- temp_facs2[my_sort]
            max_p_div <- max_p_div2[my_sort]
            my_factors <- c(my_factors, temp_facs)
            temp_facs2 <- max_p_div2 <- as.bigz(vector()); r <- r+1L
        }
    }
    my_factors[sort.list(asNumeric(my_factors))]
}

Algoritmo 2 (tamiz)

EfficientFactorList <- function(n) {
    MyFactsList <- lapply(1:n, function(x) 1)
    for (j in 2:n) {
        for (r in seq.int(j, n, j)) {MyFactsList[[r]] <- c(MyFactsList[[r]], j)}
    }; MyFactsList}

Da la factorización de cada número entre 1 y 100,000 en menos de 2 segundos. Para darle una idea de la eficiencia de este algoritmo, el tiempo para calcular el factor 1 - 100,000 usando el método de fuerza bruta toma cerca de 3 minutos.

system.time(t1 <- EfficientFactorList(10^5))
user  system elapsed 
1.04    0.00    1.05 
system.time(t2 <- sapply(1:10^5, MyFactors))
user  system elapsed 
39.21    0.00   39.23 
system.time(t3 <- sapply(1:10^5, all_divisors))
user  system elapsed 
49.03    0.02   49.05

TheTest <- sapply(1:10^5, function(x) all(t2[[x]]==t3[[x]]) && all(asNumeric(t2[[x]])==t1[[x]]) && all(asNumeric(t3[[x]])==t1[[x]]))
all(TheTest)
[1] TRUE

Pensamientos finales

El comentario original del Sr. Dontas sobre factorizar grandes números me hizo pensar, ¿qué pasa con los números realmente grandes ... como números mayores que 2 ^ 200? Verá que cualquiera que sea el algoritmo que elija en esta página, todos tardarán mucho tiempo porque la mayoría de ellos dependen de gmp::factorize que utiliza el algoritmo Pollard-Rho . A partir de esta question , este algoritmo solo es razonable para números menores de 2 ^ 70. Actualmente estoy trabajando en mi propio algoritmo de factorización que implementará el Tamiz cuadrático , que debería llevar todos estos algoritmos al siguiente nivel.

Answer 2

Actualizar

El algoritmo se ha actualizado por completo y ahora implementa múltiples polinomios, así como algunas técnicas de cribado inteligentes que eliminan millones de controles. Además de los enlaces originales, este documento junto con este post de primo fueron muy útiles para esta última etapa (muchas felicitaciones a primo). Primo hace un gran trabajo explicando las agallas del QS en un espacio relativamente corto y también escribió un algoritmo bastante sorprendente (¡¡¡factorizará el número en la parte inferior, 38! + 1, en menos de 2 segundos !! ¡¡Locura !!).

Como prometí, a continuación se muestra mi humilde implementación en R del Tamiz cuadrático . He estado trabajando en este algoritmo esporádicamente desde que lo prometí a fines de enero. No intentaré explicarlo completamente (a menos que se solicite ... también, los enlaces a continuación hacen un muy buen trabajo) ya que es muy complicado y, con suerte, mis nombres de funciones hablan por sí mismos. Este ha demostrado ser uno de los algoritmos más desafiantes que he intentado ejecutar, ya que exige tanto desde el punto de vista de un programador como matemáticamente. He leído innumerables artículos y, en última instancia, he encontrado que estos cinco son los más útiles ( QSieve1 , QSieve2 , QSieve3 , QSieve4 , QSieve5 ).

NB: Este algoritmo, en su forma actual, no sirve muy bien como un algoritmo de factorización prima general. Si se optimizó aún más, debería ir acompañado de una sección de código que elimine primos más pequeños (es decir, menos de 10 ^ 5 como lo sugiere esta publicación ), luego llame a QuadSieveAll, verifique si estos son primos, y si no, llame a QuadSieveAll sobre estos dos factores, etc. hasta que se quede con todos los números primos (todos estos pasos no son tan difíciles). Sin embargo, el punto principal de esta publicación es resaltar el corazón del Tamiz cuadrático, por lo que los ejemplos a continuación son todos semiprimos (aunque factorizará la mayoría de los números impares que no contienen un cuadrado ... Además, no he visto un ejemplo de QS que no demostró un no-semiprime). Sé que el OP estaba buscando un método para devolver todos los factores y no la factorización principal, pero este algoritmo (si se optimiza aún más) junto con uno de los algoritmos anteriores sería una fuerza a tener en cuenta como un algoritmo de factorización general (especialmente dado que el OP necesitaba algo para el Proyecto Euler , que generalmente requiere mucho más que métodos de fuerza bruta). Por cierto, la función MyIntToBit es una variación de this respuesta y PrimeSieve es de una post que el Sr. Dontas apareció hace un tiempo (Kudos en eso también).

QuadSieveMultiPolysAll <- function(MyN, fudge1=0L, fudge2=0L, LenB=0L) {
### ''MyN'' is the number to be factored; ''fudge1'' is an arbitrary number
### that is used to determine the size of your prime base for sieving;
### ''fudge2'' is used to set a threshold for sieving;
### ''LenB'' is a the size of the sieving interval. The last three
### arguments are optional (they are determined based off of the
### size of MyN if left blank)

### The first 8 functions are helper functions

    PrimeSieve <- function(n) {
        n <- as.integer(n)
        if (n > 1e9) stop("n too large")
        primes <- rep(TRUE, n)
        primes[1] <- FALSE
        last.prime <- 2L
        fsqr <- floor(sqrt(n))
        while (last.prime <= fsqr) {
            primes[seq.int(last.prime^2, n, last.prime)] <- FALSE
            sel <- which(primes[(last.prime + 1):(fsqr + 1)])
            if (any(sel)) {
                last.prime <- last.prime + min(sel)
            } else {
                last.prime <- fsqr + 1
            }
        }
        MyPs <- which(primes)
        rm(primes)
        gc()
        MyPs
    }

    MyIntToBit <- function(x, dig) {
        i <- 0L
        string <- numeric(dig)
        while (x > 0) {
            string[dig - i] <- x %% 2L
            x <- x %/% 2L
            i <- i + 1L
        }
        string
    }

    ExpBySquaringBig <- function(x, n, p) {
        if (n == 1) {
            MyAns <- mod.bigz(x,p)
        } else if (mod.bigz(n,2)==0) {
            MyAns <- ExpBySquaringBig(mod.bigz(pow.bigz(x,2),p),div.bigz(n,2),p)
        } else {
            MyAns <- mod.bigz(mul.bigz(x,ExpBySquaringBig(mod.bigz(
                pow.bigz(x,2),p), div.bigz(sub.bigz(n,1),2),p)),p)
        }
        MyAns
    }

    TonelliShanks <- function(a,p) {
        P1 <- sub.bigz(p,1); j <- 0L; s <- P1
        while (mod.bigz(s,2)==0L) {s <- s/2; j <- j+1L}
        if (j==1L) {
            MyAns1 <- ExpBySquaringBig(a,(p+1L)/4,p)
            MyAns2 <- mod.bigz(-1 * ExpBySquaringBig(a,(p+1L)/4,p),p)
        } else {
            n <- 2L
            Legendre2 <- ExpBySquaringBig(n,P1/2,p)
            while (Legendre2==1L) {n <- n+1L; Legendre2 <- ExpBySquaringBig(n,P1/2,p)}
            x <- ExpBySquaringBig(a,(s+1L)/2,p)
            b <- ExpBySquaringBig(a,s,p)
            g <- ExpBySquaringBig(n,s,p)
            r <- j; m <- 1L
            Test <- mod.bigz(b,p)
            while (!(Test==1L) && !(m==0L)) {
                m <- 0L
                Test <- mod.bigz(b,p)
                while (!(Test==1L)) {m <- m+1L; Test <- ExpBySquaringBig(b,pow.bigz(2,m),p)}
                if (!m==0) {
                    x <- mod.bigz(x * ExpBySquaringBig(g,pow.bigz(2,r-m-1L),p),p)
                    g <- ExpBySquaringBig(g,pow.bigz(2,r-m),p)
                    b <- mod.bigz(b*g,p); r <- m
                }; Test <- 0L
            }; MyAns1 <- x; MyAns2 <- mod.bigz(p-x,p)
        }
        c(MyAns1, MyAns2)
    }

    SieveLists <- function(facLim, FBase, vecLen, sieveD, MInt) {
        vLen <- ceiling(vecLen/2); SecondHalf <- (vLen+1L):vecLen
        MInt1 <- MInt[1:vLen]; MInt2 <- MInt[SecondHalf]
        tl <- vector("list",length=facLim)

        for (m in 3:facLim) {
            st1 <- mod.bigz(MInt1[1],FBase[m])
            m1 <- 1L+as.integer(mod.bigz(sieveD[[m]][1] - st1,FBase[m]))
            m2 <- 1L+as.integer(mod.bigz(sieveD[[m]][2] - st1,FBase[m]))
            sl1 <- seq.int(m1,vLen,FBase[m])
            sl2 <- seq.int(m2,vLen,FBase[m])
            tl1 <- list(sl1,sl2)
            st2 <- mod.bigz(MInt2[1],FBase[m])
            m3 <- vLen+1L+as.integer(mod.bigz(sieveD[[m]][1] - st2,FBase[m]))
            m4 <- vLen+1L+as.integer(mod.bigz(sieveD[[m]][2] - st2,FBase[m]))
            sl3 <- seq.int(m3,vecLen,FBase[m])
            sl4 <- seq.int(m4,vecLen,FBase[m])
            tl2 <- list(sl3,sl4)
            tl[[m]] <- list(tl1,tl2)
        }
        tl
    }

    SieverMod <- function(facLim, FBase, vecLen, SD, MInt, FList, LogFB, Lim, myCol) {

        MyLogs <- rep(0,nrow(SD))

        for (m in 3:facLim) {
            MyBool <- rep(FALSE,vecLen)
            MyBool[c(FList[[m]][[1]][[1]],FList[[m]][[2]][[1]])] <- TRUE
            MyBool[c(FList[[m]][[1]][[2]],FList[[m]][[2]][[2]])] <- TRUE
            temp <- which(MyBool)
            MyLogs[temp] <- MyLogs[temp] + LogFB[m]
        }

        MySieve <- which(MyLogs > Lim)
        MInt <- MInt[MySieve]; NewSD <- SD[MySieve,]
        newLen <- length(MySieve); GoForIT <- FALSE

        MyMat <- matrix(integer(0),nrow=newLen,ncol=myCol)
        MyMat[which(NewSD[,1L] < 0),1L] <- 1L; MyMat[which(NewSD[,1L] > 0),1L] <- 0L
        if ((myCol-1L) - (facLim+1L) > 0L) {MyMat[,((facLim+2L):(myCol-1L))] <- 0L}
        if (newLen==1L) {MyMat <- matrix(MyMat,nrow=1,byrow=TRUE)}

        if (newLen > 0L) {
            GoForIT <- TRUE
            for (m in 1:facLim) {
                vec <- rep(0L,newLen)
                temp <- which((NewSD[,1L]%%FBase[m])==0L)
                NewSD[temp,] <- NewSD[temp,]/FBase[m]; vec[temp] <- 1L
                test <- temp[which((NewSD[temp,]%%FBase[m])==0L)]
                while (length(test)>0L) {
                    NewSD[test,] <- NewSD[test,]/FBase[m]
                    vec[test] <- (vec[test]+1L)
                    test <- test[which((NewSD[test,]%%FBase[m])==0L)]
                }
                MyMat[,m+1L] <- vec
            }
        }

        list(MyMat,NewSD,MInt,GoForIT)
    }

    reduceMatrix <- function(mat) {
        tempMin <- 0L; n1 <- ncol(mat); n2 <- nrow(mat)
        mymax <- 1L
        for (i in 1:n1) {
            temp <- which(mat[,i]==1L)
            t <- which(temp >= mymax)
            if (length(temp)>0L && length(t)>0L) {
                MyMin <- min(temp[t])
                if (!(MyMin==mymax)) {
                    vec <- mat[MyMin,]
                    mat[MyMin,] <- mat[mymax,]
                    mat[mymax,] <- vec
                }
                t <- t[-1]; temp <- temp[t]
                for (j in temp) {mat[j,] <- (mat[j,]+mat[mymax,])%%2L}
                mymax <- mymax+1L
            }
        }

        if (mymax<n2) {simpMat <- mat[-(mymax:n2),]} else {simpMat <- mat}
        lenSimp <- nrow(simpMat)
        if (is.null(lenSimp)) {lenSimp <- 0L}
        mycols <- 1:n1

        if (lenSimp>1L) {
            ## "Diagonalizing" Matrix
            for (i in 1:lenSimp) {
                if (all(simpMat[i,]==0L)) {simpMat <- simpMat[-i,]; next}
                if (!simpMat[i,i]==1L) {
                    t <- min(which(simpMat[i,]==1L))
                    vec <- simpMat[,i]; tempCol <- mycols[i]
                    simpMat[,i] <- simpMat[,t]; mycols[i] <- mycols[t]
                    simpMat[,t] <- vec; mycols[t] <- tempCol
                }
            }

            lenSimp <- nrow(simpMat); MyList <- vector("list",length=n1)
            MyFree <- mycols[which((1:n1)>lenSimp)];  for (i in MyFree) {MyList[[i]] <- i}
            if (is.null(lenSimp)) {lenSimp <- 0L}

            if (lenSimp>1L) {
                for (i in lenSimp:1L) {
                    t <- which(simpMat[i,]==1L)
                    if (length(t)==1L) {
                        simpMat[ ,t] <- 0L
                        MyList[[mycols[i]]] <- 0L
                    } else {
                        t1 <- t[t>i]
                        if (all(t1 > lenSimp)) {
                            MyList[[mycols[i]]] <- MyList[[mycols[t1[1]]]]
                            if (length(t1)>1) {
                                for (j in 2:length(t1)) {MyList[[mycols[i]]] <- c(MyList[[mycols[i]]], MyList[[mycols[t1[j]]]])}
                            }
                        }
                        else {
                            for (j in t1) {
                                if (length(MyList[[mycols[i]]])==0L) {MyList[[mycols[i]]] <- MyList[[mycols[j]]]}
                                else {
                                    e1 <- which(MyList[[mycols[i]]]%in%MyList[[mycols[j]]])
                                    if (length(e1)==0) {
                                        MyList[[mycols[i]]] <- c(MyList[[mycols[i]]],MyList[[mycols[j]]])
                                    } else {
                                        e2 <- which(!MyList[[mycols[j]]]%in%MyList[[mycols[i]]])
                                        MyList[[mycols[i]]] <- MyList[[mycols[i]]][-e1]
                                        if (length(e2)>0L) {MyList[[mycols[i]]] <- c(MyList[[mycols[i]]], MyList[[mycols[j]]][e2])}
                                    }
                                }
                            }
                        }
                    }
                }
                TheList <- lapply(MyList, function(x) {if (length(x)==0L) {0} else {x}})
                list(TheList,MyFree)
            } else {
                list(NULL,NULL)
            }
        } else {
            list(NULL,NULL)
        }
    }

    GetFacs <- function(vec1, vec2, n) {
        x <- mod.bigz(prod.bigz(vec1),n)
        y <- mod.bigz(prod.bigz(vec2),n)
        MyAns <- c(gcd.bigz(x-y,n),gcd.bigz(x+y,n))
        MyAns[sort.list(asNumeric(MyAns))]
    }

    SolutionSearch <- function(mymat, M2, n, FB) {

        colTest <- which(apply(mymat, 2, sum) == 0)
        if (length(colTest) > 0) {solmat <- mymat[ ,-colTest]} else {solmat <- mymat}

        if (length(nrow(solmat)) > 0) {
            nullMat <- reduceMatrix(t(solmat %% 2L))
            listSol <- nullMat[[1]]; freeVar <- nullMat[[2]]; LF <- length(freeVar)
        } else {LF <- 0L}

        if (LF > 0L) {
            for (i in 2:min(10^8,(2^LF + 1L))) {
                PosAns <- MyIntToBit(i, LF)
                posVec <- sapply(listSol, function(x) {
                    t <- which(freeVar %in% x)
                    if (length(t)==0L) {
                        0
                    } else {
                        sum(PosAns[t])%%2L
                    }
                })
                ansVec <- which(posVec==1L)
                if (length(ansVec)>0) {

                    if (length(ansVec) > 1L) {
                        myY <- apply(mymat[ansVec,],2,sum)
                    } else {
                        myY <- mymat[ansVec,]
                    }

                    if (sum(myY %% 2) < 1) {
                        myY <- as.integer(myY/2)
                        myY <- pow.bigz(FB,myY[-1])
                        temp <- GetFacs(M2[ansVec], myY, n)
                        if (!(1==temp[1]) && !(1==temp[2])) {
                            return(temp)
                        }
                    }
                }
            }
        }
    }

### Below is the main portion of the Quadratic Sieve

    BegTime <- Sys.time(); MyNum <- as.bigz(MyN); DigCount <- nchar(as.character(MyN))
    P <- PrimeSieve(10^5)
    SqrtInt <- .mpfr2bigz(trunc(sqrt(mpfr(MyNum,sizeinbase(MyNum,b=2)+5L))))

    if (DigCount < 24) {
        DigSize <- c(4,10,15,20,23)
        f_Pos <- c(0.5,0.25,0.15,0.1,0.05)
        MSize <- c(5000,7000,10000,12500,15000)

        if (fudge1==0L) {
            LM1 <- lm(f_Pos ~ DigSize)
            m1 <- summary(LM1)$coefficients[2,1]
            b1 <- summary(LM1)$coefficients[1,1]
            fudge1 <- DigCount*m1 + b1
        }

        if (LenB==0L) {
            LM2 <- lm(MSize ~ DigSize)
            m2 <- summary(LM2)$coefficients[2,1]
            b2 <- summary(LM2)$coefficients[1,1]
            LenB <- ceiling(DigCount*m2 + b2)
        }

        LimB <- trunc(exp((.5+fudge1)*sqrt(log(MyNum)*log(log(MyNum)))))
        B <- P[P<=LimB]; B <- B[-1]
        facBase <- P[which(sapply(B, function(x) ExpBySquaringBig(MyNum,(x-1)/2,x)==1L))+1L]
        LenFBase <- length(facBase)+1L
    } else if (DigCount < 67) {
        ## These values were obtained from "The Multiple Polynomial
        ## Quadratic Sieve" by Robert D. Silverman
        DigSize <- c(24,30,36,42,48,54,60,66)
        FBSize <- c(100,200,400,900,1200,2000,3000,4500)
        MSize <- c(5,25,25,50,100,250,350,500)

        LM1 <- loess(FBSize ~ DigSize)
        LM2 <- loess(MSize ~ DigSize)

        if (fudge1==0L) {
            fudge1 <- -0.4
            LimB <- trunc(exp((.5+fudge1)*sqrt(log(MyNum)*log(log(MyNum)))))
            myTarget <- ceiling(predict(LM1, DigCount))

            while (LimB < myTarget) {
                LimB <- trunc(exp((.5+fudge1)*sqrt(log(MyNum)*log(log(MyNum)))))
                fudge1 <- fudge1+0.001
            }
            B <- P[P<=LimB]; B <- B[-1]
            facBase <- P[which(sapply(B, function(x) ExpBySquaringBig(MyNum,(x-1)/2,x)==1L))+1L]
            LenFBase <- length(facBase)+1L

            while (LenFBase < myTarget) {
                fudge1 <- fudge1+0.005
                LimB <- trunc(exp((.5+fudge1)*sqrt(log(MyNum)*log(log(MyNum)))))
                myind <- which(P==max(B))+1L
                myset <- tempP <- P[myind]
                while (tempP < LimB) {
                    myind <- myind + 1L
                    tempP <- P[myind]
                    myset <- c(myset, tempP)
                }

                for (p in myset) {
                    t <- ExpBySquaringBig(MyNum,(p-1)/2,p)==1L
                    if (t) {facBase <- c(facBase,p)}
                }
                B <- c(B, myset)
                LenFBase <- length(facBase)+1L
            }
        } else {
            LimB <- trunc(exp((.5+fudge1)*sqrt(log(MyNum)*log(log(MyNum)))))
            B <- P[P<=LimB]; B <- B[-1]
            facBase <- P[which(sapply(B, function(x) ExpBySquaringBig(MyNum,(x-1)/2,x)==1L))+1L]
            LenFBase <- length(facBase)+1L
        }
        if (LenB==0L) {LenB <- 1000*ceiling(predict(LM2, DigCount))}
    } else {
        return("The number you''ve entered is currently too big for this algorithm!!")
    }

    SieveDist <- lapply(facBase, function(x) TonelliShanks(MyNum,x))
    SieveDist <- c(1L,SieveDist); SieveDist[[1]] <- c(SieveDist[[1]],1L); facBase <- c(2L,facBase)
    Lower <- -LenB; Upper <- LenB; LenB2 <- 2*LenB+1L; MyInterval <- Lower:Upper
    M <- MyInterval + SqrtInt ## Set that will be tested
    SqrDiff <- matrix(sub.bigz(pow.bigz(M,2),MyNum),nrow=length(M),ncol=1L)
    maxM <- max(MyInterval)
    LnFB <- log(facBase)

    ## N.B. primo uses 0.735, as his siever
    ## is more efficient than the one employed here
    if (fudge2==0L) {
        if (DigCount < 8) {
            fudge2 <- 0
        } else if (DigCount < 12) {
            fudge2 <- .7
        } else if (DigCount < 20) {
            fudge2 <- 1.3
        } else {
            fudge2 <- 1.6
        }
    }

    TheCut <- log10(maxM*sqrt(2*asNumeric(MyNum)))*fudge2
    myPrimes <- as.bigz(facBase)

    CoolList <- SieveLists(LenFBase, facBase, LenB2, SieveDist, MyInterval)
    GetMatrix <- SieverMod(LenFBase, facBase, LenB2, SqrDiff, M, CoolList, LnFB, TheCut, LenFBase+1L)

    if (GetMatrix[[4]]) {
        newmat <- GetMatrix[[1]]; NewSD <- GetMatrix[[2]]; M <- GetMatrix[[3]]
        NonSplitFacs <- which(abs(NewSD[,1L])>1L)
        newmat <- newmat[-NonSplitFacs, ]
        M <- M[-NonSplitFacs]
        lenM <- length(M)

        if (class(newmat) == "matrix") {
            if (nrow(newmat) > 0) {
                PosAns <- SolutionSearch(newmat,M,MyNum,myPrimes)
            } else {
                PosAns <- vector()
            }
        } else {
            newmat <- matrix(newmat, nrow = 1)
            PosAns <- vector()
        }
    } else {
        newmat <- matrix(integer(0),ncol=(LenFBase+1L))
        PosAns <- vector()
    }

    Atemp <- .mpfr2bigz(trunc(sqrt(sqrt(mpfr(2*MyNum))/maxM)))
    if (Atemp < max(facBase)) {Atemp <- max(facBase)}; myPoly <- 0L

    while (length(PosAns)==0L) {LegTest <- TRUE
        while (LegTest) {
            Atemp <- nextprime(Atemp)
            Legendre <- asNumeric(ExpBySquaringBig(MyNum,(Atemp-1L)/2,Atemp))
            if (Legendre == 1) {LegTest <- FALSE}
        }

        A <- Atemp^2
        Btemp <- max(TonelliShanks(MyNum, Atemp))
        B2 <- (Btemp + (MyNum - Btemp^2) * inv.bigz(2*Btemp,Atemp))%%A
        C <- as.bigz((B2^2 - MyNum)/A)
        myPoly <- myPoly + 1L

        polySieveD <- lapply(1:LenFBase, function(x) {
            AInv <- inv.bigz(A,facBase[x])
            asNumeric(c(((SieveDist[[x]][1]-B2)*AInv)%%facBase[x],
                        ((SieveDist[[x]][2]-B2)*AInv)%%facBase[x]))
        })

        M1 <- A*MyInterval + B2
        SqrDiff <- matrix(A*pow.bigz(MyInterval,2) + 2*B2*MyInterval + C,nrow=length(M1),ncol=1L)
        CoolList <- SieveLists(LenFBase, facBase, LenB2, polySieveD, MyInterval)
        myPrimes <- c(myPrimes,Atemp)
        LenP <- length(myPrimes)
        GetMatrix <- SieverMod(LenFBase, facBase, LenB2, SqrDiff, M1, CoolList, LnFB, TheCut, LenP+1L)

        if (GetMatrix[[4]]) {
            n2mat <- GetMatrix[[1]]; N2SD <- GetMatrix[[2]]; M1 <- GetMatrix[[3]]
            n2mat[,LenP+1L] <- rep(2L,nrow(N2SD))
            if (length(N2SD) > 0) {NonSplitFacs <- which(abs(N2SD[,1L])>1L)} else {NonSplitFacs <- LenB2}
            if (length(NonSplitFacs)<2*LenB) {
                M1 <- M1[-NonSplitFacs]; lenM1 <- length(M1)
                n2mat <- n2mat[-NonSplitFacs,]
                if (lenM1==1L) {n2mat <- matrix(n2mat,nrow=1)}
                if (ncol(newmat) < (LenP+1L)) {
                    numCol <- (LenP + 1L) - ncol(newmat)
                    newmat <-     cbind(newmat,matrix(rep(0L,numCol*nrow(newmat)),ncol=numCol))
                }
                newmat <- rbind(newmat,n2mat); lenM <- lenM+lenM1; M <- c(M,M1)
                if (class(newmat) == "matrix") {
                    if (nrow(newmat) > 0) {
                        PosAns <- SolutionSearch(newmat,M,MyNum,myPrimes)
                    }
                }
            }
        }
    }

    EndTime <- Sys.time()
    TotTime <- EndTime - BegTime
    print(format(TotTime))
    return(PosAns)
}

Con antiguo algoritmo QS

> library(gmp)
> library(Rmpfr)

> n3 <- prod(nextprime(urand.bigz(2, 40, 17)))
> system.time(t5 <- QuadSieveAll(n3,0.1,myps))
  user  system elapsed 
164.72    0.77  165.63 
> system.time(t6 <- factorize(n3))
user  system elapsed 
0.1     0.0     0.1 
> all(t5[sort.list(asNumeric(t5))]==t6[sort.list(asNumeric(t6))])
[1] TRUE

Con el nuevo algoritmo Muli-polinomio QS

> QuadSieveMultiPolysAll(n3)
[1] "4.952 secs"
Big Integer (''bigz'') object of length 2:
[1] 342086446909 483830424611

> n4 <- prod(nextprime(urand.bigz(2,50,5)))
> QuadSieveMultiPolysAll(n4)   ## With old algo, it took over 4 hours
[1] "1.131717 mins"
Big Integer (''bigz'') object of length 2:
[1] 166543958545561 880194119571287

> n5 <- as.bigz("94968915845307373740134800567566911")   ## 35 digits
> QuadSieveMultiPolysAll(n5)
[1] "3.813167 mins"
Big Integer (''bigz'') object of length 2:
[1] 216366620575959221 438925910071081891

> system.time(factorize(n5))   ## It appears we are reaching the limits of factorize
   user  system elapsed 
 131.97    0.00  131.98

Nota al margen : el número n5 de arriba es un número muy interesante. Mira here

¡¡¡¡El punto de quiebre!!!!

> n6 <- factorialZ(38) + 1L   ## 45 digits
> QuadSieveMultiPolysAll(n6)
[1] "22.79092 mins"
Big Integer (''bigz'') object of length 2:
[1] 14029308060317546154181 37280713718589679646221

> system.time(factorize(n6))   ## Shut it down after 2 days of running

Último triunfo (50 dígitos)

> n9 <- prod(nextprime(urand.bigz(2,82,42)))
> QuadSieveMultiPolysAll(n9)
[1] "12.9297 hours"
Big Integer (''bigz'') object of length 2:
[1] 2128750292720207278230259 4721136619794898059404993

## Based off of some crude test, factorize(n9) would take more than a year.

Debe notarse que el QS generalmente no funciona tan bien como el algoritmo rho de Pollard en números más pequeños y la potencia del QS comienza a manifestarse a medida que los números se hacen más grandes. De todos modos, como siempre, espero que esto inspire a alguien! ¡Aclamaciones!

Answer 3

El siguiente enfoque ofrece resultados correctos, incluso en casos de números realmente grandes (que deben pasarse como cadenas). Y es realmente rápido.

# TEST
# x <- as.bigz("12345678987654321")
# all_divisors(x)
# all_divisors(x*x)

# x <- pow.bigz(2,89)-1
# all_divisors(x)

library(gmp)
  options(scipen =30)

  sort_listz <- function(z) {
  #==========================
    z <- z[order(as.numeric(z))] # sort(z)
  } # function  sort_listz  


  mult_listz <- function(x,y) {
   do.call(''c'', lapply(y, function(i) i*x)) 
  } 


  all_divisors <- function(x) {
  #==========================  
  if (abs(x)<=1) return(x) 
  else {

    factorsz <- as.bigz(factorize(as.bigz(x))) # factorize returns up to
    # e.g. x= 12345678987654321  factors: 3 3 3 3 37 37 333667 333667

    factorsz <- sort_listz(factorsz) # vector of primes, sorted

    prime_factorsz <- unique(factorsz)
    #prime_ekt <- sapply(prime_factorsz, function(i) length( factorsz [factorsz==i]))
    prime_ekt <- vapply(prime_factorsz, function(i) sum(factorsz==i), integer(1), USE.NAMES=FALSE)
    spz <- vector() # keep all divisors 
    all <-1
    n <- length(prime_factorsz)
    for (i in 1:n) {
      pr <- prime_factorsz[i]
      pe <- prime_ekt[i]
      all <- all*(pe+1) #counts all divisors 

      prz <- as.bigz(pr)
      pse <- vector(mode="raw",length=pe+1) 
      pse <- c( as.bigz(1), prz)

      if (pe>1) {
        for (k in 2:pe) {
          prz <- prz*pr
          pse[k+1] <- prz
        } # for k
      } # if pe>1

      if (i>1) {
       spz <- mult_listz (spz, pse)         
      } else {
       spz <- pse;
      } # if i>1
    } #for n
    spz <- sort_listz (spz)

    return (spz)
  }  
  } # function  factors_all_divisors  

  #====================================

Versión refinada, muy rápida. El código sigue siendo simple, legible y limpio.

PRUEBA

#Test 4 (big prime factor)
x <- pow.bigz(2,256)+1 # = 1238926361552897 * 93461639715357977769163558199606896584051237541638188580280321
 system.time(z2 <- all_divisors(x))
#   user  system elapsed 
 #  19.27    1.27   20.56


 #Test 5 (big prime factor)
x <- as.bigz("12345678987654321321") # = 3 * 19 * 216590859432531953

 system.time(x2 <- all_divisors(x^2))
#user  system elapsed 
 #25.65    0.00   25.67  

Answer 4

Para dar seguimiento a mi comentario (gracias a @Ramnath por mi error tipográfico), el método de fuerza bruta parece funcionar razonablemente bien aquí en mi máquina de 64 bits y 8 gigas:

FUN <- function(x) {
    x <- as.integer(x)
    div <- seq_len(abs(x))
    factors <- div[x %% div == 0L]
    factors <- list(neg = -factors, pos = factors)
    return(factors)
}

Algunos ejemplos:

> FUN(100)
$neg
[1]   -1   -2   -4   -5  -10  -20  -25  -50 -100

$pos
[1]   1   2   4   5  10  20  25  50 100

> FUN(-42)
$neg
[1]  -1  -2  -3  -6  -7 -14 -21 -42

$pos
[1]  1  2  3  6  7 14 21 42

#and big number

> system.time(FUN(1e8))
   user  system elapsed 
   1.95    0.18    2.14 

Answer 5

Puede obtener todos los factores de los factores primos. gmp calcula muy rápidamente.

library(gmp)
library(plyr)

get_all_factors <- function(n)
{
  prime_factor_tables <- lapply(
    setNames(n, n), 
    function(i)
    {
      if(i == 1) return(data.frame(x = 1L, freq = 1L))
      plyr::count(as.integer(gmp::factorize(i)))
    }
  )
  lapply(
    prime_factor_tables, 
    function(pft)
    {
      powers <- plyr::alply(pft, 1, function(row) row$x ^ seq.int(0L, row$freq))
      power_grid <- do.call(expand.grid, powers)
      sort(unique(apply(power_grid, 1, prod)))
    }
  )
}

get_all_factors(c(1, 7, 60, 663, 2520, 75600, 15876000, 174636000, 403409160000))

Answer 6

Mucho ha cambiado en el lenguaje R desde que se hizo esta pregunta originalmente. En la versión 0.6-3del numberspaquete, divisorsse incluyó la función que es muy útil para obtener todos los factores de un número. Satisfará las necesidades de la mayoría de los usuarios, sin embargo, si está buscando una velocidad bruta o está trabajando con números más grandes, necesitará un método alternativo. He creado dos nuevos paquetes (parcialmente inspirados en esta pregunta, que podría agregar) que contienen funciones altamente optimizadas dirigidas a problemas como este. El primero es RcppAlgosy el otro es bigIntegerAlgos.

RcppAlgos

RcppAlgoscontiene dos funciones para obtener divisores de números menores que 2^53 - 1: divisorsRcpp(una función vectorizada para obtener rápidamente la factorización completa de muchos números) y divisorsSieve(genera rápidamente la factorización completa en un rango). Primero, factorizamos muchos números aleatorios usando divisorsRcpp:

library(gmp)  ## for all_divisors by @GeorgeDontas
library(RcppAlgos)
library(numbers)
options(scipen = 999)
set.seed(42)
testSamp <- sample(10^10, 10)

## vectorized so you can pass the entire vector as an argument
testRcpp <- divisorsRcpp(testSamp)
testDontas <- lapply(testSamp, all_divisors)

identical(lapply(testDontas, as.numeric), testRcpp)
[1] TRUE

Y ahora, factoriza muchos números en un rango usando divisorsSieve:

system.time(testSieve <- divisorsSieve(10^13, 10^13 + 10^5))
 user  system elapsed 
0.586   0.014   0.602 

system.time(testDontasSieve <- lapply((10^13):(10^13 + 10^5), all_divisors))
  user  system elapsed 
54.327   0.187  54.655 

identical(lapply(testDontasSieve, asNumeric), testSieve)
[1] TRUE

Ambas divisorsRcppy divisorsSieveson funciones agradables que son flexibles y eficientes, sin embargo están limitadas a 2^53 - 1.

bigIntegerAlgos

El bigIntegerAlgospaquete presenta dos funciones, divisorsBig& quadraticSieve, que están diseñadas para números muy grandes. Se enlazan directamente a la biblioteca C de gmp . Para divisorsBig, tenemos:

library(bigIntegerAlgos)
## testSamp is defined above... N.B. divisorsBig is not quite as
## efficient as divisorsRcpp. This is so because divisorsRcpp
## can take advantage of more efficient data types..... it is
## still blazing fast!! See the benchmarks below for reference.
testBig <- divisorsBig(testSamp)

identical(testDontas, testBig)
[1] TRUE

Y aquí están los puntos de referencia definidos en mi publicación original (NB MyFactorsse reemplaza por divisorsRcppy divisorsBig).

## Case 2
library(rbenchmark)
set.seed(199)
samp <- sample(10^9, 10^5)
benchmark(RcppAlgos=divisorsRcpp(samp),
          bigIntegerAlgos=divisorsBig(samp),
          DontasDivs=lapply(samp, all_divisors),
          replications=10,
          columns = c("test", "replications", "elapsed", "relative"),
          order = "relative")

             test replications elapsed relative
1       RcppAlgos           10   8.021    1.000
2 bigIntegerAlgos           10  15.246    1.901
3      DontasDivs           10 400.284   49.905

## Case 3
set.seed(97)
samp <- sample(10^6, 10^4)
benchmark(RcppAlgos=divisorsRcpp(samp),
          bigIntegerAlgos=divisorsBig(samp),
          numbers=lapply(samp, divisors),      ## From the numbers package
          DontasDivs=lapply(samp, all_divisors),
          CottonDivs=lapply(samp, get_all_factors),
          ChaseDivs=lapply(samp, FUN),
          replications=5,
          columns = c("test", "replications", "elapsed", "relative"),
          order = "relative")

             test replications elapsed relative
1       RcppAlgos            5   0.098    1.000
2 bigIntegerAlgos            5   0.330    3.367
3         numbers            5  11.613  118.500
4      DontasDivs            5  16.093  164.214
5      CottonDivs            5  60.428  616.612
6       ChaseDivs            5 342.608 3496.000

Los siguientes puntos de referencia demuestran el verdadero poder del algoritmo subyacente en la divisorsBigfunción. El número que se factoriza es un poder de 10, por lo que el paso de factorización principal puede ignorarse casi por completo (por ejemplo, los system.time(factorize(pow.bigz(10,30)))registros 0en mi máquina). Por lo tanto, la diferencia en el tiempo se debe únicamente a la rapidez con que los factores primos se pueden combinar para producir todos los factores.

library(microbenchmark)
powTen <- pow.bigz(10,30)
microbenchmark(divisorsBig(powTen), all_divisors(powTen), unit = "relative")
Unit: relative
                expr      min      lq     mean   median       uq      max neval
 divisorsBig(powTen)  1.00000  1.0000  1.00000  1.00000  1.00000  1.00000   100
all_divisors(powTen) 32.35507 33.3054 33.28786 33.31253 32.11571 40.39236   100

## Negative numbers show an even greater increase in efficiency
negPowTen <- powTen * -1
microbenchmark(divisorsBig(negPowTen), all_divisors(negPowTen), unit = "relative")
Unit: relative
                   expr      min       lq    mean   median       uq      max neval
 divisorsBig(negPowTen)  1.00000  1.00000  1.0000  1.00000  1.00000  1.00000   100
all_divisors(negPowTen) 46.42795 46.22408 43.7964 47.93228 45.33406 26.64657   100

quadraticSieve

Os dejo con dos demostraciones de quadraticSieve.

n5 <- as.bigz("94968915845307373740134800567566911")
system.time(print(quadraticSieve(n5)))
Big Integer (''bigz'') object of length 2:
[1] 216366620575959221 438925910071081891
 user  system elapsed 
4.154   0.021   4.175   ## original time was 3.813167 mins or 228.8 seconds ~ 50x slower

n9 <- prod(nextprime(urand.bigz(2, 82, 42)))
system.time(print(quadraticSieve(n9)))
Big Integer (''bigz'') object of length 2:
[1] 2128750292720207278230259 4721136619794898059404993
    user   system  elapsed 
1010.404    2.715 1013.184   ## original time was 12.9297 hours or 46,547 seconds ~ 46x slower

Como puede ver, quadraticSievees mucho más rápido que el original QuadSieveMultiPolysAll, sin embargo, aún queda mucho trabajo por hacer. Hay investigaciones en curso para mejorar esta función con el objetivo actual de tener n9en cuenta en menos de un minuto. También hay planes de vectorización quadraticSieve, así como la integración divisorsBigcon quadraticSieve, en la actualidad, se unidas al mismo algoritmo que gmp::factorizeutiliza.