Commits

David Wilhelm  committed 179fac3

p97: implement recursive formula with non-recursive loop

Compute base**exponent by using the base 2 representation of
exponent.

To get get the necessary sequence, work backward from base**exponent
and note which of two steps is used:

Let x = exponent and let f be an empty bitstring that will hold a
flag for each step.

While x > 0:
If x is even, set x = x / 2 and prepend 0 to f. This is equivalent to
taking the square root of base**x.
If x is odd, set x = floor(x / 2) and prepend 1 to f. This is
equivalent to diving base**x by base and then taking the square root.

Then to create base**exponent, start with 1 and follow the steps
described by f reading left to right (opposite direction than written):

Let p = 1. Looping over the bits in f:
If 0, set p = p * p
If 1, set p = p * p * base

Note that the bitstring f equals the binary representation of exponent:
Squaring base**exponent is the same as doubling exponent, which is done
by bitshifting exponent left by 1, filling with 0.
Squaring base**exponent and multiplying by base is the same as doubling
exponent then adding 1, or bitshifting exponent left by 1 and filling
with 1.

  • Participants
  • Parent commits 593eb25

Comments (0)

Files changed (3)

File problem_97.py

 Find the last ten digits of this prime number.
 """
 
+import pdb
+
+from util import ops
+
+
 def main():
-    num = 28433
-    for x in xrange(7830457):
-        num <<= 1
-        num %= 10**10
-    print num + 1
+    base = 2
+    exponent = 7830457
+    modulus = 10**10
+    num = ops.exp(base, exponent, modulus)
+    print '%d ** %d == %d mod %d' % (base, exponent, num, modulus)
+    print (num * 28433 + 1) % modulus
 
 
 if __name__ == '__main__':

File util/__init__.py

 Miscellaneous utilities for Project Euler problems
 """
 
-__all__ = ['base', 'comb', 'divisors', 'primes']
+__all__ = ['base', 'comb', 'divisors', 'ops', 'primes']
+"""
+Operations on mathematical objects
+"""
+
+def exp(base, exponent, modulus):
+    """Return base**exponent as a number.
+
+    If modulus is given, then the answer and all intermediate steps are
+    mod modulus. This can greatly speed calculation.
+    """
+    trackexp = 1
+    while trackexp < exponent:
+        trackexp <<= 1
+    print 'start: trackexp=%d, exponent=%d' % (trackexp, exponent)
+    num = 1
+    while trackexp:
+        if trackexp & exponent:
+            num = num * num * base
+        else:
+            num = num * num
+        if modulus:
+            num %= modulus
+        trackexp >>= 1
+        print 'trackexp=%d, num=%d' % (trackexp, num)
+    return num
+