60

I have to write a program to calculate a**b % c where b and c are both very large numbers. If I just use a**b % c, it's really slow. Then I found that the built-in function pow() can do this really fast by calling pow(a, b, c).
I'm curious to know how does Python implement this? Or where could I find the source code file that implement this function?

2
  • 4
    The cpython source repo is at hg.python.org/cpython
    – Wooble
    Mar 9, 2011 at 14:07
  • 2
    ...under Objects/longobject.c:long_pow() (as JimB had already commented).
    – smci
    Jul 19, 2011 at 20:45

6 Answers 6

49

If a, b and c are integers, the implementation can be made more efficient by binary exponentiation and reducing modulo c in each step, including the first one (i.e. reducing a modulo c before you even start). This is what the implementation of long_pow() does indeed. The function has over two hundred lines of code, as it has to deal with reference counting, and it handles negative exponents and a whole bunch of special cases.

At its core, the idea of the algorithm is rather simple, though. Let's say we want to compute a ** b for positive integers a and b, and b has the binary digits b_i. Then we can write b as

b = b_0 + b1 * 2 + b2 * 2**2 + ... + b_k ** 2**k

ans a ** b as

a ** b = a**b0 * (a**2)**b1 * (a**2**2)**b2 * ... * (a**2**k)**b_k

Each factor in this product is of the form (a**2**i)**b_i. If b_i is zero, we can simply omit the factor. If b_i is 1, the factor is equal to a**2**i, and these powers can be computed for all i by repeatedly squaring a. Overall, we need to square and multiply k times, where k is the number of binary digits of b.

As mentioned above, for pow(a, b, c) we can reduce modulo c in each step, both after squaring and after multiplying.

4
  • 1
    Why can we reduce by modulo c in each step? Aug 30, 2015 at 22:26
  • 2
    @BenSandler: Because aa' (mod c) and bb' (mod c) imply aba'b' (mod c), or in other words, it doesn't matter whether you first reduce a and b modulo c and then multiply them, or multiply them first and then reduce modulo c. See the Wikipedia article on modular arithmetic. Aug 31, 2015 at 12:53
  • Note that long_pow is now defined at another line in that file: github.com/python/cpython/blob/master/Objects/…
    – JohanC
    Dec 8, 2019 at 21:57
  • @JohanC I've updated the link to include the commit hash, so it doesn't get out of date anymore. Dec 9, 2019 at 8:38
39

You might consider the following two implementations for computing (x ** y) % z quickly.

In Python:

def pow_mod(x, y, z):
    "Calculate (x ** y) % z efficiently."
    number = 1
    while y:
        if y & 1:
            number = number * x % z
        y >>= 1
        x = x * x % z
    return number

In C:

#include <stdio.h>

unsigned long pow_mod(unsigned short x, unsigned long y, unsigned short z)
{
    unsigned long number = 1;
    while (y)
    {
        if (y & 1)
            number = number * x % z;
        y >>= 1;
        x = (unsigned long)x * x % z;
    }
    return number;
}

int main()
{
    printf("%d\n", pow_mod(63437, 3935969939, 20628));
    return 0;
}
5
  • @Noctis, I tried running your Python implementation and got this:TypeError: ufunc 'bitwise_and' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe'' ---- As I'm learning Python right now, I thought you might have an idea about this error (a search suggests it might be a bug but I'm thinking that there's a quick workaround)
    – stackuser
    May 7, 2013 at 2:24
  • @stackuser: It appears to be working fine in the following demonstration: ideone.com/sYzqZN May 7, 2013 at 12:50
  • 5
    Can anyone explain why this solution works? I am having trouble understanding the logic behind this algorithm.
    – kilojoules
    May 24, 2015 at 3:59
  • 2
    @NoctisSkytower, what would be the benefit of this considering the native python pow() builtin function supports this as well and seems faster? >>> st_pow = 'pow(65537L, 767587L, 14971787L) >>> st_pow_mod = 'pow_mod(65537L, 767587L, 14971787L)' >>> timeit.timeit(st_pow) 4.510787010192871 >>> timeit.timeit(st_pow_mod, def_pow_mod) 10.135776996612549
    – Fabiano
    Jun 8, 2016 at 3:25
  • 7
    @Fabiano My function is not supposed to be used. It is simply an explanation of how Python works behinds the scenes without referring to its source in C. I was trying to answer wong2's question about how pow was implimented. Jun 8, 2016 at 13:18
1

I don't know about python, but if you need fast powers, you can use exponentiation by squaring:

http://en.wikipedia.org/wiki/Exponentiation_by_squaring

It's a simple recursive method that uses the commutative property of exponents.

0

Line 1426 of this file shows the Python code that implements math.pow, but basically it boils down to it calling the standard C library which probably has a highly optimized version of that function.

Python can be quite slow for intensive number-crunching, but Psyco can give you a quite speed boost, it won't be as good as C code calling the standard library though.

1
  • 7
    math.pow() does't have the modulo argument, and isn't the same function as the builtin pow(). Also FYI, Psyco is getting pretty stale, and no 64-bit support. NumPy is great for serious math.
    – JimB
    Mar 9, 2011 at 14:30
0

Python uses C math libraries for general cases and its own logic for some of its concepts (such as infinity).

0

Implement pow(x,n) in Python

def myPow(x, n):
        p = 1
        if n<0:
            x = 1/x
            n = abs(n)

        # Exponentiation by Squaring

        while n:
            if n%2:
                p*= x
            x*=x
            n//=2
        return p

Implement pow(x,n,m) in Python

def myPow(x,n,m):
            p = 1
            if n<0:
                x = 1/x
                n = abs(n)
            while n:
                if n%2:
                    p*= x%m
                x*=x%m
                n//=2
            return p

Checkout this link for explanation

Your Answer

Reminder: Answers generated by Artificial Intelligence tools are not allowed on Stack Overflow. Learn more

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.