Early termination of fractional exponent calculation?

Question

I need to write a function that takes the sixth root of something (equivalently, raises something to the 1/6 power), and checks if the answer is an integer. I want this function to be as fast and as optimized as possible, and since this function needs to run a lot, I'm thinking it might be best to not have to calculate the whole root.

How would I write a function (language agnostic, although Python/C/C++ preferred) that returns False (or 0 or something equivalent) before having to compute the entirety of the sixth root? For instance, if I was taking the 6th root of 65, then my function should, upon realizing that that the result is not an int, stop calculating and return False, instead of first computing that the 6th of 65 is 2.00517474515, then checking if 2.00517474515 is an int, and finally returning False.

Of course, I'm asking this question under the impression that it is faster to do the early termination thing than the complete computation, using something like

print(isinstance(num**(1/6), int))

Any help or ideas would be greatly appreciated. I would also be interested in answers that are generalizable to lots of fractional powers, not just x^(1/6).

Answer 1

Here are some ideas of things you can try that might help eliminate non-sixth-powers quickly. For actual sixth powers, you'll still end up eventually needing to compute the sixth root.

Check small cases

If the numbers you're given have a reasonable probability of being small (less than 12 digits, say), you could build a table of small cases and check against that. There are only 100 sixth powers smaller than 10**12. If your inputs will always be larger, then there's little value in this test, but it's still a very cheap test to make.

Eliminate small primes

Any small prime factor must appear with an exponent that's a multiple of 6. To avoid too many trial divisions, you can bundle up some of the small factors. For example, 2 * 3 * 5 * 7 * 11 * 13 * 17 * 19 * 23 = 223092870, which is small enough to fit in single 30-bit limb in Python, so a single modulo operation with that modulus should be fast.

So given a test number n, compute g = gcd(n, 223092870), and if the result is not 1, check that n is exactly divisible by g ** 6. If not, n is not a sixth power, and you're done. If n is exactly divisible by g**6, repeat with n // g**6.

Check the value modulo 124488 (for example)

If you carried out the previous step, then at this point you have a value that's not divisible by any prime smaller than 25. Now you can do a modulus test with a carefully chosen modulus: for example, any sixth power that's relatively prime to 124488 = 8 * 9 * 7 * 13 * 19 is congruent to one of the six values [1, 15625, 19657, 28729, 48385, 111385] modulo 124488. There are larger moduli that could be used, at the expense of having to check more possible residues.

Check whether it's a square

Any sixth power must be a square. Since Python (at least, Python >= 3.8) has a built-in integer square root function that's reasonably fast, it's efficient to check whether the value is a square before going for computing a full sixth root. (And if it is a square and you've already computed the square root, now you only need to extract a cube root rather than a sixth root.)

Use floating-point arithmetic

If the input is not too large, say 90 digits or smaller, and it's a sixth power then floating-point arithmetic has a reasonable chance of finding the sixth root exactly. However, Python makes no guarantees about the accuracy of a power operation, so it's worth making some additional checks to make sure that the result is within the expected range. For larger inputs, there's less chance of floating-point arithmetic getting the right result. The sixth root of (2**53 + 1)**6 is not exactly representable as a Python float (making the reasonable assumption that Python's float type matches the IEEE 754 binary64 format), and once n gets past 308 digits or so it's too large to fit into a float anyway.

Use integer arithmetic

Once you've exhausted all the cheap tricks, you're left with little choice but to compute the floor of the sixth root, then compare the sixth power of that with the original number.

Here's some Python code that puts together all of the tricks listed above. You should do your own timings targeting your particular use-case, and choose which tricks are worth keeping and which should be adjusted or thrown out. The order of the tricks will also be significant.

from math import gcd, isqrt

# Sixth powers smaller than 10**12.
SMALL_SIXTH_POWERS = {n**6 for n in range(100)}

def is_sixth_power(n):
    """
    Determine whether a positive integer n is a sixth power.

    Returns True if n is a sixth power, and False otherwise.
    """
    # Sanity check (redundant with the small cases check)
    if n <= 0:
        return n == 0

    # Check small cases
    if n < 10**12:
        return n in SMALL_SIXTH_POWERS

    # Try a floating-point check if there's a realistic chance of it working
    if n < 10**90:
        s = round(n ** (1/6.))
        if n == s**6:
            return True
        elif (s - 1) ** 6 < n < (s + 1)**6:
            return False
        # No conclusive result; fall through to the next test.

    # Eliminate small primes
    while True:
        g = gcd(n, 223092870)
        if g == 1:
            break
        n, r = divmod(n, g**6)
        if r:
            return False

    # Check modulo small primes (requires that 
    # n is relatively prime to 124488)
    if n % 124488 not in {1, 15625, 19657, 28729, 48385, 111385}:
        return False

    # Find the square root using math.isqrt, throw out non-squares
    s = isqrt(n)
    if s**2 != n:
        return False

    # Compute the floor of the cube root of s
    # (which is the same as the floor of the sixth root of n).
    # Code stolen from https://stackoverflow.com/a/35276426/270986
    a = 1 << (s.bit_length() - 1) // 3 + 1
    while True:
        d = s//a**2
        if a <= d:
            return a**3 == s
        a = (2*a + d)//3