#Nathan Fox
#Homework 3
#I give permission for this work to be posted online

#Read procedures from class
read(`C3.txt`):

#Linear algebra
with(LinearAlgebra):

#Help procedure
Help:=proc():
    print(` FindHilThreshold(eps) , inv1(p) , gf(n, q) `):
    print(` RandMat(n) , FastPer(M) `):
end:

##PROBLEM 2##

#Find number of digits such that
#abs(1-det(evalf(Hil(n)))/evalf(det(Hil(n)))) < eps;
#Uses a binary search
FindHilThreshold:=proc(n, eps) local olddig, lo, hi, ed, de, test:
    olddig:=Digits:
    lo:=1:
    hi:=infinity:
    Digits:=10:
    while lo <> hi - 1 do
        de:=Determinant(evalf(Hil(n))):
        ed:=evalf(Determinant(Hil(n))):
        test:=abs(1-de/ed):
        if test < eps then
            hi:=Digits:
            Digits:=floor((lo + hi)/2):
        else
            lo:=Digits:
            if hi = infinity then
                Digits:=Digits * 2:
            else
                Digits:=ceil((lo + hi)/2):
            fi:
        fi:
    od:
    Digits:=olddig:
    return hi:
end:
#FindHilThreshold(50, 1e-1) returns 74
#FindHilThreshold(50, 1e-5) returns 78
#FindHilThreshold(50, 1e-10) returns 83
#Hence, good agreement occurs around 75 to 80 digits
#
#seq(FindHilThreshold(10*i, 1e-10), i=1..10) returns
#22, 37, 53, 68, 83, 99, 113, 128, 142, 159
#This sequence looks rather linear, as its successive differences
#are 15, 16, 15, 15, 16, 14, 15, 14, 17

##PROBLEM 3##

#Count inversions
inv1:=proc(p) local i, j, n, s:
    s:=0:
    n:=nops(p):
    for i from 1 to n-1 do
        for j from i+1 to n do
            if p[i] > p[j] then
                s:=s+1:
            fi:
        od:
    od:
    return s:
end:

#Generating function
gf:=proc(n, q) local p:
    return add(q^inv1(p), p in permute(n)):
end:
#seq(factor(gf(n,q)), n=1..8) returns
#1,
#1+q,
#(1+q)*(q^2+q+1),
#(q^2+1)*(q^2+q+1)*(1+q)^2,
#(q^2+q+1)*(q^4+q^3+q^2+q+1)*(q^2+1)*(1+q)^2,
#(q^2-q+1)*(q^4+q^3+q^2+q+1)*(q^2+1)*(q^2+q+1)^2*(1+q)^3,
#(q^2-q+1)*(q^4+q^3+q^2+q+1)*(q^6+q^5+q^4+q^3+q^2+q+1)*(q^2+1)*(q^2+q+1)^2*(1+q)^3,
#(q^2-q+1)*(q^4+q^3+q^2+q+1)*(q^4+1)*(q^6+q^5+q^4+q^3+q^2+q+1)*(q^2+1)^2*(q^2+q+1)^2*(1+q)^4
#Conjecture: gf(n,q) = mul(cyc(n)^(n/i), i=1..n), where cyc(n) is
#the nth cyclotomic polynomial

##PROBLEM 4##

#Random nxn matrix will all entries -99..99
#Necessary because we use a nonstandard matrix format
RandMat:=proc(n) local r, i, j:
    r:=rand(-99..99):
    return [seq([seq(r(), i=1..n)], j=1..n)]:
end:

#"Fast" permanents
PerFast:=proc(M) local n, i, j, k, S:
    n:=nops(M):
    return (-1)^n*add((-1)^nops(S)*
            mul(
                add(M[i][j], j in S),
                i=1..n),
        S in powerset({seq(k, k=1..n)})):
end:
#M 8x8:
#time(Per(M)) returned 0.818
#time(PerFast(M)) returned 0.018
#
#M 9x9:
#time(Per(M)) returned 12.203
#time(PerFast(M)) returned 0.046