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

# Nathaniel Shar
# HW 3
# Experimental Mathematics

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

with(LinearAlgebra):
Help := proc(): print(`hilbertMatrix(n), inv(pi), invGF(n,q), dr(M, S), colsumprod(A), PerFast(M)`): end:
#############
# Problem 2 #
#############

hilbertMatrix := proc(n):
    return [seq([seq(1/(i+j-1), i=1..n)], j=1..n)]:
end:

# "Good agreement": ratio is 1 plus or minus 0.01
# n        Digits
# 10       14
# 20       29
# 30       44
# 40       59
# 50       75
# 60       90
# 70       105
# 80       120
# 90       135
# 100      149

# The relationship appears to be linear.

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

with(combinat):

inv := proc(pi) local n, i, j, ct:
    ct := 0:
    n := nops(pi):
    for i from 1 to n do:
        for j from i+1 to n do:
            if pi[i] > pi[j] then:
                ct := ct + 1:
            fi:
        od:
    od:
    return ct:
end:

invGF := proc(n,q) local pi:
    return add(q^inv(pi), pi in permute(n)):
end:

# equals prod((1-q^i)/(1-q), i=1..n)

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

# delete the rows indexed by the set S
dr := proc(M, S) local n, A, i:
    n := nops(M):
    A := []:
    for i from 1 to n do:
        if not i in S then:
            A := [op(A), M[i]]:
        fi:
    od:
    return A:
end:

colsumprod := proc(A) local i, j, s, R:
    s := 1:
    for i from 1 to nops(A[1]) do:
        s := s * add(j, j=map(R->R[i], A)):
    od:
end:

PerFast := proc(M) local S,n:
    n := nops(M):
    add((-1)^i*add(colsumprod(dr(M, S)), S in choose(n, i)), i=0..n-1)
end: