Python-only example of python modular exponentiation by squaring, similar to built-in pow().

mesq.py

import sys


try:
    rng = xrange
except NameError:
    rng = range  # Python 3+


def mesq(base, exp, mod, debug=False):
    """Modular exponentiation by squaring."""
    if debug:
        # Print some debugging output. Don't use with huge numbers!
        print('({}**{})%{} ='.format(base, exp, mod))
        print('(({}%{})**{})%{} ='.format(base, mod, exp, mod))
        log, cexp = '', 1
        for x in reversed('{:b}'.format(exp)):
            if x == '1':
                log = '*((({}%{})**{})%{})'.format(base, mod, cexp, mod) + log
            cexp *= 2
        print('({})%{} ='.format(log[1:], mod))
        log, nexp = '', 0
        for x in reversed('{:b}'.format(exp)):
            if x == '1':
                args = '('*nexp*2, base, mod, '**2)%{})'.format(mod)*nexp
                log = '*{}({}%{}){}'.format(*args) + log
            nexp += 1
        print('({})%{} ='.format(log[1:], mod))
    # Loop until we have the result:
    result, counter = 1, 0
    cache = (counter, base % mod)
    while exp:
        leap = 1
        if exp & 1:
            current = cache[1]
            for x in rng(counter - cache[0]):
                current = pow(current, 2) % mod  # pow() is fast!
            cache = (counter, current)
            result = (result * current) % mod
        else:
            # In case of many zeros, this will give a small speedup:
            while not exp & ((1<<(leap*2))-1):
                leap *= 2
        counter += leap
        exp >>= leap
    return result


if __name__ == '__main__':
    print(mesq(*(int(x) for x in sys.argv[1:4]), debug=False))