Adventures of a Programmer: Parser Writing Peril XXXIX

Is n a Perfect Power?

That means, is n = x^y with x,y \in \mathbb{N} . We don’t have to check all of the possible candidates, we can do a very large amount of argument reduction.
First observe that the greatest exponent y possible comes with smallest base y possible. The smallest base is 2, of course, so the largest possible exponent is y_{\text{max}} = \lfloor \log_2 x \rfloor.

By the Fundamental Theorem of Arithmetic together with the laws of exponentiation x^{a^b} = x^{ab} we can write the exponent ab as a product of prime powers ab = \prod_{i=1} p_i^{a_i} and with x^{p_1p_2\dotsm p_n} = x^{p_1}x^{p_2}\dotsm x^{p_n} and x^y = p_1^yp_2^y\dotsm p_n^y we can reduce the set of possible exponents to the set of primes; 16 in the case of the 53 bit long integers in JavaScript.
Can we reduce it further? Yes, by observing the maximum exponent for each base given a finite maximum for the number (that is, not for arbitrary large integers aka bigints simply because they have no maximum): the smaller the base the higher the maximum and vice versa. Lets take a look:

Base Factors
2 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53
3 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31
5 2, 3, 5, 7, 11, 13, 17, 19
7 2, 3, 5, 7, 11, 13, 17
11 2, 3, 5, 7, 11
13 2, 3, 5, 7, 11
17 2, 3, 5, 7, 11
19 2, 3, 5, 7
181 2, 3, 5, 7
191 2, 3, 5
1549 2, 3, 5
1553 2, 3
208057 2, 3
208067 2
< 94906266 2

The number of primes in the range is quite large \pi \Bigl( \big\lceil 2^\frac{53}{2} \big\rceil \Bigr) = 5\,484\,598 so just testing \log_b n \in \mathbb{N} won’t make much sense for more than a small handful of bases, although it is mathematically elegant.

The next simple and general test would be \sqrt[n]{x}\in \mathbb{N} . The problem here is the restricted precision of JavaScript and the quite simple implementations of Math.pow() which doesn’t care for integers. This leads to results like Math.pow(8916100448256,1/3) = 20735.999999999993 (12^{12} = 8\,916\,100\,448\,256). The usual floating point tests for equality work only for tests very close to zero.
But we can round the result of r = \sqrt[n]{x} if we test \lfloor r \rceil^n  \in \mathbb{N} . That works because the exponentiation by n \ge 2 makes the error two times larger and hence larger than the machine epsilon \epsilon_m and a check \big| x - \lfloor r \rceil^n \big| \le \epsilon_m would suffice but we can check against zero now because of the rounding..
We can even skip the test for base two because it is quite easy to check if a base 2 number is a power of two by just counting the bits: if the count comes up as zero or one it is a power of two. That makes three the smallest base which even triples the error from extracting the root.

Alternatively we could write an optimized integral root function but that is not necessary for JavaScript doubles as seen above. It will be necessary for Bigints, though.

Let’s sum it up for a positive integer x smaller than 53 bits:

  1. p_{max} \gets \log_2 x
  2. \text{if } x = 2^n \land n\in\mathbb{N}\text{ then } S_r \gets  \{2,n\}\text{and return true}
  3. sqrt is optimized in most cases, so try it first
  4. r \gets \lfloor\sqrt{x}\,\rceil
  5. \text{if } x - r^2 = 0 \text{ then } S_r \gets \{r,2 \}\text{ and return true }
  6. \text{while } p_i \le p_{\text{max}} \text{ do}
    1. r \gets \lfloor x^{\frac{1}{p_i}}\rceil
    2. \text{if } x - r^{p_i} = 0 \text{ then } S_r \gets \{r,p_i \}\text{ and return true }
  7. S_r \gets \{-1,-1 \}\text{ and return false}

Further optimization could be to tak the first values from the large table, make a small table out of it (say, 3–17) and check if it divides the input evenly e.g.: \gcd\left(x, 3^{23}\right) \not= 1.

D. J. Bernstein described in his paper “Detecting perfect powers in essentially linear time”[1] three different “levels” of such algorithms:

A perfect-power detection algorithm is an algorithm that, given an integer n > 1, figures out whether n is a perfect power. A perfect-power decomposition algorithm does more: if n is a perfect power, it finds x and k > 1 with n = x^k . A perfect-power classification algorithm does everything one could expect: it writes n in the form x^k with k maximal.
[emphasis added by D. J. Bernstein]

As we already computed the root and have the exponent, although not the largest, we can do a perfect-power decomposition algorithm easily.

Number.prototype.isPerfPow = function(nroot){
  // check if it is an positive integer omitted
  var primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53];
  var root,log2,t;
  t = Math.abs(this);
  log2 = Math.floor(Math.log(t)/Math.LN2);
  // reduces the number of tries to 11. at most
  if(this.isPow2()) return true;
  // sqrt is optimized in most cases, so try it first
  root = Math.round(Math.sqrt(t));
  if(t - (root * root) == 0){
    if(arguments.length  > 0){
      nroot[0] = root;
      nroot[1] = 2;
    }
    return true
  }
  // cube root is only in the next ECMAScript standard, sadly
  for(var i=1;primes[i] <= log2;i++){
    root = Math.round(Math.pow(t,1/primes[i]));
    if(t - (Math.pow(root,primes[i])) == 0){
      if(arguments.length  > 0){
        nroot[0] = Math.round(root);
        nroot[1] = primes[i];
      }
      return true
    }
  }
  return false;
};

To get the numbers you have to pass a two element Array to the function, but I think that is obvious 😉

To make it a perfect-power classification algorithm we need to recursively pass the root to the number until it fails and add the results up. With our example from above (12^12 = 8\,916\,100\,448\,256 = 20736^3 ) we get

  1. 20736^{\frac{1}{2}} = 144
  2. 144^{\frac{1}{2}} = 12
  3. 12^{\frac{1}{5}} \approx 1.64375

Multiplying all together 3 \times 2 \times 2 = 12 gives the highest possible exponent.
Another example (albeit base two gets done with a different algorithm)

  1. \sqrt{2^{52}} = 67108864
  2. \sqrt{67108864} = 8192
  3. 8192^\frac{1}{13} = 2

Building the product gives 2 \times 2 \times 13 = 52

The maximum number of recursion is with base two, as expected:

  1. \sqrt{2^{32}} = 65536
  2. \sqrt{65536} = 256
  3. \sqrt{256} = 16
  4. \sqrt{16} = 4
  5. \sqrt{4} = 2

I’ll leave the implementation for that to my dear readers 😉

[1]Bernstein, Daniel. “Detecting perfect powers in essentially linear time.” Mathematics of Computation of the American Mathematical Society 67.223 (1998): 1253-1283.

Advertisements

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