## Perfect Power Predicate

### March 13, 2012

We begin with the function `(ilog b n)`

that finds the largest integer *e* such that *b ^{e}* ≤

*n*; that function is part of the Standard Prelude, but we restate it here:

`(define (ilog b n)`

(let loop1 ((lo 0) (b^lo 1) (hi 1) (b^hi b))

(if (< b^hi n) (loop1 hi b^hi (* hi 2) (* b^hi b^hi))

(let loop2 ((lo lo) (b^lo b^lo) (hi hi) (b^hi b^hi))

(if (<= (- hi lo) 1) (if (= b^hi n) hi lo)

(let* ((mid (quotient (+ lo hi) 2))

(b^mid (* b^lo (expt b (- mid lo)))))

(cond ((< n b^mid) (loop2 lo b^lo mid b^mid))

((< b^mid n) (loop2 mid b^mid hi b^hi))

(else mid))))))))

`Ilog`

works by binary search. `Loop1`

counts up from 1 until it finds a *lo*/*hi* pair such that *b ^{lo}* <

*b*, and

^{hi}`loop2`

then continually bisects the interval until the result is found.We also need the `iroot`

function from the Standard Prelude. This function works by Newton’s method; if *x*_{0} is an approximation to a function *f*, then *x*_{1} = *x*_{0} − *f*(*x*_{0}) / *f* '(*x*_{0}) is a better approximation, where *f* ' is the derivative of *f*.

`(define (iroot k n)`

(let loop ((u n))

(let* ((s u)

(t (+ (* (- k 1) s)

(quotient n (expt s (- k 1)))))

(u (quotient t k)))

(if (< u s) (loop u) s))))

Now the perfect power predicate is easy; we use a `primes`

function based on the Sieve of Eratosthenes that we have seen many times previously:

`(define (perfect-power? n)`

(if (not (and (integer? n) (positive? n)))

(error 'perfect-power? "must be positive integer")

(let loop ((ps (primes (ilog 2 n))))

(if (null? ps) #f

(let ((x (iroot (car ps) n)))

(if (= (expt x (car ps)) n) x

(loop (cdr ps))))))))

Here are some examples:

`> (perfect-power? 32768)`

32

> (perfect-power? 205442259656281392806087233013)

53

> (perfect-power 213)

#f

Note that, in the case where *e* is composite, the return value will not be the fundamental base.

You can run the program at http://programmingpraxis.codepad.org/F6lXXfYo.

Here’s the smart-ass J version:

isPerfectPower=: 3 : ‘1<+./#/.~q:y'

It feels like cheating when your language has a "prime factors of" primitive.

Or you could say that you have chosen the right tool for the job.

For those of us who don’t speak J, could you provide an explanation?

Sure.

“q:” returns the prime factors of a number. “q: 216” returns “2 2 2 3 3 3”, “q: 700” returns “2 2 5 5 7”.

“#/.~” returns the counts of the unique items. “#/.~q: 216” is “3 3”. “#/.~q: 700” is “2 2 1”.

“+.” is gcd. “+./” is “foldr gcd”. So “+./#/.~q: 216” is “3”. “+/#/.~q: 700” is “1”.

And if the gcd of the counts of all the prime factors is greater than 1, it’s a perfect power. (Unless my brain is being very slow this morning.)

“foo=: 3 : ‘blah'” is just function-definition boilerplate.

The details of that inscrutable-looking “#/.~” are interesting, at least to me. “#” is count. “~” is reflexive call: “f~ x” is “x f x”. “/.” is a “keyed operation”. Given “x f/.y”, the unique items of x (the “key”) are used to break y into groups, then the function f is applied to each group. So, here we’re using the list of factors to classify itself, then we’re counting the number of values in each group. It’s an interesting higher-order function that I’ve not seen in too many places. The rather terse J documentation on it is here: http://jsoftware.com/help/dictionary/d421.htm.

Okay, one more J version, then I’m done.

ilog=: <.@^.

iroot=: : _1 p: y’

isPerfectPower=: 3 : ‘+./(primesTo 2 ilog y) isEvenRoot y’

“@” is composition, “:” is 1+, “i.” is integers-up-to.

Diving into “isEvenRoot”, we have the equivalent of:

(define (left x y) x) ; [

(define (right x y) y) ; ]

(define (fork f g h) (lambda (x y) (g (f x y) (h x y)) ; implicit

(define isEvenRoot (fork right = (fork iroot expt left))

Python code equivalent to Johanns’

Note: utils is my collection of handy functions from solving other Praxis problems

Johann and Mike: Those functions work if you can factor n. But if n is the product of two large primes, it’s unlikely that you will get an answer any time soon. You could assume that if you can’t factor n in some reasonable time, it’s not a perfect power, but if n is the square of two large primes, you’re in trouble. For instance, is 148675665359980297048795508874724049089546782584077934753925649 a perfect power?

Here is an iterative solution.

Basically, start with the smallest base and the largest exponent. If b**e is too large decrease the exponent. If b**e is too small increase the base. Repeat until b**e == n or all the exponents have been tried.

i think the best i can do is this

perfect power

All positive integers satisfy this definition of a perfect power, you need to add the condition e > 1.

Basically the same algorithm as Mike’s, except with a bunch of re-invented wheels.

OP: Factoring is difficult, but won’t the prime sieve also suffer from a time/memory issue if the input is a large prime squared? Sieving all primes up to 12193263113702594790275394159593 (to match your example (thanks wolframalpha)) would take a lot of memory and time… Once you had this table though, your root search seems to be the most efficient of anyone’s posted.

My algorithm also chokes horribly on very large input (greater than 64 bit integers). I even tried integrating perl’s arbitrary precision library, but I failed miserably

Oops, disregard me, you only need to generate the primes up to log2(n) =)

WordPress formatting ate my previous attempt to post a version that handles large numbers gracefully.

pastebin

This works for 148675665359980297048795508874724049089546782584077934753925649, saying that, yes, it’s a perfect square.

This is the least I can do for now. The function displays all b and e where b^e = n.

A version in Perl that runs much faster than the Perl version above, and handles bigints. Given that the *reason* for a is_perfect_power function is typically as an initial step in factoring or primality testing (e.g. SQUFOF, AKS), using factoring to answer the question is a non-starter.

Doing just the bigint part (the is_perfect_power_simple function) is simpler, isolated from quirks of the C library and some old Perl (pre 5.8.8) screwiness with 64-bit numbers, and still reasonably fast. The native ints/floats method is ~400x faster on my computer so worth thinking about.

For perfect square testing (e.g. for HOLF or SQUFOF), I’ve found the first one or two bloom filters from http://mersenneforum.org/showpost.php?p=110896 works better than anything else I’ve seen. E.g. for non-bigints:

His full sequence is really meant to speed up the test for bigints, trying to avoid a bigint sqrt, but it still helps a bit for native computations.