with(linalg):
 sgcc:=proc(r,s,n,m) sgcc2:=0;
 for i from 1 to n do : 
 for j from i+1 to n do: 
 for k from 1 to m do: 
 for l from 1 to k-1 do:
 sgcc2:=sgcc2+r[i,k]*s[k,i]*r[j,l]*s[l,j]:od:od:od:od:
 print(sgcc2);
 for i from 1 to n do: 
 for j from i+1 to n do: 
 for k from 1 to m do: 
 for l from 1 to k do:
 sgcc2:=sgcc2+r[i,k]*s[k,j]*r[j,l]*s[l,i]:od:od:od:od:
 print(sgcc2);
 for i from 1 to n do:
 for k from 1 to m do: 
 sgcc2:=sgcc2+r[i,k]*(r[i,k]-1)*s[k,i]^2/2:od:od:
 print(sgcc2); return(sgcc2);end;
 #this computes sgcc, the next is a test example
 #a:=matrix(2,2,[1,1,1,1]);b:=diag(1,1,1,1);sgcc(a,b,2);sgcc(b,a,2);
 #example n not m
 #a:=matrix(2,3,[1,1,0,0,0,1]);b:=matrix(3,2,[1,0,0,1,1,1]);
 #sgcc(b,a,3,2);
 a:=matrix(2,2,[1,2,3,4]);b:=matrix(2,2,[5,6,7,8]);sgcc(a,b,2,2);
 u41:=matrix(4,4,[0,0,1,1,1,0,0,0,0,1,0,0,0,0,1,0]);
 i41:=diag(1,1,1,1);r41:=evalm(u41-i41);v41:=inverse(u41);
 col:=matrix(4,3,[3,0,0,3,0,0,3,0,0,3,0,0]); 
 row:=matrix(3,4,[0,0,0,0,1,1,1,1,1,1,1,1]);
 w:=evalm(row&*v41&*col);
 z41:=matrix(3,3,[0,0,1,10,0,0,0,1,0]);
 a71:=blockmatrix(2,2,u41,col,row,z41);det(a71);
 pc:=charpoly(a71,t);fsolve(pc,t,complex);
 o43:=matrix(4,3,0);o34:=transpose(o43);i31:=diag(1,1,1);
 r71:=blockmatrix(2,2,r41,o43,o34,i31);
 s71:=evalm(a71&*inverse(r71));
 sgcc(r71,s71,7,7);evalm(r71&*s71);
 u7a:=evalm(38*a71^6/3+2*a71^5/3-1*a71^4/3-872*a71^3/3
-1108*a71^2/3-1309*a71/3-713/3);
 det(u7a);
 u7b:=evalm(842*a71^6/3+5072*a71^5/3-6847*a71^4/3
-46061*a71^3/3-34930*a71^2/3-52216*a71/3+2878/3);
 det(u7b);
 u7c:=evalm(4260971*a71^6/3-3124108*a71^5/3
+2290532*a71^4/3-99681839*a71^3/3-46221667*a71^2/3-106722952*a71/3+5811547/3);
 det(u7c);
#these matrices are algebraic integers over the field generated by a
#and are a system of fundamental units with a itself; the actual
#integer matrices giving automorphisms of the dimension group
#are a subgroup of the group they generate; the positive eigenvectors
#of a will be their eigenvectors; 
#since -I has sgcc 0 we don't need to 
#check signs on the eigenvectors; 
#also order of r,s doesn't matter since
#the shift has sgcc 0; also  we can compute sgcc with 
#negative signs and denominators 3 (same as multiplying by 9 mod 4)
i7a:=evalm(a71&*inverse(u7a));i7b:=evalm(a71&*inverse(u7b));
i7c:=evalm(a71&*inverse(u7c));
 sgcc(u7a,i7a,7,7);sgcc(u7b,i7b,7,7);sgcc(u7c,i7c,7,7);
#mod 2 all sgcc2 are 0; next check charpoly irreducible
factor(charpoly(a71,t));
#the following program can check the units are independent 
#modulo squares; the polynomials checks they are algebraic integers
 charpoly(u7a,t);charpoly(u7b,t);charpoly(u7c,t);
 with(numtheory);
 p:=11;
 for i from 1 to p-1 do: 
 p1:=i^7-23*i^4-28*i^3-33*i^2-17*i+1:
 if (p1 mod p)=0 then: print(i):
m1:=(p-(p mod 3))/3:m1:=m1*((p-1)/(p mod 3)):print(m1,(m1*3) mod p):
 x1:=m1*(38*i^6+2*i^5-i^4-872*i^3-1108*i^2-1309*i-713);
 x2:=m1*(842*i^6+5072*i^5-6847*i^4-46061*i^3-34930*i^2-52216*i+2878);
 x3:=m1*(4260971*i^6-3124108*i^5+2290532*i^4
 -99681839*i^3-46221667*i^2-106722952*i
 +5811547);
 print((x1 mod p),(x2 mod p),(x3 mod p));
 q1:=L (i,p);q2:=L (x1,p):q3:=L (x2,p):q4:=L (x3,p):q5:=L (-1,p):
  print(q1,q2,q3,q4,q5):
for j from 1 to (p-1)/2 do: 
print (j^2 mod p): od:print(0): fi: od:
 #to use this program, put in primes, m1 is 1/3 at that prime, we
#get roots mod p of charpoly, values of the 3 units, then a
#list of quadratic residues, hence we can write out the
#system composed of their Legendre symbols and test
# that it is independent; we insert the Legendre symbols of -1,t
#output data as follows: prime, value of t, units, Legendre symbols
#11  1      1,-1,1,-1,-1
#11 10     -1,-1,1,-1,-1
#17 2         1,-1,1,-1,1
#17 4         1,1,-1,1,1
#41 3         -1,-1,1,1,1
with(linalg):
m:=matrix(5,5,[0,1,0,1,0,0,0,1,0,0,0,1,0,1,1,1,1,0,1,1,1,1,0,0,0]);
 det(m);