# homology # If A and B are matrices such that A*B=0 then homology(B,A) will # compute ker(A)/im(B), i.e. it will compute the torsion coefficients, # the rank and a set of generators. `normformnew` is a package which # includes `ismithex` and should be read in first. homology:=proc(B,A) local m,k,S1,S2,L1,L2,R1,R2,r,C,R1inv,L,tors,gen,i; k:=linalg[coldim](A); m:=linalg[coldim](B); S1:=ismithex(A,'L1','R1'); r:=linalg[rank](S1); if r=k then # print(`Homology is 0`);RETURN() # inserted else # C:=linalg[submatrix](evalm(R1&*B),r+1..k,1..m); S2:=ismithex(C,'L2','R2'); R1inv:=linalg[inverse](R1); L:=linalg[matrix](k,k,0); L:=linalg[copyinto](L2,L,r+1,1); L:=linalg[copyinto](linalg[matrix](r,r,1),L,1,k-r+1); L:=evalm(R1inv&*L); tors:=[]; for i to k-r while S2[i,i]<>0 do if S2[i,i]<>1 then tors:=[op(tors),S2[i,i]] fi od; if i=k-r+1 then # print(`Homology is 0`);RETURN() # inserted else # gen:=linalg[submatrix](L,1..k,i-nops(tors)..k-r); print(`Torsion coeffs`,tors,`Rank`,k-r-i+1,`Homology generators\ are columns of`); print(gen) fi fi end: