TabProcs := module()
export MapHW, Casimir, Racah, CGC, Apply, EncodeO, DecodeO;
local Init, Gen_Tabs, Gen_HW, GenBorelOps, NullSpace,
St, St_, St_2, t_thres, status,
ApplyTable, MinorsA, MinorsB, MinorsBtoO, NumMinors, Sizes,
irrep, tensors, Wrd, piL_1r, B, M_1r, Delta_1r, R_Delta_1r, Delta_c,
R_Delta_c, Stop,
Dimension, rDim, cDim, MMM, Equal, Transpose, BinarySearch;
global Irrep, Tensors, PETE1, PETE2, PETE3;
option package;
# Written by Stephen Gliske (sgliske@umich.edu)
Dimension := LinearAlgebra:-Dimension;
rDim := LinearAlgebra:-RowDimension;
cDim := LinearAlgebra:-ColumnDimension;
MMM := LinearAlgebra:-MatrixMatrixMultiply;
Equal := LinearAlgebra:-Equal;
Transpose := LinearAlgebra:-Transpose;
BinarySearch := ListTools:-BinarySearch;
status := "initialized";
t_thres := 10;
Stop := proc(a,b) if nargs=0 or a=b then error("stop") end if end proc:
# (1.0) Creates Tableaux, Applies Operators from Tableaux, Solves Borel Condition
MapHW := proc(in1 :: {list, identical(Terse), identical(Verbose)}, in2 :: list,
in3 :: {identical(Terse), identical(Verbose)})
local verbosity, tabs, piL_, HW, `M~`, b_Ops, dim, i, j, k, c, temp, C, a, b,
alpha, N, B_, R_delta_b, mu, thres, tensors_;
St := time(); St_ := time();
St_2 := true;
status := "MapHW: 1r";
# Check inputs
if nargs > 1 then irrep := in1; tensors := in2;
else irrep := Irrep; tensors := Tensors
end if;
if nargs = 1 then verbosity := in1
elif nargs = 3 then verbosity := in3
else verbosity := 'Terse'
end if;
if verbosity = 'Terse' then verbosity := false;
elif verbosity = 'Verbose' then verbosity := true;
else error("unknown verbosity");
end if;
# Check conditions on irrep and tensors
tensors_ := select(`>`, ListTools:-Flatten(tensors), 0);
Wrd := nops(irrep);
if nops(tensors_) <> nops(irrep) then
error("Number of entries in Irrep does not equal number of non-zero entries in Tensors.")
elif convert(tensors_, `+`) <> convert(irrep, `+`) then
error("Weight of Irrep does not equal weight of Tensors.")
elif not ListTools:-Sorted( [seq(irrep[-i], i=1..Wrd)] ) then
error("Irrep does not meet dominance condition.")
end if;
# Generate Tables
status := "MapHW: init, 1r";
Init:-Gen_Tables();
NumMinors := irrep[1];
# Generate Tableau
status := "MapHW: Generating Tableau, 1r";
if ((verbosity and St_-time() > 10) or not verbosity) then
printf("\n \t %8.3f Generating Tableau", time()-St); St_ := time();
end if;
tabs := Gen_Tabs(irrep, tensors_);
if nops(tabs) = 0 then error("Multiplicity is zero %1", time()-St) end if;
if (verbosity and St_-time() > t_thres) or not verbosity then
printf("\n \t %8.3f \tNumber of Tableau is %d", time()-St, nops(tabs));
St_ := time();
end if;
# Derive the operators from all the tableau.
piL_1r := map(T->[seq(seq( ((Wrd-i)*Wrd+j) $ T[i][j]-T[i+1][j],
j=1..nops(T[i+1])), i=1..Wrd-1)], tabs);
# Generate Highest Weight Vector;
HW := Gen_HW(irrep, 0);
# Apply piL
status := "MapHW: Applying PiL, 1r";
unassign('R_Delta_1r'): unassign('Delta_1r'): k := 0:
for i from 1 to nops(piL_1r) do
c, temp := Apply( piL_1r[i], [1], [HW] ):
if c = 0 then
C[i] := Vector[row]([0]);
else
for j from 1 to nops(c) do
if not assigned(R_Delta_1r[temp[j]]) then
k := k + 1;
R_Delta_1r[temp[j]] := k;
Delta_1r[k] := temp[j];
end if;
end do:
C[i] := Vector[row](k, datatype=numeric):
for j from 1 to nops(c) do
C[i][R_Delta_1r[temp[j]]] := c[j];
end do
end if;
end do:
Delta_1r := Vector(k, Delta_1r); dim := k:
`M~` := Matrix(nops(piL_1r), k, datatype=integer, storage=sparse):
for i from 1 to nops(piL_1r) do
`M~`[i, 1..-1] := C[i];
end do:
# Apply Borel Condition
status := "MapHW: Borel, 1r";
b_Ops := GenBorelOps(tensors);
if b_Ops = [] then
B := 1; M_1r := `M~`; mu := rDim(M_1r);
else
if (verbosity and St_-time() > t_thres) or not verbosity then
printf("\n \t %8.3f Solving Borel Condition", time()-St); St_ := time();
end if;
B := 1; M_1r := Matrix(`M~`, storage=rectangular);
for i from 1 to nops(b_Ops) do
if Dimension(Delta_1r) > 5*t_thres^2 then
printf("\n \t %8.3f Solving Borel Condition %d of %d", time()-St, i, nops(b_Opts));
St_ := time();
end if;
k := 0; unassign('R_delta_b'); unassign('N');
for j from 1 to dim do
c, temp := Apply( [b_Ops[i]], [1], [Delta_1r[j]] ):
if c <> 0 then
for a from 1 to nops(c) do
if not assigned(R_delta_b[temp[a]]) then
k := k + 1;
R_delta_b[temp[a]] := k;
end if;
end do:
# Add result into Matrix Mult
for b from 1 to rDim(M_1r) do
for a from 1 to nops(c) do
alpha := R_delta_b[temp[a]];
if assigned(N[b,alpha]) then
N[b,alpha] := N[b,alpha] + M_1r[b,j]*c[a];
else
N[b,alpha] := M_1r[b,j]*c[a];
end if;
end do;
end do:
end if;
end do:
if k = 0 then error("Multiplicity is zero") end if;
if rDim(M_1r)*k > 10000 then
printf("\n\t %8.3f \t\t Solving %d by %d Matrix",
time()-St, k, rDim(M_1r));
end if;
B_ := NullSpace(Matrix(rDim(M_1r), k, N, datatype=integer)):
M_1r := LinearAlgebra:-MatrixMatrixMultiply(B_, M_1r);
if i = 1 then B := B_
else B := LinearAlgebra:-MatrixMatrixMultiply(B_, B);
end if;
end do:
mu := rDim(M_1r);
# Rescaling
if ((verbosity and St_-time() > 10) or not verbosity) then
printf("\n \t %8.3f Rescaling", time()-St, mu);
end if;
for i from 1 to mu do
B[i,1..-1] := B[i,1..-1]/igcd(seq(B[i,q], q=1..cDim(B)));
M_1r[i,1..-1] := M_1r[i,1..-1]/igcd(seq(M_1r[i,q], q=1..cDim(M_1r)));
end do:
end if;
# Complete
setattribute(M_1r, HW);
status := "MapHW: done, 1r";
if ((verbosity and St_-time() > 10) or not verbosity) then
printf("\n \t %8.3f Complete: Multiplicity is %d", time()-St, mu);
St_ := time();
end if;
# Return values if not verbosity
if not verbosity then
if nargs < 2 then
return(time()-St);
else
return(M_1r, Delta_1r, B, piL_1r);
end if;
end if;
# Print values...
thres := 4;
# Tableau
printf(" \n \n");
if nops(tabs) <= thres then print(`Tableau`)
else print(`First few Tableau`);
end if;
print(`- - - - - - - - - - - - - - - - - -`);
print( op(map(x->Matrix(x, fill=` `), tabs[1..min(thres, nops(tabs))])) ):
# Operators from Tableau
printf(" \n \n");
temp := map(X->map( a -> [ iquo(a-1, Wrd)+1, irem(a-1, Wrd)+1 ], X),
piL_1r[1..min(thres, nops(piL_1r))]);
if nops(tabs) < thres + 1 then
print(`Operators from Tableau`);
print(`- - - - - - - - - - - - - - - - - - - - - - -`);
else
print(`Operators from first few Tableau`);
print(`- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -`);
end if;
for i from 1 to nops(temp) do
print( mul(L[temp[i][j][1], temp[i][j][2]], j=1..nops(temp[i])) );
end do;
# f_max
printf(" \n \n");
print('f'['max']); print(`- - - - - - - - -`);
#print(convert(map(X->Delta^seq( i $ X[Wrd-i+1], i=1..Wrd), map(x->MinorsA[x, 1..-1], DecodeO(HW))), `*`));
print(convert(map(x->Delta^seq( i $ MinorsA[x,-i], i=1..Wrd), DecodeO(HW)), `*`));
# Basis of Minors
printf(" \n \n");
if Dimension(Delta_1r) <= 2*thres then
print(`Basis of Products of Minors of Determinants`);
else
print(`Begining of Basis of Products of Minors of Determinants`);
end if;
print(`- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -`);
print();
a := min(Dimension(Delta_1r), 2*thres);
temp := Vector([seq(convert(map(x->Delta^seq( i $ MinorsA[x,-i], i=1..Wrd), DecodeO(Delta_1r[i])), `*`),
i=1..a)]);
print(temp);
# M~, Borel Ops, M
printf(" \n \n");
if nops(b_Ops) > 1 then
print(`Matrix of Coefficients before Borel Condition`);
print(`- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -`);
print(`M~`);
printf(" \n \n");
print(`Borel Operators`);
print(`- - - - - - - - - - - - - - - - -`);
print(seq(L[op(b_Ops[i])], i=1..nops(b_Ops)));
printf(" \n \n");
print(`Borel Constraint Matrix`);
print(`- - - - - - - - - - - - - - - - - - - - -`);
print(B);
printf(" \n \n");
print(`Matrix of Coefficients after Borel`);
else
print(`Matrix of Coefficients`);
end if;
print(`- - - - - - - - - - - - - - - - - - - - - - - -`);
print(M_1r);
# Few actual vectors
printf(" \n \n ");
if mu <= thres then
if cDim(M_1r) <= Dimension(temp) then
print(`Final Vectors`);
else
print(`First Terms of final Vectors`);
end if;
else
if cDim(M_1r) <= Dimension(temp) then
print(`First several final Vectors`);
else
print(`First terms of first several final Vectors`);
end if;
end if;
print(`- - - - - - - - - - - - - - - - - - - - - - - - - -`);
for i from 1 to min(mu, thres) do
print(M_1r[i, 1..Dimension(temp)].temp);
printf(" \n ");
end do;
if nargs < 2 then
return(time()-St);
else
return(M_1r, Delta_1r, B, piL_1r);
end if;
end proc:
# (1.1) Initialization Routines
Init := module()
local BtoA, Apply_For_Table, two_powers, irrepWs, Lop_bin,
MinBase, MinBasePowers;
export Gen_Tables, MaxComOrd;
Gen_Tables := proc() local i, maxNumMin;
# Lops prepared to generate ApplyTable
two_powers := [seq(2^(i-1), i=1..Wrd)]:
Lop_bin := Matrix(Wrd,Wrd, (i,j)->two_powers[i]-two_powers[j],
shape=antisymmetric):
# determine what size minors will be needed
irrepWs := [seq(irrep[i]-irrep[i+1], i=1..nops(irrep)-1), irrep[-1]];
Sizes := zip((a,b)->if a<>0 then b end if, irrepWs,
[$1..nops(irrep)]);
Sizes := [op({op(Sizes),
op(map(op, map( Y-> zip((a,b)->if a<>0 then b end if, Y,
[$1..nops(Y)]), map( X -> [seq(X[i]-X[i+1], i=1..nops(X)-1), X[-1]],
tensors) )))})];
# Make Ordinal Lists of Minors for ApplyTable generation
MinorsB := map(x->add( two_powers[j], j in x),
ListTools:-FlattenOnce(map(x->combinat:-choose(Wrd, x), Sizes))):
MinorsA := Matrix(nops(MinorsB), Wrd, datatype=integer[1]):
for i from 1 to nops(MinorsB) do
MinorsA[i, 1..Wrd] := BtoA(MinorsB[i]);
end do:
MinorsBtoO := table( {seq(MinorsB[i]=i, i=1..nops(MinorsB)) });
# Make application table
ApplyTable := Matrix(Wrd^2, nops(MinorsB), Apply_For_Table,
storage=sparse, datatype=integer[2]);
MinBase := nops(MinorsB):
maxNumMin := max( irrep[1], seq(tensors[i][1], i=1..nops(tensors)));
MinBasePowers := Vector(maxNumMin, fill=1):
for i from 2 to maxNumMin do
MinBasePowers[-i] := MinBasePowers[-i+1]*MinBase:
end do:
MaxComOrd := MinBasePowers[1]*MinBase-1;
NULL;
end proc:
BtoA := proc(x_in) local x, l, i;
x := x_in;
l := Vector[row](Wrd, datatype=integer[1]):
for i from 1 to Wrd while x > 0 do
x := iquo(x, 2, 'l'[-i]);
end do:
l;
end proc:
# Lop a, minor b
# Uses MinorsA, MinorsB
Apply_For_Table := proc(a, b)
local Lop;
Lop := iquo(a-1, Wrd, 'Lop')+1, Lop+1;
if Lop[1] = Lop[2] and MinorsA[b, -Lop[1]] = 1 then
return(b)
elif MinorsA[b, -Lop[1]] = 0 and MinorsA[b, -Lop[2]] = 1 then
if is(add( MinorsA[b, -k], k=min(Lop)+1..max(Lop)-1 ), even) then
return(MinorsBtoO[MinorsB[b]+Lop_bin[Lop]]);
else
return(-MinorsBtoO[MinorsB[b]+Lop_bin[Lop]]);
end if;
else
return(0)
end if;
end proc:
# Encode/Decode Ordinals to Compact Ordinals
EncodeO := O -> add( MinBasePowers[i-NumMinors-1]*(sort(O, `<`)[i]-1),
i=1..NumMinors):
DecodeO := proc(x_in) local x, l, i;
x := x_in;
l := Vector[row](NumMinors):
for i from 1 to NumMinors while x > 0 do
x := iquo(x, MinBase, 'l'[-i]);
end do:
return( l+Vector[row](NumMinors, fill=1) );
end proc: end module:
# (1.2) Generate Tableau.
Gen_Tabs := proc(irrep_in, weight)
local irrep, z, dim, Rows, between, pSums, tabs, tabs_new, i, j, y;
z := ListTools:-Occurrences(0, irrep_in);
if z > 1 then irrep := irrep_in[1..-z]
else irrep := irrep_in;
end if;
dim := nops(weight);
between := (a, b) -> not (false in
{op(zip((x,y) -> evalb(x >= y), a, b)),
op(zip((x,y) -> evalb(x <= y), a[2..-1], b))}):
pSums := [seq(add(weight[j], j=1..i), i=1..dim-1)];
# Prepare Possible Rows
Rows := Vector(dim-2, x -> select( ListTools[Sorted],
convert(combinat:-composition(pSums[-x]+dim-max((z-1), x),
dim-max((z-1),x)), list)));
LinearAlgebra:-Map(X->map(x->x-[1 $ nops(x)], X), Rows);
LinearAlgebra:-Map2(select, x->x[-1] <= irrep[1], Rows):
LinearAlgebra:-Map2(select, x->x[1] >= irrep[-1], Rows):
LinearAlgebra:-Map(X->map(x->[seq(x[-q], q=1..nops(x))], X), Rows);
Rows := Vector([Rows, [[[weight[1]]]]]);
for i from 1 to dim-1 do
if nops(Rows[i][1]) <> dim-i then
Rows[i] := map(X -> [op(X), 0 $ dim-i-nops(X)], Rows[i]);
end if;
end do;
tabs := [[irrep_in]];
i := 2;
for i from 2 to dim do
tabs_new := [];
for j from 1 to nops(tabs) do
for y in Rows[i-1] do
if between(tabs[j][i-1], y) then
tabs_new := [op(tabs_new), [op(tabs[j]), y]];
end if;
end do;
end do;
tabs := tabs_new;
end do:
end proc:
# (1.3) Generate HW vector for irrep irr
Gen_HW := (irr, p) -> EncodeO( map( x->MinorsBtoO[x],
[seq( (2^i-1) $ (irr[i]-irr[i+1]), i=1..nops(irr)-1), (2^nops(irr)-1) $ irr[-1]]*2^p)):
# (1.4) Generate Operators for Borel Condition
GenBorelOps := proc(tensors) local nonzero, shifts, temp;
nonzero := map(X->nops(select(`>`, X, 0)), tensors);
shifts := [0, seq(convert(nonzero[1..i], `+`), i=1..nops(nonzero)-1)];
temp := [seq( op( map( x -> x+[shifts[i], shifts[i]],
combinat:-choose(nonzero[i], 2)) ), i=1..nops(tensors))];
map( x->((x[1]-1)*Wrd + x[2]), temp);
end proc:
# (1.5) Compute Nullspace for Borel Condition
NullSpace := proc(N)
local rd, cd, c, r, c_, r_, w, i, j, k, x, L_, B;
rd, cd := Dimension(N):
c := []: r := []:
c_ := {$1..cd}:
# Make (quasi) upper triagular with partial pivoting.
for i from 1 to rd while c_ <> {} do
# Find nonzero element for partial pivot
for w in c_ while N[i,w] = 0 do end do:
# Pivot about the point
if N[i,w] <> 0 then
c_ := c_ minus {w}:
c := [op(c), w]:
r := [op(r), i]:
N[1..-1,w] := N[1..-1,w]/igcd(seq(N[j,w], j=i..rd));
# option to set diagonal to one
# N[1..-1,w] := N[1..-1,w]/N[i,w]:
for j in c_ do
x := N[i,j];
if N[i,j] <> 0 then
for k from i to rd do N[k,j] := N[i,w]*N[k,j] - x*N[k,w] end do:
end if;
end do:
end if:
end do:
r_ := [op({$(1..rd)} minus {op(r)})]:
w := ilcm(seq(N[r[i],c[i]], i=1..nops(c)));
L_ := LinearAlgebra:-MatrixInverse(Matrix( nops(c), nops(c),
(i,j) -> N[r[i],c[j]], shape=triangular[lower]), method=subs)*w;
try:
L_ := Matrix(L_, datatype=integer, shape=triangular[lower]):
catch:
k := 1;
for i from 1 to nops(c) do
for j from 1 to i do
if not is(L_[i,j], integer) then
k := lcm(k, denom(L_[i,j]));
end if;
end do;
end do;
L_ := Matrix(k*L_, datatype=integer):
w := k*w;
end:
B := Matrix(nops(r_), rd, datatype=integer):
for j from 1 to rd do
k := ListTools:-BinarySearch(r, j);
if k = 0 then
k := ListTools:-BinarySearch(r_, j);
B[k,j] := -w;
else
for i from 1 to nops(r_) do
B[i,j] := add( N[r_[i], c[q]]*L_[q, k], q=1..nops(c));
end do:
end if;
end do:
return(B);
end proc:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# (2.0) Application of Operators
# input list of ops, list of coefs, list of encoded products of minors
Apply := proc(Lops, c_in, delta_in)
local a, i, j, k, x, W, temp, Changes, c, c_, deltaO, deltaO_, R_delta, delta_, Cas_end;
c := c_in;
deltaO := Matrix(nops(c), NumMinors):
for i from 1 to nops(c) do
deltaO[i,1..-1] := DecodeO(delta_in[i]);
end do:
Cas_end := false;
for a from 1 to nops(Lops) do
Changes := map(x->ApplyTable[Lops[a], x], deltaO);
if a = nops(Lops) and status[1..7] = "Casimir" then
Cas_end := true;
if status[-1] = "c" then
c_ := Vector[row]( Dimension(Delta_c) );
else
c_ := Vector[row]( Dimension(Delta_1r) );
end if;
else
c_ := Vector( rtable_num_elems(Changes, NonZero) );
deltaO_ := Vector( rtable_num_elems(Changes, NonZero) );
end if;
k := 0:
# Apply Chain Rule, sort like terms
# am trusting Maple to ensure W <> 0
for W in op(2, Changes) do
i, j := op(1, W):
W := op(2, W);
temp := deltaO[i, 1..-1];
temp[j] := abs(W);
x := EncodeO(convert(temp, list));
if Cas_end then
if status[-1] = "c" then
if not assigned(R_Delta_c[x]) then
error("Casimir Operators left the Product Space")
end if;
c_[R_Delta_c[x]] := c_[R_Delta_c[x]] + sign(W)*c[i];
else
if not assigned(R_Delta_1r[x]) then
error("Casimir Operators left the Product Space")
end if;
c_[R_Delta_1r[x]] := c_[R_Delta_1r[x]] + sign(W)*c[i];
end if;
elif assigned(R_delta[x]) then
c_[R_delta[x]] := c_[R_delta[x]] + sign(W)*c[i];
else
k := k + 1;
R_delta[x] := k;
c_[k] := sign(W)*c[i];
deltaO_[k] := temp;
if a=nops(Lops) then
delta_[k] := x;
end if;
end if;
end do:
# Store results, remove those with coef 0
if not Cas_end then
temp := map(lhs, op(2, c_));
c := [seq(c_[i], i in temp)];
if nops(temp) = 0 then
break;
elif a < nops(Lops) then
deltaO := Matrix(nops(temp), NumMinors, (i,j) -> deltaO_[temp[i]][j]):
c_ := 'c_': deltaO_ := 'deltaO_': R_delta := 'Rdelta':
end if;
end if;
end do:
if status = "Racah" then
if nops(c) = 0 then
return(0,0)
else
return(c, [seq(delta_[i], i in temp)] );
end if;
elif status[1..3] = "Map" then
if nops(c) = 0 then
return(0,0)
else
return(c, Vector([seq(delta_[i], i in temp)] ));
end if;
elif status[1..7] = "Casimir" then
if nops(temp) = 0 then
if status[-1] = "r" then return(Vector[row](Dimension(Delta_1r)))
else return(Vector[row](Dimension(Delta_c)))
end if;
else
return(c_);
end if;
else
error("Unknown Status %1", status)
end if;
end proc:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# (2.0) Applies Casimir Operators
Casimir := module()
export first, additional, Couple, Get_C, EV;
local UnBlock, iM_1r, lambda, P1, iP1, used_C;
# add functionality for inputing Map info later...
first := proc( in1 :: {list, set}, in2 :: posint );
St := time(); St_ := time();
if status = "initialized" then
error("Must run MapHW first")
else
status := "Casimir: first, 1r";
end if;
try:
# tables already set, since so is M, Delta...
iM_1r := LinearAlgebra:-MatrixInverse( Matrix(M_1r, datatype=algebraic),
outputoptions=[datatype=rational]);
used_C := Vector(1);
printf("\n \t %8.3f Applying Operator", time()-St);
used_C[1] := MMM(M_1r, MMM( Get_C(args), iM_1r));
printf("\n \t %8.3f Matrix Found--determining eigenvalues", time()-St);
lambda, P1 := LinearAlgebra:-Eigenvectors(used_C[1]);
P1 := simplify(P1, factor);
lambda := Matrix([lambda]);
try:
iP1 := LinearAlgebra:-MatrixInverse(P1);
catch "singular":
print("Singular Matrix--Do not attempt to Simultaniously Diagnolize");
end:
printf("\n \t %8.3f Complete", time()-St);
return(lambda, P1)
finally:
St_ := time()-St;
status := "first Casimir complete":
end:
end proc:
additional := proc( in1 :: {list, set}, in2 :: posint )
local C2, i, temp, lambda2, P2, iPt, Pt;
if not (status = "first Casimir complete" or status = "additional Casimir complete") then
error("Must run Casimir[first] prior to running Casimir[additional]");
else
status := "Casimir, additional, 1r";
end if;
try:
St := time()-St_;
printf("\n\t %8.3f Applying Operator", time()-St);
C2 := MMM(M_1r, MMM( Get_C(args), iM_1r));
# check comutation
for i from 1 to Dimension(used_C) do
if not Equal(MMM(used_C[i], C2), MMM(C2, used_C[i])) then
error("Casimir does not commute with previous Casimir #%1", i);
end if;
end do;
# simultaniously diagnolize
printf("\n\t %8.3f Simultainiously Diagnolizing", time()-St);
# attemp diagnolizing on previous eigenvectors
temp := simplify(MMM(MMM(iP1, C2), P1), factor);
if Equal(temp, Matrix(temp, shape=diagonal)) then
printf("\n \t %8.3f \t Operator diagonalizes on same eigenvectors.", time()-St);
lambda2 := LinearAlgebra:-Diagonal(temp);
for i from 1 to cDim(lambda) do
if Equal(lambda[1..-1,i], lambda2) then
printf("\n \t %8.3f \t Eigenvalues duplicate operator #%d--No new information.",
time()-St, i);
return();
end if;
end do;
lambda := Matrix([lambda, lambda2]);
used_C := Vector([used_C, [C2]]);
else
# Further diagnolize
lambda2, P2 := LinearAlgebra:-Eigenvectors(temp);
if false in convert(map(type, lambda2, realcons), set) then
printf("\n \t %8.3f \t ERROR: Complex eigenvalues", time()-St);\
return(lambda2, P2, temp);
end if;
Pt := P1.P2;
iPt := LinearAlgebra:-MatrixInverse(Pt);
temp := simplify( MMM(MMM(iPt, used_C[1]), Pt), factor);
if not Equal(temp, Matrix(temp, shape=diagonal)) then
printf("\n \t %8.3f \t ERROR: Commutes, yet cannot simultaneously diagnolize",\
time()-St);
return(C2);
end if;
# reorder eigenvectors.
printf("\n \t %8.3f Correcting Ordering of Eigenvectors", time()-St);\
temp := Transpose(LinearAlgebra[LUDecomposition](P2, output='`P`'));\
Pt := MMM(Pt, temp);
iPt := MMM(LinearAlgebra:-MatrixInverse(temp), iPt);
temp := simplify(MMM(MMM(iPt, used_C[1]), Pt), factor);
if not Equal(LinearAlgebra:-Diagonal(temp), lambda[1..-1, 1]) then\
printf("\n \t %8.3f \t ERROR: Cannot order Eigen vectors", time()-St);\
return(LinearAlgebra:-Diagonal(temp), lambda, C2);
end if;
lambda2 := simplify(LinearAlgebra:-Diagonal(MMM(MMM(iPt, C2), Pt)), factor);\
lambda := Matrix([lambda, lambda2]);
P1 := simplify(Pt, factor);
iP1 := simplify(iPt, factor);
used_C := Vector([used_C, [C2]]);
end if;
printf("\n \t %8.3f Complete", time()-St);
return(lambda, P1);
finally:
St_ := time()-St;
status := "additional Casimir complete";
end:
end proc:
# Computes Operator acting spanning set Delta
Get_C := proc( in1 :: {list, set}, in2 :: posint )
local C_op, cop, C_, i, d;
if whattype( in1 ) = list then C_op := UnBlock( in1 )
else C_op := Couple( in1, in2 )
end if;
if status[-1] = "c" then
C_ := Matrix( Dimension(Delta_c) );
d := nops(C_op[1]);
for i from 1 to rDim(C_) do
if time()-St_ > t_thres then
printf("\n \t %8.3f \t Basis element %d of %d",
time()-St, i, rDim(C_)); St_ := time();
end if;
for cop in C_op do
C_[i,1..-1] := C_[i,1..-1] +
Apply( [seq(cop[-q], q=1..d)], [1], [Delta_c[i]] );
end do:
end do;
elif status[-2..-1] = "1r" then
C_ := Matrix( Dimension(Delta_1r) );
d := nops(C_op[1]);
for i from 1 to rDim(C_) do
if time()-St_ > t_thres then
printf("\n \t %8.3f \t Basis element %d of %d",
time()-St, i, rDim(C_)); St_ := time();
end if;
for cop in C_op do
C_[i,1..-1] := C_[i,1..-1] +
Apply( [seq(cop[-q], q=1..d)], [1], [Delta_1r[i]] );
end do:
end do;
else
error("Unknown status: %1", status);
end if;
return(C_);
end proc:
# Converts L_ops from Irrep space indices to variable indices
UnBlock := proc( Lops_in :: list )
local nonzero, shifts, limits, Lops, Lops_out, X, temp, i;
if status[1..8] = "coupling" then
Lops := Lops_in;
else
if type( Lops_in, list(integer)) then Lops := [[Lops_in]]
elif type( Lops_in, list(list(integer))) then Lops := [Lops_in]
elif type( Lops_in, list(list(list(integer)))) then Lops := Lops_in\
else error("Lops must of type list(list(list(integer)))");
end if;
if false in { seq(seq( evalb(Lops[i][j][2] = Lops[i][j+1][1]),
j=1..nops(Lops[i])-1), i=1..nops(Lops)) } or false in
{seq( evalb(Lops[i][1][1] = Lops[i][-1][2]), i=1..nops(Lops))} then\
error(" Invalid Casimir Operator ");
end if;
end if;
nonzero := map(X->nops(select(`>`, X, 0)), tensors);
shifts := [0, seq(convert(nonzero[1..i], `+`), i=1..nops(nonzero))];
limits := [seq(shifts[i]+1..shifts[i+1], i=1..nops(shifts)-1)];
Lops_out := [];
for i from 1 to nops(Lops) do
X := combinat:-cartprod( map(x -> {$ limits[x]},
[seq(Lops[i][j][1], j=1..nops(Lops[i]))]));
while not X[finished] do
temp := X[nextvalue]();
Lops_out := [op(Lops_out), zip( (i,j) -> Wrd*(i-1) + j, temp[1..-1],\
[op(temp[2..-1]), temp[1]]) ];
end do;
end do;
return(Lops_out);
end proc:
# Returns operator representing full coupling
Couple := proc( in1::{list({set(integer), integer}), set(integer)},in2::integer )\
local spaces, power, d, A, T, i;
if nargs = 2 then spaces := in1; power := in2;
elif nargs = 1 and whattype(in1) = set then error("Need set and power")\
elif nargs = 1 then spaces, power := op(in1);
else error("One or Two inputs required")
end if;
d := nops(spaces)^power;
A := Vector(d);
T := combinat[cartprod]( [ spaces $ power ] ):
for i from 1 to d while not T[finished] do
A[i] := T[nextvalue]();
end do;
LinearAlgebra:-Map(x->[seq([x[i], x[i+1]], i=1..nops(x)-1), [x[-1], x[1]]],\
A, inplace);
if status[1..13] = "Casimir (CGC)" then
A := convert(map(X->map(x -> Wrd*(x[1]-1) + x[2], X), A), list);
elif status[1..7] = "Casimir" then
try:
status := "coupling" || status;
A := UnBlock(convert(A, list));
finally:
status := status[9..-1];
end:
end if;
return(A);
end proc:
EV := proc(k :: posint, m :: list) local n, l, all, Gamma;
n := nops(m);
l := [seq(m[i]+(n-i), i=1..n)];
all := {seq(i, i=1..n)};
Gamma := [seq( mul(1-1/(l[i]-l[j]), j=1..i-1) * mul(1-1/(l[i]-l[j]),
j=i+1..n), i=1..n)];
return add(Gamma[i]*(l[i])^k, i=1..n);
end proc:
end module:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #\
#
# (3.0) Finds Racah Coefficients
Racah := proc ( P1 :: Matrix, P2 :: Matrix)
local Adjoint, piL_, HW, i, j, temp, Q, iP1, iP2, X, Y, R;
St := time();
if status[1..7] = "initial" then
error("Must run MapHW first")
else
status := "Racah";
end if;
HW := attributes(M_1r);
# Take adjoint
Adjoint := [seq( seq( (j-1)*Wrd + i, j=1..Wrd), i=1..Wrd)];
piL_ := map(X->[seq(Adjoint[X[-i]], i=1..nops(X))], piL_1r);
# Apply Adjoint
printf("\n\t %8.3f Applying Hermition of Operators", time()-St); St_ := time();\
Q := Matrix(nops(piL_1r), Dimension(Delta_1r));
for j from 1 to Dimension(Delta_1r) do
if time()-St_ > t_thres then
printf("\n \t %8.3f \t Basis Element %d of %d", time()-St,
j, nops(Delta_1r_)); St_ := time();
end if;
for i from 1 to nops(piL_) do
temp := Apply( piL_[i], [1], [Delta_1r[j]] );
if temp[1] <> 0 then
if op(temp[2]) = HW then
Q[i,j] := op(temp[1]);
else
error("Adjoint left space");
end if;
end if;
end do:
end do;
Q := MMM(MMM(B, Q), Transpose(M_1r));
printf("\n\t %8.3f Finding Racah Coefficients", time()-St);
iP1 := map(simplify, LinearAlgebra:-MatrixInverse(P1), factor);
iP2 := map(simplify, LinearAlgebra:-MatrixInverse(P2), factor);
temp := simplify(MMM(iP1, Q), factor);
R := MMM(temp, Transpose(iP2), outputoptions=[datatype=realcons]);
X := Vector(rDim(temp), i -> simplify(LinearAlgebra:-DotProduct(temp[i,1..-1],\
iP1[i,1..-1])), datatype=realcons);
Y := Vector(rDim(temp), i -> simplify(LinearAlgebra:-DotProduct(MMM(iP2, Q)[i,1..-1],\
iP2[i,1..-1])), datatype=realcons);
try:
R := Matrix(rDim(temp), (i,j) -> R[i,j]*(X[i])^(-1/2)*(Y[j])^(-1/2),
datatype=realcons);
catch:
PETE1 := R, X, Y, Q;
print("ERROR: Division by zero");
return(NULL);
end:
R := map(x->combine(simplify(x)), R);
printf("\n\t %8.3f Racah Coefficients found:", time()-St);
return(R);
end proc:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # (4.0) Computes Clebsch-Gordan Coefficients\
CGC := module()
export Compute, Rotate;
local Check_Tabs, MapVec, BreakTables, Break, norm_sq, cDET, DET, ri, ci,
h1, h2, Delta_2r, Delta_1c, Delta_2c, cZ_1, piZ_1, cZ_2, piZ_2,
cZ_2_, piZ_2_, R_Z_1, Da, Norm1, Norm2;
global IrrepTableau, ProductTableau;
# (4.0.0) Computes Clebsch-Gordan Coefficients given sufficent Tableau
Compute := proc(P)
local N, p, irrep_same, product_same, i, j, temp, iP, d1, d2, D_1, D_2;
St := time();
if status = "initialized" then error("Must run MapHW first");
else status := "CGC";
end if;
# Check Tableau
Check_Tabs[do_it]();
# Generate Delta_1c if necessary
status := "CGC: 1c";
if assigned(Delta_1c) and attributes(Delta_1c) = IrrepTableau then
printf("\n\t %8.3f Using previous basis element labeled by IrrepTableau",
time()-St);
irrep_same := true;
else
printf("\n\t %8.3f Determining basis element labeled by IrrepTableau",
time()-St);
St_ := time();
irrep_same := false;
h1, Delta_1c := MapVec(IrrepTableau);
setattribute(Delta_1c, IrrepTableau);
end if;
# Generate Delta_2r (is always Highest Weight)
# Shifted?
status := "CGC: 2r";
Delta_2r := Vector(nops(tensors)): p := 0;
for i from 1 to nops(tensors) do
NumMinors := tensors[i][1];
Delta_2r[i] := <Gen_HW(tensors[i], p)>;
p := p + nops(tensors[i])-ListTools:-Occurrences(0, tensors[i]);
end do:
# Generate Delta_2c
status := "CGC: 2c";
temp := nops(tensors);
if attributes(Delta_2c) = NULL then
product_same := [false $ temp];
h2 := Vector(temp); Delta_2c := Vector(temp);
else
product_same := zip( (a,b) -> evalb( a = b ), ProductTableau,
attributes(Delta_2c));
end if;
if nops(product_same) < temp then
product_same := [op(product_same), false $ (temp - nops(product_same))];
end if;
# Could make more intellegent in using past information
for i from 1 to temp do
if product_same[i] then
printf("\n\t %8.3f Space %d of Tensor Product:", time()-St, i);
printf("\n\t\t\t using previous Vector");
else
printf("\n\t %8.3f Space %d of Tensor Product:", time()-St, i);
printf("\n\t\t\t Generating Vector corresponding to Tableau");
St_ := time();
h2[i], Delta_2c[i] := MapVec(ProductTableau[i]);
end if;
end do;
setattribute(Delta_2c, ProductTableau);
# Finished mapping
#print(Delta_1r, Delta_1c, h1, Delta_2r, Delta_2c, h2);
if not assigned(cDET) then BreakTables() end if;
# Break apart Delta_1r, Delta_1c;
status := "CGC: Breaking #1";
if not irrep_same then
printf("\n\t %8.3f Breaking minors for mapped Irrep Space Vector", time()-St); St_ := time();
cZ_1, piZ_1 := Break(0);
#print(cZ_1, piZ_1);
end if;
# Break apart Delta_2r, Delta_2c
status := "CGC: Breaking #2";
printf("\n\t %8.3f Breaking & multiplying minors in Product Space", time()-St);
St_ := time(); St_2 := true;
if (not assigned(cZ_2)) or Dimension(cZ_2)<>nops(tensors) then
cZ_2 := Vector(nops(tensors));
piZ_2 := Vector(nops(tensors));
end if;
for i from 1 to nops(product_same) do
if not product_same[i] then
if time()-St_ > t_thres then
printf("\n\t %8.3f Breaking Space %d of %d", time()-St, i, temp);
St_ := time(); St_2 := false;
end if;
cZ_2[i], piZ_2[i] := Break(i);
#print(cZ_2[i], piZ_2[i]);
end if;
end do;
status := "CGC: finishing";
if false in product_same then
# Multiply Vectors in Product
if not St_2 then
printf("\n \t %8.3f Multiplying Vectors of Product Space", time()-St);\
St_ := time();
end if;
piZ_2_ := piZ_2[1];
cZ_2_ := cZ_2[1];
for i from 2 to nops(tensors) do
d1 := Dimension(cZ_2_);
d2 := Dimension(cZ_2[i]);
cZ_2_ := Vector(d1*d2, k->cZ_2_[irem(k-1,d1,'temp')+1]*cZ_2[i][temp+1]);
piZ_2_ := Vector(d1*d2,k->piZ_2_[irem(k-1,d1,'temp')+1]*piZ_2[i][temp+1]);
end do:
end if;
# Compute Normalizations
if time()-St_ > t_thres then
printf("\n \t %8.3f Normalizing Vectors", time()-St); St_ := time();
end if;
D_1 := map(norm_sq, piZ_1);
D_2 := map(norm_sq, piZ_2_);
Norm1:= MMM(MMM(M_1r, MMM(cZ_1, MMM(Matrix(D_1, shape=diagonal),
Transpose(cZ_1)))), Transpose(M_1r));
Norm2 := add(cZ_2_[k]^2*D_2[k], k=1..Dimension(D_2))^(1/2);
# Compute Overlap
printf("\n\t %8.3f Computing Overlap", time()-St); St_ := time();
Da := Vector(Dimension(piZ_1));
for i from 1 to Dimension(cZ_2_) do
if assigned(R_Z_1[piZ_2_[i]]) then
Da[R_Z_1[piZ_2_[i]]] := D_1[R_Z_1[piZ_2_[i]]]*cZ_2_[i];
end if;
end do;
# evaluate CGC
printf("\n\t %8.3f Evaluating Coefficient", time()-St); St_ := time();
if Norm2 = 0 then
error("%1 ERROR: Product Space Vectors have norm of zero", time()-St);
end if;
Da := M_1r.(cZ_1.Da);
if nargs = 0 then
temp := zip( (x,y) -> x / (y^(1/2)*Norm2), Da,
LinearAlgebra:-Diagonal(Norm1) );
else
iP := P^(-1);
temp := zip( (x,y) -> x / (y^(1/2)*Norm2), iP.Da,
LinearAlgebra:-Diagonal(MMM(MMM(iP, Norm1), Transpose(iP))));
end if;
printf("\n\t %8.3f Complete", time()-St);
status := "CGC: complete";
return(map(combine@simplify, simplify(temp, factor)));
end proc:
# (4.1) Checks tableau meet required conditions
Check_Tabs := module()
export do_it;
local between, check;
do_it := proc()
local temp;
if not assigned(IrrepTableau) then
error("Must define IrrepTableau")
elif not assigned(ProductTableau) then
error("Must define ProductTableau")
elif whattype(IrrepTableau) <> list then
error("IrrepTableau must be of type list")
elif whattype(ProductTableau) <> list then
error("ProductTableau must be of type list")
end if;
# Check betweeness
temp := check(IrrepTableau);
if temp then
error("Number of entries in each row of IrrepTableau incorrect")
elif not temp then
error("IrrepTableau does not meet betweeness and dominance conditions.")\
end if;
temp := map(check, ProductTableau);
if true in temp then
error("Number of entries in rows of ProductTableau incorrect");
elif false in temp then
error("ProductTableau do not meet betweeness and dominance conditions");\
end if;
# Check in same spaces as Mapped in MapHW();
if irrep[1..nops(IrrepTableau[1])] <> IrrepTableau[1] then
error("IrrepTableau not in same space as mapped irrep.")
elif nops(tensors) <> nops(ProductTableau) then
error("Not same number of Product Tableau as there are irreps in the Tensor Product")\
elif false in zip( (X, Y) -> evalb(X = Y[1]), tensors, ProductTableau) then\
error("Not all Product Tableau in correct irrep space.")
end if;
# Check if CGC zeros
temp := map(x->convert(ListTools[Flatten](x), `+`),
ListTools:-Transpose(ProductTableau));
if temp <> map(convert, IrrepTableau, `+`) then
printf("\n Weights of Irrep Tableau do not match that of ProductTableau.");\
error("CGC identically zero");
end if;
end proc:
check := proc(gzTab);
if map(nops, gzTab) <> [seq(nops(gzTab[1])-i+1, i=1..nops(gzTab[1]))] then\
return(true)
elif false in zip(between, gzTab[1..-2], gzTab[2..-1]) then
return(false)
else
return(FAIL)
end if;
end proc:
between := (a, b) -> not (false in
{op(zip((x,y) -> evalb(x >= y), a, b)),
op(zip((x,y) -> evalb(x <= y), a[2..-1], b))}):
end module:
# (4.2)
MapVec := proc( gzTab :: list )
local weight, temp, tabs, piL, HW, M_c, i, j, k, c, tabs_, ROW, cSet, ev, iM,\
h, C, lambda, lambda_, P;
status := cat("MapVec: ", status[-2..-1]);
temp := map(convert, gzTab, `+`);
weight := [temp[-1], seq(temp[-i-1]-temp[-i], i=1..nops(temp)-1)];
NumMinors := gzTab[1][1];
# if highest weight vector, then just return it.
if weight = gzTab[1] then
return( Vector[row]([1]), Vector[column]([[ Gen_HW(gzTab[1], 0) ]]) )\
end if;
# Get tableau & Operators from tableau
tabs := Gen_Tabs(gzTab[1], weight);
if nops(tabs)=0 then error("Could not recreate tableau %1",time()-St) end if;\
piL := map(T->[seq(seq( ((nops(T[1])-i)*Wrd+j) $ T[i][j]-T[i+1][j],
j=1..nops(T[i+1])), i=1..nops(T[1])-1)], tabs);
# Generate Highest Weight Vector & tuples;
HW := Gen_HW(irrep, 0);
# Apply piL
unassign('R_Delta_c'): unassign('Delta_c'): k := 0:
for i from 1 to nops(piL) do
c, temp := Apply( piL[i], [1], [HW] ):
if c = 0 then
C[i] := Vector[row]([0]);
else
for j from 1 to nops(c) do
if not assigned(R_Delta_c[temp[j]]) then
k := k + 1;
R_Delta_c[temp[j]] := k;
Delta_c[k] := temp[j];
end if;
end do:
C[i] := Vector[row](k, datatype=numeric):
for j from 1 to nops(c) do
C[i][R_Delta_c[temp[j]]] := c[j];
end do
end if;
end do:
Delta_c := Vector(k, Delta_c):
M_c := Matrix(nops(piL), k, datatype=integer, storage=sparse):
for i from 1 to nops(piL) do
M_c[i, 1..-1] := C[i];
end do:
# If only one Vector with same weight, return values
if nops(piL) = 1 then
return(convert(M_c, Vector[row]), Delta_c);
end if;
status := "Casimir (CGC) " || (status[-2..-1]);
# Determine Casimir Operators Needed to Break Degeneracy
tabs_ := map(x->subsop(1=NULL, -1=NULL, x), tabs);
cSet := Vector(nops(tabs_[1]), fill=1):
for i from nops(tabs_[1]) to 1 by -1 while nops(tabs_) > 1 do
ROW := convert({seq(tabs_[j][i], j=1..nops(tabs_))}, list);
if nops(ROW) = 1 then next end if;
j := 2; ev := [0, 0];
while nops(ListTools:-MakeUnique(ev)) < nops(ev) do
if j > nops(tabs_) - i + 1 then
error("Casimir didn't break degeneracy") end if;
ev := map2(Casimir[EV], j, ROW);
cSet[i] := cSet[i] + 1;
j := j + 1;
end do;
end do:
iM := LinearAlgebra:-MatrixInverse(Matrix(M_c, datatype=anything),
method=pseudo);
h := 1;
for i from Dimension(cSet) to 1 by -1 do
for j from 2 to cSet[i] do
C := MMM(MMM( M_c, Casimir[Get_C]( {seq(k, k=1..nops(tabs_[1])+2-i)}, j ) ), iM);\
lambda_, P := LinearAlgebra:-Eigenvectors(C);
P := LinearAlgebra:-MatrixInverse(P);
lambda := Casimir[EV]( j, gzTab[i+1] );
temp := zip( (x,i) -> if x = lambda then i end if,
convert(lambda_, list), [ $ 1..Dimension(lambda_) ] );
P := P[temp, 1..-1];
if h = 1 then h := P
else h := MMM(P, h);
end if;
M_c := MMM(P, M_c);
iM := LinearAlgebra:-MatrixInverse(Matrix(M_c, datatype=anything),\
method=pseudo);
end do;
end do;
if rDim(M_c) <> 1 then error("Degeneracy not broken") end if;
# Remove entries with zero coef
temp := zip( (a,b) -> if b=0 then NULL else a end if, [$1..cDim(M_c)],
[seq(M_c[1,q], q=1..cDim(M_c))]);
if temp = [$1..cDim(M_c)] then
M_c := M_c[1, 1..-1];
else
M_c := Vector[row]([seq(M_c[1,q], q in temp)]);
Delta_c := Vector([seq(Delta_c[q], q in temp)]);
end if;
status := "CGC: mapped vector " || (status[-2..-1]);
return(M_c, Delta_c);
end proc:
# (4.5)
BreakTables := proc()
local par, s, a, p;
par := table(antisymmetric):
for s in Sizes do
a := [$ 1..s];
p := combinat:-permute(a):
cDET[s] := Vector(s!, i-> signum(par[op(p[i])] / par[op(a)]),
datatype=integer[1]);
DET[s] := Vector(s!, i-> mul(ithprime((ri[a[k]]-1)*nops(IrrepTableau[1])\
+ ci[p[i][k]]), k=1..s));
end do:
end proc:
# (4.6)
Break := proc(min_k :: integer)
local spots, Rows, Cols, h, temp, divs, i, j, k, l, s, ris, cis, ctemp, ztemp, i2, j2,
c_z, z, d1, d2, cZ, piZ, R_piZ, perm_cis, n_fact, ctemp_a, ztemp_a, ctemp_b, ztemp_b, pcis, cB, B;
spots := ListTools:-PartialSums(map2(combinat:-numbcomb, Wrd, Sizes));
if min_k = 0 then
NumMinors := irrep[1];
h := h1;
Rows := map(DecodeO, Delta_1r);
Cols := map(DecodeO, Delta_1c);
else
NumMinors := tensors[min_k][1];
h := h2[min_k];
Rows := map(DecodeO, Delta_2r[min_k]);
Cols := map(DecodeO, Delta_2c[min_k]);
end if;
# divide into sizes
temp := convert(Rows[1], list);
divs := [0,seq(nops(select(x-> x <= spots[i], temp)), i=1..nops(spots))];
temp := convert(Cols[1], list);
if divs <> [0,seq(nops(select(x-> x <= spots[i], temp)), i=1..nops(spots))] then
PETE1 := Rows, Cols; error("Size error, %1", status)
end if;
l := 0;
for i from 1 to Dimension(Rows) do
for j from 1 to Dimension(Cols) do
unassign('c_z'): unassign('z'):
for k from 1 to nops(divs)-1 do #each size minor
if divs[k]=divs[k+1] then next end if; # doesn't include any of that size minor
ris := convert(Rows[i][divs[k]+1..divs[k+1]], list);
cis := convert(Cols[j][divs[k]+1..divs[k+1]], list);
s := add(MinorsA[ris[1],q], q=1..Wrd);
if nops(ris) = 1 then
# temp converts from binary to index list
temp := [seq( i $ MinorsA[ris[1], -i], i=1..Wrd)], [seq( i $ MinorsA[cis[1], -i], i=1..Wrd)];
ctemp := subs( ri=temp[1], ci=temp[2], cDET[s]);
ztemp := subs( ri=temp[1], ci=temp[2], DET[s]);
else
# Can improve this ruitine greatly
if nops(ris)<>nops(cis) then error("Bizarre Error in Break\n") end if;
ctemp := [];
ztemp := [];
# Consider all permutations and fraction
perm_cis := combinat:-permute(cis);
n_fact := nops(cis)!;
# for each permutation, multiply together, and sum over permutations
for pcis in perm_cis do
# break apart first pair
temp := [seq( i $ MinorsA[ris[1], -i], i=1..Wrd)], [seq( i $ MinorsA[pcis[1], -i], i=1..Wrd)];
ctemp_a := subs( ri=temp[1], ci=temp[2], cDET[s]);
ztemp_a := subs( ri=temp[1], ci=temp[2], DET[s]);
# now multiply by each next pair.
for i2 from 2 to nops(ris) do
# convert next pair
temp := [seq( i $ MinorsA[ris[ i2 ], -i], i=1..Wrd)], [seq( i $ MinorsA[pcis[ i2 ], -i],
i=1..Wrd)];
ctemp_b := subs( ri=temp[1], ci=temp[2], cDET[s]);
ztemp_b := subs( ri=temp[1], ci=temp[2], DET[s]);
# multiply
d1 := Dimension(ctemp_a);
d2 := Dimension(ctemp_b);
ctemp_a := Vector(d1*d2, k->ctemp_a[irem(k-1,d1,'temp')+1]*ctemp_b[temp+1]);
ztemp_a := Vector(d1*d2, k->ztemp_a[irem(k-1,d1,'temp')+1]*ztemp_b[temp+1]);
end do;
# Add results
if ctemp = [] then
ctemp := (1/n_fact)*ctemp_a;
ztemp := convert(ztemp_a, list);
else
# Need to create unique basis to add
B := sort(ListTools:-MakeUnique([ op(convert(ztemp_a, list)), op(ztemp) ]));
cB := Vector(nops(B)):
for i2 from 1 to Dimension( ctemp ) do
j2 := ListTools:-BinarySearch(B, ztemp[i2]):
if j2=0 then error("Data loss") end if:
cB[j2] := cB[j2] + ctemp[i2]:
end do:
for i2 from 1 to Dimension( ctemp_a ) do
j2 := ListTools:-BinarySearch(B, ztemp_a[i2]):
if j2=0 then error("Data loss") end if:
cB[j2] := cB[j2] + ctemp_a[i2] / n_fact:
end do:
# Check for zeros?
ctemp := cB;
ztemp := B;
end if;
end do;
ztemp := Vector(ztemp);
end if;
# Multiply each size minor.
if not assigned(c_z) then
c_z := ctemp;
z := ztemp;
else
d1 := Dimension(c_z);
d2 := Dimension(ctemp);
c_z := Vector(d1*d2, k->c_z[irem(k-1,d1,'temp')+1]*ctemp[temp+1]);
z := Vector(d1*d2, k-> z[irem( k-1, d1, 'temp')+1]*ztemp[temp+1]);
end if;
end do; # k
for k from 1 to Dimension(z) do
if assigned(R_piZ[z[k]]) then
cZ[i,R_piZ[z[k]]] := cZ[i,R_piZ[z[k]]] + h[j]*c_z[k];
else
l := l + 1;
cZ[i,l] := h[j]*c_z[k];
piZ[l] := eval(z[k]);
R_piZ[z[k]] := l;
end if;
end do;
for k from 1 to l do
cZ[i+1,k] := 0;
end do:
end do;
end do;
if status[-1] = "1" then
R_Z_1 := R_piZ;
return(Matrix(i-1, l, cZ), Vector(l, [seq(piZ[i], i=1..l)]));
else
return(Vector(l, [seq(cZ[1,i],i=1..l)]), Vector(l, [seq(piZ[i],i=1..l)]));
end if;
end proc:
# (4.6)
norm_sq := proc(x) local w;
w := [op(ifactor(x))];
w := [seq([op(w[i])], i=1..nops(w))];
convert(map( A-> if nops(A) > 1 then A[2] end if, w), `*`);
end proc:
# (4.7) Apply different coupling scheme
Rotate := proc(P) local iP, temp;
iP := P^(-1);
temp := zip( (x,y) -> x / (y^(1/2)*Norm2), iP.Da,
LinearAlgebra:-Diagonal(MMM(MMM(iP, Norm1), Transpose(iP))));
return(map(combine@simplify, simplify(temp, factor)));
end proc:
end module:
end module:
with(TabProcs);
# Do datatype=integer[]
# changes to paper: ordering of Minors, Prime Numbers
# exp > 2
# Irrep := [3,3,0]; Tensors := [[4,0,0],[2,0,0]];
# IrrepTableau := [[4,1],[2]];
# ProductTableau := [ [[1,0],[0]], [[1,0],[1]], [[1,0],[0]], [[1,0],[0]], [[1,0],[1]] ];