Is the following divide and conquer recursive algorithm for the exponentiation more efficient than the iterative one for large numbers?

Question

I have the two following algorithms. My analysis says that both of them are \$\mathcal O(m^24^n)\$ i.e they are equivalent for big numbers (more than 32 bits). Is this right? Note that m and n are the bits numbers for x and y

def pow1(x, y):
    if y == 0:
        return 1
    temp = x
    while y > 1:
        y -= 1
        temp *= x
    return temp

def pow2(x, y):
    if y == 0:
        return 1
    temp = pow2(x, y//2)
    if y & 1:    return temp * temp * x
    return temp * temp

The complexity of the first

y-1 iterations and in each iteration, a subtraction taking \$\mathcal O (\lg (y-i))\$, and a multiplication taking \$\mathcal O (im\cdot m)\$, hence the total work takes

\$T=\mathcal O(\sum\limits_{i=1}^ym^2i)=\mathcal O\left(m^2\frac{y(y+1)}2\right)=\mathcal O\left(m^24^n\right)\$

Of the second

We have \$n=\lg y\$ calls and for each, we have multiplications taking \$\mathcal O (2^im\cdot 2^im)\$, hence the total work takes \$ T=\mathcal O\left(\sum\limits_{i=1}^{n}\left(4^im^2\right)\right)=\mathcal O\left({4^n}{m^2}\right)\tag{for large $n$}\$

Answer 1

I think your analysis is correct if multiplication of k-digit and ℓ-digit numbers takes Θ(kℓ) time, but more efficient multiplication algorithms are known. The most efficient algorithm in wide use (implemented in GMP) is Schönhage-Strassen, which is O(k log k log log k) for numbers of equal length. I don't know the complexity for unequal lengths, but I suspect it's O(ℓ log k log log k) for k < ℓ. Using that algorithm, or any sub-kℓ algorithm, you should find that the divide-and-conquer approach is faster.

Answer 2

As benrg answered, it looks correct if multiplication of k-digit and ℓ-digit numbers takes Θ(kℓ) time. We can also somewhat verify it experimentally. Let's write a number class that keeps track of the bit operations:

class Int(int):
    def __mul__(self, other):
        global bitops
        bitops += self.bit_length() * other.bit_length()
        return Int(int(self) * other)

And now let's test that a bit, first increasing n by 1:

m  n  pow1     pow2     pow1 / pow2        pow1 / prev_pow1
10 10 52272170 34951501 1.4955629516454816 None
10 11 209388450 139788522 1.4978944408611745 4.005734791572648
10 12 838148190 559136896 1.4990035463515539 4.002838695257546
10 13 3353781770 2236448811 1.4996014008925151 4.001418615483737
10 14 13417505370 8945532982 1.4999112291015417 4.0007091367784495

Here pow1 is the number of bit operations with pow1 and likewise for pow2. The pow1 / pow2 column shows that pow1 takes about a constant 1.5 times as many bit operations as pow2. And the last column shows that increasing n by 1 quadruples pow1 as predicted by your analysis saying \$O(4^nm^2)\$.

Now let's instead repeatedly double m:

m  n  pow1     pow2     pow1 / pow2        pow1 / prev_pow1
10 10 52272170 34951501 1.4955629516454816 None
20 10 209101200 139806021 1.495652322441821 4.000239515596923
40 10 836404800 559224041 1.4956524374459073 4.0
80 10 3345619200 2236896081 1.4956524929420716 4.0
160 10 13382476800 8947584161 1.4956525201886839 4.0

We see that pow1 and pow2 again differ only by constant factor 1.5 and that doubling m quadruples the bit operations as expected from \$O(4^nm^2)\$.

Whole code:

class Int(int):
    def __mul__(self, other):
        global bitops
        bitops += self.bit_length() * other.bit_length()
        return Int(int(self) * other)

def pow1(x, y):
    if y == 0:
        return Int(1)
    temp = x
    while y > 1:
        y -= 1
        temp *= x
    return temp

def pow2(x, y):
    if y == 0:
        return Int(1)
    temp = pow2(x, y//2)
    if y & 1:    return temp * temp * x
    return temp * temp

m = 10
n = 10
prev_bitops1 = None
for _ in range(5):
    x = Int(2**m - 1)
    y = 2**n - 1
    bitops = 0; pow1(x, y); bitops1 = bitops
    bitops = 0; pow2(x, y); bitops2 = bitops
    print(m, n,
          bitops1, bitops2,
          bitops1 / bitops2,
          prev_bitops1 and bitops1 / prev_bitops1)
    prev_bitops1 = bitops1
    # n += 1
    m *= 2