Re: [S] Prime divisors

Paul A Tukey (
Wed, 16 Sep 1998 18:27:15 -0400 (EDT)

This discussion has been fun.

Seems to me we've been gradually heading toward
computing the prime factorization of a number.
That is, a collection of prime numbers (with possible
duplication) whose product is the given number.

A recursive layer on top of Frank Harrell's factorize() does
it -- but the code below uses a slightly shortened version
of factorize() that only returns the smallest divisor > 1.

fac <- function(n) {
p <- n/(z <- 1:floor(sqrt(n)))
z <- z[trunc(p) == p]
c(z, rev(n/z))[2]

pfac <- function(n, nn = 0)
if(nn == 0)
pfac(n, fac(n))
else if(n <= nn)
else c(nn, pfac(n/nn))

Note that prod(pfac(n)) == n.

Now I'm sure someone can write a more elegant version.
Also, recursion is probably neither memory-efficient
nor CPU-efficient in Splus.

-- Paul Tukey

This message was distributed by To unsubscribe
send e-mail to with the BODY of the
message: unsubscribe s-news