##################################################################################### ##################################################################################### ## intHaar ## ## A Maple package for integrals over compact lie groups ## ## Written by Alejandro Ginory and Jongwon Kim ## ## May 29th, 2016 ## ##################################################################################### ##################################################################################### withSF(); with(ListTools); with(combinat); with(group); with(Statistics); with(GraphTheory); jacks(2); print(`intHaar`): print(`Following functions are supported for the integrals`): print(`IntU(polynomial function in g[i,j] and conjugate(g[i,j]))`): print(`IntO(polynomial function in g[i,j])`): print(`IntSp(polynomial function in g[i,j], n :: dimension of Sp)`): ##################################################################################### ## Helper Functions with self explanatory names ## ##################################################################################### nonConjugate := proc (x) not has(x, conjugate) end proc; isPerm := proc (list) type(list, permlist) end proc; PermToPart := proc (P, n) local pt, df; if evalb(nops(P) <> 0) then pt := [seq(nops(P[i]), i = 1 .. nops(P))] else pt := [`\$`(1, n)] end if; df := n-add(pt); sort([op(pt), seq(1, i = 1 .. df)], `>`); end proc; PartToPowerSum := proc (part::list) local partl, var; partl := [seq(cat(p, i), i in part)]; var := mul(i, i in partl) end proc; CoefZ := proc (P, e, vlist::list) local g, res; if e = 1 then res := p else g := freeze(e); res := coeff(subs(e = g, P), g) end if; eval(res, [seq(v = 0, v in vlist)]) end proc; zl := proc (lambda) mul(i^numboccur(lambda, i)*factorial(numboccur(lambda, i)), i = 1 .. max(lambda)); end proc; fLambda := proc (lambda, n) Chi(Reverse(lambda), [`\$`(1, n)]) end proc; CosetType := proc (sigma, n) local G, perm, edges; G := Graph({seq({2*i, 2*i-1}, i = 1 .. n)}); perm := convert(sigma, permlist, 2*n); edges := {seq({perm[2*i], perm[2*i-1]}, i = 1 .. n)}; G := ConnectedComponents(AddEdge(G, edges)); [(1/2)*seq(nops(G[i]), i = 1 .. nops(G))] end proc; mInterpreter := proc (m) convert(Flatten(m), disjcyc) end proc; ##################################################################################### ## Unitary Base ## ##################################################################################### Schur1 := proc (lambda) mul(mul(d+j-i, j = 1 .. lambda[i]), i = 1 .. nops(lambda)) end proc; WeingartenU := proc (sigma, n) simplify(add(Chi(Q, [`\$`(1, n)])*Chi(Q, Reverse(sigma))/Schur1(Reverse(Q)), Q in partition(n))/factorial(n)); end proc; GetPermU := proc (I1, I2) local check, path, T, n, P; check := evalb(sort(I1) = sort(I2)); if check then path := [seq([SearchAll(I1[i], I2)], i = 1 .. nops(I1))]; T := cartprod(path); n := 0; while not T[finished] do n := n+1; P[n] := `~`[op](T[nextvalue]()) end do; P := map(convert, select(isPerm, convert(P, list)), disjcyc); else 0 end if; end proc; IUL := proc (I1, J1, I2, J2) local m, pI, pJ, int, i, j, P; m := nops(I1); pI := GetPermU(I1, I2); pJ := GetPermU(J1, J2); int := 0; if type(pI, numeric) or type(pJ, numeric) then 0; else for i to nops(pI) do for j to nops(pJ) do int := int+WeingartenU(PermToPart(mulperms(pJ[j], invperm(pI[i])), m), m); end do; end do; end if; simplify(int); end proc; ##################################################################################### ## Orthogonal Base ## ##################################################################################### Zonal1 := proc (lambda) mul(mul(d+2*j-i-1, j = 1 .. lambda[i]), i = 1 .. nops(lambda)) end proc; CoefNormalizerO := proc (rho, n) zl(2*rho)/(2^n*factorial(n)) end proc; OLS := proc (lambda, sigma) local n, j, P, vb, cf; n := add(lambda); j := J[seq(i, i in lambda)]; P := subs(a = 2, top(j)); vb := PartToPowerSum(sigma); cf := CoefNormalizerO(sigma, n)*CoefZ(P, vb, [seq(cat(p, i), i = 1 .. n)]); end proc; WeingartenO := proc (p1, p2, n) 2^n*factorial(n)*add(fLambda(2*lambda, 2*n)*OLS(lambda, CosetType(mulperms(p1, p2), n))/Zonal1(lambda), lambda in map(Reverse, partition(n)))/factorial(2*n); end proc; GetPermO := proc (lst) local de, df, dp, T, n, P; de := MakeUnique(lst); df := [seq([SearchAll(i, lst)], i in de)]; if convert(map(type, map(nops, df), even), `and`) then dp := map(setpartition, df, 2); T := combinat[cartprod](dp); n := 0; while not T[finished] do n := n+1; P[n] := `~`[op](T[nextvalue]()) end do; P := map(`@`(mInterpreter, sort), convert(P, list)); else 0 end if; end proc; IOL := proc (Il, Jl) local m, pI, pJ, int, i, j; m := nops(Il); pI := GetPermO(Il); pJ := GetPermO(Jl); int := 0; if type(pI, numeric) or type(pJ, numeric) then 0 else for i to nops(pI) do for j to nops(pJ) do int := int+WeingartenO(pJ[j], invperm(pI[i]), (1/2)*m); end do; end do; end if; simplify(int); end proc; ##################################################################################### ## Symplectic Base ## ##################################################################################### Bijections := proc (A, B) map[3](zip, `[]`, A, combinat:-permute(B)) end proc; GetPermSp := proc (lst, m) local k, de, dI, dJ, check, dp, T, n, P; k := nops(lst); de := remove(proc (x) options operator, arrow; m < x end proc, MakeUnique(lst)); dI := [seq([SearchAll(i, lst)], i = 1 .. m)]; dJ := [seq([SearchAll(i+m, lst)], i = 1 .. m)]; check := convert([seq(evalb(nops(dI[i]) = nops(dJ[i])), i = 1 .. m)], `and`); if check then dp := [seq(Bijections(dI[i], dJ[i]), i = 1 .. nops(dI))]; T := combinat[cartprod](dp); n := 0; while not T[finished] do n := n+1; P[n] := `~`[op](T[nextvalue]()) end do; P := map(`@`(mInterpreter, sort), convert(P, list)); else 0 end if; end proc; WeingartenSp := proc (p1, p2, n) simplify(subs(d = -2*d, WeingartenO(p1, p2, n)*(-1)^n*parity(mulperms(p1, p2)))); end proc; WgSign := proc (p1, p2, n) sign(mul([seq((p1[2*i]-p1[2*i-1])*(p2[2*i]-p2[2*i-1]), i = 1 .. n)])); end proc; ISpL := proc (Il, Jl, n) local m, pI, pJ, int, i, j; m := (1/2)*nops(Il); pI := GetPermSp(Il, n); pJ := GetPermSp(Jl, n); int := 0; if type(pI, numeric) or type(pJ, numeric) then 0 else for i to nops(pI) do for j to nops(pJ) do int := int+WgSign(convert(pJ[j], permlist, 2*m), convert(pI[i], permlist, 2*m), m)*WeingartenSp(pJ[j], invperm(pI[i]), m); end do; end do; end if; subs(d = n, simplify(int)); end proc; ##################################################################################### ## IntHaar ## ##################################################################################### GetCoeff := proc (exp) local parts; parts := [op(exp)]; if has(exp, g) then if type(exp, `^`) then 1 else mul(select(not has, parts, g)); end if; else exp end if; end proc; OneVar := proc (mn) local parts, IJ, i; parts := [op(mn)]; IJ := []; if type(mn, 'indexed') then parts else for i to parts[2] do IJ := [op(IJ), op(parts[1])]; end do; end if; end proc; MonomialToList := proc (mn) if type(mn, `*`) then Flatten(map(OneVar, [op(mn)])) else OneVar(mn) end if; end proc; IOM := proc (mn) local cf, vr, IJ, Il, Jl; cf := GetCoeff(mn); vr := mn/cf; if type(vr, 'numeric') then cf*vr else IJ := MonomialToList(vr); Il := [seq(IJ[2*i-1], i = 1 .. (1/2)*nops(IJ))]; Jl := [seq(IJ[2*i], i = 1 .. (1/2)*nops(IJ))]; cf*IOL(Il, Jl); end if; end proc; IntO := proc (pn) if type(pn, `+`) then add(map(IOM, [op(pn)])) else IOM(pn) end if; end proc; ISpM := proc (mn, n) local cf, vr, IJ, Il, Jl; cf := GetCoeff(mn); vr := mn/cf; if type(vr, 'numeric') then cf*vr else IJ := MonomialToList(vr); Il := [seq(IJ[2*i-1], i = 1 .. (1/2)*nops(IJ))]; Jl := [seq(IJ[2*i], i = 1 .. (1/2)*nops(IJ))]; cf*ISpL(Il, Jl, n) end if; end proc; IntSp := proc (pn, n) if type(pn, `+`) then add(map(ISpM, [op(pn)], n)) else ISpM(pn, n) end if; end proc; hasAbs := proc (x) has(x, abs) end proc; absSelector := proc (x) if type(x, `^`) and has(x, abs) then x else select(hasAbs, x) end if; end proc; breaker := proc (x) if type(x, `^`) then breakerH(x) else mul(map(breakerH, [op(x)])) end if; end proc; breakerH := proc (x) if type(op(x)[2], even) then (op(op(x)[1])*conjugate(op(op(x)[1])))^((1/2)*op(x)[2]) else 0 end if; end proc; IUM := proc (mn) local cf, vr, as, nc, cc, IJ, Il, Jl, KL, Kl, Ll; cf := GetCoeff(mn); vr := mn/cf; as := absSelector(vr); if type(as, numeric) then 0 else vr := vr*breaker(as)/as end if; if type(vr, 'numeric') then cf*vr else nc := select(nonConjugate, vr); cc := conjugate(vr/nc); IJ := MonomialToList(nc); Il := [seq(IJ[2*i-1], i = 1 .. (1/2)*nops(IJ))]; Jl := [seq(IJ[2*i], i = 1 .. (1/2)*nops(IJ))]; KL := MonomialToList(cc); Kl := [seq(KL[2*i-1], i = 1 .. (1/2)*nops(KL))]; Ll := [seq(KL[2*i], i = 1 .. (1/2)*nops(KL))]; cf*IUL(Il, Jl, Kl, Ll) end if; end proc; IntU := proc (pn, n) if type(pn, `+`) then add(map(IUM, [op(pn)])) else IUM(pn) end if; end proc;