#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