Can integer division of numer power be optimized?

Question

In Python there is a builtin function pow which optimizes calculation of a**b%c. Why isn't there a function that calculates a**b//c?

Answer 1

I'll try to conclude here why calculation of a ** b // c can't be optimized, unlikely to a ** b % c.

Prerequisites: calculation of a ** b % c

Let's take an illustrative example: say a=2 and b=11. If we assume it's quite slow, then we are able to deduce that b = 1 + 2 + 8 = 2**0 + 2**1 + 0*2**2 + 2**3. After that this deduction can be used as a rule for multiplication of results a, a**2, a**4, a**8. Each result is assigned after squaring a previous one. Finally, a**11 = a*(a**2)*(a**8) and this process requires only 3 squarings.

If we generalize this process, this can be done like so:

a, b, r = 2, 11 , []

while b>0:
    if b % 2: r.append(a)
    b = b//2
    a = a*a
else:
    if b % 2: r.append(a)

print(r)

An output is r=[2, 4, 256]. Next, we need to multiply these multipliers. It can be done using from functools import reduce with a command reduce(lambda x,y: x*y, r).

Finally, if multipliers gets quite huge, multiplication get's very slow, so we need to replace each multiplier m with its modulus m%c and do the same in reduce function. Finally, we have:

from functools import reduce
def pow(a, b, c):
    # returns a ** b % c
    r = []
    while b > 0:
        if b % 2: r.append(a)
        b = b // 2
        a = (a * a) % c
    else:
        if b % 2: r.append(a)

    return reduce(lambda x, y: (x * y) % c, r)

Output is 4 since 2 ** 11 % 7 is 4.

I've tested also result 2046457 ** 1103207 % 71872 on my computer. An output was 18249 and it took 9 seconds to calculate while pow(2046457, 1103207, 71872) gave the same result instantly.

Update: plugging a ** b // c into calculation

Following abovementioned ideas, I'll try to achieve similar optimisation for calculation of a**b // c. I assume that squaring process remains the same and the main difference here is that we need to consider both integral and residual parts while taking squares (previous was easy because integral part was not important). If x is an integral part and y is residual one, we have a relation:

We also need to introduce similar calculation for two distinct multipliers:

My script looks like this now:

from functools import reduce
def pow(a, b, c):
    #returns tuple (a ** b // c,  a ** b % c)
    print(f'calculating: {a}**{b} = ', end='')
    r = []
    ir = (a//c, a%c) # we keep integral and residual part of a instead of a
    while b > 0:
        if b % 2: r.append(ir)
        b = b // 2
        ir = (ir[0]*ir[0]*c + 2*ir[0]*ir[1]+ (ir[1]*ir[1])//c, (ir[1]*ir[1]) % c)
    else:
        if b % 2: r.append(ir)
    out = reduce(lambda x, y: (c*x[0]*y[0] + x[0]*y[1] + x[1]*y[0] + (x[1] * y[1])//c, (x[1] * y[1]) % c), [(2, 2)]+[r[-1]])

    print(' * '.join(str(n[0]*c+n[1]) for n in r), end=' = ')
    print(' * '.join(str(n) for n in r),'=', out)
    return out

pow(2,7,3)

Output

calculating: 2**7 = 2 * 4 * 16 = (0, 2) * (1, 1) * (5, 1) = (42, 2)

Notes

Why it's not optimization still? We can see that second term in each factor remains always small but this is not a rule for the first terms like in this example of pow(26,31,35):

calculating: 26**31 = 26 * 676 * 456976 * 208827064576 * 43608742899428874059776 = 
(0, 26) * (19, 11) * (13056, 16) * (5966487559, 11) * (1245964082840824973136, 16) = 
(89709413964539398065824, 32)

In this case we can't avoid exponential growth of a%b // c. This is why no existent builtin function for a%b // c seems reasonable for me.