#!/usr/local/bin/maple
# -*- maplev -*-

# Nathaniel Shar
# HW 4 - 9 Feb 2014
# Experimental Mathematics

# It is okay to link to this assignment on the course webpage.

Help := proc(): print(`Tra(A), MatAdd(A,B), MatMult(A,B), PadZeros(A), MatMultR(A,B), Tra(A), RandMat(m,n,R), isSquare(A), submatrix(A, rows, cols), MatMultRG(A, B), Strassen(A, B), StrassenR(A, B), RandMat(n,R)`): end:

#########
# From c4.txt:

#C4.txt: Feb. 3 2014. Distance-learning class (due to University
# closing due to snow storm)


#Tra(A): the transpose of the matrix A
Tra:=proc(A) local i,j:
    [seq([seq(A[i][j],i=1..nops(A))],j=1..nops(A[1]))]:
end:


# Commented out due to conflict with function defined below.

##RandMat(m,n,R): a random m by n matrix with entries between -R and R
##given as a list of lists
#RandMat:=proc(m,n,R) local ra,i,j:
#    ra:=rand(-R..R):
#    [seq([seq(ra(),j=1..n)],i=1..m)]:
#end:


#MatAdd(A,B): adds two matrices of the same size given as lists of
# lists
MatAdd:=proc(A,B) local i,j:

    if nops(A)<>nops(B) then
        print(`Not the same number of rows`):
        RETURN(FAIL):
    fi:

    if nops({seq(nops(A[i]),i=1..nops(A))})<>1 then
        print(`Rows do not all the same lengths, hence`, A, `is not a matrix `):
        RETURN(FAIL):
    fi:


    if nops({seq(nops(B[i]),i=1..nops(B))})<>1 then
        print(`Rows do not all the same lengths, hence`, B, `is not a matrix `):
        RETURN(FAIL):
    fi:

    if nops(A)>0 and nops(A[1])<>nops(B[1]) then
        print(`Not the same number of columns`):
        RETURN(FAIL):
    fi:


    [seq([seq(A[i][j]+B[i][j],j=1..nops(A[i]))],i=1..nops(A))]

end:


#MatMult(A,B): multiplies two matrices A and B
MatMult:=proc(A,B) local i,j,k:


    if nops({seq(nops(A[i]),i=1..nops(A))})<>1 then
        print(`Rows do not all the same lengths, hence`, A, `is not a matrix `):
        RETURN(FAIL):
    fi:


    if nops({seq(nops(B[i]),i=1..nops(B))})<>1 then
        print(`Rows do not all the same lengths, hence`, B, `is not a matrix `):
        RETURN(FAIL):
    fi:

    if nops(A[1])<>nops(B) then
        print(`number of columns of`, A, `is not the same as the number of rows of`, B):
        print(`hence can't multiply them.`):
        RETURN(FAIL):
    fi:


    [seq([seq( add(A[i][k]*B[k][j],k=1..nops(A[i])),
               j=1..nops(B[1]))],i=1..nops(A))]:

end:

#PadZeros(A): given  a square matrix A pads it with zero to make its
#dimension a powr of 2
PadZeros:=proc(A) local i,n,n1:

    if nops(A)<>nops(A[1]) then
        print(`not a square matrix`):
        RETURN(FAIL):
    fi:

    n:=nops(A):

    if type(log[2](n),integer) then
        A:
    else
        n1:=2^ceil(log[2](n)):
        [seq([op(A[i]),0$(n1-n)],i=1..nops(A)),[0$n1]$(n1-n)]:
    fi:

end:


#MatMultR(A,B): multiplies two square matrices A and B whose dimension
#is a power of 2
MatMultR:=proc(A,B) local n,i,A11,A12,A21,A22,B11,B12,B21,B22,
    M1,M2,M3,M4,M5,M6,M7,M8,C11,C12,C21,C22:
    n:=nops(A):

    if nops(B)<>n then
        RETURN(FAIL):
    fi:

    if not type(log[2](n),integer) then
        RETURN(FAIL):
    fi:

    if n=1 then
        RETURN([[A[1][1]*B[1][1]]]):
    fi:

    A11:=[seq([op(1..n/2,A[i])],i=1..n/2)]:
    A12:=[seq([op(n/2+1..n,A[i])],i=1..n/2)]:
    A21:=[seq([op(1..n/2,A[i])],i=n/2+1..n)]:
    A22:=[seq([op(n/2+1..n,A[i])],i=n/2+1..n)]:



    B11:=[seq([op(1..n/2,B[i])],i=1..n/2)]:
    B12:=[seq([op(n/2+1..n,B[i])],i=1..n/2)]:
    B21:=[seq([op(1..n/2,B[i])],i=n/2+1..n)]:
    B22:=[seq([op(n/2+1..n,B[i])],i=n/2+1..n)]:



    M1:=MatMultR(A11,B11):
    M2:=MatMultR(A12,B21):
    M3:=MatMultR(A11,B12):
    M4:=MatMultR(A12,B22):


    M5:=MatMultR(A21,B11):
    M6:=MatMultR(A22,B21):
    M7:=MatMultR(A21,B12):
    M8:=MatMultR(A22,B22):

    C11:=MatAdd(M1,M2):
    C12:=MatAdd(M3,M4):
    C21:=MatAdd(M5,M6):
    C22:=MatAdd(M7,M8):

    [seq([op(C11[i]),op(C12[i])],i=1..nops(C11)),
     seq([op(C21[i]),op(C22[i])],i=1..nops(C21))]:

end:






#Tra(A): the transpose of the matrix A
Tra:=proc(A) local i,j:
    [seq([seq(A[i][j],i=1..nops(A))],j=1..nops(A[1]))]:
end:


#RandMat(m,n,R): a random m by n matrix with entries between -R and R
#given as a list of lists
RandMat:=proc(m,n,R) local ra,i,j:
    ra:=rand(-R..R):
    [seq([seq(ra(),j=1..n)],i=1..m)]:
end:

#############
# Problem 3 #
#############

isSquare := proc(A):
    if nops(A) <> nops(A[1]) then:
        return false:
    else:
        return true:
    fi:
end:

submatrix := proc(A, rows, cols) local i, j:
    return [seq([seq(A[i][j], j in cols)], i in rows)]:
end:

MatMultRG := proc(A, B) local n:
    if not(isSquare(A) and isSquare(B) and evalb(nops(A) = nops(B))) then:
        return FAIL:
    fi:
    n := nops(A):
    return submatrix(MatMultR(PadZeros(A), PadZeros(B)), {seq(i, i=1..n)}, {seq(i, i=1..n)}):
end:

#############
# Problem 4 #
#############

Strassen := proc(A, B) local n:
    if not(isSquare(A) and isSquare(B) and evalb(nops(A) = nops(B))) then:
        return FAIL:
    fi:
    n := nops(A):
    return submatrix(StrassenR(PadZeros(A), PadZeros(B)), {seq(i, i=1..n)}, {seq(i, i=1..n)}):
end:

StrassenR := proc(A, B) local n, A11, A12, A21, A22, B11, B12, B21, B22, C11, C12, C21, C22, M1, M2, M3, M4, M5, M6, M7:
    n := nops(A):

    if n = 1 then:
        return MatMultR(A, B):
    fi:

    A11 := submatrix(A, 1..n/2, 1..n/2):
    A12 := submatrix(A, 1..n/2, n/2+1..n):
    A21 := submatrix(A, n/2+1..n, 1..n/2):
    A22 := submatrix(A, n/2+1..n, n/2+1..n):

    B11 := submatrix(B, 1..n/2, 1..n/2):
    B12 := submatrix(B, 1..n/2, n/2+1..n):
    B21 := submatrix(B, n/2+1..n, 1..n/2):
    B22 := submatrix(B, n/2+1..n, n/2+1..n):


    M1 := StrassenR(A11 + A22, B11 + B22):
    M2 := StrassenR(A21 + A22, B11):
    M3 := StrassenR(A11, B12 - B22):
    M4 := StrassenR(A22, B21 - B11):
    M5 := StrassenR(A11 + A12, B22):
    M6 := StrassenR(A21 - A11, B11 + B12):
    M7 := StrassenR(A12 - A22, B21 + B22):

    C11 := M1 + M4 - M5 + M7:
    C12 := M3 + M5:
    C21 := M2 + M4:
    C22 := M1 - M2 + M3 + M6:

    return [seq([op(C11[i]), op(C12[i])], i=1..nops(C11)),
            seq([op(C21[i]), op(C22[i])], i=1..nops(C21))]:
end:

#############
# Problem 5 #
#############

RandMat := proc(n,R) local r, i, j:
    r := rand(0..R):
    return [seq([seq(r(), j=1..n)], i=1..n)]:
end:

#############
# Problem 6 #
#############

# Strassen does run significantly faster.

#############
# Problem 7 #
#############

# det(H_n) = 2^n/((2n)!n!) * C_0 * C_1 * ... * C_n, where C_n is the
# nth Catalan number.

# I don't know if this counts as a closed form, but it's the best I've
# got.