modular exponentiation python program

This is a simple – not efficient – but doable way to do modular exponentiation

#!/usr/bin/python

def modexp(base, pow, mod):
  exponent = 1
  i = 0
  while i < pow:
    exponent = (exponent * base) % mod
    i += 1
  return exponent

if __name__ == "__main__":
  import sys
  print modexp(int(sys.argv[1]), int(sys.argv[2]), int(sys.argv[3]))

Example usage

$ modexponent.py 2 1234 789
481

For large numbers this method isn’t practical. The below program implements the improved python algorithms for modular exponents.

#!/usr/bin/python
#Let x be an integer written in binary, y and n be integers
#I'm sure this isn't the best, but this is only called once per program
def binconvert(n):
  barray = []
  if n < 0: 
    raise ValueError, "must be positive"
  if n == 0:
    return 0
  while n > 0:
    #barray = n%2 + barray[:]
    barray.append(n%2)
    n = n >> 1
  barray.reverse()
  return barray
#this is the intuitive but slow way
def modexp(base, pow, mod):
  exponent = 1
  i = 0
  while i < pow:
    exponent = (exponent * base) % mod  
    i += 1
  return exponent
#y**x mod n
def modexp1(y, x, n):
  #convert x to a binary list
  x = binconvert(x)
  
  s = [1]
  r = x[:]
  for k in range (0, len(x)):
    if x[k] == 1:
      r[k] = (s[k] * y) % n
    else:
      r[k] = s[k]
    s.append ((r[k]**2)%n)
  #print s
  #print r
  return r[-1]
#similar to modexp1 except uses the bits in reverse order
def modexp2(y, x, n):
  a = x
  b = 1
  c = y
  while a != 0:
    if a % 2 == 0:
      a = a/2
      c = (c**2) % n
    else:
      a = a -1
      b = (b * c) % n
  return b
if __name__ == "__main__":
import sys
  print modexp(int(sys.argv[1]), int(sys.argv[2]), int(sys.argv[3]))
  print modexp1(int(sys.argv[1]), int(sys.argv[2]), int(sys.argv[3]))
  print modexp2(int(sys.argv[1]), int(sys.argv[2]), int(sys.argv[3]))

Using the above program, some interesting things are easy to demonstrate. For example, here is a test that shows for 100,000 tries, the 2 algorithms get the correct answer.

import prob23
import random
import sys

BASEMAX = 60000
EXPMAX = 60000
MODMAX = 100000
ITERATIONS = 100000

print "beginning program"
print "This will take awhile... you need", ITERATIONS, "of these dots to complete."
p
for i in range (0, ITERATIONS):
  if i % 10 == 0:
    sys.stdout.write(".")
    sys.stdout.flush()
  r_base = random.randint(1, BASEMAX)
  r_exp = random.randint(1, EXPMAX)
  r_mod = random.randint(max(r_base, r_exp), MODMAX)
  trueval = prob23.modexp(r_base, r_exp, r_mod)
  testval1 = prob23.modexp1(r_base, r_exp, r_mod)
  testval2 = prob23.modexp2(r_base, r_exp, r_mod)
  if trueval != testval1 or trueval != testval2:
    print "One of these exponent things doesn't work"
    #sys.exit(1)
print "Done, all of these answers were the same."

Running this program several times with various parameters verifies the correctness of the algorithm (like how an Engineer verifies things).

Also, it’s interesting to time these algorithms when compared to the brute force method (multiplying every number). All three algorithms were timed to find a solution. Below is some example code that accomplishes this (this particular code is for the brute force version, but it is almost identical to the other algorithms – the only change is the modexponent call).

import modexponent

for i in range(1, 65538):
  thistry = modexponent.modexp(3, i, 65537)
  if thistry == 2:
    print i
    break
  if i % 15000 == 0:
    print "trying", i

Running all three algorithms, the results are shown below.

me@Opteron:~/math/chp3$ time python bruteforce.py
55296
real 11m28.131s
user 11m18.710s
sys 0m1.170s
me@Opteron:~/math/chp3$ time python algorithm1.py
55296
real 0m3.285s
user 0m3.290s
sys 0m0.000s
me@Opteron:~/math/chp3$ time python algorithm2.py
55296
real 0m2.621s
user 0m2.610s
sys 0m0.010s

This is a huge speedup – from 11 minutes in the case of brute forcing to under 4 seconds for the two algorithms.