Implement the function fast modular exponentiation

Question

I am trying to implement the function fast modular exponentiation(b, k, m) which computes:
b^(2k) mod m using only around 2k modular multiplications.

I tried this method:

def FastModularExponentiation(b, k, m):
    res = 1
    b = b % m
    while (k > 0):
        if ((k & 1) == 1):
            res = (res * b) % m
        k = k >> 1
        b = (b * b) % m
    return res

but I am still stuck in same problem which is if I try b = 2, k = 1, m = 10, my code returns 22. However, the correct answer is:

2^(2^1) mod 10 = 2^2 mod 10 = 4

and I cannot find the reason why.

Answer 1

Update: I finally understood that you do not want regular modular exponentiation (i.e., b^k mod m), but b^(2^k) mod m (as you plainly stated).

Using the regular built-in Python function pow this would be:

def FastModularExponentiation(b, k, m):
    return pow(b, pow(2, k), m)

Or, without using pow:

def FastModularExponentiation(b, k, m):
    b %= m
    for _ in range(k):
        b = b ** 2 % m
    return b

If you know r = phi(m) (Euler's totient function), you could reduce the exponent first: exp = pow(2, k, r) and then calculate pow(b, exp, m). Depending on the input values, this might speed things up.

---

(This was the original answer when I thought you wanted, b^k mod m)

This is what works for me:

def fast_mod_exp(b, exp, m):
    res = 1
    while exp > 1:
        if exp & 1:
            res = (res * b) % m
        b = b ** 2 % m
        exp >>= 1
    return (b * res) % m

The only significant differences I spot is in the last line: return (b * res) % m and that my while loop terminates earlier: while exp > 1 (which should be the same thing you do - except it saves an unnecessary squaring operation).

---

Also note that the built-in function pow will do all that for free (if you supply a third argument):

pow(4, 13, 497)
# 445