#OK to post homework
#Joseph Koutsoutis, 03-09-2025, Assignment 13

read `C13.txt`:

with(combinat):

#1

FindAllTriples := proc() local i,j,k,r,triples,M:
  triples := {}:
  for i from 2 to 7 do:
    for j from 2 to 7 do:
      for k from 2 to 7 do:
        M := simplify(Mul(Mul(C1QG()[i],C1QG()[j]),C1QG()[k])):
        for r from 1 to 7 do:
          if M = simplify(C1QG()[r]) then: 
            triples := triples union {[i,j,k,r]}: 
          fi:
        od:
      od:
    od:
  od:
  triples:
end:

#FindAllTriples outputted 57 values. I'm not positive this is all of them though since I had to
#add simplify to get equality sometimes.


#2

FindAllQuadruples := proc() local i,j,k,r,c,quads,M,mul_entries_possible,c1qg_entries_possible,possible_c:
  c1qg_entries_possible := {}:
  mul_entries_possible := {}:
  quads := {}:
  for i from 2 to 7 do:
    c1qg_entries_possible := c1qg_entries_possible union {op(C1QG()[i][1]), op(C1QG()[i][2])}:
    for j from 2 to 7 do:
      for k from 2 to 7 do:
        M := simplify(Mul(Mul(C1QG()[i],C1QG()[j]),C1QG()[k])):
        mul_entries_possible := mul_entries_possible union {op(M[1]), op(M[2])}:
      od:
    od:
  od:

  mul_entries_possible := mul_entries_possible minus {0}:
  c1qg_entries_possible := c1qg_entries_possible minus {0}:
  possible_c := {}:
  for i from 1 to nops(mul_entries_possible) do:
    for j from 1 to nops(c1qg_entries_possible) do:
      if simplify(abs(c1qg_entries_possible[j] / mul_entries_possible[i])) = 1 then:
        possible_c := possible_c union {simplify(c1qg_entries_possible[j] / mul_entries_possible[i])}:
      fi:
    od:
  od:

  quads := {}:
  for i from 2 to 7 do:
    for j from 2 to 7 do:
      for k from 2 to 7 do:
        M := simplify(expand(Mul(Mul(C1QG()[i],C1QG()[j]),C1QG()[k]))):
        for r from 1 to 7 do:
          for c in possible_c do:
            if M = simplify(expand(c * C1QG()[r])) then: 
              quads := quads union {[i,j,k,r,c]}: 
            fi:
          od:
        od:
      od:
    od:
  od:
  quads:
end:


#FindAllQuadruples() outputted 80. I'm not positive this is all of them though since I had to
#add simplify/expand to various places to get equality sometimes.


#3

ApplyT := proc(T,X,n,u) local T1, V, v, val:
  T1, V := op(T):
  val := 0:
  for v in V do:
    val += coeff(u, v, 1) * T1[v]:
  od:
  simplify(val):
end:

#ApplyT(T,X,2,X[[0,0]]+I*X[[0,1]]+2*I*X[[1,0]]+5*X[[1,1]]) outputted
#(((5 - I)*X[[1, 1]] - (5 + I)*X[[0, 1]] + 
#  (1 - 2*I)*X[[1, 0]] + (1 + 2*I)*X[[0, 0]])*sqrt(2))/2


#5

RM := proc() local ra:
  ra := rand(1..10):
  [[ra(),ra()],[ra(),ra()]]:
end:

RMs := proc(n) local i:
  [seq(RM(), i=1..n)]:
end:

ExperimentS3 := proc() local A, pi, eqs, M, var, i:
  var := {seq(c[pi], pi in permute(3))}:
  eqs := {}:

  for i from 1 to 2 do:
    A := RMs(3):
    M := expand(add(c[pi]*Mul(Mul(A[pi[1]],A[pi[2]]),A[pi[3]]), pi in permute(3))):
    eqs := eqs union {op(M[1]), op(M[2])}:
  end:
  solve(eqs, var):
end:


#The only solution found was c[pi] = 0 for pi in permute(3).


ExperimentS4 := proc() local A, pi, eqs, M, var, i, var1, v, T:
  var := {seq(c[pi], pi in permute(4))}:
  eqs := {}:

  for i from 1 to 6 do:
    A := RMs(4):
    M := expand(add(c[pi]*Mul(Mul(A[pi[1]],A[pi[2]]),Mul(A[pi[3]],A[pi[4]])), pi in permute(4))):
    eqs := eqs union {op(M[1]), op(M[2])}:
  end:
  var := solve(eqs, var):

  var1:={}:
  for v in var do
    if op(1,v)=op(2,v) then
      var1:=var1 union {op(1,v)}:
    fi:
  od:

  for v in var do:
    T[op(1, v)] := subs({var1[1]=1,seq(var1[i]=0,i=2..nops(var1))},op(2,v)):
  od:
  op(T):
end:


#This outputted
#table([c[[4, 1, 2, 3]] = -1, c[[2, 1, 3, 4]] = -1, c[[3, 4, 2, 1]] = -1, c[[1, 2, 4, 3]] = -1, 
#       c[[1, 4, 3, 2]] = -1, c[[2, 1, 4, 3]] = 1, c[[3, 4, 1, 2]] = 1, c[[4, 2, 1, 3]] = 1, 
#       c[[2, 3, 1, 4]] = 1, c[[4, 1, 3, 2]] = 1, c[[2, 4, 1, 3]] = -1, c[[2, 3, 4, 1]] = -1, 
#       c[[1, 4, 2, 3]] = 1, c[[4, 3, 2, 1]] = 1, c[[2, 4, 3, 1]] = 1, c[[3, 1, 2, 4]] = 1, 
#       c[[4, 2, 3, 1]] = -1, c[[4, 3, 1, 2]] = -1, c[[1, 3, 4, 2]] = 1, c[[1, 2, 3, 4]] = 1, 
#       c[[3, 1, 4, 2]] = -1, c[[1, 3, 2, 4]] = -1, c[[3, 2, 1, 4]] = -1, c[[3, 2, 4, 1]] = 1])


TestConjecture := proc(T) local A, pi:
  A := [[[a[1,1,1], a[1,1,2]], [a[1,2,1], a[1,2,2]]],
        [[a[2,1,1], a[2,1,2]], [a[2,2,1], a[2,2,2]]],
        [[a[3,1,1], a[3,1,2]], [a[3,2,1], a[3,2,2]]],
        [[a[4,1,1], a[4,1,2]], [a[4,2,1], a[4,2,2]]]]:
  expand(add(T[c[pi]]*Mul(Mul(A[pi[1]],A[pi[2]]),Mul(A[pi[3]],A[pi[4]])), pi in permute(4))):
end:

#TestConjecture(T) outputs [[0,0],[0,0]] with the output of ExperimentS4() from above.