r1 := [[1, 2, 3, 4, 5, 6, 7, 8]];
r2 := [[9, 10, 11, 12, 13, 14, 15, 16]];
r3 := [[17, 18, 19, 20, 21, 22, 23, 24]];
r4 := [[25, 26, 27, 28, 29, 30, 31, 32]];
f5 := [[1,28],[2,27],[3,26],[4,25],[9,20],[10,19],[11,18],[12,17]];
f4 := [[2,29],[3,28],[4,27],[5,26],[10,21],[11,22],[12,23],[13,24]];
f3 := [[3,30],[4,29],[5,28],[6,27],[11,22],[12,21],[13,20],[14,19]];
f2 := [[4,31],[5,30],[6,29],[7,28],[12,23],[13,22],[14,21],[15,20]];
f1 := [[5,32],[6,31],[7,30],[8,29],[13,24],[14,23],[15,22],[16,21]];
f8 := [[6,25],[7,32],[8,31],[1,30],[14,17],[15,24],[16,23],[9,22]];
f7 := [[7,26],[8,25],[1,32],[2,31],[15,18],[16,17],[9,24],[10,23]];
f6 := [[8,27],[1,26],[2,25],[3,32],[16,19],[9,18],[10,17],[11,24]];

order_of_group := 437763136697395052544000000;
#divide this by 2^20 to get the number of distinct positions

perm_ball 
:= permgroup(32,{
r1=[[1, 2, 3, 4, 5, 6, 7, 8]], 
r2=[[9, 10, 11, 12, 13, 14, 15, 16]], 
r3=[[17, 18, 19, 20, 21, 22, 23, 24]], 
r4=[[25, 26, 27, 28, 29, 30, 31,32]], 
f5=[[1,28],[2,27],[3,26],[4,25],[9,20],[10,19],[11,18],[12,17]], 
f4=[[2,29],[3,28],[4,27],[5,26],[10,21],[11,22],[12,23],[13,24]], 
f3=[[3,30],[4,29],[5,28],[6,27],[11,22],[12,21],[13,20],[14,19]], 
f2=[[4,31],[5,30],[6,29],[7,28],[12,23],[13,22],[14,21],[15,20]], 
f1=[[5,32],[6,31],[7,30],[8,29],[13,24],[14,23],[15,22],[16,21]], 
f8=[[6,25],[7,32],[8,31],[1,30],[14,17],[15,24],[16,23],[9,22]], 
f7=[[7,26],[8,25],[1,32],[2,31],[15,18],[16,17],[9,24],[10,23]], 
f6=[[8,27],[1,26],[2,25],[3,32],[16,19],[9,18],[10,17],[11,24]]});
#this is "the group" of the masterball - it does not act freely

flip_poles := [[1, 25], [2, 26], [3, 27], [4, 28], [5, 29], 
[6, 30], [7, 31], [8, 32]];
#this group element flips the facelets from the north pole
#to the south pole, leaving the middle facelets intact
#It *does* belong to perm_group (and obviously acts trivially
#on the solved masterball). It also belongs to the derived subgroup.

flip_wide_bands := [[9, 17], [10, 18], [11, 19], [12, 20], 
[13, 21], [14, 22], [15, 23], [16, 24]];
#this group element flips the facelets from the upper middle band
#to the lower middle band, leaving the poles intact
#It *does* belong to perm_group (and obviously acts trivially
#on the masterball). It also belongs to the derived subgroup.

long2cycle:=
[f1,f3,r1^(-1),f3,f1,r4^(-1),r4^(-1),f1,f3,r1,f3,f1,r4^(-1),f1,
f3,r1^(-1),f3,f1,r4^(-1),f1,f3,r1^(-1),f3,f1,r4^(-1),f1,
f3,r1,r1,f3,f1,r4^(-1),f1,f3,r1^(-1),f3,f1,r4,f1,
f3,r1,f3,f1,r4^(-1),r4^(-1),r4^(-1),f1,f3,r1^(-1),f3,f1,r4^(-1),f1,
f3,r1^(-1),f3,f1,r4^(-1),f1,f3,r1,r1,f3,f1,r4^(-1),r4^(-1),r4^(-1),
f1,f3,r1^(-1),f3,f1,r4^(-1),f1,f3,r1^(-1),f3,f1,
r4^(-1),f1,f3,r1,f3,f1,r4^(-1),f1,f3,r1,f3,f1,r4,
f1,f3,r1^(-1),r1^(-1),f3,f1,r4,f1,f3,r1,f3,f1,r4,f1,
f3,r1,f3,f1,r4,f1,f3,r1^(-1),f3,f1,r4,f1,
f3,r1,f3,f1,r4,f1,f3,r1^(-1),r1^(-1),f3,f1,r4,f1,
f3,r1,f3,f1,r4,f1,f3,r1,f3,f1,r4,f1,
f3,r1^(-1),f3,f1,r4,r4,f1,f3,r1,f3,f1,
r4^(-1),r4^(-1),r4^(-1),f1,f3,r1^(-1),f3,f1,r4^(-1),f1,f3,r1^(-1),
f3,f1,r4^(-1),f1,f3,r1,r1,f3,f1,r4^(-1),r4^(-1),r4^(-1),
f1,f3,r1^(-1),f3,f1,r4^(-1),f1,f3,r1^(-1),f3,
f1,r4^(-1),f1,f3,r1,f3,f1,r4^(-1),f1,f3,r1,
f3,f1,r4,f1,f3,r1^(-1),r1^(-1),f3,f1,r4,f1,f3,
r1,f3,f1,r4,f1,f3,r1,f3,f1,r4,f1,f3,
r1^(-1),f3,f1,r4,r4,f1,f3,r1^(-1),f3,f1,
r4,f1,f3,r1,f3,f1,r4,f1,f3,r1,f3,f1,
r4,r4,r4,f1,f3,r1^(-1),r1^(-1),f3,f1,r4,f1,f3,r1,f3,
f1,r4,f1,f3,r1,f3,f1,r4,r4,r4,f1,f3,r1^(-1),
f3,f1,r4^(-1),f1,f3,r1,f3,f1,r4,f1,f3,
r1^(-1),r1^(-1),f3,f1,r4,f1,f3,r1,f3,f1,r4,f1,
f3,r1,f3,f1,r4,f1,f3,r1^(-1),f3,f1,
r4,r4,f1,f3,r1,f3,f1,r4^(-1),r4^(-1),f1,f3,r1^(-1),
f3,f1,r4,f1,f3,r1,f3,f1,r4^(-1),r4^(-1),f1,
f3,r1^(-1),f3,f1,r4^(-1),f1,f3,r1^(-1),f3,f1,
r4^(-1),f1,f3,r1,r1,f3,f1,r4^(-1),r4^(-1),r4^(-1),f1,f3,r1^(-1),
f3,f1,r4^(-1),f1,f3,r1^(-1),f3,f1,r4^(-1),
f1,f3,r1,r1,f3,f1,r4^(-1),f1,f3,r1^(-1),f3,
f1,r4^(-1),f1,f3,r1,f3,f1,r4,r4];
#this manuever evaluates to the 2-cycle [[1,2]]

der_perm_ball := permgroup(32,{
[[1, 26, 30], [2, 6, 25], [3, 7, 32], [8, 27, 31], 
[9, 18, 22], [10, 14, 17], [11, 15, 24], 
[16, 19, 23]], 
[], 
[[1, 5, 30, 7, 32], [6, 31, 8, 25, 29], [9, 13, 22, 15, 24], 
[14, 23, 16, 17, 21]], 
[[9, 22, 23, 24, 17, 18, 14, 15, 16]], 
[[1, 30, 5], [4, 8, 31], [6, 29, 25], [7, 28, 32], [9, 22, 13], 
[12, 16, 23], [14, 21, 17], [15, 20, 24]], 
[[1, 30, 31, 32, 25, 26, 6, 7, 8]]});
#this is the derived group of perm_ball

order_of_der_perm_ball := 109440784174348763136000000;
#this is 1/4th order_of_group

der_2_perm_ball := permgroup(32,{
[], [[1, 6, 7, 8, 30, 2, 3], [25, 31, 32, 26, 27]], 
[[1, 26], [2, 25], [3, 32], [4, 31], [5, 30], [6, 29], 
[7, 28], [8, 27], [9, 18], [10, 17], [11, 24], [12, 23], 
[13, 22], [14, 21], [15, 20], [16, 19]], 
[[1, 30, 7, 32, 5, 26, 3], [2, 27, 25, 6, 31, 8, 29], 
[9, 22, 15, 24, 13, 18, 11], [10, 19, 17, 14, 23, 16, 21]], 
[[9, 14, 15, 16, 22, 10, 11], [17, 23, 24, 18, 19]]});
#this is the derived group of the derived group of perm_ball

order_of_der_2_perm_ball := 109440784174348763136000000;
#its order is the same as der_perm_ball

mod8 := proc (j) 
local k, l; 
k := `mod`(j,8); if k = 0 then l := 8 else l := k
fi; RETURN(l) end;

mod84_row := proc (i, j, n) 
local temp, tempo; 
if i = 1 then temp := 4 else if
i = 2 then temp := 3 else if i = 3 then 
temp := 2 else if i = 4 then temp := 1
fi fi fi fi; 
if j = n or j = mod8(n+1) or j = mod8(n+2) or j = mod8(n+3) then
tempo := temp else temp := i fi; 
RETURN(temp) 
end;

mod84_col := proc (i, j, n) 
local temp; 
if j = n then temp := mod8(n+3) else 
if j = mod8(n+1) then temp := mod8(n+2) else if j = mod8(n+2) then 
temp := mod8(n+1) else if j = mod8(n+3) then temp := n else 
temp := j 
fi fi fi fi; 
RETURN(temp) 
end;

move_long := proc (A, i) 
local j, k, AA; 
AA:=array(1..4,1..8);
#print(`MOVE_LONG`);
for k to 8 do 
for j to 4 do 
AA[j,k] := A[mod84_row(j,k,i),mod84_col(j,k,i)] 
od 
od; 
#print(AA,`AAMATRIX_IN_MOVE_LONG`);
RETURN(AA) 
end;

color1 := red;
color2 := blue;
color3 := coral;
color4 := plum;
color5 := green;
color6 := yellow;
color7 := cyan;
color8 := COLOR(RGB,.5607,.7372,.5607);

master_squares := proc () 
local i; 
for i from 0 to 3 do 
point(A1.i,0,i), point(B1.i,1,i), point(C1.i,1,1+i), point(D1.i,0,1+i); 
square(Sq.1.i,[A1.i, B1.i, C1.i, D1.i]); 
point(A2.i,1,i), point(B2.i,2,i), point(C2.i,2,1+i), point(D2.i,1,1+i); 
square(Sq.2.i,[A2.i, B2.i, C2.i, D2.i]); 
point(A3.i,2,i), point(B3.i,3,i), point(C3.i,3,1+i), point(D3.i,2,1+i); 
square(Sq.3.i,[A3.i, B3.i, C3.i, D3.i]); 
point(A4.i,3,i), point(B4.i,4,i), point(C4.i,4,1+i), point(D4.i,3,1+i);
square(Sq.4.i,[A4.i, B4.i, C4.i, D4.i]); 
point(A5.i,4,i), point(B5.i,5,i), point(C5.i,5,1+i), point(D5.i,4,1+i); 
square(Sq.5.i,[A5.i, B5.i, C5.i, D5.i]);
point(A6.i,5,i), point(B6.i,6,i), point(C6.i,6,1+i), point(D6.i,5,1+i); 
square(Sq.6.i,[A6.i, B6.i, C6.i, D6.i]); 
point(A7.i,6,i), point(B7.i,7,i), point(C7.i,7,1+i), point(D7.i,6,1+i); 
square(Sq.7.i,[A7.i, B7.i, C7.i, D7.i]); 
point(A8.i,7,i), point(B8.i,8,i), point(C8.i,8,1+i), point(D8.i,7,1+i); 
square(Sq.8.i,[A8.i, B8.i, C8.i, D8.i]) od 
end;

square_draw := proc (A) 
local L, i, j; 
master_squares(); L := []; for j to 8 
do for i from 0 to 3 do 
L := [op(L), Sq.j.i(color = color.(A[i+1,j]))] od od;
draw(L,filled = true,printtext = false,scaling = constrained,axes = none)
end;

sphere_draw := proc (A::array) 
local i, j, p1, p2, p3, p4, B11, B12, B13, B14,
B15, B16, B17, B18, B21, B22, B23, B24, B25, B26, 
B27, B28, B31, B32, B33, B34, B35, B36, B37, B38, 
B41, B42, B43, B44, B45, B46, B47, B48, L; 
for i to 8 do
p1.i := [8*cos(1/4*i*Pi), 8*sin(1/4*i*Pi), 10/sqrt(2)] od; 
for i to 8 do p2.i:= [10*cos(1/4*i*Pi), 10*sin(1/4*i*Pi), 0] od; 
for i to 8 do p4.i := [8*cos(1/4*i*Pi), 8*sin(1/4*i*Pi), -10/sqrt(2)] od;
B11:=polygonplot3d([p11,p12,[0,0,10]],style=patch,color=color.(A[1,1]));
B12:=polygonplot3d([p12,p13,[0,0,10]],style=patch,color=color.(A[1,2]));
B13:=polygonplot3d([p13,p14,[0,0,10]],style=patch,color=color.(A[1,3]));
B14:=polygonplot3d([p14,p15,[0,0,10]],style=patch,color=color.(A[1,4]));
B15:=polygonplot3d([p15,p16,[0,0,10]],style=patch,color=color.(A[1,5])); 
B16:=polygonplot3d([p16,p17,[0,0,10]],style=patch,color=color.(A[1,6]));
B17 := polygonplot3d([p17,p18,[0,0,10]],style=patch,color=color.(A[1,7]));
B18:=polygonplot3d([p18, p11, [0, 0, 10]],style = patch,color = 
color.(A[1,8])); 
B21 := polygonplot3d([p11, p12, p22, p21],style = patch,color
= color.(A[2,1])); 
B22 := polygonplot3d([p12, p13, p23, p22],style = patch,
color = color.(A[2,2])); 
B23 := polygonplot3d([p13, p14, p24, p23],style = patch,
color = color.(A[2,3])); 
B24 := polygonplot3d([p14, p15, p25, p24],style= patch,
color = color.(A[2,4])); 
B25 := polygonplot3d([p15, p16, p26, p25],style = patch,
color = color.(A[2,5])); 
B26 := polygonplot3d([p16, p17, p27, p26],style = patch,
color = color.(A[2,6])); 
B27 := polygonplot3d([p17, p18, p28, p27],style = patch,
color = color.(A[2,7])); 
B28 := polygonplot3d([p18, p11, p21, p28],style = patch,
color = color.(A[2,8])); 
B31 := polygonplot3d([p41, p42, p22, p21],style = patch,
color = color.(A[3,1])); 
B32 := polygonplot3d([p42, p43, p23, p22],style = patch,
color = color.(A[3,2])); 
B33:= polygonplot3d([p43, p44, p24, p23],style = patch,
color = color.(A[3,3])); 
B34 := polygonplot3d([p44, p45, p25, p24],style = patch,
color = color.(A[3,4])); 
B35 := polygonplot3d([p45, p46, p26, p25],style = patch,
color = color.(A[3,5])); 
B36 := polygonplot3d([p46, p47, p27, p26],style = patch,
color = color.(A[3,6])); 
B37 := polygonplot3d([p47, p48, p28, p27],style = patch,
color = color.(A[3,7])); 
B38 := polygonplot3d([p48, p41, p21, p28],style = patch,
color = color.(A[3,8])); 
B41 := polygonplot3d([p41, p42, [0, 0, -10]],style = patch,
color = color.(A[4,1])); 
B42 := polygonplot3d([p42, p43, [0, 0, -10]],style =patch,
color = color.(A[4,2])); 
B43 := polygonplot3d([p43, p44, [0, 0, -10]],style = patch,
color = color.(A[4,3])); 
B44 := polygonplot3d([p44, p45, [0, 0,-10]],style = patch,
color = color.(A[4,4])); 
B45 := polygonplot3d([p45, p46, [0, 0, -10]],style = patch,
color = color.(A[4,5])); 
B46 := polygonplot3d([p46,p47, [0, 0, -10]],style = patch,
color = color.(A[4,6])); 
B47 := polygonplot3d([p47, p48, [0, 0, -10]],style = patch,
color = color.(A[4,7])); 
B48 := polygonplot3d([p48, p41, [0, 0, -10]],style = patch,
color = color.(A[4,8])); 
L:= [B21, B22, B23, B24, B25, B26, B27, B28, 
B11, B12, B13, B14, B15, B16, B17,B18, 
B31, B32, B33, B34, B35, B36, B37, B38, 
B41, B42, B43, B44, B45, B46, B47, B48]; 
RETURN(L) 
end;

move_lat_pow := proc (A::array,i,pow) 
local j, k, AA; 
#print(`MOVE_LAT_POW`,i,pow);
AA:=array(1..4,1..8);
for k to 8 do 
for j to 4 do if j = i then AA[j,k] := A[j,mod8(k+pow)] 
else AA[j,k] := A[j,k] fi od od; 
#print(AA,`AAMATRIX_IN_MOVE_LAT_POW`);
RETURN(AA) 
end;

move_lat_power := proc (A::array,i,n) 
local nn, j, BB; 
#print(`MOVE_LAT_POWER`);
BB:=array(1..4,1..8);
if n=0 then RETURN(A); fi;  
nn := mod8(n); 
###print(`POWER`,nn);
BB := move_lat_pow(A,i,nn); 
#print(BB,`BBMATRIX_IN_MOVE_LAT_POWER`); 
RETURN(BB) 
end;

power_move:= proc(g,A::array)
#g should be a power of a generator or 1
#A should be a 4x8 matrix
local i,j,col,row,power,gen,
AA, generators1, generators2, generators; 
AA:=array(1..4,1..8);
generators1 := {r1, r4, r3, r2}; 
generators2 := {f1, f2, f3, f4, f5, f6, f7, f8}; 
generators := `union`(generators1,generators2); 
gen:= op(1,g); 
power := op(2,g);
for i from 1 to 4 do 
if r.i=gen then
row := i; 
AA:=move_lat_pow(A,row,`mod`(power,8)); 
fi;
od; 
for j to 8 do if gen = f.j then 
col := j; 
if `mod`(power,2)<>0 then 
AA := move_long(A,col) else AA:=A; fi; 
fi 
od; 
if g=1 then AA:=A; fi;
RETURN(AA);
end;

move := proc (L, A::array) 
#L should be a list or 1
#A should be a 4x8 matrix
local i,j,k, col,row,power,gen,
DD,EE, length_L,generators1, generators2, generators; 
DD:=array(1..4,1..8);
EE:=array(1..4,1..8);
generators1 := {r1, r4, r3, r2}; 
generators2 := {f1, f2, f3, f4, f5, f6, f7, f8}; 
generators := `union`(generators1,generators2); 
if type(L,`^`) then 
EE:=power_move(L,A); 
fi;
if member(L,generators1) then 
for j from 1 to 4 do 
if L=r.j then 
RETURN(move_lat_pow(A,j,1));
fi; 
od; ###########  for j 
fi;
if member(L,generators2) then 
for k from 1 to 8 do 
if L=f.k then 
RETURN(move_long(A,k)); 
fi;
od;  ########### for k 
fi;
if L=1 then RETURN(A); fi;
if type(L,list) then
length_L := nops(L); 
DD:=A;  #########initializing DD
for i to length_L do 
EE:=eval(DD,1);  ################ new, was EE:=DD;
if type(op(i,L),`^`) then 
DD:=move(op(i,L),EE); ############ new line
fi; #########if type(op(i,L),`^`) then 
if member(op(i,L),generators1) then 
for j from 1 to 4 do 
if op(i,L)=r.j then 
DD := move_lat_pow(EE,j,1);
fi;
od;
fi;
if member(op(i,L),generators2) then 
for k from 1 to 8 do 
if op(i,L)=f.k then 
DD := move_long(EE,k);
fi; #############    if op(i,L)=f.k then 
od; ######    for k 
fi; ###          if member(op(i,L),generators2) then 
EE:=eval(DD,1); ################ new, was EE:=DD;
od;  ################   for i 
fi;
RETURN(EE) 
end;

color_count := proc (A::array) 
local count, rows, cols, i, j, count_okay, k; 
for i to 8 do count[i] := 0 od; 
count_okay := 'true'; rows := rowdim(A); cols:= coldim(A); 
if rows <> 4 or cols <> 8 then 
ERROR(`The input matrix must be 4x8`,rows,cols) fi; 
for i to 4 do for j to 8 do 
if not (type(A[i,j],integer) and 0 < A[i,j] and A[i,j] < 9) then 
ERROR(`The input matrix must be integral with entries in [1,8]`,A[i,j]) fi
od od; 
for i to 4 do for j to 8 do for k to 8 do if A[i,j] =
k then count[k] := 1+count[k] fi od od od; 
for i to 8 do if count[i] <> 4 then
count_okay := 'false' fi od; 
RETURN(count_okay) 
end;

matrix_to_perm1 := proc (A::array) 
#this procedure will not work properly unless the
#cycles are disjoint
local L, i, j, k; 
if color_count(A) <> true then 
ERROR(`Must have 8 colors each occurring 4 times`) fi; 
L := []; for i to 4 do for j to 8 do 
if i = 1 and A[i,j] <> j and not member([j, A[i,j]],L) 
then L := [op(L), [A[i,j], j]] fi; 
if i =2 and A[i,j] <> j and not member([j+8, A[i,j]+8],L) 
then L := [op(L), [A[i,j]+8, j+8]] fi; 
if i = 3 and A[i,j] <> j and not member([j+16, A[i,j]+16],L) then
L := [op(L), [A[i,j]+16, j+16]] fi; 
if i = 4 and A[i,j] <> j and not member([j+24, A[i,j]+24],L) 
then L := [op(L), [A[i,j]+24, j+24]] fi od od; 
RETURN(L) 
end;

color_position := proc (i, k, A) 
local L, j; L := []; for j to 8 do if A[i,j]
= k then L := [op(L), j] fi od; 
RETURN(L) 
end;

double_colors := proc (i, A) 
local j, k, L; L := []; 
for k to 8 do if 1 < nops(color_position(i,k,A)) 
then L := [op(L), k] fi od; 
RETURN(L) 
end;

undouble_colors := proc (MM::array) 
local i, j, k, A_mat, d, col1, col2, col3,
col4, dd1, dd4, dd2, dd3, L, num_double_colors1, 
num_double_colors2; 
A_mat :=matrix(4,8); 
for i to 4 do for j to 8 do A_mat[i,j] := MM[i,j] od od; 
L := [];
num_double_colors1 := nops(double_colors(1,MM)); 
num_double_colors2 := nops(double_colors(2,MM)); 
if nops(double_colors(1,MM)) <> nops(double_colors(4,MM)) 
then 
ERROR(`Colors are not possible for a rainbow matrix`)
fi; 
if nops(double_colors(2,MM)) <> nops(double_colors(3,MM)) then 
ERROR(`Colors are not possible for a rainbow matrix`)
fi; 
if 0 < num_double_colors1 then for d to nops(double_colors(1,MM)) do 
dd1:= op(d,double_colors(1,MM)); 
dd4 := op(d,double_colors(4,MM)); 
col1 := op(1,color_position(1,dd1,MM)); 
col4 := op(1,color_position(4,dd4,MM)); 
L := [op(L), [col1, 24+col4]]; A_mat[1,col1] := MM[4,col4]; 
A_mat[4,col4] := MM[1,col1] od fi; 
if 0 < num_double_colors2 then for d to nops(double_colors(2,MM)) do 
dd2 := op(d,double_colors(2,MM)); 
dd3 := op(d,double_colors(3,MM)); 
col2 := op(1,color_position(2,dd2,MM)); 
col3 := op(1,color_position(3,dd3,MM)); 
L := [op(L), [8+col2, 16+col3]]; A_mat[2,col2] := MM[3,col3]; 
A_mat[3,col3] := MM[2,col2] od fi; 
RETURN(A_mat) 
end;

undouble_colors_perm := proc (MM::array) 
local i, j, k, d, col1, col2, col3, 
col4, dd1, dd4, dd2, dd3, L, num_double_colors1, 
num_double_colors2;
 L := [];
num_double_colors1 := nops(double_colors(1,MM)); 
num_double_colors2 := nops(double_colors(2,MM)); 
if nops(double_colors(1,MM)) <> nops(double_colors(4,MM)) then 
ERROR(`Colors are not possible for a rainbow matrix`)
fi; 
if nops(double_colors(2,MM)) <> nops(double_colors(3,MM)) then 
ERROR(`Colors are not possible for a rainbow matrix`)
fi; 
if 0 < num_double_colors1 then for d to nops(double_colors(1,MM)) do 
dd1:= op(d,double_colors(1,MM)); dd4 := op(d,double_colors(4,MM)); 
col1 := op(1,color_position(1,dd1,MM)); 
col4 := op(1,color_position(4,dd4,MM)); 
L := [op(L), [col1, 24+col4]] od fi; if 0 < num_double_colors2 then 
for d to nops(double_colors(2,MM)) do 
dd2 := op(d,double_colors(2,MM)); dd3 := op(d,double_colors(3,MM)); 
col2 := op(1,color_position(2,dd2,MM)); 
col3 := op(1,color_position(3,dd3,MM)); 
L := [op(L), [8+col2, 16+col3]] od fi;
RETURN(L) 
end;

matrix_to_perm := proc (MM::array) 
local L, rows, cols; 
rows := rowdim(MM); 
cols := coldim(MM); 
if rows <> 4 or cols <> 8 then 
ERROR(`The input matrix must be 4x8`,rows,cols) fi; 
L := [op(undouble_colors_perm(MM)), 
op(matrix_to_perm1(undouble_colors(MM)))]; 
RETURN(L) 
end;

rb := proc (L, A) local AA, LL; 
AA := move(L,A); LL := sphere_draw(AA); 
display3d(LL) end;