# 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)<Leading1(B): return -1 # remove leading zeros and make sure len(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()