Any pointers on this troublesome algorithm?

[t3xc0d3]

Your algorithm is one form of fast binary exponentiation, known as exponentiation by squaring. It is a common algorithm for asymmetric cryptography since it is much more efficient than naive multiplication. It is, for instance, an important core of RSA and algorithms based on discrete logarithms. Variations do also exist for other domains, such as the double-and-add algorithm for elliptic curves.

For x ** e mod n, you calculate x * (x ** 2) ** ((n-1)/2) if n is odd (square and multiply), otherwise (x**2)**(n/2) (square).

To exemplify:

x=17
e=33
n=73

33 is the same as 0b100001.

1. square -> 17 * 17 = 70 mod 73 = x ** {0b10}
2. square -> 70 * 70 = 9 mod 73 = x ** {0b100}
3. square -> 9 * 9 = 8 mod 73 = x ** {0b1000}
4. square -> 8 * 8 = 64 mod 73 = x ** {0b10000}
5. square -> 64 * 64 = 8 mod 73 = x ** {0b100000}
multiply -> 8 * 17 = 63 mod 73 = x ** {0b100001}

In (not so great) python code it looks as

def exp(a,e,p):
        binary = bin(e).lstrip("0b")

        res=a
        for i in xrange(1, len(binary)):
                if binary[i] == "1":
                        res = (res ** 2) % p
                        res= (res * a) % p
                else:
                        res = (res ** 2) % p
        return res