Determine powers of a prime and other helper functions

The first steps of the AKS primality are used to detect prime-powers with exponent bigger than one. This step can be done in many ways, and almost any algorithm is efficient enough to be used into the AKS test.
The algorithm used here wont impact the time of the AKS test on primes, since all the time is taken in the last loop that checks the polynomial congruence.

In my implementation of the AKS test I used the simplest algorithm I could think of. The number $n$ can be, at most, a power of exponent $\log_2{n}$, thus we can loop over the exponents from $2$ to $\log_2{n}$ and for every exponent $e$ do a bisection search to find the base. If there exist an integer base $b$ such that $b^e = n$ then $n$ is a power of a prime, hence not a prime. Otherwise, if no such $b$ exist for every exponent, the number is not a power of a prime.

The resultant code is this:

import math


def find_int_power(n):
    """Return two integers m > 0, k > 1 such that n = m^k. (*n* > 1)

    If two such integers do not exist return None.
    """
    
    if n < 2:
        raise ValueError('This function requires n > 1.')
    
    h = int(math.log(n, 2)) + 1
    
    for k in xrange(2, h + 1):
        a = 0
        b = n
        
        while b - a > 1:
            m = (a + b) // 2
            
            res = m ** k - n
            
            if not res:
                return m, k
            elif res > 0:
                b = m
            else:
                a = m
        
        if m ** k == n:
            return m, k
    
    return None  

Most versions of the AKS test require to compute Euler's $\varphi$ function on the exponent of the polynomial. This again can be done in the naive way, without impacting much the algorithm efficiency:

def factorize(n):
    """Yield the factors of *n* with the given multiplicity."""
    
    exp = 0
    while not n & 1:
        exp += 1
        n //= 2
    
    if exp:
        yield (2, exp)
    
    for div in xrange(3, int(n**.5 + 1.5), 2):
        exp = 0
        while not n % div:
            exp += 1
            n //= div
        if exp:
            yield (div,  exp)
    
    if n > 1:
        yield (n, 1)


def phi(n):
    """Calculates the euler's function at value *n*."""
    
    tot = 1
    for div, exp in factorize(n):
        tot *= (div**(exp - 1)) * (div - 1)
    return tot

Later, in an other post,  we will see that, the timings of these functions do not affect at all the total running time of the algorithm.