[S] Summary: Factors (divisors) of an inreger

mona kanaan (M.N.Kanaan@open.ac.uk)
Fri, 11 Sep 1998 08:52:59 +0100 (BST)


Thanks a lot, for everybody who replied to my query.
Here is a summary of what was passed on.
The first two codes due to Bill Venables and Bill Dunlap give the Prime
divisors of an integer(this is what i was actually looking for), the last
code gives all divisors but is not efficient for "large"
integers (this is what i was trying to avoid).

Thanks again
Mona

-----------------------------------------------------------------

Bill Venables

This turns out to be a pretty little programming exercise.
Here's a vectorized version, even, although it only returns the
*prime* divisors, not all the devisors. That a supplmentary
exercise...

-----------------------------------------------------------------
> factorize <- function(n) {
if(!is.numeric(n))
stop("cannot factorize non-numeric arguments")
if(length(n) > 1) {
l <- list()
for(i in seq(along = n))
l[[i]] <- Recall(n[i])
return(l)
}
if(n != round(n) || n < 2)
return(n)
tab <- 2:n
fac <- numeric(0)
while(n > 1) {
while(n %% tab[1] == 0) {
fac <- c(fac, tab[1])
n <- n/tab[1]
}
tab <- tab[tab <= n]
omit <- tab[1] * c(1, tab[tab <= n/tab[1]])
tab <- tab[ - match(omit, tab, nomatch = 0)]
}
fac
}
> factorize(6)
[1] 2 3
> factorize(4:8)
[[1]]:
[1] 2 2

[[2]]:
[1] 5

[[3]]:
[1] 2 3

[[4]]:
[1] 7

[[5]]:
[1] 2 2 2
-----------------------------------------------------------------
-----------------------------------------------------------------
-----------------------------------------------------------------
Bill Dunlap

I use the following factors(), which uses the enclosed primes():

> factors( round(gamma(13:14)))
$"479001600":
[1] 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 5 5 7 11

$"6227020800":
[1] 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 5 5 7 11 13

You can use table() to collect repeated factors

> lapply( factors( round(gamma(13:14))), table)
$"479001600":
2 3 5 7 11
10 5 2 1 1

$"6227020800":
2 3 5 7 11 13
10 5 2 1 1 1

"factors"<-
function(x)
{
factor1 <- function(y, max.factor, .Primes)
{
if(missing(.Primes))
.Primes <- primes(max.factor)
else .Primes <- primes(max.factor, .Primes)
f <- numeric(0)
while(y > 1) {
# note 1 has no factors according to this
which <- y %% .Primes == 0
if(sum(which) == 0) {
f <- c(f, y)
break
}
else f <- c(f, .Primes[which])
y <- y/prod(.Primes[which])
}
val <- sort(f)
if(length(val) && any(big <- val > max.factor^2)) {
if(sum(big) != 1)
stop("internal error: sum(big)!=1")
val <- sort(c(val[!big], Recall(val[big], min(ceiling(
sqrt(val[big])), max.factor^2), .Primes)))
}
val
}
val <- lapply(x, factor1, 43)
names(val) <- as.character(x)
val
}
"primes"<-
function(n, .Primes = c(2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43))
{
if(max(.Primes) < n) {
# compute longer .Primes by sieve
.Primes <- seq(from = 2, to = n)
for(i in 1:length(.Primes)) {
composite <- .Primes %% .Primes[i] == 0
composite[i] <- F
if(all(!composite))
break
.Primes <- .Primes[!composite]
if(i >= length(.Primes))
break
}
}
.Primes[.Primes <= n]
}

factors.simple() is easier to understand and is faster on small numbers
but can work very slowly on large numbers with lots of small factors
(like numbers arising in combinatorics).

"factors.simple"<-
function(x)
{
factor1 <- function(y, .Primes)
{
f <- numeric(0)
while(y > 1) {
# note 1 has no factors according to this
which <- y %% .Primes == 0
if(sum(which) == 0) {
f <- c(f, y)
break
}
else f <- c(f, .Primes[which])
y <- y/prod(.Primes[which])
}
sort(f)
}
val <- lapply(x, factor1, primes(ceiling(sqrt(max(x)))))
names(val) <- as.character(x)
val
}
----------------------------------------------------------------------------
----------------------------------------------------------------------------
----------------------------------------------------------------------------

Guido Schwarzer ,Gardar Johannesson, Remy vande Ven, Henrik Aalborg-Nielsen

DIV <- function(N){
N.seq <- 1:N
N.seq[(N %% N.seq) == 0]
}

_/_/ _/_/ _/_/_/_/ _/_/ _/ _/_/_/_/ email:M.N.Kanaan@open.ac.uk
_/ _/ _/ _/ _/ _/ _/ _/ _/ _/ _/ Statistics Dept.
_/ _/ _/ _/ _/ _/ _/ _/ _/ _/_/_/_/ The Open University
_/ _/_/ _/ _/ _/ _/ _/ _/ _/ _/ Milton Keynes
_/ _/ _/ _/_/_/_/ _/ _/_/ _/ _/ U.K. MK7 6AA

-----------------------------------------------------------------------
This message was distributed by s-news@wubios.wustl.edu. To unsubscribe
send e-mail to s-news-request@wubios.wustl.edu with the BODY of the
message: unsubscribe s-news