# Alfonso Cevallos. March 17 2015.
# Please report any mistakes.
# Rules of the game: do not use Python integers, only bit lists!
# We start by defining the same functions as last time.
# Some were modified to be less readable, but more efficient.
# Most functions use ITERATION, as opposed to the RECURSION seen in class. Compare them!
# The "Leading1" function now also removes any leading zeros from the bit list.
import math
import os
import random
import time
def Leading1(A): #Traverse the list starting from most significant bit, until we find the first 1.
for i in xrange(len(A) - 1, -1, -1):
if A[i] != 0: return i
else: A.pop() #Remove unnecessary zeros on the fly
return -1
def Intify(A): # Convert a bit list into a Python integer
s = 0 # We use this only for checking correctness of our algorithms
for b in reversed(A):
s = 2*s + b
return s
def Rand(n): # Generate a random bit list with n entries, the most significant bit being 1
A = [] #Here n is assumed to be an integer, not a list.
for i in range(n - 1):
A.append(random.randint(0, 1))
A.append(1)
return A
def Add(A, B): # Add two bit lists
if Leading1(A) < Leading1(B): A,B = B,A #swaps the pointers, to ensure len(A)>=len(B)
Out = []
carry = 0
for i in range(len(A)):
if i < len(B): h = carry + A[i] + B[i]
else: h = carry + A[i]
Out.append(h % 2)
if h > 1: carry = 1
else: carry = 0
if carry == 1: Out.append(1)
return Out
def Sub(A, B): # Compute A-B
if Leading1(A)=len(B), else report error
Out = []
carry = 0
for i in range(len(A)):
if i < len(B): h = A[i] - B[i] - carry
else: h = A[i] - carry
Out.append(h % 2)
if h < 0: carry = 1
else: carry = 0
if carry == 1: return -1 # if the last carry is 1, L2 was bigger than L1. Report an error
Leading1(Out) # this line removes any leading zeros
return Out
def SimpleMult(A, B): # Standard multiplication A * B
if Leading1(A) < Leading1(B): A,B = B,A # remove leading zeros and ensure len(L1)>=len(L2)
Out = []
for b in reversed(B):
Out.insert(0,0)
if b==1: Out = Add(Out,A)
return Out
def Mult(A, B):
# Remark: In this assignment we use this hybrid multiplication algorithm by default.
# It is the most efficient multiplication algorithm (see the first assignment).
# We compute A*B using the Karatsuba algorithm; except if input is smaller than the limit.
limit = 60
if Leading1(A) < Leading1(B): A,B = B,A # swap the pointers, to ensure len(a)>=len(b)
if len(A) < limit: return SimpleMult(A,B) #Use simple multiplication
k = len(A) // 2 # split the lists at this point
A0 = A[:k]
A1 = A[k:]
B0 = B[:k]
B1 = B[k:] # we have A*B = (A1*2^k + A0) * (B1*2^k + B0)
S2 = Mult(A1, B1)
S0 = Mult(A0, B0)
SS = Mult(Add(A1, A0), Add(B1, B0))
S1 = Sub(SS, Add(S2, S0)) # we have A*B = (S2*2^k + S1)*2^k + S0
return Add(([0]*k + Add(([0]*k + S2), S1)), S0)
# Exercise 3: We define division with remainder first, and then fast modular exponentiation (FME).
def Div(A, B): #Computes division with remainder. Outputs Q,R, with A = BQ + R, R < B.
if Leading1(B) == -1: return -1 # If B = 0, report an error.
Q,R = [],[]
for b in reversed(A):
R.insert(0,b)
difference = -1 #if R < B, difference is -1 (code for error), else difference = R - B
if len(R) >= len(B): difference = Sub(R,B)
if difference != -1: # Case when R >= B
R = difference
Q.insert(0,1)
else: # Case when R is still smaller than B
if Q!=[]: Q.insert(0,0)
return Q,R
def FME(A, E, N): # Computes A^E mod N, using fast modular exponentiation.
# Starting from Out = 1, we just need to perform two operations:
# squaring the number, and multiplying it by A.
# The order at which we perform these operations is encoded in the bits of E.
# We also reduce mod N after each iteration.
Out = [1]
for b in reversed(E):
Out = Mult(Out,Out) # Square it
if b==1: Out = Mult(Out, A) # Multiply it by A
Q, Out = Div(Out, N) # Reduction mod N
return Out
def ExampleDivisionAndFME(): #some examples of division with remainder and fast mod. exp.
A = [1, 1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1]
B = [0, 0, 1, 0, 0, 0, 0, 0, 1]
Q,R = Div(A, B)
print "Example of division with remainder \n"
print "%d = %d * %d + %d \n" % (Intify(A), Intify(B), Intify(Q), Intify(R))
N = [1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1]
A = [1, 0, 1, 0, 1, 0, 1]
E = [1, 1, 0, 1]
F = FME(A, E, N)
print "Example of fast modular exponentiation \n"
print "%d ^ %d mod %d = %d \n" % (Intify(A), Intify(E), Intify(N), Intify(F))
# Exercise 5: Fibonacci functions.
# Remark. We treat the input N as a bit list.
def Fib1(N): # Computes the N-th fibonacci number using the trivial algorithm
if Leading1(N) < 1: return N
return Add(Fib1(Sub(N,[1])), Fib1(Sub(N,[0,1])))
def MatrixMult(A, B, C, D, E, F, G, H):
# This matrix multiplication computes (P,Q,R,S) where
# |P Q| = |A B| * |E F|
# |R S| |C D| |G H|
# Notice we perform 8 integer multiplications. Can this be improved?
P = Add(Mult(A,E),Mult(B,G))
Q = Add(Mult(A,F),Mult(B,H))
R = Add(Mult(C,E),Mult(D,G))
S = Add(Mult(C,F),Mult(D,H))
return P,Q,R,S
def Fib2a(N): #Here we do slow exponentiation
P,Q,R,S = [1],[0],[0],[1] #Is a power of matrix A. Start with identity matrix
for i in range(Intify(N)):
P,Q,R,S = MatrixMult(P,Q,R,S,[1],[1],[1],[0]) #Go from A^(i-1) to A^i
return Q
def Fib2b(N): #Here we do fast exponentiation. Similar to FME except for modular reduction
P,Q,R,S = [1],[0],[0],[1] #Is a power of matrix A. Start with identity matrix
for b in reversed(N):
P,Q,R,S = MatrixMult(P,Q,R,S,P,Q,R,S) #Square the matrix
if b==1: P,Q,R,S = MatrixMult(P,Q,R,S,[1],[1],[1],[0]) #Multiply by A
return Q
def Fib2c(N):
# Besides using fast exponentiation, here we optimize for the special case
# of the powers of matrix A. Any power of A is of the form
# |P+Q P|, and its square is |R+S R|
# |P Q| |R S|,
# where R = P*(P+2Q) and S = R-Q*(2P-Q) (check it!).
# Notice we only need to perform two integer multiplications.
# (Moreover, multiplying by A only requires one sum.)
# We gain a multiplicative constant factor compared to Fib2b.
P,Q = [],[1] # P,Q a power of matrix A. Start with identity matrix
for b in reversed(N):
if P != []: # Square the matrix (redundant if P,Q is identity)
R = Mult(P, Add(P, [0]+Q))
S = Sub(R, Mult(Q, Sub([0]+P, Q)))
P,Q = R,S
if b==1: P, Q = Add(P,Q), P # Multiply matrix by A
return P
def Test_Compare_Fib():
# Fibonacci tests: We compute Fibonacci numbers using the 4 algorithms
# Fib1 is way slower than the rest starting from N = 2^4
# For N = 2^5 Fib1 will crash.
print "Compare the 4 Fibonacci algorithms: \n"
for i in range(2,6):
N = Rand(i)
time0 = time.clock()
fib1 = Intify(Fib1(N))
time1 = time.clock()
fib2a = Intify(Fib2a(N))
time2 = time.clock()
fib2b = Intify(Fib2b(N))
time3 = time.clock()
fib2c = Intify(Fib2c(N))
time4 = time.clock()
if fib1 != fib2a or fib1 != fib2b or fib1 != fib2c: print "Error in the computation of Fibonacci numbers"
else: print "The %d-th Fibonacci number is %d" % (Intify(N), fib1)
print "Running time of Fib1 : %f" % (time1-time0)
print "Running time of Fib2a: %f" % (time2-time1)
print "Running time of Fib2b: %f" % (time3-time2)
print "Running time of Fib2c: %f \n" % (time4-time3)
def Test_Fib2():
# Now we compare the 3 Fib2 algorithms.
# Fib2a becomes significantly slower for N = 2^9
# Meanwhile, Fib2c seems to be around 4 times faster than Fib2b
print "Test for the 3 Fib2 algorithms: \n"
for i in range(6,10):
N = Rand(i)
time1 = time.clock()
fib2a = Intify(Fib2a(N))
time2 = time.clock()
fib2b = Intify(Fib2b(N))
time3 = time.clock()
fib2c = Intify(Fib2c(N))
time4 = time.clock()
if fib2a != fib2b or fib2a != fib2c: print "Error in the computation of Fibonacci numbers"
else: print "The %d-th Fibonacci number is %d" % (Intify(N), fib2a)
print "Running time of Fib2a: %f" % (time2-time1)
print "Running time of Fib2b: %f" % (time3-time2)
print "Running time of Fib2c: %f \n" % (time4-time3)
def Test_Fast_Fib():
# Finally we study the capacity of the fastest algorithm.
# It allows us to go up to N = 2^13 in reasonable time.
print "Test for the fast Fibonacci algorithm: \n"
for i in range(10,14):
N = Rand(i)
time3 = time.clock()
fib2c = Fib2c(N)
time4 = time.clock()
print "The %d-th Fibonacci number is %d" % (Intify(N), Intify(fib2c))
print "Running time of Fib2c: %f \n" % (time4-time3)
ExampleDivisionAndFME()
Test_Compare_Fib()
Test_Fib2()
Test_Fast_Fib()