"""
Compute powers b^e mod m using the
method of repeated squaring.  Two
implementations.  One using a loop,
one using recursion.

Apply this to Fermat's test for composite
numbers, and to the generation of "industrial-grade"
primes.
"""

######################################################
#
# Powers mod m

def modpower(b,e,m):
   """ modpower(b,e,m) = b^e mod m """
   P = 1
   S = b
   while e > 0:
      r = e%2
      e = e/2
      if r == 1:
         P = P*S % m
      S = S*S % m
   return P

def modpower2(b,e,m):
   """ modpower(b,e,m) = b^e mod m, recursive implementation """
   P = modpower(b,e/2,m)
   if e % 2 == 0:
      return (P*P) % m
   else:
      return (P*P*b) % m

#######################################################
#
# The Fermat test, and an application to
# generating industrial-grade primes

def ft(n):
   """   Fermat test   """
   if modpower(2,n-1,n) > 1:
      return 0
   else:
      return 1

def fpp(n):
   """ fpp(n) == first probable prime >= n """
   p = n
   if p == 2:
     return p
   if p % 2 == 0:
     p = p + 1
   while ft(p) == 0:
      p = p + 2
   return p

"""
10^1000 + 453 is an industrial-grade prime
"""