Perfect Power Predicate

March 13, 2012

A number n is a perfect power if there exists a b and e for which be = n. For instance 216 = 63 = 23 · 33 is a perfect power, but 72 = 23 · 32 is not. Testing for perfect powers is similar to other powering predicates we have seen, and is useful in some factoring algorithms.

The trick to determining if a number is a perfect power is to know that, if the number is a perfect power, then the exponent e must be less than log2 n, because if e is greater then 2e will be greater than n. Further, it is only necessary to test prime es, because if a number is a perfect power to a composite exponent it will also be a perfect power to the prime factors of the composite component; for instance, 215 = 32768 = 323 = 85 is a perfect cube root and also a perfect fifth root.

Your task is to write a function to determine if a number is a perfect power. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

About these ads

Pages: 1 2

14 Responses to “Perfect Power Predicate”

  1. Johann Hibschman said

    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.

  2. programmingpraxis said

    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?

  3. Johann Hibschman said

    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.

  4. Johann Hibschman said

    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))

  5. Mike said

    Python code equivalent to Johanns’

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

    from collections import Counter
    from utils import factor, gcd
    
    def is_perfect_power(n):
        return reduce(gcd, Counter(factor(n)).values()) > 1
    
    
  6. programmingpraxis said

    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?

  7. Mike said

    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.

    from math import log
    from utils import primes_to
    
    def is_perfect_power(n):
        b = 2
        for e in reversed(primes_to(int(log(n,2)))):
            pwr = b**e
            while pwr < n:
                b += 1
                pwr = b**e
    
            if pwr == n:
                return True
    
        return False
    
    
  8. phillip said

    i think the best i can do is this
    perfect power

  9. ardnew said

    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.

    use strict;
    use warnings;
    
    #
    # naive base 2 logarithm
    #
    sub log2
    {
      my ($n, $c) = (int(shift), 0);
      return 0 if $n <= 1; # hehe
      ++$c while $n >>= 1;
      return $c;
    }
    
    #
    # Modular exponentiation (using binary method)
    #
    sub mod_power
    {
      my ($b, $e, $m) = @_;
      my $r = 1;
      
      while ($e > 0)
      {
        $r = ($r * $b) % $m if $e & 1;
        $e >>= 1;
        $b = ($b * $b) % $m;
      }  
      return $r;
    }
    
    #
    # Miller-Rabin primality test
    #
    sub miller_rabin
    {
      my ($n, $k) = @{$_[0]}; 
      
      $k //= 5; # default accuracy
      
      # don't test unless $n > 3 and odd
      if ($n > 3 and $n & 1)
      {
        my ($s, $d, $a, $x);    
        
        # LSB(n) == n & -n, so LSB(n) - 1 effectively flips all 
        # trailing unset bits of n; then use unpack to count set bits.
        $s = unpack("%32b*", pack("L", (($n - 1) & (1 - $n)) - 1));
        $d = $n >> $s;
        
        # try to be a little more deterministic
        $k = $n - 3 if $n - 3 < $k;    
        
        # generate k distinct random integers to prevent duplicate bases 
        # and thus increase accuracy
        my %base = ();
        
        for (1 .. $k)
        {
          $a = int(rand($n - 3)) + 2; # interval [2, n - 2]
          exists $base{$a} and redo or $base{$a} = undef;
        }
        
        # repeatedly interrogate the witnesses            
        OUT: while ($a = each %base)
        {
          $x = mod_power($a, $d, $n);     
    
          next if ($x == 1 or $x == $n - 1);
    
          foreach my $r (1 .. $s - 1)
          {
            $x = mod_power($x, 2, $n);
            return 0 if $x == 1;
            next OUT if $x == $n - 1;
          }     
          return 0;
        }    
        return 1;
      }
      else
      {
        return ($n == 2 or $n == 3);
      }
    }
    
    #
    # MAIN PROCEDURE
    #
    # Accepts as input a positive integer n and prints all 
    # integer pairs b and e (prime) such that b^e = n
    #
    sub perfect_power
    {
      my $n = shift;
      
      for (my $e = log2($n); $e > 0; --$e)
      {
        next unless miller_rabin([$e]);        
        
        for (my ($b, $p) = (2, 0); ($p = $b ** $e) <= $n; ++$b)
        {
          printf("%d^%d = %.f\n", $b, $e, $n) if $p == $n;
        }
      }
    }
    
    scalar @ARGV or 
      die "\nusage:\n\tperl $0 <integer>\n";
      
    perfect_power($ARGV[0]);
    
  10. ardnew said

    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

  11. ardnew said

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

  12. Johann Hibschman said

    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.

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

    import math
    def powerPredicate(n):
    	e_limit  = int(math.log(n,2))
    	b_limit = int(math.sqrt(n))
    	li = [(b,e) for b in xrange(2,b_limit+1) for e in xrange(2,e_limit+1) if pow(b,e) == n]
    	
    	if li:
    		print li
    
  14. danaj said

    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.

    #!/usr/bin/env perl
    use strict; use warnings;
    use Math::BigInt;
    use Math::Prime::Util qw/-nobigint next_prime/;
    
    die "Usage: $0 <positive integer>\n" unless scalar @ARGV > 0;
    my $n = $ARGV[0];
    die "$n is not a positive integer" if $n =~ tr/0123456789//c;
    $n = Math::BigInt->new("$n") if $n > ~0;
    print is_perfect_power($n), "\n";
    
    sub is_perfect_power {
      my $n = shift;
      return 0 if $n <= 3 || $n != int($n);
      return 1 if ($n & ($n-1)) == 0;                       # Power of 2
      # double pow doesn't give all bits of accuracy after this
      $n = Math::BigInt->new("$n") if ref($n) ne 'Math::BigInt' && $n > 2**52;
      my $log2n = 0; { my $num = $n; $log2n++ while $num >>= 1; }
      if (ref($n) eq 'Math::BigInt') {
        for (my $e = 2; $e <= $log2n; $e = next_prime($e)) {
          return 1 if $n->copy()->broot($e)->bpow($e) == $n;
        }
      } else {
        my $b = int(sqrt($n)+0.5); return 1 if int($b*$b) == $n; # perfect square
        for (my $e = 3; $e <= $log2n; $e = next_prime($e)) {
          my $root = int($n ** (1.0 / $e) + .5);
          return 1 if $root > 2 && int($root ** $e + .5) == $n;
        }
      }
      0;
    }
    
    # Simple version that always uses BigInt.
    sub is_perfect_power_simple {
      my $n = shift;
      return 0 if $n <= 3 || $n != int($n);
      $n = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt';
      my $log2n = 0; { my $num = $n; $log2n++ while $num >>= 1; }
      for (my $e = 2; $e <= $log2n; $e = next_prime($e)) {
        return 1 if $n->copy()->broot($e)->bpow($e) == $n;
      }
      0;
    }
    
    $ perl is-perfect-power.pl 148675665359980297048795508874724049089546782584077934753925649
    1
    
    $ factor.pl 148675665359980297048795508874724049089546782584077934753925649
    148675665359980297048795508874724049089546782584077934753925649: 1234567890123493 1234567890123493 9876543210987701 9876543210987701
    

    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:

    /* returns 1 and sets sqrtn to isqrt(n) if a perfect square.  Returns 0 if not. */
    static int is_perfect_square(UV n, UV* sqrtn)
    {
      UV m;  /* lm */
      m = n & 127;
      if ((m*0x8bc40d7d) & (m*0xa1e2f5d1) & 0x14020a)  return 0;
      m = n % 63;
      if ((m*0x3d491df7) & (m*0xc824a9f9) & 0x10f14008) return 0;
      m = sqrt(n);
      if (n != (m*m))
        return 0;
    
      if (sqrtn != 0) *sqrtn = m;
      return 1;
    }
    

    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.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 632 other followers

%d bloggers like this: