# 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: